PyMultiNest icon indicating copy to clipboard operation
PyMultiNest copied to clipboard

dump_callback etc.

Open jtlz2 opened this issue 10 years ago • 16 comments

Hi - I'm trying to pass the current sample values to a python wrapper to plot them using pyplot as MultiNEST runs, rather than using show as in your demo script, because I hope it's more efficient than using show -> pdf.

I have a function to do the plot update, but I can't figure out how to extract the current sample values. Should I use dump_callback, or progressPrinter, or go via an ascii file..? How should I do that if so?

Thanks v much.

jtlz2 avatar Aug 07 '13 09:08 jtlz2

PS: I suppose I could add a routine to the prior or likelihood functions to do this - is that the best way forward?

jtlz2 avatar Aug 07 '13 10:08 jtlz2

Yes, the dump_callback would be the right way to do it. However, since I have not used it before, the implementation is not really tested. You will receive nsamples, nlive, n, and then what is in the postdist file as a 2-dimensional double array (nsamples x n). I would be happy if you could test it.

The benefit of writing to a pdf file and watching it with a viewer is that the output is permanent (does not go away if your code crashed), and that you are independent from graphical displays.

The second option is progressPrinter, which just reads the MultiNest output every so often. Note that there can be a partial output and so it might fail once in some read-outs, but it should be stable against that and just try again. You can build on e.g. ProgressPrinter, and read self.live using numpy.loadtxt

The third option is as you say, just output on every evaluation. This will probably make your code very slow, because also points that will be rejected are used then.

JohannesBuchner avatar Aug 07 '13 10:08 JohannesBuchner

Thanks for the speedy reply. Happy to test - I want to try dump_callback - could you give a simple example of the syntax? Something like:

pymultinest.run(..., dump_callback=func)

func(???): print ??? return

What would be an example of the syntax for ProgressPrinter (sounds a bit more complicated).

jtlz2 avatar Aug 07 '13 10:08 jtlz2

You can always start with

def func(*args): 
   print args 

but in this case, using ctypes to convert the array

def func(nsamples, nlive, n, postdist):
   print nsamples, nlive, n, postdist
   arr_type = ((ctypes.c_float * nsamples) * n)
   values = arr_type(postdist)
   print nsamples, nlive, n, values

JohannesBuchner avatar Aug 07 '13 10:08 JohannesBuchner

something along those lines anyways, you will have to try around a bit

JohannesBuchner avatar Aug 07 '13 10:08 JohannesBuchner

Thanks - so I get an error with dump_wrapper:

evidence_tolerance=0.5,mode_tolerance=-1e90,dump_callback=func)

File "/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/pymultinest-0.1-py2.7.egg/pymultinest/run.py", line 187, in run lib.set_dumper(dumper_type(dump_wrapper)) NameError: global name 'dump_wrapper' is not defined

I tried setting the name of func to dump_wrapper, but get the same.

jtlz2 avatar Aug 07 '13 10:08 jtlz2

That should be dump_callback instead of dump_wrapper

JohannesBuchner avatar Aug 07 '13 10:08 JohannesBuchner

(bug in https://github.com/JohannesBuchner/PyMultiNest/blob/master/pymultinest/run.py#L199 )

JohannesBuchner avatar Aug 07 '13 10:08 JohannesBuchner

Now I get (1000 live points):

Acceptance Rate: 0.847836 Replacements: 1900 Total Samples: 2241 ln(Z): -88079.423159 Acceptance Rate: 0.838710 Replacements: 1950 Total Samples: 2325 ln(Z): -77446.214381 converting for access ... --> DumperConverter called with: 2000 1000 4 converting for access ... done. dumper set to 0x1002fff60 calling client. calling 0x1002fff60 Illegal instruction

Could this take a while to debug? :-S

jtlz2 avatar Aug 07 '13 11:08 jtlz2

Try as a starting point

def func(nsamples, nlive, n, postdist):
   print 'python dumper callback called!'
   sys.stdout.flush()
   print nsamples, nlive, n
   sys.stdout.flush()

JohannesBuchner avatar Aug 07 '13 13:08 JohannesBuchner

Here you go (n_iter_before_update=10):

generating live points live points generated, starting sampling Acceptance Rate: 0.999049 Replacements: 1050 Total Samples: 1051 ln(Z): ************** converting for access ... --> DumperConverter called with: 1100 1000 4 converting for access ... done. dumper set to 0x100598f60 calling client. calling 0x100598f60 Segmentation fault

jtlz2 avatar Aug 07 '13 13:08 jtlz2

In the multinest bridge try simplifying

#define DUMPERTYPE(f) void (f)(int nsamples, int nlive, int npar, \
double ** physlive, double ** posterior, \
double *mean, double *std, double *best, double *map, \
double maxloglike, double logz, double logzerr)

to

#define DUMPERTYPE(f) void (f)(int nsamples, int nlive, int npar, \
double ** physlive)

and

p.Dumper(*nsamples, *nlive, n,
(double**)postdist, (double**)pLivePts,
mean, std, best, map, *maxloglike, *logz, *logzerr);

to

p.Dumper(*nsamples, *nlive, n, (double**)postdist);

JohannesBuchner avatar Aug 07 '13 13:08 JohannesBuchner

Same problem:


Starting MultiNest generating live points live points generated, starting sampling Acceptance Rate: 0.998099 Replacements: 1050 Total Samples: 1052 ln(Z): ************** converting for access ... --> DumperConverter called with: 1100 1000 4 converting for access ... done. dumper set to 0x100598f60 calling client. calling 0x100598f60 Segmentation fault

jtlz2 avatar Aug 07 '13 15:08 jtlz2

Try putting a few printf in before the p.Dumper() call to check whether *nsamples, *nlive, n, and postdist contain reasonable values (e.g. not NULL). You can also try replacing

dumper_type = CFUNCTYPE(c_void_p,
c_int, c_int, c_int, POINTER(c_double))

with

dumper_type = CFUNCTYPE(c_void_p, 
c_int, c_int, c_int, c_void_p)

in run.py

JohannesBuchner avatar Aug 07 '13 17:08 JohannesBuchner

*nsamples, *nlive and n are fine, so the problem must be with postdist (and hence with posterior)? It might be an overflow, but I'm not really a C programmer(!)...

jtlz2 avatar Aug 07 '13 19:08 jtlz2

You can check by replacing postdist in the p.Dumper call by NULL for the moment. printf("%p\n", postdist) would tell you the address.

JohannesBuchner avatar Aug 07 '13 19:08 JohannesBuchner