Skip to content

GitLab

Commit 542d87fd authored by Johan Richard's avatar Johan Richard
Browse files

Added double shapemodel parameters in Bayes Optimisation

parent f77c5a73
Hide whitespace changes
Inline Side-by-side
src/bayesapp.c
View file @ 542d87fd
......@@ -784,6 +784,21 @@ char* getParName(int ipx, char* name, int type)
case(SFLUX):
strcpy( name, "mag" );
break;
case(SA2):
strcpy( name, "a2 (arcsec)" );
break;
case(SEPS2):
strcpy( name, "eps2" );
break;
case(STHETA2):
strcpy( name, "theta2 (deg)" );
break;
case(SINDEX2):
strcpy( name, "index2" );
break;
case(SFLUX2):
strcpy( name, "mag2" );
break;
case(VFCX):
strcpy( name, "velocity x (arcsec)");
break;
......@@ -1161,11 +1176,12 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
extern struct galaxie source[NFMAX];
for( i = 0; i < S.ns; i++ )
for( ipx =SCX; ipx <= SFLUX; ipx++ )
for( ipx =SCX; ipx <= SFLUX2; ipx++ )
if( sblock[i][ipx] != 0 )
{
val = o_get_source(&source[i], ipx);
if( ipx == STHETA ) val *= RTD;
if( ipx == STHETA2 ) val *= RTD;
fprintf(bayes, "%f ", val);
addStat(val, UserCommon, iparam++, samp);
}
......
src/checkpar.c
View file @ 542d87fd
......@@ -73,6 +73,32 @@ void checkpar()
source[i].E.b=source[i].E.a;
}
}
if((source[i].eps2>0)||(source[i].E2.a>0))
{
if(source[i].type!=6)
{
fprintf(stderr, "ERROR: wrong type (should be 6) for 2nd Sersic profile in source %d\n", i);
exit(-1);
}
if(source[i].eps2>=0)
{
// second size has been defined in shapemodel
if(source[i].E2.a>0)
source[i].E2.b = source[i].E2.a * (1. - source[i].eps2) / (1. + source[i].eps2);
else
{
fprintf(stderr, "ERROR: 2nd size not defined in source %d type 6\n", i);
exit(-1);
}
}
else // only size defined
{
fprintf(stderr, "WARNING: second ellipticity for source %d not defined. Assuming 0\n", i);
source[i].eps2 = 0;
source[i].E2.b=source[i].E2.a;
}
}
}
return;
......
src/o_global.c
View file @ 542d87fd
......@@ -941,7 +941,7 @@ int getNParameters()
// source parameters
for ( i = 0; i < S.ns; i++ )
for (ipx = SCX; ipx <= SFLUX; ipx++ )
for (ipx = SCX; ipx <= SFLUX2; ipx++ )
if ( sblock[i][ipx] != 0 )
parameters++;
......
src/o_print_res.c
View file @ 542d87fd
......@@ -285,6 +285,7 @@ void o_print_res(double chi0, double evidence)
// SHAPE MODEL
fprintf(best, "shapemodel\n");
fprintf(best, "\tid %s\n", source[i].n);
fprintf(best, "\ttype %d\n", source[i].type);
fprintf(best, "\ts_center_x %.6lf\n", source[i].C.x);
fprintf(best, "\ts_center_y %.6lf\n", source[i].C.y);
fprintf(best, "\ts_angle %.6lf\n", source[i].E.theta * RTD);
......@@ -292,8 +293,16 @@ void o_print_res(double chi0, double evidence)
fprintf(best, "\ts_sigy %.6lf\n", source[i].E.b);
fprintf(best, "\ts_eps %.6lf\n", source[i].eps);
fprintf(best, "\tmag %.6lf\n", source[i].mag);
fprintf(best, "\ttype %d\n", source[i].type);
fprintf(best, "\tindex %.6lf\n", source[i].var1);
if(source[i].type==6) // double Sersic
{
fprintf(best, "\ts_angle2 %.6lf\n", source[i].E2.theta * RTD);
fprintf(best, "\ts_sigx2 %.6lf\n", source[i].E2.a);
fprintf(best, "\ts_sigy2 %.6lf\n", source[i].E2.b);
fprintf(best, "\ts_eps2 %.6lf\n", source[i].eps2);
fprintf(best, "\tmag2 %.6lf\n", source[i].mag2);
fprintf(best, "\tindex2 %.6lf\n", source[i].var2);
}
fprintf(best, "\tend\n");
}
......@@ -572,16 +581,25 @@ void o_print_res(double chi0, double evidence)
// SHAPE MODEL
fprintf(besto, "shapemodel\n");
fprintf(besto, "\tid %s\n", source[i].n);
fprintf(besto, "\ttype %d\n", source[i].type);
fprintf(besto, "\ts_center_x %.6lf\n", source[i].C.x);
fprintf(besto, "\ts_center_y %.6lf\n", source[i].C.y);
fprintf(besto, "\ts_angle %.6lf\n", source[i].E.theta * RTD);
fprintf(besto, "\ts_sigx %.6lf\n", source[i].E.a);
fprintf(besto, "\ts_sigy %.6lf\n", source[i].E.b);
fprintf(besto, "\ts_eps %.6lf\n", source[i].eps);
fprintf(besto, "\tmag %.6lf\n", source[i].mag);
fprintf(besto, "\ttype %d\n", source[i].type);
fprintf(besto, "\ts_mag %.6lf\n", source[i].mag);
fprintf(besto, "\tindex %.6lf\n", source[i].var1);
fprintf(besto, "\tend\n");
if(source[i].type==6) //double Sersic
{
fprintf(besto, "\ts_angle2 %.6lf\n", source[i].E2.theta * RTD);
fprintf(besto, "\ts_sigx2 %.6lf\n", source[i].E2.a);
fprintf(besto, "\ts_sigy2 %.6lf\n", source[i].E2.b);
fprintf(besto, "\ts_eps2 %.6lf\n", source[i].eps2);
fprintf(besto, "\ts_mag2 %.6lf\n", source[i].mag2);
fprintf(besto, "\tindex2 %.6lf\n", source[i].var2);
fprintf(besto, "\tend\n");
}
// SHAPE LIMITS
fprintf(besto, "shapelimit\n");
......@@ -599,6 +617,19 @@ void o_print_res(double chi0, double evidence)
fprintf(besto, "\tmag %d %.6lf %.6lf\n", sblock[i][SFLUX], smin[i].mag, smax[i].mag);
if( sblock[i][SINDEX] )
fprintf(besto, "\tindex %d %.6lf %.6lf\n", sblock[i][SINDEX], smin[i].var1, smax[i].var1);
if(source[i].type==6) //double Sersic
{
if( sblock[i][STHETA2] )
fprintf(besto, "\ts_angle2 %d %.6lf %.6lf\n", sblock[i][STHETA2], smin[i].E2.theta * RTD, smax[i].E2.theta * RTD);
if( sblock[i][SA2] )
fprintf(besto, "\ts_sigx2 %d %.6lf %.6lf\n", sblock[i][SA2], smin[i].E2.a, smax[i].E2.a);
if( sblock[i][SEPS2] )
fprintf(besto, "\ts_eps2 %d %.6lf %.6lf\n", sblock[i][SEPS2], smin[i].eps2, smax[i].eps2);
if( sblock[i][SFLUX2] )
fprintf(besto, "\ts_mag2 %d %.6lf %.6lf\n", sblock[i][SFLUX2], smin[i].mag2, smax[i].mag2);
if( sblock[i][SINDEX2] )
fprintf(besto, "\tindex2 %d %.6lf %.6lf\n", sblock[i][SINDEX2], smin[i].var2, smax[i].var2);
}
fprintf(besto, "\tend\n");
}
......
src/o_rescale.c
View file @ 542d87fd
......@@ -273,6 +273,21 @@ int rescaleParam(int id, int type, int ipx, double *pval)
case(SFLUX):
val = prior(sblock[id][SFLUX], val, smin[id].mag, smax[id].mag);
break;
case(SA2):
val = prior(sblock[id][SA2], val, smin[id].E2.a, smax[id].E2.a);
break;
case(SEPS2):
val = prior(sblock[id][SEPS2], val, smin[id].eps2, smax[id].eps2);
break;
case(STHETA2):
val = prior(sblock[id][STHETA2], val, smin[id].E2.theta, smax[id].E2.theta);
break;
case(SINDEX2):
val = prior(sblock[id][SINDEX2], val, smin[id].var2, smax[id].var2);
break;
case(SFLUX2):
val = prior(sblock[id][SFLUX2], val, smin[id].mag2, smax[id].mag2);
break;
}
}
......@@ -545,7 +560,7 @@ int rescaleCube_1Atom(double *cube, int npar)
// Rescale the source parameters
if ( M.iclean == 2 )
for ( k = 0; k < S.ns; k++ )
for ( ipx = SCX; ipx <= SFLUX; ipx++ )
for ( ipx = SCX; ipx <= SFLUX2; ipx++ )
if ( sblock[k][ipx] != 0 )
valid *= rescaleParam(k, SOURCES, ipx, &cube[ipar++]);
......
src/o_set_lens_bayes.c
View file @ 542d87fd
......@@ -146,7 +146,7 @@ void o_set_limit_bayes(int nParam, long int nVal, double **array,
extern struct g_source S;
for ( i = 0; i < S.ns; i++ )
for ( ipx = SCX; ipx <= SFLUX; ipx++ )
for ( ipx = SCX; ipx <= SFLUX2; ipx++ )
if ( sblock[i][ipx] != 0 )
{
sblock[i][ipx] = prior;
......
src/o_set_lmax.c
View file @ 542d87fd
......@@ -102,6 +102,21 @@ void o_set_lmax(int i, int ipx, double x)
case(SINDEX):
smax[i].var1 = x;
break;
case(SA2):
smax[i].E2.a = x;
break;
case(SEPS2):
smax[i].eps2 = x;
break;
case(STHETA2):
smax[i].E2.theta = x;
break;
case(SFLUX2):
smax[i].mag2 = x;
break;
case(SINDEX2):
smax[i].var2 = x;
break;
case(VFCX):
vfmax.C.x = x;
break;
......
src/o_set_lmin.c
View file @ 542d87fd
......@@ -103,6 +103,21 @@ void o_set_lmin(int i, int ipx, double x)
case(SINDEX):
smin[i].var1 = x;
break;
case(SA2):
smin[i].E2.a = x;
break;
case(SEPS2):
smin[i].eps2 = x;
break;
case(STHETA2):
smin[i].E2.theta = x;
break;
case(SFLUX2):
smin[i].mag2 = x;
break;
case(SINDEX2):
smin[i].var2 = x;
break;
case(VFCX):
vfmin.C.x = x;
break;
......
src/o_set_source.c
View file @ 542d87fd
......@@ -38,6 +38,21 @@ void o_set_source(struct galaxie *source, int ipx, double val)
case(SFLUX):
source->mag = val;
break;
case(SA2):
source->E2.a = val;
break;
case(SEPS2):
source->eps2 = val;
break;
case(STHETA2):
source->E2.theta = val;
break;
case(SINDEX2):
source->var2 = val;
break;
case(SFLUX2):
source->mag2 = val;
break;
}
}
......@@ -68,6 +83,21 @@ double o_get_source(struct galaxie *source, int ipx)
case(SFLUX):
val = source->mag;
break;
case(SA2):
val = source->E2.a;
break;
case(SEPS2):
val = source->eps2;
break;
case(STHETA2):
val = source->E2.theta;
break;
case(SINDEX2):
val = source->var2;
break;
case(SFLUX2):
val = source->mag2;
break;
}
return val;
......
src/r_shapemodel.c
View file @ 542d87fd
......@@ -124,7 +124,8 @@ void r_shapemodel(FILE *IN,FILE *OUT, long int i)
fprintf(OUT,"\t%s\t%lf\n",second,sshape->E.theta);
fprintf(OUT,"\t%s\t%lf\n",second,sshape->E2.theta);
}
else if (!strcmp(second,"s_mag2") )
else if (!strcmp(second,"s_mag2") ||
!strcmp(second, "mag2"))
{
sscanf(third,"%lf",&sshape->mag2);
fprintf(OUT,"\t%s\t%lf\n",second,sshape->mag2);
......
src/readBayesModels.c
View file @ 542d87fd
......@@ -290,12 +290,14 @@ void setBayesModel( long int iVal, long int nVal, double **array)
if( pch != NULL )
source[i].C.x = source[i].C.y = 0.;
for( ipx = SCX ; ipx <= SFLUX; ipx++ )
for( ipx = SCX ; ipx <= SFLUX2; ipx++ )
if( sblock[i][ipx] != 0 )
o_set_source( &source[i], ipx, method(nVal, iParam++, array));
if((sblock[i][SEPS] != 0 )|| (sblock[i][SA] !=0))
source[i].E.b = source[i].E.a * (1. - source[i].eps) / (1. + source[i].eps);
if((sblock[i][SEPS2] != 0 )|| (sblock[i][SA2] !=0))
source[i].E2.b = source[i].E2.a * (1. - source[i].eps2) / (1. + source[i].eps2);
}
}
......@@ -843,11 +845,13 @@ static void convertArray(long int nVal, double **array )
// Sourcelimit parameters
if( M.iclean == 2 )
for ( i = 0; i < S.ns; i++ )
for( ipx = SCX; ipx <= SFLUX; ipx ++ )
for( ipx = SCX; ipx <= SFLUX2; ipx ++ )
if( sblock[i][ipx] != 0 )
{
if( ipx == STHETA )
array[nParam][iVal] = array[nParam][iVal] * DTR;
if( ipx == STHETA2 )
array[nParam][iVal] = array[nParam][iVal] * DTR;
nParam ++;
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment