Commit 02398346 authored by Eric Jullo's avatar Eric Jullo
Browse files

Externalize bayesapp.c UserMonitor functions

parent b028a4b6
......@@ -28,5 +28,10 @@ typedef struct // INDIVIDUAL PARAMETERS
.....
} UserObjectStr;
*/
#ifdef __cplusplus
extern "C" FILE* open_mcmc_files(double cool);
extern "C" void write_bayes_line(FILE *bayes, double o_lhood_rez, double last_col, long int samp, UserCommonStr* UserCommon);
#endif
#endif
......@@ -2,6 +2,10 @@
#define LT_H_
#include <stdio.h>
#ifdef __cplusplus
extern "C" {
#endif
/* Functions declaration */
double *** rdf_cube_fits(char *filename,int *nx,int *ny,int *nz, char **pheader,double *lmin, double*lmax );
//double **rdf_fits(char *filename,int *nx,int *ny,char **pheader);
......@@ -66,6 +70,10 @@ void wrf_ipx(char *file,long int data,int dim,int *size,char *type,char *mode,ch
double xmin,double xmax,double ymin,double ymax);
void printerror( int status);
#ifdef __cplusplus
};
#endif
#define median_WIRTH(n,a) kth_smallest(a,n,(((n)&1)?((n)/2):(((n)/2)-1)))
#endif /*LT_H_*/
......@@ -3,12 +3,11 @@
# see http://dan.iel.fm/emcee/current/api/#the-affine-invariant-ensemble-sampler
# Parameters
nwalkers = 1000 # number of walkers
nwalkers = 100 # number of walkers
nburnin = 1000 # number of iterations for Burn-in phase
nsample = 1000 # number of iterations for Sample phase
nthin = 10 # write each nthin from sample during Sample phase
nthin_bayes = 10 # write each nthin from sample during burin phase
nthin = 100 # write each nthin from sample during Sample phase
url = "ipc:///tmp/sm_lenstool.test"
......@@ -28,12 +27,13 @@ sm.req_init_files();
# First, define the probability distribution that we would like to sample.
# Use a flat prior
def lnprob(x):
return sm.req_calc_lnprob(x);
lnp = sm.req_calc_lnprob(x);
return lnp
# Choose an initial set of positions for the walkers.
p0 = np.random.uniform(low=0, high=1.0, size=(nwalkers, ndim))
for i in xrange(nwalkers):
for i in range(nwalkers):
p0[i] = np.random.rand(ndim)
while (lnprob(p0[i]) == -np.inf):
p0[i] = np.random.rand(ndim)
......@@ -44,21 +44,21 @@ sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob)
# BURN-IN phase
iteration = 0
for pos, prob, state in sampler.sample(p0, iterations=nburnin, storechain=False):
print "burning: ", iteration;
print "Mean acceptance fraction:", np.mean(sampler.acceptance_fraction)
print "min/max prob=", np.amin(prob), np.amax(prob)
for pos, prob, state in sampler.sample(p0, iterations=nburnin, store=False):
print("burning: ", iteration)
print("Mean acceptance fraction:", np.mean(sampler.acceptance_fraction))
print("min/max prob=", np.amin(prob), np.amax(prob))
# ask server to write in burin.dat
if (iteration % nthin_bayes == 0):
print "Ask server to write in burnin.dat"
if (iteration % nthin == 0):
print("Ask server to write in burnin.dat")
for i, p in enumerate(prob):
sm.req_write_burnin(iteration, p, pos[i])
iteration = iteration + 1
print "\n"
print "Start sampling:"
print("\n")
print("Start sampling:")
# Reset the chain to remove the burn-in samples.
sampler.reset()
......@@ -66,14 +66,14 @@ sampler.reset()
# Sample phase
iteration = 0
for pos, prob, state in sampler.sample(pos, prob, state, iterations=nsample, storechain=False):
print "sampling: ", iteration;
print "Mean acceptance fraction:", np.mean(sampler.acceptance_fraction)
print "min/max prob=", np.amin(prob), np.amax(prob)
for pos, prob, state in sampler.sample(pos, prob, state, iterations=nsample, store=False):
print("sampling: ", iteration)
print("Mean acceptance fraction:", np.mean(sampler.acceptance_fraction))
print("min/max prob=", np.amin(prob), np.amax(prob))
# ask server to write in bayes.dat
if (iteration % nthin == 0):
print "Ask server to write in bayes.dat"
print("Ask server to write in bayes.dat")
for i, p in enumerate(prob):
sm.req_write_bayes(iteration, p, pos[i])
iteration = iteration + 1
......
......@@ -39,7 +39,7 @@ def logp(x):
# Choose an initial set of positions for the walkers.
p0 = np.random.uniform(low=0, high=1.0, size=(ntemps, nwalkers, ndim))
for t in range(ntemps):
p0[t] = [np.random.rand(ndim) for i in xrange(nwalkers)]
p0[t] = [np.random.rand(ndim) for i in range(nwalkers)]
for i in range(nwalkers):
while (lnprob(p0[t,i]) == -np.inf):
p0[t,i] = np.random.rand(ndim)
......@@ -52,19 +52,19 @@ sampler = emcee.PTSampler(ntemps, nwalkers, ndim, lnprob, logp)
# BURN-IN phase
iteration = 0
for pos, prob, state in sampler.sample(p0, iterations=nburnin, storechain=False):
print "burning: ", iteration;
print "Mean acceptance fraction:", np.mean(sampler.acceptance_fraction)
print "min/max prob=", np.amin(prob[0]), np.amax(prob[0])
print("burning: ", iteration)
print("Mean acceptance fraction:", np.mean(sampler.acceptance_fraction))
print("min/max prob=", np.amin(prob[0]), np.amax(prob[0]))
# ask server to write in burin.dat (only from first temperature level)
if (iteration % nthin_bayes == 0):
print "Ask server to write in burnin.dat"
print("Ask server to write in burnin.dat")
for i, p in enumerate(prob[0]):
sm.req_write_burnin(iteration, p, pos[0,i])
iteration = iteration + 1
print "\n"
print "Start sampling:"
print("\n")
print("Start sampling:")
# Reset the chain to remove the burn-in samples.
sampler.reset()
......@@ -73,13 +73,13 @@ sampler.reset()
# Sample phase
iteration = 0
for pos, prob, state in sampler.sample(pos, prob, state, iterations=nsample, storechain=False):
print "sampling: ", iteration;
print "Mean acceptance fraction:", np.mean(sampler.acceptance_fraction)
print "min/max prob=", np.amin(prob[0]), np.amax(prob[0])
print("sampling: ", iteration)
print("Mean acceptance fraction:", np.mean(sampler.acceptance_fraction))
print("min/max prob=", np.amin(prob[0]), np.amax(prob[0]))
# ask server to write in bayes.dat (only from first temperature level)
if (iteration % nthin == 0):
print "Ask server to write in bayes.dat"
print("Ask server to write in bayes.dat")
for i, p in enumerate(prob[0]):
sm.req_write_bayes(iteration, p, pos[0,i])
iteration = iteration + 1
......
......@@ -3,6 +3,7 @@ runmode
shearfield 0 1.24 nfw_shear.txt 1
verbose 0
inverse 6 ipc:///tmp/sm_lenstool.test
# inverse 3 0.1 10
mass 0 200 0.5 mass.fits
end
image
......
......@@ -58,7 +58,7 @@ class client:
def req_write_wtag(self, idx, lnprob, vx, tag):
# vx could be numpy array, so we must not use "+"
to_send = [idx] + [lnprob];
to_send = [float(idx)] + [lnprob];
for x in vx:
to_send.append(x)
self.send(tag, to_send)
......
This diff is collapsed.
......@@ -76,23 +76,36 @@ static void cpa2d_free(struct cpa2d* c)
free(c);
}
//
inline void cpa2d_set(struct cpa2d* c, int i, int j, double re, double im)
#ifndef DEBUG
inline
#endif
void cpa2d_set(struct cpa2d* c, int i, int j, double re, double im)
{
c->v[(i * c->ny + j) * 2] = re;
c->v[(i * c->ny + j) * 2 + 1] = im;
}
//
inline void cpa2d_set_re(struct cpa2d* c, int i, int j, double re)
#ifndef DEBUG
inline
#endif
void cpa2d_set_re(struct cpa2d* c, int i, int j, double re)
{
c->v[(i * c->ny + j) * 2] = re;
}
//
inline double cpa2d_get_re(struct cpa2d* c, int i, int j)
#ifndef DEBUG
inline
#endif
double cpa2d_get_re(struct cpa2d* c, int i, int j)
{
return c->v[(i * c->ny + j) * 2];
}
//
inline double cpa2d_get_im(struct cpa2d* c, int i, int j)
#ifndef DEBUG
inline
#endif
double cpa2d_get_im(struct cpa2d* c, int i, int j)
{
return c->v[(i * c->ny + j) * 2 + 1];
}
......
......@@ -242,21 +242,6 @@ int init_grille(char *infile, int noedit)
exit(E_PARFILE);
}
// Check that we do not pass MPMAX redshift planes
j = 0; double zcur = lens[0].z;
for( i = 0; i < G.nlens; i++ )
if ( lens[i].z > zcur )
{
zcur = lens[i].z;
j++;
}
if( j >= MPMAX )
{
fprintf(stderr, "ERROR: Too many redshift planes (%d). Increase MPMAX variable.\n", j+1);
exit(E_PARFILE);
}
NPRINTF(stderr, "Total number of clumps :%ld ", G.nlens);
NPRINTF(stderr, "(no_lens:%ld, ", G.no_lens);
for( i = 0; i < G.npot; i++ )
......
......@@ -282,6 +282,15 @@ static void scanInverse(char * third)
if ( pch != NULL && sscanf(pch, "%lf", &inverse1) == 1 )
M.itmax = (int) inverse1;
}
else if (M.inverse == 6) //server mode
{
pch = strtok( NULL, " ");
if (pch == NULL || !sscanf(pch, "%s", &M.servermode_url))
{
printf("Error! You should specify url of endpoint in service mode (inverse 6)\n");
exit(EXIT_FAILURE);
}
}
else
{
M.itmax = 1000; // Default number of posterior samples
......
......@@ -3,6 +3,8 @@
#include "sm_writeline.h"
#include <fonction.h>
#include <iostream>
#include "userstr.h"
using namespace std;
extern "C" {
......@@ -48,10 +50,10 @@ void run_servermode()
//clear bayes.dat and write header
bayesHeader();
//open burnin.dat and bayes.dat (we don't close them to speed-up output)
f_burnin = fopen("burnin.dat", "w");
f_bayes = fopen("bayes.dat", "a");
f_burnin = open_mcmc_files(0.);
f_bayes = open_mcmc_files(1.);
server.send(0, vector<double>()); //send ok to the client
}
else if (tag == 11) //print bayes line in burnin.dat
......
......@@ -30,7 +30,8 @@ void sm_server::recv(int& tag, vector<double>& v)
iss>>tag;
v.clear();
double d;
while (iss>>d)
while (iss>>d) {
v.push_back(d);
}
}
//
......@@ -16,61 +16,7 @@
#include "fonction.h"
#include "constant.h"
#include "lt.h"
/* Convert input Lenstool values to the corresponding values in <bayes.dat> file.
* bayes2lt() function in readBayesModels.c
*/
static double lt2bayes(int i, int ipx, double val)
{
extern struct pot lens[];
extern struct g_cosmo C;
switch (ipx)
{
case(THETA):
val *= RTD;
break;
case(B0):
// copied from o_print_res()
if ( lens[i].type <= 1 )
val = sqrt(val / 4. / pia_c2);
else if ( lens[i].type == 13 )
val = val * 1e4 * (cH0_4piG * C.h / distcosmo1(lens[i].z)); // in 10^8 Msol/kpc^2
else
val = sqrt(val / 6. / pia_c2);
break;
/* Everything in solar masses
case(PMASS):
if( lens[i].type == 12 )
val /= 1e14; // rhos in 10^14 Msol -> Msol
break;
case(MASSE):
if( lens[i].type == 12 )
val /= 1e14; // masse in 10^14 Msol -> Msol
break;
*/
default:
break;
}
return val;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: getParVal
//
// Purpose: Return the parameter value corresponding to the ilens and ipx parameters
//=============================================================================
static double getParVal(int i, int ipx)
{
double x;
x = o_get_lens(i, ipx);
return lt2bayes(i, ipx, x);
}
#include "userstr.h"
//
//reduced version of UserMonitor from bayesapp.c
......@@ -78,251 +24,39 @@ static double getParVal(int i, int ipx)
//index, lhood, parameters
void sm_writeline(FILE* bayes, const std::vector<double>& data, bool is_recalc_chi2)
{
int N = getNParameters();
if (data.size() != N + 2)
{
fprintf(stderr, "ERROR: sm_writeline: Wrong size of data\n");
exit(EXIT_FAILURE);
}
int idx = data[0];
double lhood = data[1];
std::vector<double> data2(data.begin() + 2, data.end());
rescaleCube_1Atom(data2.data(), N);
//if we need chi2, we have to recalculate o_chi_lhood0
double chi2 = 0 ;
double o_lhood_rez = lhood;
if (is_recalc_chi2)
{
double lhood0;
int error = o_chi_lhood0(&chi2, &lhood0, NULL); // TODO: catch error
o_lhood_rez = (-0.5*(chi2 + lhood0));
if (fabs(o_lhood_rez - lhood) > fabs(1e-10 * lhood))
printf("WARN: sm_writeline: lhood and o_lhood_rez are different %g %g %g\n",
o_lhood_rez, lhood, fabs(o_lhood_rez - lhood));
}
fprintf( bayes, "%ld %g ", idx, o_lhood_rez);
// SAVE: save lens parameters
int iparam = 0;
// long int samp = UserCommon->Nsample * Common->ENSEMBLE + k + 1;
extern int block[][NPAMAX];
double val = 0.;
//bad C-like style from UserMonitor
int i, j, ipx;
extern struct g_mode M;
extern struct g_grille G;
for ( i = 0; i < G.nmsgrid; i++ )
{
for ( ipx = CX; ipx <= PMASS; ipx++ )
if ( block[i][ipx] != 0 )
{
val = getParVal(i, ipx);
fprintf(bayes, "%lf ", val );
// addStat(val, UserCommon, iparam++, samp);
}
}
int nzeros= 0;
for ( i = 0; i < G.nlens - G.nmsgrid; i++ )
{
const extern double *np_b0;
val = np_b0[i];
if ( fabs(val) > 1E-6 )
{
if( nzeros > 0 )
fprintf(bayes, "%dx0.0 %lf ", nzeros, val);
else
fprintf(bayes, "%lf ", val);
nzeros = 0;
}
else
nzeros++;
// addStat(val, UserCommon, iparam++, samp);
}
if ( nzeros > 0 )
fprintf(bayes, "%dx0.0 ", nzeros);
/* // SAVE: source positions for grid model (for DEBUG purposes for now)
extern struct g_image I;
if ( Common->Valency && I.n_mult != 0)
{
double *Ps = (double *)calloc(I.n_mult * 2, sizeof(double));
long int ilens;
for ( i = 0; i < Objects[k].Natoms; i++ )
{
ilens = (int)((G.nlens - G.nmsgrid + I.n_mult) * Cubes[i][Common->Ndim - 1]);
if( ilens >= G.nlens - G.nmsgrid )
{
int isrc = ilens - (G.nlens - G.nmsgrid);
Ps[isrc] += Cubes[i][Common->Ndim + 1];
Ps[isrc + I.n_mult] += Cubes[i][Common->Ndim + 2];
}
}
for ( i = 0; i < I.n_mult * 2; i++ )
fprintf(bayes, "%lf ", Ps[i]);
free(Ps);
}*/
// SAVE: source parameters for M.iclean == 2
if ( M.iclean == 2 )
{
extern int sblock[NFMAX][NPAMAX];
extern struct g_source S;
extern struct galaxie source[NFMAX];
for( i = 0; i < S.ns; i++ )
for( ipx =SCX; ipx <= SFLUX; ipx++ )
if( sblock[i][ipx] != 0 )
{
val = o_get_source(&source[i], ipx);
if( ipx == STHETA ) val *= RTD;
fprintf(bayes, "%f ", val);
// addStat(val, UserCommon, iparam++, samp);
}
}
// SAVE: save cosmological parameters
extern int cblock[NPAMAX];
for ( ipx = OMEGAM; ipx <= WA; ipx++ )
if ( cblock[ipx] != 0 )
{
val = getParVal(0, ipx);
fprintf(bayes, "%lf ", val);
// addStat(val, UserCommon, iparam++, samp);
}
// SAVE: save zmlimit parameters
extern struct g_image I;
extern struct z_lim zlim[];
extern struct galaxie multi[NFMAX][NIMAX];
char limages[ZMBOUND][IDSIZE];
int nimages;
for ( ipx = 0; ipx < I.nzlim; ipx++ )
if ( zlim[ipx].bk != 0 )
{
// look for the images families corresponding to the zmlimit to optimize
nimages = splitzmlimit(zlim[ipx].n, limages);
i = 0;
while ( indexCmp( multi[i][0].n, limages[0] ) ) i++;
val = multi[i][0].z;
fprintf( bayes, "%lf ", val );
// addStat(val, UserCommon, iparam++, samp);
/* i = 0;
while( indexCmp( multi[i][0].n, zlim[ipx].n ) ) i++;
for( k = 0; k < I.mult[i]; k++ )
fprintf( bayes, "%lf ", multi[i][k].z );
*/
}
// SAVE: save z_a_limit parameter
extern struct z_lim zalim;
if( zalim.bk != 0 )
{
val = I.zarclet;
fprintf( bayes, "%lf ", val);
// addStat(val, UserCommon, iparam++, samp);
}
// SAVE: save velocity field parameters
extern int vfblock[NPAMAX];
for ( ipx = VFCX; ipx <= VFSIGMA; ipx++ )
if ( vfblock[ipx] != 0 )
{
val = getParVal(0, ipx);
if( ipx == VFTHETA ) val *= RTD;
if( ipx == VFI ) val *= RTD;
fprintf(bayes, "%lf ", val);
// addStat(val, UserCommon, iparam++, samp);
}
// SAVE : save pot parameters
extern struct g_pot P[NPOTFILE];
for( i = 0; i < G.npot; i++ )
{
struct g_pot *pot = &P[i];
if ( pot->ftype != 0 )
{
if ( pot->ircut != 0 )
{
fprintf( bayes, "%lf ", pot->cut );
// addStat(pot->cut, UserCommon, iparam++, samp);
}
if ( pot->isigma != 0 )
{
fprintf( bayes, "%lf ", pot->sigma );
// addStat(pot->sigma, UserCommon, iparam++, samp);
}
if ( pot->islope != 0 )
{
fprintf( bayes, "%lf ", pot->slope );
// addStat(pot->slope, UserCommon, iparam++, samp);
}
if ( pot->ivdslope != 0 )
{
fprintf( bayes, "%lf ", pot->vdslope );
// addStat(pot->vdslope, UserCommon, iparam++, samp);
}
if ( pot->ivdscat != 0 )
{
fprintf( bayes, "%lf ", pot->vdscat );
// addStat(pot->vdscat, UserCommon, iparam++, samp);
}
if ( pot->ircutscat != 0 )
{
fprintf( bayes, "%lf ", pot->rcutscat );
// addStat(pot->rcutscat, UserCommon, iparam++, samp);
}
if ( pot->ia != 0 )
{
fprintf( bayes, "%lf ", pot->a );
// addStat(pot->rcutscat, UserCommon, iparam++, samp);
}
if ( pot->ib != 0 )
{
fprintf( bayes, "%lf ", pot->b );
// addStat(pot->rcutscat, UserCommon, iparam++, samp);
}
}
}
// SAVE: the nuisance parameters
extern struct sigposStr sigposAs;
if ( sigposAs.bk != 0 )
{
for ( i = 0; i < I.n_mult; i++ )
for ( j = 0; j < I.mult[i]; j++ )
{
val = sqrt(I.sig2pos[i][j]);
fprintf( bayes, "%lf ", val );
// addStat(val, UserCommon, iparam++, samp);
}
}
if ( I.dsigell != -1. )
{
val = sqrt(I.sig2ell);
fprintf( bayes, "%lf ", val );
// addStat(val, UserCommon, iparam++, samp);
}