dfm-io / dfm.io Goto Github PK

Maybe time to switch to hugo or something?

Hi Dan -

I really appreciated the applied example of the exoplanet to showcase how and why someone might want to plug external C++ into Stan.
When adding better support for this feature in cmdstanpy I adapted the code in the blog, both to switch over from pystan and to update the C++ used to calculate the derivatives to work with more recent versions of Stan's math library.

The code is here: https://github.com/WardBrian/notebooks/blob/main/DFM-Exoplanets.ipynb

I'm not sure if you want to do anything with it, but at the very least I thought someone looking for an updated version might look at these issues

@dfm I was using your excellent notebook to learn about how to call C++ from PyStan. The features are in the develop branch of PyStan and gets picked up by conda, which is great.

I came across one line
if (abs((E - E0) / E) <= 1.234e-10) return E;
in the C++ code in cell 4 and cell 17 that caused a problem for my gcc/g++ stack (GCC 7.4).

Changing this line to read
if (fabs((E - E0) / E) <= 1.234e-10) return E;
fixes the problem.

Posting this, in case others run into it.

The origin of the issue is discussed in Stack Overflow here.
The use of abs() is strictly correct C++, so this isn't really a bug - but it may be useful to others to know this trick.

The C++ extension is excellent.

The post on Python modules in C is one of the best introductory post which I have ever seen on that topic. It will be great if it could be updated to work with Python3. To make it work, we need to do the following changes:

• Use PyModuleDef to define the module in _chi2.c

static PyModuleDef chi2_module = {
    PyModuleDef_HEAD_INIT,
    .m_name = "_chi2",
    .m_doc = module_docstring,
    .m_size = -1,
    .m_methods = module_methods,
};

• Use PyInit__chi2 instead of init_chi2, update the function to return NULL when the module does not get initialized and PyModule_Create instead of Py_InitModule3. The updated version:

PyMODINIT_FUNC PyInit__chi2(void)
{ 
  PyObject *m = PyModule_Create(&chi2_module);
  if (m == NULL) {
    return NULL; 
  }
  
  /* Load numpy functionality */
  import_array();

  return m;
}

Thank you for the wonderful post.