Commit 6d0e9ee6 authored by John Douglas Veitch's avatar John Douglas Veitch
Browse files

Merge branch 'printinjection_hdf5' into 'master'

Store injection sample in HDF5 file

See merge request !193
parents b657ecbb 944ddd77
Pipeline #25972 passed with stages
in 155 minutes and 32 seconds
......@@ -1414,6 +1414,13 @@ void LALInferenceWriteMCMCSamples(LALInferenceRunState *runState) {
}
/* TODO: Write metadata */
LALInferenceVariables *injParams = NULL;
if ( (injParams=LALInferencePrintInjectionSample(runState)) )
{
LALInferenceH5VariablesArrayToDataset(group, &injParams, 1, "injection_params");
LALInferenceClearVariables(injParams);
XLALFree(injParams);
}
}
XLALH5FileClose(group);
......
......@@ -98,6 +98,7 @@ src/lalinference_bench
src/lalinference_burst
src/lalinference_datadump
src/lalinference_nest
src/lalinference_injectedlike
src/lalinference_version
src/stamp-h1
src/stamp-h2
......
/*
* LALInferenceInjectedLike.c: Util bin to create an *injection file with the true parameters and the injected logL/P
*
* Copyright (C) 2018 salvatore vitale
*
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with with program; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
* MA 02111-1307 USA
*/
#include <stdio.h>
#include <lal/Date.h>
#include <lal/GenerateInspiral.h>
#include <lal/LALInference.h>
#include <lal/FrequencySeries.h>
#include <lal/Units.h>
#include <lal/StringInput.h>
#include <lal/LIGOLwXMLInspiralRead.h>
#include <lal/TimeSeries.h>
#include <lal/LALInferencePrior.h>
#include <lal/LALInferenceReadData.h>
#include <lal/LALInferenceLikelihood.h>
#include <lal/LALInferenceTemplate.h>
#include <lal/LALInferenceProposal.h>
#include <lal/LALInferenceInit.h>
#include <lal/LALInferenceCalibrationErrors.h>
#include <lal/LALInferenceVCSInfo.h>
/*************** MAIN **********************/
int main(int argc, char *argv[]){
int helpflag=0;
char help[]="\
LALInferenceInjectedLike:\n\
Print the injected values of parameters and (delta)LogL/P \n";
LALInferenceRunState *state;
ProcessParamsTable *procParams=NULL;
LALInferenceIFOData *data = NULL;
/* Read command line and parse */
procParams=LALInferenceParseCommandLine(argc,argv);
if(LALInferenceGetProcParamVal(procParams,"--help"))
{
helpflag=1;
fprintf(stdout,"%s",help);
}
/* write down git information */
fprintf(stdout,"\n\nLALInference version:%s,%s,%s,%s,%s\n\n", lalInferenceVCSInfo.vcsId,lalInferenceVCSInfo.vcsDate,lalInferenceVCSInfo.vcsBranch,lalInferenceVCSInfo.vcsAuthor,lalInferenceVCSInfo.vcsStatus);
/* initialise runstate based on command line */
/* This includes reading in the data */
/* And performing any injections specified */
/* And allocating memory */
state = LALInferenceInitRunState(procParams);
/* Create header */
if (state!=NULL && !helpflag){
ProcessParamsTable *ppt=NULL;
ppt=LALInferenceGetProcParamVal(state->commandLine,"--outfile");
if(!ppt){
ppt=LALInferenceGetProcParamVal(state->commandLine,"--outhdf");
if(!ppt){
fprintf(stderr,"Must specify --outfile <filename.dat> or --outhdf <filename.h5>\n");
exit(1);
}
}
}
if (state == NULL) {
if (!helpflag) {
fprintf(stderr, "run state not allocated (%s, line %d).\n",
__FILE__, __LINE__);
}
} else {
data = state->data;
}
/* Perform injections if data successful read or created */
if (state&&!helpflag){
LALInferenceInjectInspiralSignal(data, state->commandLine);
}
/* Simulate calibration errors.
* NOTE: this must be called after both ReadData and (if relevant)
* injectInspiralTD/FD are called! */
LALInferenceApplyCalibrationErrors(data, procParams);
if (!helpflag && LALInferenceGetProcParamVal(state->commandLine, "--roqtime_steps")){
LALInferenceSetupROQdata(state->data, state->commandLine);
fprintf(stderr, "done LALInferenceSetupROQdata\n");
}
/* Set up the threads */
LALInferenceInitCBCThreads(state,1);
/* Init the prior */
LALInferenceInitCBCPrior(state);
/* Choose the likelihood and set some auxiliary variables */
LALInferenceInitLikelihood(state);
/* Exit since we printed all command line arguments */
if(state == NULL || LALInferenceGetProcParamVal(state->commandLine,"--help"))
{
exit(0);
}
/* write injection with noise evidence information from algorithm */
LALInferencePrintInjectionSample(state);
/* end */
return(0);
}
......@@ -485,7 +485,6 @@ static REAL8 LALInferenceFusedFreqDomainLogLikelihood(LALInferenceVariables *cur
INT4 SKY_FRAME=0;
if(LALInferenceCheckVariable(currentParams,"SKY_FRAME"))
SKY_FRAME=*(INT4 *)LALInferenceGetVariable(currentParams,"SKY_FRAME");
if(SKY_FRAME==0){
/* determine source's sky location & orientation parameters: */
ra = *(REAL8*) LALInferenceGetVariable(currentParams, "rightascension"); /* radian */
......
......@@ -21,6 +21,7 @@
#include <lal/LALInferencePrior.h>
#include <lal/LALInferenceLikelihood.h>
#include <lal/LALInferenceProposal.h>
#include <lal/LALInferenceReadData.h>
#include <lal/LALInferenceHDF5.h>
#include <lal/LALInferencePriorVolumes.h>
......@@ -1209,6 +1210,13 @@ void LALInferenceNestedSamplingAlgorithm(LALInferenceRunState *runState)
XLALH5FileAddScalarAttribute(groupPtr, "cpu_time", &execution_time, LAL_D_TYPE_CODE);
}
LALInferenceVariables *injParams = NULL;
if ( (injParams=LALInferencePrintInjectionSample(runState)) )
{
LALInferenceH5VariablesArrayToDataset(groupPtr, &injParams, 1, "injection_params");
LALInferenceClearVariables(injParams);
XLALFree(injParams);
}
XLALH5FileClose(h5file);
}
......
......@@ -1805,6 +1805,7 @@ void InjectFD(LALInferenceIFOData *IFOdata, SimInspiralTable *inj_table, Process
char SNRpath[FILENAME_MAX+50];
ProcessParamsTable *ppt=NULL;
int flipped_masses=0;
LALInferenceIFOData *dataPtr;
ppt = LALInferenceGetProcParamVal(commandLine,"--outfile");
if (ppt)
......@@ -1895,6 +1896,8 @@ void InjectFD(LALInferenceIFOData *IFOdata, SimInspiralTable *inj_table, Process
*/
}
/* FIXME: One also need to code the same manipulations done to f_max in LALInferenceTemplateXLALSimInspiralChooseWaveform line 856 circa*/
/* Set up LAL dictionary */
LALDict* LALpars=XLALCreateDict();
......@@ -1931,7 +1934,17 @@ void InjectFD(LALInferenceIFOData *IFOdata, SimInspiralTable *inj_table, Process
if(LALInferenceGetProcParamVal(commandLine,"--inj-fref")) {
fref = atoi(LALInferenceGetProcParamVal(commandLine,"--inj-fref")->value);
}
if (approximant == TaylorF2)
f_max = 0.0; /* this will stop at ISCO */
else{
/* get the max f_max as done in Template. Has to do this since I cannot access LALInferenceModel here*/
dataPtr=IFOdata;
while(dataPtr){
if(dataPtr->fHigh>f_max) f_max=dataPtr->fHigh;
dataPtr=dataPtr->next;
}
}
/* Print a line with information about approximant, amp_order, phaseorder, tide order and spin order */
fprintf(stdout,"\n\n---\t\t ---\n");
fprintf(stdout,"Injection will run using Approximant %i (%s), phase order %i, amp order %i, spin order %i, tidal order %i, in the frequency domain.\n",approximant,XLALSimInspiralGetStringFromApproximant(approximant),phase_order,amp_order,(int) spinO,(int) tideO);
......@@ -1959,7 +1972,6 @@ void InjectFD(LALInferenceIFOData *IFOdata, SimInspiralTable *inj_table, Process
exit(1);
}
LALInferenceIFOData *dataPtr;
REAL8 Fplus, Fcross;
REAL8 plainTemplateReal, plainTemplateImag;
REAL8 templateReal, templateImag;
......@@ -2099,7 +2111,7 @@ static void PrintSNRsToFile(LALInferenceIFOData *IFOdata , char SNRpath[] ){
*/
void LALInferenceInjectionToVariables(SimInspiralTable *theEventTable, LALInferenceVariables *vars)
{
UINT4 spinCheck=0;
//UINT4 spinCheck=0;
if(!vars) {
XLALPrintError("Encountered NULL variables pointer");
XLAL_ERROR_VOID(XLAL_EINVAL);
......@@ -2108,46 +2120,13 @@ void LALInferenceInjectionToVariables(SimInspiralTable *theEventTable, LALInfere
REAL8 q = theEventTable->mass2 / theEventTable->mass1;
if (q > 1.0) q = 1.0/q;
REAL8 sx = theEventTable->spin1x;
REAL8 sy = theEventTable->spin1y;
REAL8 s1z = theEventTable->spin1z;
REAL8 a_spin1 = sqrt(sx*sx + sy*sy + s1z*s1z);
REAL8 theta_spin1, phi_spin1;
if (a_spin1 == 0.0) {
theta_spin1 = 0.0;
phi_spin1 = 0.0;
} else {
theta_spin1 = acos(s1z / a_spin1);
phi_spin1 = atan2(sy, sx);
if (phi_spin1 < 0.0) phi_spin1 += 2.0*M_PI;
}
sx = theEventTable->spin2x;
sy = theEventTable->spin2y;
REAL8 s2z = theEventTable->spin2z;
REAL8 a_spin2 = sqrt(sx*sx + sy*sy + s2z*s2z), theta_spin2, phi_spin2;
if (a_spin2 == 0.0) {
theta_spin2 = 0.0;
phi_spin2 = 0.0;
} else {
theta_spin2 = acos(s2z / a_spin2);
phi_spin2 = atan2(sy, sx);
if (phi_spin2 < 0.0) phi_spin2 += 2.0*M_PI;
}
/* Check for presence of spin in the injection */
if(a_spin1!=0.0 || a_spin2!=0.0) spinCheck=1;
REAL8 psi = theEventTable->polarization;
if (psi>=M_PI) psi -= M_PI;
REAL8 injGPSTime = XLALGPSGetREAL8(&(theEventTable->geocent_end_time));
REAL8 dist = theEventTable->distance;
REAL8 cosinclination = cos(theEventTable->inclination);
//REAL8 cosinclination = cos(theEventTable->inclination);
REAL8 phase = theEventTable->coa_phase;
REAL8 dec = theEventTable->latitude;
REAL8 ra = theEventTable->longitude;
......@@ -2158,38 +2137,60 @@ void LALInferenceInjectionToVariables(SimInspiralTable *theEventTable, LALInfere
REAL8 m1=theEventTable->mass1;
REAL8 m2=theEventTable->mass2;
REAL8 chirpmass = theEventTable->mchirp;
LALInferenceAddVariable(vars, "mass1", &m1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "mass2", &m2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "chirpmass", &chirpmass, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "q", &q, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "time", &injGPSTime, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "distance", &dist, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "costheta_jn", &cosinclination, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
if (LALInferenceCheckVariable(vars,"distance"))
LALInferenceAddVariable(vars, "distance", &dist, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
else if (LALInferenceCheckVariable(vars,"logdistance")){
REAL8 logdistance=log(dist);
LALInferenceAddVariable(vars, "logdistance", &logdistance, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
}
LALInferenceAddVariable(vars, "polarisation", &(psi), LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "phase", &phase, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
/* Those will work even if the user is working with the detector-frame variables because SKY_FRAME is set
to zero while calculating the injected logL in LALInferencePrintInjectionSample */
LALInferenceAddVariable(vars, "declination", &dec, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "rightascension", &ra, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "time", &injGPSTime, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "LAL_APPROXIMANT", &injapprox, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "LAL_PNORDER",&order,LALINFERENCE_INT4_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "LAL_AMPORDER", &(theEventTable->amp_order), LALINFERENCE_INT4_t, LALINFERENCE_PARAM_FIXED);
if(spinCheck){
if (theEventTable->spin1x==0 && theEventTable->spin1y==0.0 && theEventTable->spin2x==0.0 && theEventTable->spin2y==0.0){
LALInferenceAddVariable(vars, "a_spin1", &s1z, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "a_spin2", &s2z, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
}
else{
LALInferenceAddVariable(vars, "a_spin1", &a_spin1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "a_spin2", &a_spin2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "theta_spin1", &theta_spin1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "theta_spin2", &theta_spin2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "phi_spin1", &phi_spin1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "phi_spin2", &phi_spin2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
}
}
REAL8 thetaJN,phiJL,theta1,theta2,phi12,chi1,chi2;
/* Convert cartesian spin coordinates to system-frame variables*/
/* This fref is --inj-fref, which has been set previously by LALInferencePrintInjectionSample
I don't call it inj_fref since LALInferenceTemplate looks for fref, and that is what will be called to calculate
the logL at injval
*/
REAL8 fref=100.0;
if (LALInferenceCheckVariable(vars,"f_ref"))
fref= *(REAL8*) LALInferenceGetVariable(vars,"f_ref");
XLALSimInspiralTransformPrecessingWvf2PE(&thetaJN,&phiJL,&theta1,&theta2,&phi12,&chi1,&chi2,theEventTable->inclination,theEventTable->spin1x,theEventTable->spin1y,theEventTable->spin1z, theEventTable->spin2x, theEventTable->spin2y, theEventTable->spin2z,m1,m2,fref,phase);
if (LALInferenceCheckVariable(vars,"a_spin1"))
LALInferenceAddVariable(vars,"a_spin1", &chi1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
if (LALInferenceCheckVariable(vars,"a_spin2"))
LALInferenceAddVariable(vars,"a_spin2", &chi2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
if (LALInferenceCheckVariable(vars,"tilt_spin1"))
LALInferenceAddVariable(vars,"tilt_spin1", &theta1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
if (LALInferenceCheckVariable(vars,"tilt_spin2"))
LALInferenceAddVariable(vars,"tilt_spin2", &theta2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
if (LALInferenceCheckVariable(vars,"phi_jl"))
LALInferenceAddVariable(vars,"phi_jl", &phiJL, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
if (LALInferenceCheckVariable(vars,"phi12"))
LALInferenceAddVariable(vars,"phi12", &phi12, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
REAL8 costhetajn=cos(thetaJN);
LALInferenceAddVariable(vars, "costheta_jn", &costhetajn, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
}
void LALInferencePrintInjectionSample(LALInferenceRunState *runState) {
LALInferenceVariables *LALInferencePrintInjectionSample(LALInferenceRunState *runState) {
int errnum=0;
char *fname=NULL;
char defaultname[]="injection_params.dat";
......@@ -2208,7 +2209,7 @@ void LALInferencePrintInjectionSample(LALInferenceRunState *runState) {
ProcessParamsTable *ppt = LALInferenceGetProcParamVal(runState->commandLine,"--inj");
if (!ppt)
return;
return(NULL);
SimInspiralTableFromLIGOLw(&injTable, ppt->value, 0, 0);
......@@ -2239,8 +2240,25 @@ void LALInferencePrintInjectionSample(LALInferenceRunState *runState) {
if (!(approx && order)){
fprintf(stdout,"Unable to print injection sample: No approximant/PN order set\n");
return;
return(NULL);
}
REAL8 fref = 100.;
if(LALInferenceGetProcParamVal(runState->commandLine,"--inj-fref")) {
fref = atoi(LALInferenceGetProcParamVal(runState->commandLine,"--inj-fref")->value);
}
LALInferenceAddVariable(injparams,"f_ref",(void *)&fref,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
UINT4 azero=0;
LALInferenceAddVariable(injparams,"SKY_FRAME",(void *)&azero,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
/* remove eventual SKY FRAME vars since they will contain garbage*/
if (LALInferenceCheckVariable(injparams,"t0"))
LALInferenceRemoveVariable(injparams,"t0");
if (LALInferenceCheckVariable(injparams,"cosalpha"))
LALInferenceRemoveVariable(injparams,"cosalpha");
if (LALInferenceCheckVariable(injparams,"azimuth"))
LALInferenceRemoveVariable(injparams,"azimuth");
/* Fill named variables */
LALInferenceInjectionToVariables(theEventTable, injparams);
......@@ -2254,10 +2272,10 @@ void LALInferencePrintInjectionSample(LALInferenceRunState *runState) {
}
}
LALInferenceAddVariable(injparams, "logL", (void *)&injL,LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
if (LALInferenceCheckVariable(runState->algorithmParams, "logZnoise")){
REAL8 tmp=injL-*(REAL8 *)LALInferenceGetVariable(runState->algorithmParams,"logZnoise");
LALInferenceAddVariable(injparams,"deltalogL",(void *)&tmp,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
}
REAL8 logZnoise=LALInferenceNullLogLikelihood(runState->data);
REAL8 tmp2=injL-logZnoise;
LALInferenceAddVariable(injparams,"deltalogL",(void *)&tmp2,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
LALInferenceIFOData *data=runState->data;
while(data) {
char tmpName[320];
......@@ -2266,6 +2284,7 @@ void LALInferencePrintInjectionSample(LALInferenceRunState *runState) {
LALInferenceAddVariable(injparams, tmpName, &tmp, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
data=data->next;
}
/* Save to file */
outfile=fopen(fname,"w");
if(!outfile) {fprintf(stderr,"ERROR: Unable to open file %s for injection saving\n",fname); exit(1);}
......@@ -2275,9 +2294,8 @@ void LALInferencePrintInjectionSample(LALInferenceRunState *runState) {
LALInferencePrintSample(outfile, injparams);
fclose(outfile);
LALInferenceClearVariables(injparams);
XLALFree(injparams);
return;
//LALInferenceClearVariables(injparams);
return(injparams);
}
int enforce_m1_larger_m2(SimInspiralTable* injEvent){
......
......@@ -66,7 +66,7 @@ void LALInferenceInjectionToVariables(SimInspiralTable *theEventTable, LALInfere
* \brief Function to output a sample with logL values etc for the injection, if one is made.
* Requires --inj, --outfile and optionally --event (if not 0).
*/
void LALInferencePrintInjectionSample(LALInferenceRunState *runState);
LALInferenceVariables *LALInferencePrintInjectionSample(LALInferenceRunState *runState);
void LALInferenceInjectFromMDC(ProcessParamsTable *commandLine, LALInferenceIFOData *IFOdata);
/*@}*/
......
......@@ -91,6 +91,7 @@ LDADD = liblalinference.la
bin_PROGRAMS = \
lalinference_nest \
lalinference_injectedlike \
lalinference_burst \
lalinference_datadump \
lalinference_bench \
......@@ -98,6 +99,7 @@ bin_PROGRAMS = \
$(END_OF_LIST)
lalinference_nest_SOURCES = LALInferenceNest.c
lalinference_injectedlike_SOURCES = LALInferenceInjectedLike.c
lalinference_burst_SOURCES = LALInferenceBurst.c
lalinference_datadump_SOURCES = LALInferenceDataDump.c
lalinference_bench_SOURCES = LALInferenceBench.c
......
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