Skip to content

GitLab

Commit 3f6b7ec8 authored by Eric Jullo's avatar Eric Jullo
Browse files

This shm version works with one time creation/join of slave threads. It takes... 

This shm version works with one time creation/join of slave threads. It takes twice as long as the std shm version with as many creation/join as operations to process.
parent 3c40a406
Expand all Hide whitespace changes
Inline Side-by-side
Makefile.am
View file @ 3f6b7ec8
......@@ -4,71 +4,49 @@ noinst_LIBRARIES = liblenstool.a
AM_CFLAGS = @CFITSIO_HDR@ @WCS_HDR@ -I@top_srcdir@/include \
-I@top_srcdir@/liblt
liblenstool_a_SOURCES = \
al_sq_point.c amplif.c amplif_mat.c amplif_matinv.c \
bayesapp.c bayesys3.c bicubic_coef.c bicubic_int.c \
box_caustic.c brent.c chi_invim.c chsigne.c \
classer.c cleanlens.c comp_chi_os.c comp_chi_osv.c \
comp_dchi_osv.c copyright.c cor_seeing.c cor_shear.c \
cosmoratio.c cp_errdiff.c cp_images.c cp_im.c \
crea_filtre.c critic_an.c criticinv.c criticnew.c \
cv_cpsf.c d_binning.c d_bruiter.c d_gauss.c \
diag.c diff_mag.c d_integrer.c dist.c \
dist_min.c distor.c do_itos.c d_poisson.c \
d_profil.c d_rndschechter.c d_rndtype.c d_rndz.c \
d_seeing.c e_amp.c ecrire_r.c e_dpl.c \
e_g2cpx.c e_gcpx.c e_giant.c e_grad2.c \
e_grad.c e_im_prec.c e_inthere.c e_lens.c \
e_lensing.c e_lens_P.c e_mag.c e_mass.c \
e_nfw.c e_nfwgt.c e_pcpx.c e_pixel.c \
e_pot.c err_invim.c e_sersic.c e_tau.c \
e_test.c e_testg.c e_test_P.c e_time.c \
e_unlens.c e_unlens_fast.c e_unlensgrid.c e_unmag.c \
e_zeroamp.c fft.c fftcc_im.c fftc_im.c \
fmin_ell.c follow.c followi.c formeli.c \
frprmn.c fr_sq_point.c f_shape2.c f_shape3.c \
f_shape4.c f_shape_abs.c f_shape.c f_shmap.c \
f_source.c fz_dlsds.c g_ampli.c g_amplif.c \
g_curv.c g_dpl.c g_grid.c g_mass.c \
g_poten.c g_profil.c g_prop.c g_radial.c \
grid.c gridp.c gs_ecrire.c gs_ecrire_t.c \
g_shear.c g_shearf.c g_time.c hilbert.c \
ic_product.c im_add_const.c i_marker.c im_bin.c \
im_convolve.c im_stat.c imtosou.c inconvexe.c \
init_grille.c inside.c interpol.c inverse.c \
isoima.c iter.c keep_cl.c lire.c \
local.c lt_math.c min_parabol.c \
min_slo_par.c multiscale_grid.c nimage.c norm_angle.c \
o_big_slope.c o_chi.c o_chi_flux.c o_chi_pos.c \
o_chires.c o_dpl.c o_flux.c o_global.c \
o_invim.c o_keep_min.c o_keepz_min.c o_kp_sp_min.c \
o_mag.c o_mag_m.c o_min_loc.c o_min_slope.c \
o_prep.c o_prep_mult.c o_print.c o_print_res.c \
opt_source.c o_rescale.c o_run1.c o_run2.c \
o_run_bayes.c o_run.c o_run_mc0.c o_run_mc.c \
o_runpot1.c o_runpot2.c o_runz0.c o_runz2.c \
o_runz.c o_runz_mc0.c o_runz_mc.c o_set_err.c \
o_set_exc.c o_set_lens_bayes.c o_set_lens.c o_set_lmax.c \
o_set_lmin.c o_set_map.c o_set_map_mc.c o_set_map_z.c \
o_set_start.c o_shape.c o_slope_sp.c o_stat.c \
o_step.c o_swi_big.c polxy.c pro_arclet.c \
pro_arclet_m.c ptau.c random.c random_lt.c \
r_cleanlens.c r_cline.c r_cosmolimit.c r_cosmologie.c \
readBayesModels.c read.c readimage.c read_lenstable.c \
r_frame.c r_grille.c r_image.c r_large.c \
r_limit.c r_msgrid.c r_observ.c rotation.c \
rotmatrix.c r_potentiel.c r_potfile2.c r_potfile.c \
r_runmode.c r_source.c s_compmag.c set_default.c \
set_lens.c set_lens_par.c set_potfile2.c set_potfile.c \
set_res_par.c sig2posS4j.c sig2posS.c sig2posSe.c \
sig2posSj.c signe.c slope.c sort.c \
sortf.c spline.c split_image.c sp_main.c \
sp_set.c s_sof.c s_sourcebox.c s_source.c \
st_opt.c study.c tirer.c tracepot.c \
unlens1.c unlens_bc.c updatecut.c update_emass.c \
update_epot.c w_critic.c weight_baryc.c w_gianti.c \
w_prop.c wr_arclet.c wr_mass.c wr_pot.c \
w_sicat.c zero.c zero_t.c zonemult.c
liblenstool_a_SOURCES = \
al_sq_point.c amplif.c amplif_mat.c amplif_matinv.c bayesapp.c \
bayesys3.c bicubic_coef.c bicubic_int.c brent.c chi_invim.c chsigne.c \
classer.c cleanlens.c copyright.c cor_seeing.c cor_shear.c \
cosmoratio.c cp_errdiff.c cp_images.c crea_filtre.c critic_an.c \
criticinv.c criticnew.c d_binning.c d_bruiter.c Debug d_gauss.c diag.c \
diff_mag.c d_integrer.c dist.c dist_min.c distor.c do_itos.c \
d_poisson.c d_profil.c d_rndschechter.c d_rndtype.c d_rndz.c \
d_seeing.c e_amp.c ecrire_r.c e_dpl.c e_g2cpx.c e_gcpx.c e_giant.c \
e_grad2.c e_grad.c e_im_prec.c e_inthere.c e_lens.c e_lensing.c \
e_lens_P.c e_mag.c e_mass.c e_nfw.c e_nfwgt.c e_pcpx.c e_pixel.c \
e_pot.c err_invim.c e_sersic.c e_tau.c e_test.c e_testg.c e_test_P.c \
e_time.c e_unlens.c e_unlens_fast.c e_unlensgrid.c e_unmag.c \
e_zeroamp.c fft.c fftcc_im.c fftc_im.c fmin_ell.c follow.c followi.c \
formeli.c frprmn.c fr_sq_point.c f_shape2.c f_shape3.c f_shape4.c \
f_shape_abs.c f_shape.c f_shmap.c fz_dlsds.c g_ampli.c g_amplif.c \
g_curv.c g_dpl.c g_grid.c g_mass.c g_poten.c g_profil.c g_prop.c \
g_radial.c grid.c gridp.c gs_ecrire.c gs_ecrire_t.c g_shear.c \
g_shearf.c g_time.c hilbert.c ic_product.c im_add_const.c i_marker.c \
im_bin.c im_convolve.c im_stat.c imtosou.c inconvexe.c init_grille.c \
inside.c interpol.c inverse.c isoima.c iter.c keep_cl.c liblenstool.a \
lire.c local.c lt_math.c main.c main.o Makefile Makefile.am \
Makefile.in min_parabol.c min_slo_par.c multiscale_grid.c nimage.c \
norm_angle.c *.o o_big_slope.c o_chi.c o_chi_flux.c o_chires.c o_dpl.c \
o_flux.c o_global.c o_invim.c o_keep_min.c o_keepz_min.c o_kp_sp_min.c \
o_mag.c o_mag_m.c o_min_loc.c o_min_slope.c o_prep.c o_prep_mult.c \
o_print.c o_print_res.c opt_source.c o_rescale.c o_run1.c o_run2.c \
o_run_bayes.c o_run.c o_run_mc0.c o_run_mc.c o_runpot1.c o_runpot2.c \
o_runz_mc0.c o_runz_mc.c o_set_err.c o_set_exc.c o_set_lens_bayes.c \
o_set_lens.c o_set_lmax.c o_set_lmin.c o_set_map.c o_set_map_mc.c \
o_set_map_z.c o_set_start.c o_shape.c o_slope_sp.c o_stat.c o_step.c \
o_swi_big.c polxy.c pro_arclet.c pro_arclet_m.c ptau.c random.c \
random_lt.c r_cleanlens.c r_cline.c r_cosmolimit.c r_cosmologie.c \
readBayesModels.c read.c readimage.c read_lenstable.c r_frame.c \
r_grille.c r_image.c r_large.c r_limit.c r_msgrid.c r_observ.c \
rotation.c rotmatrix.c r_potentiel.c r_potfile.c r_runmode.c \
r_source.c s_compmag.c set_default.c set_lens.c set_lens_par.c \
set_potfile.c set_res_par.c sig2posS4j.c sig2posS.c sig2posSe.c \
sig2posSj.c signe.c slope.c sort.c sortf.c spline.c split_image.c \
sp_main.c sp_set.c s_sof.c s_sourcebox.c s_source.c st_opt.c study.c \
tirer.c tracepot.c unlens1.c unlens_bc.c updatecut.c update_emass.c \
update_epot.c weight_baryc.c w_gianti.c wr_arclet.c wr_mass.c wr_pot.c \
w_sicat.c zero.c zero_t.c
bin_PROGRAMS = lenstool
......@@ -76,4 +54,4 @@ lenstool_SOURCES = main.c
lenstool_LDFLAGS = @CFITSIO_LIB@ @WCS_LIB@
lenstool_LDADD = liblenstool.a ../liblt/liblt.a -lcfitsio -lwcs -lm
lenstool_LDADD = liblenstool.a ../liblt/liblt.a -lcfitsio -lwcs -lpthread -lm
Makefile.in
View file @ 3f6b7ec8
This diff is collapsed.
amplif.c
View file @ 3f6b7ec8
#include<stdio.h>
#include<math.h>
#include<fonction.h>
#include<constant.h>
#include<dimension.h>
#include<structure.h>
#include "fonction.h"
#include "constant.h"
#include"dimension.h"
#include "structure.h"
/****************************************************************/
/* nom: amplif */
......@@ -22,28 +22,21 @@
* - in e_amp() : G, lens, lens_table
* - in distcosmo1() : C
*/
void amplif()
void amplif(int model, double amplifi[NFMAX][NIMAX])
{
extern struct g_image I;
extern struct galaxie multi[NIMAX][NIMAX];
extern double amplifi[NIMAX][NIMAX];
extern struct galaxie multi[NFMAX][NIMAX];
double dlsds,A;
register int i,j;
int n;
/*For each image*/
for(i=0;i<I.n_mult;i++)
{
n=I.mult[i];
dlsds=multi[i][0].dr;
/* NPRINTF(stderr,"n=%d dlsds=%.3lf\n",n,dlsds);*/
for(j=0;j<n;j++)
dlsds=multi[i][0].dr[model];
for(j=0;j<I.mult[i];j++)
{
A=1./e_amp_gal(&multi[i][j],dlsds);
A=1./e_amp_gal(model,&multi[i][j],dlsds);
amplifi[i][j]=A;
/* NPRINTF(stderr," A[%d][%d]=%.3lf\n ",i,j,amplifi[i][j]); */
};
};
}
}
}
amplif_mat.c
View file @ 3f6b7ec8
#include<stdio.h>
#include<math.h>
#include<fonction.h>
#include<constant.h>
#include<dimension.h>
#include<structure.h>
#include "fonction.h"
#include "constant.h"
#include"dimension.h"
#include "structure.h"
/****************************************************************/
/* nom: amplif_mat */
......@@ -19,33 +19,27 @@
* - in distcosmo1() : C
* - in e_amp() : G, lens, lens_table
*/
void amplif_mat()
void amplif_mat(int model, struct matrix amplifi_mat[NFMAX][NIMAX] )
{
extern struct g_image I;
extern struct galaxie multi[NIMAX][NIMAX];
extern struct matrix amplifi_mat[NIMAX][NIMAX];
extern struct galaxie multi[NFMAX][NIMAX];
struct matrix Mat;
double dlsds,A;
register int i,j;
int n;
/* NPRINTF(stderr,"Coucou.......\n"); */
for(i=0;i<I.n_mult;i++)
{
n=I.mult[i];
dlsds=multi[i][0].dr;
/* NPRINTF(stderr,"n=%d dlsds=%.3lf\n",n,dlsds);*/
for(j=0;j<n;j++)
dlsds=multi[i][0].dr[model];
for(j=0;j<I.mult[i];j++)
{
Mat=e_grad2_gal(&multi[i][j]);
Mat=e_grad2_gal(model,&multi[i][j]);
Mat.a*=dlsds;
Mat.b*=dlsds;
Mat.c*=dlsds;
Mat.d*=dlsds;
A=e_amp_gal(&multi[i][j],dlsds);
A=e_amp_gal(model,&multi[i][j],dlsds);
amplifi_mat[i][j].a=(1.-Mat.c)/A;
amplifi_mat[i][j].b=-Mat.b/A;
amplifi_mat[i][j].c=(1.-Mat.a)/A;
......@@ -55,6 +49,6 @@ void amplif_mat()
NPRINTF(stderr," A[%d][%d].b=%.3lf\n ",i,j,amplifi_mat[i][j].b);
NPRINTF(stderr," A[%d][%d].c=%.3lf\n ",i,j,amplifi_mat[i][j].c);
NPRINTF(stderr," A[%d][%d].d=%.3lf\n ",i,j,amplifi_mat[i][j].d); */
};
};
}
}
}
amplif_matinv.c
View file @ 3f6b7ec8
#include<stdio.h>
#include<math.h>
#include<fonction.h>
#include<constant.h>
#include<dimension.h>
#include<structure.h>
#include "fonction.h"
#include "constant.h"
#include"dimension.h"
#include "structure.h"
/****************************************************************/
/* nom: amplif_matinv */
......@@ -19,27 +19,21 @@
* - in distcosmo1() : C
* - in e_amp() : G, lens, lens_table
*/
void amplif_matinv()
void amplif_matinv(int model, struct matrix amplifi_matinv[NFMAX][NIMAX])
{
extern struct g_image I;
extern struct galaxie multi[NIMAX][NIMAX];
extern struct matrix amplifi_matinv[NIMAX][NIMAX];
extern struct galaxie multi[NFMAX][NIMAX];
struct matrix Mat;
double dlsds;
register int i,j;
int n;
/* NPRINTF(stderr,"Coucou.......\n"); */
for(i=0;i<I.n_mult;i++)
{
n=I.mult[i];
dlsds=multi[i][0].dr;
/* NPRINTF(stderr,"n=%d dlsds=%.3lf\n",n,dlsds);*/
for(j=0;j<n;j++)
dlsds=multi[i][0].dr[0];
for(j=0;j<I.mult[i];j++)
{
Mat=e_grad2_gal(&multi[i][j]);
Mat=e_grad2_gal(0,&multi[i][j]);
Mat.a*=dlsds;
Mat.b*=dlsds;
Mat.c*=dlsds;
......
bayesapp.c
View file @ 3f6b7ec8
......@@ -34,12 +34,10 @@
#include "fonction.h"
#include "constant.h"
double getLhood0(); // defined in o_chi.c
static int minusone; // count chi2 errors in image plane
static double nchi2; // count number of computed chi2
static time_t start,end; // computation time
static int UserBuild(double*, CommonStr*, ObjectStr*, int, int);
static int UserBuild(double*, CommonStr*, ObjectStr*, int, int, int);
//=============================================================================
// SIMPLE CODE
......@@ -58,9 +56,10 @@ static int UserBuild(double*, CommonStr*, ObjectStr*, int, int);
int UserEmpty( // O >=0, or -ve error code
double* Lhood, // O loglikelihood
CommonStr* Common, // I O general information
ObjectStr* Object) // I O sample object, output new Lhood
ObjectStr* Object, // I O sample object, output new Lhood
int slave) // I ID of the processor to work with
{
UserBuild(Lhood, Common, Object, 0, 1);
UserBuild(Lhood, Common, Object, 0, 1,slave);
return 1;
}
......@@ -83,7 +82,8 @@ ObjectStr* Object) // I O sample object, output new Lhood
int UserTry1( // O +ve = OK, 0 = DO NOT USE, -ve = error
double* dLtry, // O trial d(logLikelihood) value
CommonStr* Common, // I general information
ObjectStr* Object) // I sample object (DO NOT UPDATE Lhood)
ObjectStr* Object, // I sample object (DO NOT UPDATE Lhood)
int slave) // I ID of the processor to work with
{
int Natoms = Object->Natoms;
double Lhood = Object->Lhood; // existing Lhood
......@@ -91,7 +91,7 @@ ObjectStr* Object) // I sample object (DO NOT UPDATE Lhood)
int OK;
*dLtry = 0.0;
OK = UserBuild(&Ltry, Common, Object, Natoms+1, 0); // trial
OK = UserBuild(&Ltry, Common, Object, Natoms+1, 0,slave); // trial
if( OK > 0 )
*dLtry = Ltry - Lhood; // increment
return OK; // OK?
......@@ -123,7 +123,8 @@ int UserTry2( // O +ve = OK, 0 = DO NOT USE, -ve = error
double* dLtry1, // O trial d(logLikelihood) value for 1st atom
double* dLtry2, // O trial d(logLikelihood) value for both atoms
CommonStr* Common, // I general information
ObjectStr* Object) // I sample object (DO NOT UPDATE Lhood)
ObjectStr* Object, // I sample object (DO NOT UPDATE Lhood)
int slave) // I ID of the processor to work with
{
int Natoms = Object->Natoms;
double Lhood = Object->Lhood; // existing Lhood
......@@ -131,11 +132,11 @@ ObjectStr* Object) // I sample object (DO NOT UPDATE Lhood)
int OK;
*dLtry1 = *dLtry2 = 0.0;
OK = UserBuild(&Ltry1, Common, Object, Natoms+1, 0); //trial for 1 more
OK = UserBuild(&Ltry1, Common, Object, Natoms+1, 0,slave); //trial for 1 more
if( OK > 0 )
{
*dLtry1 = Ltry1 - Lhood; // increment for 1
OK = UserBuild(&Ltry2, Common, Object, Natoms+2, 0); //trial for 2 more
OK = UserBuild(&Ltry2, Common, Object, Natoms+2, 0,slave); //trial for 2 more
if( OK > 0 )
*dLtry2 = Ltry2 - Lhood; // increment for 2
}
......@@ -164,13 +165,14 @@ ObjectStr* Object) // I sample object (DO NOT UPDATE Lhood)
int UserInsert1( // O >=0, or -ve error code
double* dL, // O d(loglikelihood)
CommonStr* Common, // I general information
ObjectStr* Object) // I O sample object
ObjectStr* Object, // I O sample object
int slave) // I ID of the processor to work with
{
int Natoms = Object->Natoms; // new number
double Lold = Object->Lhood; // not yet updated
double Lnew;
UserBuild(&Lnew, Common, Object, Natoms, 1); // new updated state
UserBuild(&Lnew, Common, Object, Natoms, 1,slave); // new updated state
*dL = Lnew - Lold; // I will update Object->Lhood
return 1;
}
......@@ -199,13 +201,14 @@ ObjectStr* Object) // I O sample object
int UserInsert2( // O >=0, or -ve error code
double* dL, // O d(loglikelihood)
CommonStr* Common, // I general information
ObjectStr* Object) // I O sample object
ObjectStr* Object, // I O sample object
int slave) // I ID of the processor to work with
{
int Natoms = Object->Natoms; // new number
double Lold = Object->Lhood; // not yet updated
double Lnew;
UserBuild(&Lnew, Common, Object, Natoms, 1); // new updated state
UserBuild(&Lnew, Common, Object, Natoms, 1,slave); // new updated state
*dL = Lnew - Lold; // I will update Object->Lhood
return 1;
}
......@@ -226,13 +229,14 @@ ObjectStr* Object) // I O sample object
int UserDelete1( // O >=0, or -ve error code
double* dL, // O d(loglikelihood)
CommonStr* Common, // I general information
ObjectStr* Object) // I O sample object
ObjectStr* Object, // I O sample object
int slave) // I ID of the processor to work with
{
int Natoms = Object->Natoms; // new number
double Lold = Object->Lhood; // not yet updated
double Lnew;
UserBuild(&Lnew, Common, Object, Natoms, 1); // new updated state
UserBuild(&Lnew, Common, Object, Natoms, 1, slave); // new updated state
*dL = Lnew - Lold; // I will update Object->Lhood
return 1;
}
......@@ -253,10 +257,12 @@ double* Lhood, // O loglikelihood
CommonStr* Common, // I General information
ObjectStr* Object, // I(O) Cubes in, perhaps Mock out
int Natoms, // I # atoms
int output) // I +ve if Mock to be written out, -ve if not
int output, // I +ve if Mock to be written out, -ve if not
int slave) // I ID of the processor to work with
{
double** Cubes = Object->Cubes; // I Cubes in [0,1) [Natoms][Ndim]
double* Cube=Cubes[0];
double chi2,lhood0;
int valid;
UserCommonStr* UserCommon = Common->UserCommon;
......@@ -266,21 +272,21 @@ int output) // I +ve if Mock to be written out, -ve if not
if( Natoms == 1 )
{
/* Set the clumps parameters from the cube and eventually reject the whole object*/
valid = rescaleCube_1Atom(Cube,Common->Ndim);
valid = rescaleCube_1Atom(slave-1,Cube,Common->Ndim);
/* Compute the likelyhood for this object */
if( valid )
{
UserCommon->Nchi2++;
*Lhood = o_chi();
chi2 = o_chi(slave-1,&lhood0); // slave starts from 1
nchi2++;
if( *Lhood == -1 )
if( chi2 == -1 )
{
minusone++;
return 0;
}
*Lhood = -0.5 * (*Lhood);
*Lhood += -0.5 * getLhood0();
*Lhood = -0.5 * chi2 - lhood0;
// printf("%lf %lf %lf\n", Cube[0], chi2, lhood0);
}
}
return valid;
......@@ -374,6 +380,10 @@ char* getParName(int ipx, char* name, int type)
case(WX):
strcpy( name, "wX" );
break;
case(WA):
strcpy( name, "wa" );
break;
};
return name;
}
......@@ -382,9 +392,8 @@ char* getParName(int ipx, char* name, int type)
/* 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)
static double lt2bayes(int model,struct pot *ilens, int ipx, double val)
{
extern struct pot lens[];
extern struct g_cosmo C;
switch(ipx)
......@@ -394,21 +403,21 @@ static double lt2bayes(int i, int ipx, double val)
break;
case(B0):
// copied from o_print_res()
if ( lens[i].type <= 1 )
if ( ilens->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 if( ilens->type == 13 )
val=val*1e4*(cH0_4piG*C.h/distcosmo1(model,ilens->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 )
if( ilens->type == 12 )
val /= 1e14; // rhos in 10^14 Msol -> Msol
break;
case(MASSE):
if( lens[i].type == 12 )
if( ilens->type == 12 )
val /= 1e14; // masse in 10^14 Msol -> Msol
break;
*/
......@@ -425,14 +434,13 @@ static double lt2bayes(int i, int ipx, double val)
//
// Purpose: Return the parameter value corresponding to the ilens and ipx parameters
//=============================================================================
static double getParVal(int i, int ipx)
static double getParVal(int model,struct pot *ilens, int ipx)
{
extern struct pot lens[];
double x;
x = o_get_lens(i, ipx);
x = o_get_lens(model,ilens, ipx);
return lt2bayes(i, ipx, x);
return lt2bayes(model,ilens, ipx, x);
}
......@@ -468,7 +476,7 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
extern int block[][NPAMAX];
extern int cblock[NPAMAX];
extern struct z_lim zlim[];
extern struct galaxie multi[NIMAX][NIMAX];
extern struct galaxie multi[NFMAX][NIMAX];
extern int optInterrupt; // Global variable modified in o_run_bayes.c/signalReset()
extern struct sigposStr sigposAs;
......@@ -484,12 +492,13 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
// local
ObjectStr* Object; // individual object
int k, i, j; // counters
int ilens, ipx;
double chi2;
int ipx;
double chi2,lhood0;
FILE *bayes;
int nimages;
double chi2rate;
char limages[ZMBOUND][IDSIZE];
int model0 = 0; // in UserMonitor always with model0
// FILE *debug;
//.............................................................................
......@@ -512,8 +521,8 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
if( Common->cool >= 1.0 && M.inverse == 3 )
Common->cool = 1.0;
rescaleCube_1Atom(Objects[0].Cubes[0], Common->Ndim);
chi2 = o_chi();
rescaleCube_1Atom(model0,Objects[0].Cubes[0], Common->Ndim);
chi2 = o_chi(model0,&lhood0);
// if( getLhood0() != 0. ) chi2 = -0.5 * ( chi2 + getLhood0() );
if( Common->cool < 1.0 )
......@@ -557,25 +566,26 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
for( j = 0; j < Common->ENSEMBLE; j++ )
{
// rescale the cube of object_j
rescaleCube_1Atom(Objects[j].Cubes[0], Common->Ndim);
chi2 = o_chi();
if( getLhood0() != 0. ) chi2 = -0.5 * ( chi2 + getLhood0() );
rescaleCube_1Atom(model0,Objects[j].Cubes[0], Common->Ndim);
chi2 = o_chi(model0,&lhood0);
if( sigposAs.bk != 0 || I.dsigell != -1 )
chi2 = -0.5 * chi2 - lhood0;
fprintf( bayes, "%d %lf ", UserCommon->Nsample+1, chi2);
// SAVE: save lens parameters
for( ilens = 0; ilens < G.no_lens; ilens++ )
for( i = 0; i < G.no_lens; i++ )
for( ipx = CX; ipx <= PMASS; ipx++ )
if( block[ilens][ipx] != 0 )
fprintf(bayes, "%lf ", getParVal(ilens, ipx));
if( block[i][ipx] != 0 )
fprintf(bayes, "%lf ", getParVal(model0,&lens[i], ipx));
for( ilens = G.nmsgrid; ilens < G.nlens; ilens++ )
fprintf(bayes, "%lf ", getParVal(ilens, B0));
for( i = G.nmsgrid; i < G.nlens; i++ )
fprintf(bayes, "%lf ", getParVal(model0,&lens[i], B0));
// SAVE: save cosmological parameters
for( ipx = OMEGAM; ipx <= WX; ipx++ )
for( ipx = OMEGAM; ipx <= WA; ipx++ )
if( cblock[ipx] != 0 )
fprintf(bayes, "%lf ", getParVal(0, ipx));
fprintf(bayes, "%lf ", getParVal(model0,NULL, ipx));
// SAVE: save zmlimit parameters
for( ipx = 0; ipx < I.nzlim; ipx++ )
......@@ -585,7 +595,7 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
nimages = splitzmlimit(zlim[ipx].n, limages);
i = 0;
while( indexCmp( multi[i][0].n, limages[0] ) ) i++;
fprintf( bayes, "%lf ", multi[i][0].z );
fprintf( bayes, "%lf ", multi[i][0].z[model0] );
/* i = 0;
while( indexCmp( multi[i][0].n, zlim[ipx].n ) ) i++;
for( k = 0; k < I.mult[i]; k++ )
......@@ -597,27 +607,27 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
if( P.ftype != 0 )
{
if( P.ircut != 0 )
fprintf( bayes, "%lf ", P.cut );
fprintf( bayes, "%lf ", P.cut[model0] );
if( P.isigma != 0 )
fprintf( bayes, "%lf ", P.sigma );
fprintf( bayes, "%lf ", P.sigma[model0] );
if( P.islope != 0 )
fprintf( bayes, "%lf ", P.slope );
fprintf( bayes, "%lf ", P.slope[model0] );
if( P.ivdslope != 0 )
fprintf( bayes, "%lf ", P.vdslope );
fprintf( bayes, "%lf ", P.vdslope[model0] );
if( P.ivdscat != 0 )
fprintf( bayes, "%lf ", P.vdscat );
fprintf( bayes, "%lf ", P.vdscat[model0] );
if( P.ircutscat != 0 )
fprintf( bayes, "%lf ", P.rcutscat );
fprintf( bayes, "%lf ", P.rcutscat[model0] );
}
// SAVE: the nuisance parameters
if( sigposAs.bk != 0 )
{
for( i = 0; i < I.n_mult; i++ )
for( k = 0; k < I.mult[i]; k++ )
fprintf( bayes, "%lf ", sqrt(I.sig2pos[i][k]) );
fprintf( bayes, "%lf ", sqrt(multi[i][k].sig2pos[model0]) );
}
if( I.dsigell != -1. )
fprintf( bayes, "%lf ", sqrt(I.sig2ell) );
fprintf( bayes, "%lf ", sqrt(I.sig2ell[model0]) );
/*
// SAVE : save the gradx and grady values of the multiple images
......
bayesys3.c
View file @ 3f6b7ec8
This diff is collapsed.