Commit 6747bf30 authored by Miroslav Shaltev's avatar Miroslav Shaltev
Browse files

another FStatNomad cache writer fix, AfterMath compute some stats

parent f72097f6
......@@ -68,6 +68,8 @@ int main(int argc, char *argv[]) {
// compute distance for each point to each point
MAFTERMath->ComputeAMS();
// compare 2F for each point to each point
......@@ -166,3 +168,14 @@ AFTERMath::AFTERMath(int argc, char *argv[]) {
AFTERMath::~AFTERMath() {}
INT4 AFTERMath::ComputeAMS() {
for (int i=0; i < indata.size() - 1; i++) {
for (int j = i + 1; j < indata.size(); j++) {
REAL8 dst = compute_metric_distance(uvar->sType, metric, indata.at(i), indata.at(j));
REAL8 guess_stat = guess_2F(dst,gsl_vector_get(indata.at(i),DP_STAT));
REAL8 mess_stat = gsl_vector_get(indata.at(j),DP_STAT);
printf("Distance point %d to point %d: %e, measured stat: %e, guess: %e, perc: %f\n",i,j,dst,mess_stat,guess_stat,mess_stat/guess_stat*100);
}
}
return 0;
}
......@@ -259,6 +259,8 @@ public:
return indata.at(v);
};
INT4 ComputeAMS();
private:
CHAR *VCSInfoString;
BOOLEAN var_BinarySearch;
......@@ -273,6 +275,7 @@ private:
std::vector<gsl_vector*> indata;
INT4 compute_metric(gsl_matrix *g, LALSegList *seglist, INT4 SType, ExPoint *p, INT4 EType, INT4 MAType);
REAL8 compute_metric_distance(INT4 SType, gsl_matrix *g,gsl_vector *vs,gsl_vector *vc);
bool BinarySearch(){
if (uvar->sType > STYPE_BNS_ALPHA_DELTA_FREQ_ASINI_TASC_OMEGA_KAPPA_ETA) {
......@@ -283,7 +286,84 @@ private:
}
}
REAL8 n3x_equ(REAL8 Alpha, REAL8 Delta) {
return cos(Alpha) * cos(Delta);
};
REAL8 n3y_equ(REAL8 Alpha, REAL8 Delta) {
return sin(Alpha) * cos(Delta);
};
REAL8 n3z_equ(REAL8 Alpha, REAL8 Delta) {
return sin(Delta);
};
REAL8 bnsstdecc(REAL8 k, REAL8 argp) {
REAL8 cargp = cos(argp);
if ( cargp > 0 || cargp < 0 ) {
return k/cargp;
}
else {
return 0;
}
}
REAL8 bnsstdargp(REAL8 eta, REAL8 k) {
if (k > 0 or k < 0) {
return atan(eta/k);
}
else {
return 0;
}
}
REAL8 bnsstdtp(REAL8 tasc,REAL8 omega, REAL8 argp) {
if ( omega > 0 or omega < 0) {
return tasc + argp/omega;
}
else {
return 0;
}
}
REAL8 bnsstdperiod(REAL8 omega) {
if ( omega > 0 or omega < 0 ) {
return 2*LAL_PI/omega;
}
else {
return 0;
}
}
REAL8 bnsleltasc(REAL8 tp, REAL8 omega, REAL8 argp) {
if ( omega > 0 || omega < 0) {
return tp - argp / omega ;
}
else {
return 0;
}
};
REAL8 bnslelomega(REAL8 P) {
if ( P > 0 ) {
return LAL_TWOPI / P;
}
else {
return 0;
}
};
REAL8 bnskappa(REAL8 e, REAL8 w) {
return e * cos(w);
}
REAL8 bnseta(REAL8 e, REAL8 w) {
return e * sin(w);
}
REAL8 guess_2F(REAL8 mu, REAL8 twoF) {
return twoF - mu * ( twoF - 4 );
}
};
#endif /* Double-include protection. */
......@@ -317,5 +317,219 @@ INT4 AFTERMath::ComputeMetric(gsl_vector *v){
p.period = gsl_vector_get(v,DP_PERIOD);
p.stat = gsl_vector_get(v,DP_STAT);
return compute_metric(metric, segList, uvar->sType, &p, USETYPE_EXTENDS, uvar->metricType);
return compute_metric(metric, segList, uvar->sType, &p, 0, uvar->metricType);
}
REAL8 AFTERMath::compute_metric_distance(INT4 SType, gsl_matrix *g,gsl_vector *vs,gsl_vector *vc) {
REAL8 sumg = 0;
gsl_vector *c;
gsl_vector *s;
switch (SType) {
case STYPE_ISO_ALPHA_DELTA:
s = gsl_vector_calloc(3);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,n3x_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,1,n3y_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,2,n3z_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
c = gsl_vector_calloc(3);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,n3x_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,1,n3y_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,2,n3z_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
break;
case STYPE_ISO_ALPHA_DELTA_FREQ_F1DOT:
s = gsl_vector_calloc(5);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,n3x_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,1,n3y_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,2,n3z_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,3,gsl_vector_get(vs,DP_FREQ));
gsl_vector_set(s,4,gsl_vector_get(vs,DP_F1DOT));
c = gsl_vector_calloc(5);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,n3x_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,1,n3y_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,2,n3z_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,3,gsl_vector_get(vc,DP_FREQ));
gsl_vector_set(c,4,gsl_vector_get(vc,DP_F1DOT));
break;
case STYPE_ISO_ALPHA_DELTA_FREQ_F1DOT_F2DOT:
s = gsl_vector_calloc(6);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,n3x_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,1,n3y_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,2,n3z_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,3,gsl_vector_get(vs,DP_FREQ));
gsl_vector_set(s,4,gsl_vector_get(vs,DP_F1DOT));
gsl_vector_set(s,5,gsl_vector_get(vs,DP_F2DOT));
c = gsl_vector_calloc(6);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,n3x_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,1,n3y_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,2,n3z_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,3,gsl_vector_get(vc,DP_FREQ));
gsl_vector_set(c,4,gsl_vector_get(vc,DP_F1DOT));
gsl_vector_set(c,5,gsl_vector_get(vc,DP_F2DOT));
break;
case STYPE_ISO_ALPHA_DELTA_FREQ_F1DOT_F2DOT_F3DOT:
s = gsl_vector_calloc(7);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,n3x_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,1,n3y_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,2,n3z_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,3,gsl_vector_get(vs,DP_FREQ));
gsl_vector_set(s,4,gsl_vector_get(vs,DP_F1DOT));
gsl_vector_set(s,5,gsl_vector_get(vs,DP_F2DOT));
gsl_vector_set(s,6,gsl_vector_get(vs,DP_F3DOT));
c = gsl_vector_calloc(7);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,n3x_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,1,n3y_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,2,n3z_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,3,gsl_vector_get(vc,DP_FREQ));
gsl_vector_set(c,4,gsl_vector_get(vc,DP_F1DOT));
gsl_vector_set(c,5,gsl_vector_get(vc,DP_F2DOT));
gsl_vector_set(c,6,gsl_vector_get(vc,DP_F3DOT));
break;
case STYPE_ISO_FREQ_F1DOT:
s = gsl_vector_calloc(2);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,gsl_vector_get(vs,DP_FREQ));
gsl_vector_set(s,1,gsl_vector_get(vs,DP_F1DOT));
c = gsl_vector_calloc(2);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,gsl_vector_get(vc,DP_FREQ));
gsl_vector_set(c,1,gsl_vector_get(vc,DP_F1DOT));
break;
case STYPE_ISO_FREQ_F1DOT_F2DOT:
s = gsl_vector_calloc(3);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,gsl_vector_get(vs,DP_FREQ));
gsl_vector_set(s,1,gsl_vector_get(vs,DP_F1DOT));
gsl_vector_set(s,2,gsl_vector_get(vs,DP_F2DOT));
c = gsl_vector_calloc(3);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,gsl_vector_get(vc,DP_FREQ));
gsl_vector_set(c,1,gsl_vector_get(vc,DP_F1DOT));
gsl_vector_set(c,2,gsl_vector_get(vc,DP_F2DOT));
break;
case STYPE_ISO_FREQ_F1DOT_F2DOT_F3DOT:
s = gsl_vector_calloc(4);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,gsl_vector_get(vs,DP_FREQ));
gsl_vector_set(s,1,gsl_vector_get(vs,DP_F1DOT));
gsl_vector_set(s,2,gsl_vector_get(vs,DP_F2DOT));
gsl_vector_set(s,3,gsl_vector_get(vs,DP_F3DOT));
c = gsl_vector_calloc(4);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,gsl_vector_get(vc,DP_FREQ));
gsl_vector_set(c,1,gsl_vector_get(vc,DP_F1DOT));
gsl_vector_set(c,2,gsl_vector_get(vc,DP_F2DOT));
gsl_vector_set(c,3,gsl_vector_get(vc,DP_F3DOT));
break;
case STYPE_BNS_ALPHA_DELTA_FREQ_ASINI_TASC_OMEGA_KAPPA_ETA:
s = gsl_vector_calloc(9);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,n3x_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,1,n3y_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,2,n3z_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,0,gsl_vector_get(vs,DP_FREQ));
gsl_vector_set(s,4,gsl_vector_get(vs,DP_ASINI));
gsl_vector_set(s,5,bnsleltasc(gsl_vector_get(vs,DP_TPSSB),bnslelomega(gsl_vector_get(vs,DP_PERIOD)),gsl_vector_get(vs,DP_ARGP)));
gsl_vector_set(s,6,bnslelomega(gsl_vector_get(vs,DP_PERIOD)));
gsl_vector_set(s,7,bnskappa(gsl_vector_get(vs,DP_ECC),gsl_vector_get(vs,DP_ARGP)));
gsl_vector_set(s,8,bnseta(gsl_vector_get(vs,DP_ECC),gsl_vector_get(vs,DP_ARGP)));
c = gsl_vector_calloc(9);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,n3x_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,1,n3y_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,2,n3z_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,0,gsl_vector_get(vc,DP_FREQ));
gsl_vector_set(c,4,gsl_vector_get(vc,DP_ASINI));
gsl_vector_set(c,5,bnsleltasc(gsl_vector_get(vc,DP_TPSSB),bnslelomega(gsl_vector_get(vc,DP_PERIOD)),gsl_vector_get(vc,DP_ARGP)));
gsl_vector_set(c,6,bnslelomega(gsl_vector_get(vc,DP_PERIOD)));
gsl_vector_set(c,7,bnskappa(gsl_vector_get(vc,DP_ECC),gsl_vector_get(vc,DP_ARGP)));
gsl_vector_set(c,8,bnseta(gsl_vector_get(vc,DP_ECC),gsl_vector_get(vc,DP_ARGP)));
break;
case STYPE_BNS_ALPHA_DELTA_FREQ_F1DOT_ASINI_TASC_OMEGA_KAPPA_ETA:
s = gsl_vector_calloc(10);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,n3x_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,1,n3y_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,2,n3z_equ(gsl_vector_get(vs,DP_ALPHA),gsl_vector_get(vs,DP_DELTA)));
gsl_vector_set(s,0,gsl_vector_get(vs,DP_FREQ));
gsl_vector_set(s,1,gsl_vector_get(vs,DP_F1DOT));
gsl_vector_set(s,5,gsl_vector_get(vs,DP_ASINI));
gsl_vector_set(s,6,bnsleltasc(gsl_vector_get(vs,DP_TPSSB),bnslelomega(gsl_vector_get(vs,DP_PERIOD)),gsl_vector_get(vs,DP_ARGP)));
gsl_vector_set(s,7,bnslelomega(gsl_vector_get(vs,DP_PERIOD)));
gsl_vector_set(s,8,bnskappa(gsl_vector_get(vs,DP_ECC),gsl_vector_get(vs,DP_ARGP)));
gsl_vector_set(s,9,bnseta(gsl_vector_get(vs,DP_ECC),gsl_vector_get(vs,DP_ARGP)));
c = gsl_vector_calloc(10);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,n3x_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,1,n3y_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,2,n3z_equ(gsl_vector_get(vc,DP_ALPHA),gsl_vector_get(vc,DP_DELTA)));
gsl_vector_set(c,0,gsl_vector_get(vc,DP_FREQ));
gsl_vector_set(c,1,gsl_vector_get(vc,DP_F1DOT));
gsl_vector_set(c,5,gsl_vector_get(vc,DP_ASINI));
gsl_vector_set(c,6,bnsleltasc(gsl_vector_get(vc,DP_TPSSB),bnslelomega(gsl_vector_get(vc,DP_PERIOD)),gsl_vector_get(vc,DP_ARGP)));
gsl_vector_set(c,7,bnslelomega(gsl_vector_get(vc,DP_PERIOD)));
gsl_vector_set(c,8,bnskappa(gsl_vector_get(vc,DP_ECC),gsl_vector_get(vc,DP_ARGP)));
gsl_vector_set(c,9,bnseta(gsl_vector_get(vc,DP_ECC),gsl_vector_get(vc,DP_ARGP)));
break;
case STYPE_BNS_FREQ_ASINI_TASC_OMEGA_KAPPA_ETA:
s = gsl_vector_calloc(6);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,gsl_vector_get(vs,DP_FREQ));
gsl_vector_set(s,1,gsl_vector_get(vs,DP_ASINI));
gsl_vector_set(s,2,bnsleltasc(gsl_vector_get(vs,DP_TPSSB),bnslelomega(gsl_vector_get(vs,DP_PERIOD)),gsl_vector_get(vs,DP_ARGP)));
gsl_vector_set(s,3,bnslelomega(gsl_vector_get(vs,DP_PERIOD)));
gsl_vector_set(s,4,bnskappa(gsl_vector_get(vs,DP_ECC),gsl_vector_get(vs,DP_ARGP)));
gsl_vector_set(s,5,bnseta(gsl_vector_get(vs,DP_ECC),gsl_vector_get(vs,DP_ARGP)));
c = gsl_vector_calloc(6);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,gsl_vector_get(vc,DP_FREQ));
gsl_vector_set(c,1,gsl_vector_get(vc,DP_ASINI));
gsl_vector_set(c,2,bnsleltasc(gsl_vector_get(vc,DP_TPSSB),bnslelomega(gsl_vector_get(vc,DP_PERIOD)),gsl_vector_get(vc,DP_ARGP)));
gsl_vector_set(c,3,bnslelomega(gsl_vector_get(vc,DP_PERIOD)));
gsl_vector_set(c,4,bnskappa(gsl_vector_get(vc,DP_ECC),gsl_vector_get(vc,DP_ARGP)));
gsl_vector_set(c,5,bnseta(gsl_vector_get(vc,DP_ECC),gsl_vector_get(vc,DP_ARGP)));
break;
case STYPE_BNS_FREQ_F1DOT_ASINI_TASC_OMEGA_KAPPA_ETA:
s = gsl_vector_calloc(7);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,gsl_vector_get(vs,DP_FREQ));
gsl_vector_set(s,1,gsl_vector_get(vs,DP_F1DOT));
gsl_vector_set(s,2,gsl_vector_get(vs,DP_ASINI));
gsl_vector_set(s,3,bnsleltasc(gsl_vector_get(vs,DP_TPSSB),bnslelomega(gsl_vector_get(vs,DP_PERIOD)),gsl_vector_get(vs,DP_ARGP)));
gsl_vector_set(s,4,bnslelomega(gsl_vector_get(vs,DP_PERIOD)));
gsl_vector_set(s,5,bnskappa(gsl_vector_get(vs,DP_ECC),gsl_vector_get(vs,DP_ARGP)));
gsl_vector_set(s,6,bnseta(gsl_vector_get(vs,DP_ECC),gsl_vector_get(vs,DP_ARGP)));
c = gsl_vector_calloc(7);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,gsl_vector_get(vc,DP_FREQ));
gsl_vector_set(c,1,gsl_vector_get(vc,DP_F1DOT));
gsl_vector_set(c,2,gsl_vector_get(vc,DP_ASINI));
gsl_vector_set(c,3,bnsleltasc(gsl_vector_get(vc,DP_TPSSB),bnslelomega(gsl_vector_get(vc,DP_PERIOD)),gsl_vector_get(vc,DP_ARGP)));
gsl_vector_set(c,4,bnslelomega(gsl_vector_get(vc,DP_PERIOD)));
gsl_vector_set(c,5,bnskappa(gsl_vector_get(vc,DP_ECC),gsl_vector_get(vc,DP_ARGP)));
gsl_vector_set(c,6,bnseta(gsl_vector_get(vc,DP_ECC),gsl_vector_get(vc,DP_ARGP)));
break;
default:
LogPrintf(LOG_CRITICAL,"Unknown STYPE: %d! Abort!\n",SType);
exit(-1);
break;
}
for (INT4 i=0; i < g->size1; i++) {
for (INT4 j=0; j < g->size2; j++) {
sumg += gsl_matrix_get(g,i,j)*(gsl_vector_get(c,i)-gsl_vector_get(s,i))*(gsl_vector_get(c,j)-gsl_vector_get(s,j));
}
}
gsl_vector_free(s);
gsl_vector_free(c);
return sumg;
}
......@@ -484,7 +484,7 @@ INT4 FCSearch::write_cache(NOMAD::Cache *cache, INT4 cnt, INT4 type) {
LogPrintf(LOG_NORMAL,"Write cache file: %s\n",CacheFile(cnt,type));
}
for (INT4 i = 0; i < eCache.size(); i++) {
fprintf(fpCACHE,"%.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n",eCache.at(i).alpha,eCache.at(i).delta,eCache.at(i).freq,eCache.at(i).f1dot,eCache.at(i).f2dot,eCache.at(i).f3dot,eCache.at(i).TpSSB,eCache.at(i).argp,eCache.at(i).asini,eCache.at(i).ecc,eCache.at(i).period,eCache.at(i).stat);
fprintf(fpCACHE,"%.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n",eCache.at(i).alpha,eCache.at(i).delta,eCache.at(i).freq,eCache.at(i).f1dot,eCache.at(i).f2dot,eCache.at(i).f3dot,XLALGPSGetREAL8(&eCache.at(i).TpSSB),eCache.at(i).argp,eCache.at(i).asini,eCache.at(i).ecc,eCache.at(i).period,eCache.at(i).stat);
}
fclose(fpCACHE);
}
......
......@@ -3198,6 +3198,13 @@ void FStatNomad::CacheIn(NOMAD::Cache *cache) {
n++;
}
}
else {
p.asini = 0;
p.argp = 0;
p.ecc = 0;
XLALGPSSetREAL8(&p.TpSSB,0);
p.period = 0;
}
p.stat = -x->get_f().value();
eCache.push_back(p);
}
......
Supports Markdown
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