Commit 811a814e authored by John Douglas Veitch's avatar John Douglas Veitch

Merge branch 'master' into dist_marg_bayestar

parents 21991f58 f366366f
......@@ -609,6 +609,283 @@ int XLALAdaptiveRungeKutta4NoInterpolate(LALAdaptiveRungeKuttaIntegrator * integ
return outputlength;
}
int XLALAdaptiveRungeKuttaDenseandSparseOutput(LALAdaptiveRungeKuttaIntegrator * integrator,
void * params, REAL8 * yinit, REAL8 tinit, REAL8 tend, REAL8 deltat,
REAL8Array ** sparse_output, REAL8Array ** dense_output)
{
/* Error-checking variables used throughout */
int errnum = 0;
int status;
/* Integration and interpolation variables */
size_t dim = integrator->sys->dimension;
UINT4 sparse_outputlength = 0;
UINT4 dense_outputlength = 0;
size_t sparse_bufferlength, dense_bufferlength, retries;
REAL8 t = tinit;
REAL8 tnew;
REAL8 h0 = deltat;
REAL8Array *sparse_buffers = NULL;
REAL8Array *dense_buffers = NULL;
REAL8 *temp = NULL, *y, *y0, *dydt_in, *dydt_in0, *dydt_out, *yerr;
REAL8 interp_t = tinit + h0;
/* For speed, this replaces the single CALLGSL wrapper applied before each GSL call */
XLAL_BEGINGSL;
/* Allocate buffers!
* Note: REAL8Array has a field dimLength (UINT4Vector) with dimensions, and a field data that points to a single memory block;
* dimLength itself has fields length and data */
sparse_bufferlength = (int)((tend - tinit) / h0) + 2; /* allow for the initial value and possibly a final semi-step */
dense_bufferlength = (int)((tend - tinit) / h0) + 2;
const UINT4 dimn = dim + 1;/* Time is not included in input dimesions, but is included in ouput arrays. */
sparse_buffers = XLALCreateREAL8ArrayL(2, dimn, sparse_bufferlength);
dense_buffers = XLALCreateREAL8ArrayL(2, dimn, dense_bufferlength);
temp = LALCalloc(6 * dim, sizeof(REAL8));
if (!sparse_buffers || !dense_buffers || !temp) {
errnum = XLAL_ENOMEM;
goto bail_out;
}
/* Aliases */
y = temp;
y0 = temp + dim;
dydt_in = temp + 2 * dim;
dydt_in0 = temp + 3 * dim;
dydt_out = temp + 4 * dim;
yerr = temp + 5 * dim;
/* Integrator set up */
integrator->sys->params = params;
integrator->returncode = 0;
retries = integrator->retries;
memcpy(y, yinit, dim * sizeof(REAL8));
/* Store the first data point. */
sparse_buffers->data[0] = t;
dense_buffers->data[0] = t;
for (UINT4 i = 1; i <= dim; i++){
UINT4 iminus1 = i - 1;
sparse_buffers->data[i * sparse_bufferlength] = y[iminus1];
dense_buffers->data[i * dense_bufferlength] = y[iminus1];
}
/* Compute derivatives at the initial time (dydt_in); bail out if impossible. */
if ((status = integrator->dydt(t, y, dydt_in, params)) != GSL_SUCCESS) {
integrator->returncode = status;
errnum = XLAL_EFAILED;
goto bail_out;
}
while (1) {
if (!integrator->stopontestonly && t >= tend) {
break;
}
if (integrator->stop) {
if ((status = integrator->stop(t, y, dydt_in, params)) != GSL_SUCCESS) {
integrator->returncode = status;
break;
}
}
/* Try stepping! */
try_step:
/* If we would be stepping beyond the final time, stop there instead. */
if (!integrator->stopontestonly && t + h0 > tend)
h0 = tend - t;
/* Save y to y0 and dydt_in to dydt_in0. */
memcpy(y0, y, dim * sizeof(REAL8));
memcpy(dydt_in0, dydt_in, dim * sizeof(REAL8));
/* Call the GSL stepper function. */
status = gsl_odeiv_step_apply(integrator->step, t, h0, y, yerr, dydt_in, dydt_out, integrator->sys);
/* Note: If the user-supplied functions defined in the system dydt return a status other than GSL_SUCCESS,
* the step will be aborted. In this case, the elements of y will be restored to their pre-step values,
* and the error code from the user-supplied function will be returned. */
/* Did the stepper report a derivative-evaluation error? */
if (status != GSL_SUCCESS) {
if (retries--) {
/* If we have singularity retries left, reduce the timestep and try again... */
h0 = h0 / 10.0;
goto try_step;
} else {
integrator->returncode = status;
/* ...otherwise exit the loop. */
break;
}
} else {
/* We stepped successfully; reset the singularity retries. */
retries = integrator->retries;
}
tnew = t + h0;
/* Call the GSL error-checking function. */
status = gsl_odeiv_control_hadjust(integrator->control, integrator->step, y, yerr, dydt_out, &h0);
/* Did the error-checker reduce the stepsize?
* Note: other possible return codes are GSL_ODEIV_HADJ_INC if it was increased;
* GSL_ODEIV_HADJ_NIL if it was unchanged. */
if (status == GSL_ODEIV_HADJ_DEC) {
/* If so, undo the step, and try again */
memcpy(y, y0, dim * sizeof(REAL8));
memcpy(dydt_in, dydt_in0, dim * sizeof(REAL8));
goto try_step;
}
if (2*dense_outputlength >= dense_bufferlength) {
REAL8Array *rebuffers;
/* Sadly, we cannot use ResizeREAL8Array because it would only work if we extended the first array dimension;
* thus we copy everything and switch the buffers. */
if (!(rebuffers = XLALCreateREAL8ArrayL(2, dimn, 2 * dense_bufferlength))) {
errnum = XLAL_ENOMEM;
goto bail_out;
} else {
for (UINT4 k = 0; k <= dim; k++)
memcpy(&rebuffers->data[k * 2 * dense_bufferlength], &dense_buffers->data[k * dense_bufferlength], dense_outputlength * sizeof(REAL8));
XLALDestroyREAL8Array(dense_buffers);
dense_buffers = rebuffers;
dense_bufferlength *= 2;
}
}
{
const REAL8 interp_t_old = interp_t;
const UINT4 dense_outputlength_old = dense_outputlength;
while (interp_t < tnew) {
dense_outputlength++;
dense_buffers->data[dense_outputlength] = interp_t;
interp_t += deltat;
}
interp_t = interp_t_old;
dense_outputlength = dense_outputlength_old;
}
{
const REAL8 h = tnew - t;
const REAL8 h_inv = 1.0/h;
for (UINT4 i = 0; i < dim; i++) {
REAL8 interp_t_old = interp_t;
UINT4 dense_outputlength_old = dense_outputlength;
REAL8 y0i = y0[i];
REAL8 yi = y[i];
while (interp_t < tnew) {
dense_outputlength++;
const REAL8 theta = (interp_t - t)*h_inv;
dense_buffers->data[(i+1)*dense_bufferlength + dense_outputlength] =
(1.0 - theta)*y0i + theta*yi + theta*(theta-1.0)*( (1.0 - 2.0*theta)*(yi - y0i) + h*( (theta-1.0)*dydt_in[i] + theta*dydt_out[i]));
interp_t += deltat;
}
if(i<(dim-1)){
interp_t = interp_t_old;
dense_outputlength = dense_outputlength_old;
}
}
}
/* Update the current time and input derivatives. */
t = tnew;
memcpy(dydt_in, dydt_out, dim * sizeof(REAL8));
sparse_outputlength++;
/* Check if interpolation buffers need to be extended. */
if (sparse_outputlength >= sparse_bufferlength) {
REAL8Array *rebuffers;
/* Sadly, we cannot use ResizeREAL8Array because it would only work if we extended the first array dimension;
* thus we copy everything and switch the buffers. */
if (!(rebuffers = XLALCreateREAL8ArrayL(2, dimn, 2 * sparse_bufferlength))) {
errnum = XLAL_ENOMEM;
goto bail_out;
} else {
for (UINT4 i = 0; i <= dim; i++)
memcpy(&rebuffers->data[i * 2 * sparse_bufferlength], &sparse_buffers->data[i * sparse_bufferlength], sparse_outputlength * sizeof(REAL8));
XLALDestroyREAL8Array(sparse_buffers);
sparse_buffers = rebuffers;
sparse_bufferlength *= 2;
}
}
/* Copy time and state into buffers. */
sparse_buffers->data[sparse_outputlength] = t;
for (UINT4 i = 1; i <= dim; i++)
sparse_buffers->data[i * sparse_bufferlength + sparse_outputlength] = y[i - 1];
}
/* Copy the final state into yinit and dense_output. */
memcpy(yinit, y, dim * sizeof(REAL8));
dense_outputlength++;
dense_buffers->data[dense_outputlength] = t;
for (UINT4 i = 1; i <= dim; i++)
dense_buffers->data[i * dense_bufferlength + dense_outputlength] = y[i - 1];
if (sparse_outputlength == 0 || dense_outputlength == 0)
goto bail_out;
sparse_outputlength++;
(*sparse_output) = XLALCreateREAL8ArrayL(2, dim+1, sparse_outputlength);
(*dense_output) = XLALCreateREAL8ArrayL(2, dim+1, dense_outputlength);
if (!(*sparse_output) || !(*dense_output)) {
errnum = XLAL_ENOMEM; /* ouch again, ran out of memory */
if (*sparse_output){
XLALDestroyREAL8Array(*sparse_output);
sparse_outputlength = 0;
}
if (*dense_output){
XLALDestroyREAL8Array(*dense_output);
dense_outputlength = 0;
}
goto bail_out;
}
for(UINT4 j = 0; j < sparse_outputlength; j++) {
(*sparse_output)->data[j] = sparse_buffers->data[j];
for(UINT4 i = 1; i <= dim; i++) {
(*sparse_output)->data[i * sparse_outputlength + j] = sparse_buffers->data[i * sparse_bufferlength + j];
}
}
for(UINT4 j = 0; j < dense_outputlength; j++) {
(*dense_output)->data[j] = dense_buffers->data[j];
for(UINT4 i = 1; i <= dim; i++) {
(*dense_output)->data[i * dense_outputlength + j] = dense_buffers->data[i * dense_bufferlength + j];
}
}
/* Deallocate memory and return sparse_outputlength. */
bail_out:
XLAL_ENDGSL;
if (sparse_buffers)
XLALDestroyREAL8Array(sparse_buffers);
if (dense_buffers)
XLALDestroyREAL8Array(dense_buffers);
if (temp)
LALFree(temp);
if (errnum)
XLAL_ERROR(errnum);
return sparse_outputlength;
}
int XLALAdaptiveRungeKutta4(LALAdaptiveRungeKuttaIntegrator * integrator,
void *params, REAL8 * yinit, REAL8 tinit, REAL8 tend, REAL8 deltat, REAL8Array ** yout)
......
......@@ -103,6 +103,9 @@ int XLALAdaptiveRungeKutta4( LALAdaptiveRungeKuttaIntegrator *integrator,
int XLALAdaptiveRungeKutta4NoInterpolate(LALAdaptiveRungeKuttaIntegrator * integrator,
void * params, REAL8 * yinit, REAL8 tinit, REAL8 tend, REAL8 deltat_or_h0,
REAL8Array ** t_and_yout,INT4 EOBversion);
int XLALAdaptiveRungeKuttaDenseandSparseOutput(LALAdaptiveRungeKuttaIntegrator * integrator,
void * params, REAL8 * yinit, REAL8 tinit, REAL8 tend, REAL8 deltat,
REAL8Array ** sparse_output, REAL8Array ** dense_output);
/* END OPTIMIZED */
int XLALAdaptiveRungeKutta4Hermite( LALAdaptiveRungeKuttaIntegrator *integrator,
......
......@@ -348,7 +348,7 @@ int compareFstatCandidates ( const void *candA, const void *candB );
int compareFstatCandidates_BSGL ( const void *candA, const void *candB );
int WriteFstatLog ( const CHAR *log_fname, const CHAR *logstr );
CHAR *XLALGetLogString ( const ConfigVariables *cfg );
CHAR *XLALGetLogString ( const ConfigVariables *cfg, const BOOLEAN verbose );
int write_TimingInfo ( const CHAR *timingFile, const timingInfo_t *ti, const ConfigVariables *cfg );
......@@ -411,7 +411,8 @@ int main(int argc,char *argv[])
XLAL_CHECK_MAIN ( InitFstat( &GV, &uvar) == XLAL_SUCCESS, XLAL_EFUNC );
/* ----- produce a log-string describing the specific run setup ----- */
XLAL_CHECK_MAIN ( (GV.logstring = XLALGetLogString ( &GV )) != NULL, XLAL_EFUNC );
BOOLEAN logstrVerbose = TRUE;
XLAL_CHECK_MAIN ( (GV.logstring = XLALGetLogString ( &GV, logstrVerbose )) != NULL, XLAL_EFUNC );
LogPrintfVerbatim( LOG_NORMAL, "%s", GV.logstring );
/* keep a log-file recording all relevant parameters of this search-run */
......@@ -797,6 +798,8 @@ int main(int argc,char *argv[])
timing.tauFstat /= num_templates;
timing.tauTemplate /= num_templates;
timing.tauF0 = timing.tauFstat / timing.NSFTs;
timing.tauTransFstatMap /= num_templates;
timing.tauTransMarg /= num_templates;
XLAL_CHECK_MAIN ( write_TimingInfo ( uvar.outputTiming, &timing, &GV ) == XLAL_SUCCESS, XLAL_EFUNC );
......@@ -1685,7 +1688,7 @@ InitFstat ( ConfigVariables *cfg, const UserInput_t *uvar )
* Produce a log-string describing the present run-setup
*/
CHAR *
XLALGetLogString ( const ConfigVariables *cfg )
XLALGetLogString ( const ConfigVariables *cfg, const BOOLEAN verbose )
{
XLAL_CHECK_NULL ( cfg != NULL, XLAL_EINVAL );
......@@ -1701,6 +1704,11 @@ XLALGetLogString ( const ConfigVariables *cfg )
XLAL_CHECK_NULL ( (logstr = XLALStringAppend ( logstr, "\n" )) != NULL, XLAL_EFUNC );
XLAL_CHECK_NULL ( (logstr = XLALStringAppend ( logstr, cfg->VCSInfoString )) != NULL, XLAL_EFUNC );
if ( !verbose ) {
/* don't include search details */
return logstr;
}
XLAL_CHECK_NULL ( snprintf ( buf, BUFLEN, "%%%% FstatMethod used: '%s'\n", XLALGetFstatInputMethodName ( cfg->Fstat_in ) ) < BUFLEN, XLAL_EBADLEN );
XLAL_CHECK_NULL ( (logstr = XLALStringAppend ( logstr, buf )) != NULL, XLAL_EFUNC );
......@@ -2330,6 +2338,10 @@ write_TimingInfo ( const CHAR *fname, const timingInfo_t *ti, const ConfigVariab
{
XLAL_CHECK ( (fname != NULL) && (ti != NULL), XLAL_EINVAL );
CHAR *logstr_short;
BOOLEAN logstrVerbose = FALSE;
XLAL_CHECK ( (logstr_short = XLALGetLogString ( cfg, logstrVerbose )) != NULL, XLAL_EFUNC );
FILE *fp;
if ( (fp = fopen(fname,"rb" )) != NULL )
{
......@@ -2339,14 +2351,20 @@ write_TimingInfo ( const CHAR *fname, const timingInfo_t *ti, const ConfigVariab
else
{
XLAL_CHECK ( (fp = fopen( fname, "wb" ) ), XLAL_ESYS, "Failed to open new timing-file '%s' for writing\n", fname );
fprintf ( fp, "%2s%6s %10s %10s %10s %10s %10s\n",
"%%", "NSFTs", "NFreq", "tauF[s]", "tauFEff[s]", "tauF0[s]", "FstatMethod" );
fprintf ( fp, "%s", logstr_short );
fprintf ( fp, "%2s%6s %10s %10s %10s %10s %10s %10s %10s %10s %10s %16s %12s\n",
"%%", "NSFTs", "NFreq", "tauF[s]", "tauFEff[s]", "tauF0[s]", "FstatMethod",
"tauMin", "tauMax", "NStart", "NTau", "tauTransFstatMap", "tauTransMarg"
);
}
fprintf ( fp, "%8d %10d %10.1e %10.1e %10.1e %10s\n",
ti->NSFTs, ti->NFreq, ti->tauFstat, ti->tauTemplate, ti->tauF0, XLALGetFstatInputMethodName(cfg->Fstat_in) );
fprintf ( fp, "%8d %10d %10.1e %10.1e %10.1e %10s %10d %10d %10d %10d %16.1e %12.1e\n",
ti->NSFTs, ti->NFreq, ti->tauFstat, ti->tauTemplate, ti->tauF0, XLALGetFstatInputMethodName(cfg->Fstat_in),
ti->tauMin, ti->tauMax, ti->NStart, ti->NTau, ti->tauTransFstatMap, ti->tauTransMarg
);
fclose ( fp );
XLALFree ( logstr_short );
return XLAL_SUCCESS;
} /* write_TimingInfo() */
......
......@@ -32,6 +32,7 @@
#standard library imports
import sys
import os
import socket
from math import ceil,floor
import cPickle as pickle
......@@ -84,68 +85,20 @@ __author__="Ben Aylott <benjamin.aylott@ligo.org>, Ben Farr <bfarr@u.northwester
__version__= "git id %s"%git_version.id
__date__= git_version.date
from lalinference.lalinference_pipe_utils import guess_url
def email_notify(address,path):
import smtplib
import subprocess
import socket
import os
USER = os.environ('USER')
HOST = socket.getfqdn()
address=address.split(',')
SERVER="localhost"
USER=os.environ['USER']
HOST=socket.getfqdn()
FROM=USER+'@'+HOST
SUBJECT="LALInference result is ready at "+HOST+"!"
# Guess the web space path for the clusters
fslocation=os.path.abspath(path)
webpath='posplots.html'
if 'public_html' in fslocation:
k='public_html/'
elif 'WWW' in fslocation:
k='WWW/'
elif 'www_html' in fslocation:
k='www_html/'
else:
k=None
if k is not None:
(a,b)=fslocation.split(k)
webpath=os.path.join('~%s'%(USER),b,webpath)
onweb=True
else:
(c,d)=outpath.split(os.environ['USER'])
for k in ['public_html','WWW','www_html']:
trypath=c+os.environ['USER']+'/'+k+d
#Follow symlinks
if os.path.realpath(trypath)==os.path.normpath(outpath):
(a,b)=trypath.split(k)
webpath=os.path.join('~%s'%(USER),b,webpath)
onweb=True
break
else:
webpath=os.path.join(fslocation,'posplots.html')
onweb=False
if 'atlas' in HOST:
url="https://atlas1.atlas.aei.uni-hannover.de/"
elif 'cit' in HOST or 'caltech' in HOST:
url="https://ldas-jobs.ligo.caltech.edu/"
elif 'ligo-wa' in HOST:
url="https://ldas-jobs.ligo-wa.caltech.edu/"
elif 'ligo-la' in HOST:
url="https://ldas-jobs.ligo-la.caltech.edu/"
elif 'uwm' in HOST or 'nemo' in HOST:
url="https://ldas-jobs.phys.uwm.edu/"
elif 'phy.syr.edu' in HOST:
url="https://sugar-jobs.phy.syr.edu/"
elif 'arcca.cf.ac.uk' in HOST:
url="https://geo2.arcca.cf.ac.uk/"
elif 'vulcan' in HOST:
url="https://galahad.aei.mpg.de/"
else:
if onweb:
url="http://%s/"%(HOST)
else:
url=HOST+':'
url=url+webpath
url = guess_url(os.path.join(fslocation,webpath))
TEXT="Hi "+USER+",\nYou have a new parameter estimation result on "+HOST+".\nYou can view the result at "+url+"\n"
cmd='echo "%s" | mail -s "%s" "%s"'%(TEXT,SUBJECT,', '.join(address))
proc = subprocess.Popen(cmd, stdout=subprocess.PIPE,stderr=subprocess.STDOUT, shell=True)
......
......@@ -23,6 +23,62 @@ import numpy as np
# type of job. Each class has inputs and outputs, which are used to
# join together types of jobs into a DAG.
def guess_url(fslocation):
"""
Try to work out the web address of a given path
"""
SERVER="localhost"
USER=os.environ['USER']
HOST=socket.getfqdn()
if 'public_html' in fslocation:
k='public_html/'
elif 'WWW' in fslocation:
k='WWW/'
elif 'www_html' in fslocation:
k='www_html/'
else:
k=None
if k is not None:
(a,b)=fslocation.split(k)
webpath=os.path.join('~%s'%(USER),b)
onweb=True
else:
(c,d)=fslocation.split(USER)
for k in ['public_html','WWW','www_html']:
trypath=c+os.environ['USER']+'/'+k+d
#Follow symlinks
if os.path.realpath(trypath)==os.path.normpath(fslocation):
(a,b)=trypath.split(k)
webpath=os.path.join('~%s'%(USER),b,webpath)
onweb=True
break
else:
webpath=fslocation
onweb=False
if 'atlas' in HOST:
url="https://atlas1.atlas.aei.uni-hannover.de/"
elif 'cit' in HOST or 'caltech' in HOST:
url="https://ldas-jobs.ligo.caltech.edu/"
elif 'ligo-wa' in HOST:
url="https://ldas-jobs.ligo-wa.caltech.edu/"
elif 'ligo-la' in HOST:
url="https://ldas-jobs.ligo-la.caltech.edu/"
elif 'uwm' in HOST or 'nemo' in HOST:
url="https://ldas-jobs.phys.uwm.edu/"
elif 'phy.syr.edu' in HOST:
url="https://sugar-jobs.phy.syr.edu/"
elif 'arcca.cf.ac.uk' in HOST:
url="https://geo2.arcca.cf.ac.uk/"
elif 'vulcan' in HOST:
url="https://galahad.aei.mpg.de/"
else:
if onweb:
url="http://%s/"%(HOST)
else:
url=HOST+':'
url=url+webpath
return(url)
class Event():
"""
Represents a unique event to run on
......@@ -587,7 +643,7 @@ def Query_ROQ_Bounds_Type(path, roq_paths):
return roq_bounds
class LALInferencePipelineDAG(pipeline.CondorDAG):
def __init__(self,cp,dax=False,first_dag=True,previous_dag=None,site='local'):
def __init__(self,cp,dax=False,site='local'):
self.subfiles=[]
self.config=cp
self.engine=get_engine_name(cp)
......@@ -599,22 +655,15 @@ class LALInferencePipelineDAG(pipeline.CondorDAG):
self.basepath=os.getcwd()
print 'No basepath specified, using current directory: %s'%(self.basepath)
mkdirs(self.basepath)
print("Generating LALInference DAG in {0}".format(self.basepath))
if dax:
os.chdir(self.basepath)
self.posteriorpath=os.path.join(self.basepath,'posterior_samples')
mkdirs(self.posteriorpath)
if first_dag:
daglogdir=cp.get('paths','daglogdir')