aps_120426.tar

README.txt0000664000437200043720000005751111744047515012741 0ustar  danielsfdanielsfThis document will provide a guided tour of the code provided in this package. 

Section I will point out the requirements and dependencies of the code.

Section II will give a qualitative description of the source code files. 

Section III will describe aps_interface.cpp the user interface code for the 
APS alogrithm.

Section IV will provide a description of the subroutines in
likelihoodinterface.cpp, which is where most of the work implementing the
APS algorithm is done.  

Section V will describe the output produced by the code.

Readers only interested in how to interface their particular data and likelihood
function with this code should read Section III and the discussion of the
likelihood::call_likelihood(double *v) subroutine in Section IV.

You are welcome to use, modify, and distribute this code as you see fit.  If you
end up using it, please cite

B. Bryan, J. Schneider, C. J. Miller, C. Genovese, and L. Wasserman, Apj 665:225
(2007)

and

S.F. Daniel, A. Connolly, and J. Schneider [wherever we publish] (2012)

Throughout this document ``the text'' will refer to Daniel, Connolly, and
Schneider (2012).

////////////////////////

I. Required software

In order to run this software, you must have a version of the Fortran linear
algebra libraries BLAS (http://www.netlib.org/blas/faq.html) and LAPACK
(http://www.netlib.org/lapack/).

Two versions of the code are presented: one (current configuration) uses a 
cartoon likelihood function in two-dimensional parameter space with two islands 
of high likelihood; the other uses the WMAP 7 year likelihood function in 6 
dimensional cosmological parameter space.  If you are interested in running the 
CMB likelihood function, change the Makefile to uncomment the CMB compile 
command for likelihoodinterface.o and comment out the cartoon version.
Make a similar change to the subroutine call_likelihood(double*) in
likelihoodinterface.cpp (i.e. uncomment the WMAP 7 likelihood function and
comment out the cartoon likelihood).

The CMB version of the code interfaces with the January 2011 release of 
CAMB (http://camb.info/readme.html) and the WMAP 7 likeihood code
(http://lambda.gsfc.nasa.gov/product/map/current/likelihood_info.cfm) to
calculate constraints on cosmological parameters.  In order to run this code,
you will also need to download and compile a cfitsio library
(http://heasarc.nasa.gov/fitsio/fitsio.html).  This is actually a requirement of
WMAP7 and CAMB, not APS.

To compile the code, edit the Makefile so that LAPACK, BLAS, CAMB, CFITSIO,
WMAP, and WMAPLIBS point to the appropriate directories.  

Also alter gg and ff so that they refer to your native C++ and Fortran 
compilers.  

FORTRAN should refer to whatever libraries are necessary to make your C++
compiler talk to your Fortran compiler.  This is necessary if LAPACK and BLAS
are compiled in Fortran.  In our case, we compiled LAPACK and BLAS using the
gfotran compiler, so our FORTRAN flag had to be -lgfortran.

You can uncomment the openmp flags if you want to parallelize the code using
openmp.

To compile the CMB code, simply type "make aps".  This will compile an 
executable ./aps which can be called with a parameter input file (for an
example, see the included files params.sav and params_CMB.sav*) to run the APS
algorithm on the likelihood code and data that has been hard-coded into
likelihoodinterface.cpp (see Sections III and IV below).

*params.sav is written to run on the cartoon likelihood function. 
params_CMB.sav is written for the WMAP 7 likelihood function.

Running this code will produce two output files.  One, whose name is specified
in the parameter file (in this case, test_out.sav) will record all of the points
sampled by the algorithm, along with several statistics which will allow you to
check on the accuracy of the Gaussian process.  Another file, timingfile.sav 
will log the progress of the code, including statistics indicating how much 
time is being spent on the likelihood function and how much time is being 
spent on the Gaussian process (so that the user can set parameters such as the 
number of Gaussian process candidates and nearest neighbors to optimize 
resource usage).

After the algorithm has run for the specified number of iterations, it will
write the minimum and maximum acceptable (chi-squared <= chi_squared_lim) values
found for each parameter to stdout.

/////////////////////////

II. Overview of provided source code

goto_tools.cpp

This is a file containing some general purpose routines (most of them taken 
from Numerical Recipes).  Most notably, there is a polynomial interpolater and
wrapper as well as a merge sort routine.  goto_tools.h encodes the uniform
random number class Ran, also taken from Numerical Recipes.

eigen_wrapper.cpp

This contains a matrix inversion routine based on functions in the LAPACK
library.  It is used to invert the covariogram matrices in the Gaussian process.

kd.cpp

This contains code implementing a kd-tree.  This tree is used to store the set
of points in parameter space sampled by the code.  They are organized to
facilitate efficient nearest-neighbor searches in the Gaussian process.

likelihoodinterface.cpp

This file contains all of the subroutines unique to the APS algorithm. 
They exist within a class of object called 'likelihood'.  Section IV will
detail the use of each routine.

aps_interface.cpp

This is the user interface that ultimately initializes and calls likelihood.cpp.
It will read in the parameter file (e.g. params.sav) and use the information
therein to set the customizable parameters controlling how APS is run.

camb_wrapper_wmap.F90

This contains the code which calls CAMB to convert cosmological parameters into
a CMB anisotropy spectrum.

wmap_wrapper.F90

This contains the code that calls the WMAP 7 likelihood function on the output
from camb_wrapper_wmap.F90

////////////////////////

III. aps_interface.cpp

After typing

    make aps
    
you will have an executable that can be run in the form

    ./aps params.sav

params.sav is a list of keywords and parameter values.  There should be a space
separating the keyword and the corresponding value.  The keyword/parameter pairs
can appear in any order except that 'dim' (setting the number of parameters in
your parameter space) must come first.  The acceptable parameters are

#dim N -- set the dimensionality of your parameter space to N (again,
this must be the first line in params.sav
		
#param i name min max -- set the name of the ith parameter to 'name'; set the
minimum and maximum allowed values of the ith parameter to min and max.  Note:
parameters are zero-indexed, so i goes from 0 to N-1 for N-dimensional space

#Ns S -- start APS by randomly sampling S points uniformly

#Nk K -- after randomly sampling the points above, sample another K points
uniformly and use them to set the Kriging parameter

#Ng G -- use G nearest neighbors in the Gaussian process

#Nc D -- use D random candidates at each iteration of the Gaussian process

#end E -- end APS after it has sampled E points

#chitarget C -- set the target value of chi square for your desired confidence
limit to C

#chimin M -- set the expected minimum value of chi squared (used for the
gradient descent in the lingering algorithm) to M

#linger J -- turn on the lingering algorithm.  Each time lingering is
invoked, linger for J iterations

#outputfile filename.txt -- output the points sampled to filename.txt.  Note, if
filename.txt already exists, the code will just append its output to the end of
the file.

#resumefile oldfilename.txt -- read in oldfilename.txt (which is an old iteration of
the output above) and resume running APS (writing the output to the file
specified by the output keyword)

#write_every W -- write the output file every W sampled points.  This also
governs how often the code assesses the appropriateness of its Kriging
parameter.  W is set to 100 by default.

#guess (list of N numbers) -- if you are starting a new iteration of APS, you
can use guesses to make sure that the initial 'random' sample of S points
contains certain points of interest.  This could be useful for cases where you
know where the maximum likelihood point is in parameter space.  You do not need
to supply any guesses.  It is purely optional.


///////////////////////

IV. likelihoodinterface.cpp

In this section, we will outline the behavior of the member subroutines of the
class likelihood, which do all of the work associated with the APS
algorithm.  We will start out of order with

double likelihood::call_likelihood(double *v)

This is the only subroutine that user absolutely must change in order to
interface with her particular data and likelihood function.  The subroutine
accepts a point in parameter space v[] and uses that point (and calls to the
user's own external functions) to calculate a value of chi-squared.  The
subroutine returns chi-squared. As presented, it uses a cartoon likelihood 
function in two-dimensional parameter space.  There is also commented code 
which uses CAMB to constrain cosmological parameters against the 7 year 
release WMAP CMB data.

likelihood::likelihood(int nn, double *mns, double *mxs)

This class constructor only alots space in certain arrays so that the 
algorithm can do its work.  nn is the number of dimensions in your parameter 
space.  mns[] and mxs[] are the upper and lower bounds allowed
for each of your parameters.  This routine will set default values of kk (the
number of nearest neighbors used in the Gaussian process), nsamples 
(the number of candidates called for the Guassian process), and npts (the 
number of random points with which to start your algorithm).  You can assign 
different values for these parameters (as happens in aps_interface.cpp) 
before calling likelihood::initialize() with no ill effect.

The arrays allotted by this algorithm are as follows:

CLmin[] and CLmax[] -- once you start finding points that are within your
confidence limit (i.e. chi-squared <= target), these will store the minimum and
maximum value of each parameter in parameter space that satisfies your
confidence limit.

mins[] and maxs[] -- these will copy mns[] and mxs[] for use inside the
algorithm.


kptr[] -- points sampled by the algorithm are stored in kd trees to facilitate
rapid nearest neighbor searches.  In order to parallelize the Gaussian process 
portion of the code, each separate openmp thread must have its own kd tree to 
prevent simultaneous nearest neighbor searches from interfering with each other.
Note: each kd tree object has an associated integer 'diagnostic'.  The integer
is set to unity if the tree is self-consistent (all of the nodes are arranged
and related to each other the way they should be).  It is set to zero if they
are not.  If you ever doubt that the kd tree code is working properly, call
kptr[i]->check_tree(-1).  This will inspect the tree and then set 'diagnostic' 
to the appropriate value (the argument -1 tells the routine to check the entire
tree).  The subroutine likelihood::write_pts(char*) does this automatically. 
Depending on how often you are calling write_pts(), that may slow your code
down, so you can elect to comment that line out without ill-effect.

void likelihood::initialize(double **guesses, int nguess)

This routine will sample the initial random set of points and arrange them into
the kd trees to prepare for Kriging.  It will add the nguess points stored in
guesses[][] to that initial sample.

void likelihood::resume(char *inname)

If a run of the code was interrupted for some reason, you can call this routine
with inname referring to the name of the old output file.  This will read in
the data stored in that file and then resume the code where you left off.  Note:
reading in the data does not automatically set parameters such as kk, nsamples
or lingering, so aps_interface still sets those ``by hand''.

void likelihood::add_pt(double *v, double *dd, double fin, double sin, \
double chitrue)

After the Gaussian process has chosen a point to evaluate (this is done in
likelihood::sample_pts(int)), and likelihood::call_likelihood(double*) has
evaluated chi-squared for that point, this routine will add it to the kd trees 
storing the data.  It will also update the array storing the values of chi-
squared and other run-time data. The arguments are

v[] -- the point in parameter space to be added
dd[] -- the list of distances to the Gaussian process nearest neighbors
fin -- the value of chi-squared predicted by the Gaussian process (mu in text)
sin -- the value of sigma yielded by the Gaussina process (uncertainty on mu)
chitrue -- the value of chi-squared yielded by likelihood::call_likelihood(v)

In the event that there is not enough memory in the data arrays to accommodate
the new point (i.e. npts==room), this routine will also alot new memory for
those arrays with appropriate calls of delete [] and new [].

The data arrays that are updated by this routine are as follows

kpa[] -- this stores the value of the Kriging parameter at the time
likelihood::add_pt() was called

fpredicted[] -- this stores the value of chi-squared predicted by Kriging for
each point when it was sampled (\mu in the text)

sigma[] -- this stores the Gaussian process uncertainty associated with
fpredicted[] for each point.  Note: it actually stores 
\sigma/sqrt(Kriging parameter) so that the Kriging parameter can be
heuristically reset in the future (this will hopefully make more sense after we
have discussed the likelihood::set_kp() routine)

chisq[] -- stores the values of chi-squared associated with each pooint

ddmin[], ddmax[], ddavg[] -- store the minimum, maximum, and average distances
to the kk nearest neighbors used for each point during the Gaussian process

suboptime[] -- this stores the number of points chosen suboptimally (i.e.
because of large \sigma rather than proximity to chi-squared==target) at the
point when likelihood::add_pt() was called.

suboptct[] -- this stores the number of points sampled while the code was not
lingering, so that suboptime[]/suboptct[] gives a fraction of points that have
been chosen suboptimally so far.

If, after this algorithm has been called, the number of sampled points is an
integer multiple of writevery (W in section III), the algorithm will call 
write_pts() and write  the output.  It will also assess whether or not the 
Kriging parameter should be adjusted so that more or fewer points are chosen 
suboptimally.  

void likelihood::sample_pts(int delswit)

This routine will generate kk random points in parameter space, use a Gaussian
process to predict the associated values of chi-squared, choose (using the 
statistical test described in equation (2) of the text) one to evaluate, 
actually evaluate that point with likelihood::call_likelihood(v) and call l
ikelihood::add_pt() to add that point to the kd trees and other data arrays.  
The argument 'delswit' is just a switch.  If delswit>0, the routine behaves as 
just described.  If delswith<0, the routine will delete several static arrays 
that are necessary for its functioning and return without sampling anything.

This is also the routine where the code governing the linger option exists.

This routine also keeps track of how many points are chosen subopitmally (for
large values of Gaussian process \sigma, rather than proximity to 
chi-squared==target) and keeps track of the average time (both wall time and 
cpu ticks) spent on Kriging and the likelihood function.

void likelihood::write_pts()

This algorithm will write the data recorded so far to the file whose name is
specified in the parameter file (and stored in char masteroutname[]).  
It will also write the data related to timing the code to timingfile.sav.

There is a for loop at the beginning of this algorithm which checks whether or
not the kd trees are still properly constructed.  This can be excised for
efficiency if the user anticipates calling likelihood::write_pts() frequently.

double likelihood::covariogram(double *v1, double *v2, int coords)

Given two points in parameter space v1[] and v2[], this algorithm returns the
covariogram between them.  It is used in the Gaussian process.

void likelihood::predict(double **old, double *ffn, double *q, double *dd, \
int pts, int coords, double *mu, double *sig, int delswit)

This routine (also called in likelihood::sample_pts()) takes a test point q[], a
set of known points old[][] and uses a Gaussian process to predict chi-squared 
at q[].  The predicted value of chi-squared is stored in *mu.  The uncertainty 
on that prediction is stored in *sig.  ffn[] stores the chi-squared values 
associated with the points in old[][].  dd[] stores the distances in parameter 
space between q[] and the points in old[][] (this is given for free by the 
nearest neighbor search algorithm called in likelihood::sample_pts()).  
'delswit' behaves the same as in likelihood::sample_pts() i.e. if delswit>0 
the routine does what it is supposed to do; if delswit<0, the routine deletes 
several static arrays and returns.

double likelihood::test_kp(int pts, char *name)

This subroutine provides one means of testing the code's current choice of
Kriging parameter.  When called, it will generate pts random points in parameter
space, use a Gaussian process to predict the corresponding values of 
chi-squared, and then actually call the likelihood function to find their true 
values of chi-squared. The routine returns the fraction of these random points 
for which the true value of chi-squared was within 1-\sigma of the 
Gaussian-process-predicted value (where \sigma is the Gaussian process 
uncertainty as in the text).  If this value is too small, you can
increase the Kriging parameter.  If the value is too large, you can decrease 
the Kriging parameter.  

name[] specifies a file in which the results of this test (the number of points
generated, the number and fraction of points trapped, the number of data points 
sampled so far, the value of the Kriging parameter, and the number of points 
chosen suboptimally thus far) is stored.  Note: this routine actually uses the
parameter rkp in place of kriging_parameter, just in case it is called while the
code is lingering (and kriging_parameter is set to a very small number).

Because the routine actually spends the time to call the likelihood function, 
it calls likelihood::add_pt() to add those points to the list of known points.
``Waste not, want not.''

void likelihood::set_kp(int ignore)

This routine provides a way to heuristically set the Kriging parameter after the
algorithm has run for a while.  When called, this routine will consider all of
the data points currently stored, ignoring the first 'ignore' points, and then
set the Kriging parameter to whatever value is necessary such that 68% of those
points satisfy the condition 

|chisquared-mu|<=sigma*sqrt(Kriging parameter)

where mu is the value predicted by the Gaussian process and sigma is the 
uncertainty on that value (in the case that the Kriging parameter = 1).  This 
is why the array  sigma[] ignores the present value of the Kriging parameter 
(see the discussion of the arrays updated by likelihood::add_pt())

///////////////////////

V. Output produced by this software

This code will produce two main output files. 

filename.txt -- the file specified in the params.sav file

This is the principal file produced by likelihood::write_pts().
Each line of this file corresponds to a different point in parameter space
sampled by the code.  The first N columns are the values of the N-dimensions of
parameter space (they will be labeled with the names stored in the character
array likelihood.paramnames[][]).  The subsequent columns are as follows:

'chisq' -- chi-squared for the point (the actual chi-squared)

'predicted' -- the value of chi-squared predicted by Kriging for that point
(\mu in the text)

'sig' -- the uncertainty on 'predicted' (\sigma in the text).  Note: this is
the value stored in the sigma[] array multiplied by sqrt(Kriging parameter at
the time the point was sampled). See, for example, the discussion about the 
arrays updated by likelihood::add_pt() in Section IV.  Points sampled under
lingering will have sig=0.

'diffpct' -- the difference (as a ratio) between 'chisq' and 'predicted' i.e.
|(chisq - predicted)/chisq|

'ddmin', 'ddmax', 'ddavg' -- the minimum, maximum, and average distance from
the point to its kk nearest neighbors when the Gaussian process was used on the 
point.  Though you may not see a correlation between these numbers and 
'diffpct', you should see an obvious correlation wth 'sig' (if the 'dds' are 
small, 'sig' should be small and vice-versa)

'goodpct' -- the fraction of points sampled that are within the desired
confidence limit, i.e. that satisfy chi-squared <= target

'triedpct' and 'wayoff' -- the fraction of points for which 'predicted' was
within 10% of target or differed from target by more than 90%.  These are
provided as another gauge of the Kriging parameter.  If the Kriging parameter
is too small, 'triedpct' will near 1.0 as the code focuses only on finding
points that are near chi-squared == target.  If the Kriging parameter is too
large, way off will be near 1.0 as the opposite behavior will be emphasized. 
Note: points sampled during lingering do not count towards these ratios.

'suboptimal' -- the number of points sampled suboptimally (because of large
\sigma rather than proximity to chi-squared == target)

'suboptpct' -- the fraction of points (which are not sampled under lingering)
that are sampled suboptimally

'trappedpct' -- the fraction of points for which

     |chisquared-predicted|<=sqrt(Kriging parameter)*sig

Note: this statistic assesses each point against the Kriging parameter at the
time that point was sampled.  It does not go back and reevaluate whether old 
points would have been trapped using new values of the Kriging parameter.

'kp' -- the Kriging parameter at the time the point was evaluated.

Note: if you alter likelihood::write_pts() it is important that you do not
change any of the current labels and that 'kp' is always the last column. 
likelihood::resume() requires these conditions in order to properly parse the
data.

timingfile.sav

This file contains information about each run as an average.  

The first column is the name of the root_out.sav file to which the line is 
referring.  

The second column is the total number of points sampled.  

The third column ('ling') is the number of iterations spent lignering

The fourth column ('gd') is the number of points found within your desired 
confidence limit.  

'avgkrigtimewall' refers to the average number of seconds (since the last time 
likelihood::write_pts() was called) spent on sampling and Gaussian process
candidate points.

'avgaddtimewall' refers to the average number of seconds spent calling the 
likelihood function and adding points to the data arrays.  

'ns' refers to the number of random samples chosen for the Gaussian process.

'kk' is the number of nearest neighbors chosen for the Gaussian process.

'avgkrigtimecpu' and 'avgaddtimecpu' are the same as 'avgkrigtimewall' and
'avgaddtimewall' except that they refere to the number of elapsed cpu clock 
ticks (the output of the C++ subtoutine clock() )

'suboptpct' refers to the percentage of points chosen suboptimally. 

'ct' is the number of points chosen without lingering (i.e. the denominator of 
'suboptpct').

'kp' lists the current value of the Kriging parameter. 

The 'kd#' columns display the value of 'diagnostic' for the kd trees (they 
should be unity if everything is well-behaved). 

'target' is the chi-squared value you are looking for.

test_kp_file.txt -- the char* argument to likelihood::test_kp(int,char*)

This will record statistics whenever you use likelihood::test_kp() to evaluate
the current choice of Kriging parameter.  It will record: the number of random 
points used to evaluate the Kriging parameter; the number and percentage of 
those points for which the true chi-squared value was within 1-\sigma of the 
Kriging prediction, the total number of data points stored in the data arrays, 
the current value of the Kriging parameter, and the total number of points 
chosen suboptimally so far.

If errors are encountered in the Gaussian process, they are output to a file
named gp_error.sav
Makefile0000644000437200043720000000517511746310221012665 0ustar  danielsfdanielsfLAPACK =-L/astro/users/danielsf/lapack-3.1.1/ -llapack

BLAS =-L/astro/users/danielsf/BLAS/ -lblas

FORTRAN = -lgfortran

CAMB = -I/astro/users/danielsf/CAMB_110419/ \
-L/astro/users/danielsf/CAMB_110419/ -lcamb
#CAMB = -I/astro/users/danielsf/camb/ -L/astro/users/danielsf/camb -lcamb

CFITSIO = /astro/users/danielsf/cfitsio

WMAP = -I/astro/users/danielsf/WMAP7likelihood/ \
-L/astro/users/danielsf/WMAP7likelihood -lwmap7

WMAPLIBS = -L$(CFITSIO) -lcfitsio -L/astro/users/danielsf/lapack-3.1.1 \
-llapack -L/astro/users/danielsf/BLAS -lblas -lpthread

#C++ compiler
gg = g++ #-fopenmp -D OPENMP

#Fortran compiler
ff = gfortran #-fopenmp

goto_tools.o: goto_tools.h goto_tools.cpp
	$(gg) -c goto_tools.cpp $(LAPACK) $(BLAS) $(FORTRAN)

eigen_wrapper.o: goto_tools.o eigen_wrapper.cpp eigen_wrapper.h
	$(gg) -c eigen_wrapper.cpp goto_tools.o $(LAPACK) $(BLAS) \
	$(FORTRAN) -Wno-deprecated

camb_wrapper_wmap.o: camb_wrapper_wmap.F90
	$(ff) -c camb_wrapper_wmap.F90 $(CAMB)
	
wmap_wrapper.o: wmap_wrapper.F90
	$(ff) -c wmap_wrapper.F90 -I/astro/users/danielsf/WMAP7likelihood/ \
	$(WMAP) $(FORTRAN) $(WMAPLIBS) -lgfortran

kd.o: kd.cpp kd.h goto_tools.o
	$(gg) -c kd.cpp goto_tools.o $(LAPACK) $(BLAS) $(FORTRAN) \
         -Wno-deprecated

#use the compile command below if you want to use the CAMB and WMAP 7 
#likelihood function
#
#make sure that likelihood::call_likelihood() in likelihoodinterface.cpp
#uses the proper likelihood function
#
#likelihoodinterface.o: camb_wrapper_wmap.o likelihoodinterface.h \
likelihoodinterface.cpp kd.o goto_tools.o eigen_wrapper.o wmap_wrapper.o
#	$(gg) -D CMB -c likelihoodinterface.cpp goto_tools.o \
	eigen_wrapper.o camb_wrapper_wmap.o kd.o wmap_wrapper.o \
	$(LAPACK) $(BLAS) $(FORTRAN) -Wno-deprecated

#use the compile command below if you only want to use the cartoon likelihood 
#function

likelihoodinterface.o: likelihoodinterface.h \
likelihoodinterface.cpp kd.o goto_tools.o eigen_wrapper.o
	$(gg) -c likelihoodinterface.cpp goto_tools.o \
	eigen_wrapper.o kd.o \
	$(LAPACK) $(BLAS) $(FORTRAN) -Wno-deprecated


#CMB compilation statement
#aps: aps_interface.cpp goto_tools.o eigen_wrapper.o \
likelihoodinterface.o kd.o camb_wrapper_wmap.o wmap_wrapper.o
#	$(gg) -o aps aps_interface.cpp goto_tools.o eigen_wrapper.o kd.o \
	 likelihoodinterface.o camb_wrapper_wmap.o wmap_wrapper.o \
	 $(LAPACK) $(BLAS) $(FORTRAN) $(CAMB) $(WMAP) $(WMAPLIBS) \
	 -Wno-deprecated


#cartoon compilation statement
aps: aps_interface.cpp goto_tools.o eigen_wrapper.o \
likelihoodinterface.o kd.o 
	$(gg) -o aps aps_interface.cpp goto_tools.o eigen_wrapper.o kd.o \
	 likelihoodinterface.o \
	 $(LAPACK) $(BLAS) $(FORTRAN) \
	 -Wno-deprecated


goto_tools.cpp0000644000437200043720000001140611732143647014126 0ustar  danielsfdanielsf#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "goto_tools.h"

double normal_deviate(Ran *chaos, double mu, double sig){
 
 int i;
 double x1,x2,y1,y2;
 
 x1=chaos->doub();
 x2=chaos->doub();
 
 y1=sqrt(-2.0*log(x1))*cos(2.0*pi*x2);
 y2=sqrt(-2.0*log(x1))*sin(2.0*pi*x2);
 
 
 return mu+y1*sig;
 
}

int compare_char(char *s1, char *s2){
 //are two character strings the same?
 //if so, return 1
 //if not, return 0

 int i;
 //printf("\ncomparing %s %s\n",s1,s2);
 for(i=0;i<letters && (s1[i]!=0 || s2[i]!=0);i++){
  if(s1[i]!=s2[i])return 0;
 }
 return 1;
 
}

double power(double arg,int raised){

  //return arg raised to the integer power 'raised'

  int n;
  double ans;

  if(raised==0)return 1.0;
  else{ans=1.0;
  for(n=0;n<raised;n++){
    ans=ans*arg;
  }

  return ans;
  }

}

void kill(char *words){
 double junk;
 FILE *output;
 
 //write the character string words[] to the terminal
 //then hang, waiting for an input double
 
 printf("%s\n",words);
 scanf("%lf",&junk);
}

/*just going to lift this out of Numerical Recipes p 109*/
void polint(double *xa, double *ya, int n, double x, double *y, double *dy){
/*recall: because this is from NR, it has arrays start at element unity*/
	int i,m,ns=1,isfd;
	double den,dif,dift,ho,hp,w;
	double c[n+1],d[n+1];
	char scream[100];
	dif=fabs(x-xa[1]);
	for(i=1;i<=n;i++){
		if((dift=fabs(x-xa[i]))<dif){
			ns=i;
			dif=dift;
		}
		c[i]=ya[i];
		d[i]=ya[i];
	}
	*y=ya[ns--];
	for(m=1;m<n;m++){
		for(i=1;i<=n-m;i++){
			ho=xa[i]-x;
			hp=xa[i+m]-x;
			w=c[i+1]-d[i];
			if((den=ho-hp)==0.0){printf("Error in routine polint ");
				for(isfd=1;isfd<=n;isfd++)printf(" (%e, %e) ",xa[isfd],ya[isfd]);
				printf(" want %e \n",x); 
				sprintf(scream,"stop");
				kill(scream);
			
			
			}
			/*This error can occur only if two input xas re (to within roundoff) identical*/
			den=w/den;
			d[i]=hp*den;
			c[i]=ho*den;
		}
		*y+=(*dy=(2*ns<(n-m)?c[ns+1]:d[ns--]));
	}
}


double interpolate(double *x, double *y, double target, int el){

  //this is a wrapper for the Numerical Recipes polynomial interpolation
  //routine polint()
  
  //x[] is the array of x values
  
  //y[] is the array of y values
  
  //'target' is the x value for which you want to interpolate a y value
  
  //'el' is the number of elements in x[] and y[]
  
  //x[] must be monotonic, but it can be either ascending or descending

 double *xt,*yt;
 double xint[5],yint[5],err,ans;
 int i,n,min;
 
 
 if(x[0]>x[1]){
  xt=new double[el];
  yt=new double[el];
  for(i=0;i<el;i++){
   xt[i]=x[el-1-i];
   yt[i]=y[el-1-i];
  
  }
 }
 else{
  xt=x;
  yt=y;
 }

 
 for(i=0;xt[i]<target && i<el;i++);
 if(xt[i]==target){ans=yt[i];}
 else{
  if(i<2)min=0;
  else if(i>el-2)min=el-4;
  else min=i-2;
  for(n=1;n<5;n++){
   xint[n]=xt[min+n-1];
   yint[n]=yt[min+n-1];
  }
  //printf("min was %d %e i %d %e target %e lim %e %e\n",min,xt[min],i,xt[i],target,xt[0],xt[el-1]);
  polint(xint,yint,4,target,&ans,&err);
 }
 
 if(x[0]>x[1]){
  delete [] xt;
  delete [] yt;
 }
 
 //if(!(ans>0.0))scanf("%lf",&junk);
 return ans;

}


//the routines below are just the merge sort routine from Numerical Recipes;

int scanner(double *m, int *indices, int index, int el){
/*this will take the matrix m and put everything in it with value
greater than element m[index] to the right of that element
and everything less than m[index] to the left; it then returns
the new index of that anchored element (which you now *know* is in
the right spot*/

	double toobig,toosmall;
	int i,j,n,k,pp;
	int itoobig,itoosmall;
        FILE *output;

	itoobig=0;
	itoosmall=0;
	for(i=0;i<el;i++){
	  if(m[i]<=m[index])itoosmall++;
	  
	} 
	
	toobig=m[index];
	itoobig=indices[index];
	
	indices[index]=indices[itoosmall-1];
	m[index]=m[itoosmall-1];
	
	index=itoosmall-1;
	m[index]=toobig;
	indices[index]=itoobig;

	i=0;
	j=el-1;
	n=0;
	pp=1;
	while(index<el-1 && n<el-1){
	 for(;m[i]<=m[index] && i<index;i++,n++);
	 itoobig=indices[i];
	 toobig=m[i];
	 
	 for(;m[j]>m[index] && j>index;j--,n++){
	  
	 }
	 itoosmall=indices[j];
	 toosmall=m[j];
	
	 
	 
	 if(toosmall<toobig){
	
	 //in case it ran to the end and i found m[index]
	  m[i]=toosmall;
	  indices[i]=itoosmall;
	 
	  m[j]=toobig;
	  indices[j]=itoobig;
	 }
	 
	}

	return index;
}

void sort(double *m, int *indices, int el){
	double *newm,junk;
	int k,i,*newi,ii,diff;
	
	if(el==2){
	 if(m[0]>m[1]){
	  
	   junk=m[1];
	   m[1]=m[0];
	   m[0]=junk;
	   
	   k=indices[1];
	   indices[1]=indices[0];
	   indices[0]=k;
	  
	 }
	}
	else if(el>2){
	
	diff=0; //just in case all elements are identical
	for(ii=1;ii<el;ii++)if(m[ii]!=m[0])diff=1;
	  
	if(diff==1){

	   i=scanner(m,indices,el/2,el);

	   if(i+1<el){
	     newm=&m[i+1];
	    sort(m,indices,i+1);
	     newi=&indices[i+1];
	    sort(newm,newi,el-i-1);
	  }
	  else{
	    sort(m,indices,el-1);
	  }
	
	 }
	
	}
}

goto_tools.h0000644000437200043720000000147511732141364013572 0ustar  danielsfdanielsf#define pi 3.141592654
#define letters 100

void kill(char*);


int compare_char(char*,char*);

double power(double,int);

struct Ran{

unsigned long long u,v,w;
Ran(unsigned long long j) : v(4101842887655102017LL), w(1) {
	u = j^v; int64();
	v=u; int64();
	w=v; int64();
}
inline unsigned long long int64(){
	u = u * 2862933555777941757LL + 7046029254386353087LL;
	v ^= v >> 17; v ^= v<<31; v ^= v>>8;
	w = 4294957665U*(w & 0xffffffff) + (w>>32);
	unsigned long long x = u^(u<<21); x^=x>>35;x^=x<<4;
	return (x+v)^w;
}
inline double doub(){return 5.42101086242752217e-20 * int64();}
inline unsigned int int32(){return (unsigned int)int64();}
};

void polint(double*,double*,int,double,double*,double*);

double interpolate(double*,double*,double,int);

void sort(double*,int*,int);

double normal_deviate(Ran*,double,double);


eigen_wrapper.h0000644000437200043720000000025011732141364014217 0ustar  danielsfdanielsfvoid matrix_multiply(double**, int, int, double**, int, int, double**);

double check_inversion(double**,double**,int);

void invert_lapack(double**,double**,int,int);
eigen_wrapper.cpp0000644000437200043720000000460711732141360014560 0ustar  danielsfdanielsf#include <stdio.h>
#include <math.h>
#include "eigen_wrapper.h"
#include "goto_tools.h"

extern "C" void dgetrf_(int*,int*,double*,int*,int*,int*);

extern "C" void dgetri_(int*,double*,int*,int*,double*,int*,int*);

void matrix_multiply(double **a, int ar, int ac, double **b, int br, int bc, \
double **p){

	int i,j,m,n;
	
	if(ac!=br)printf(\
	"WARNING dimensions are wrong  %d x %d  times %d x %d\n", \
	ar,ac,br,bc);
	for(i=0;i<ar;i++){
	for(j=0;j<bc;j++){
		p[i][j]=0.0;
		for(m=0;m<ac;m++){
			p[i][j]+=a[i][m]*b[m][j];
		}
	
	}}

}


double check_inversion(double **m, double **min, int el){
	double ans,**i;
	int j,k;
	
	
	i=new double*[el];
	for(j=0;j<el;j++){
	 i[j]=new double[el];
	}

	
	matrix_multiply(m,el,el,min,el,el,i);
	
	ans=-1.0;
	for(j=0;j<el;j++){
	for(k=0;k<el;k++){
		if(j==k){
			if(fabs(1.0-i[j][k])>ans){ans=fabs(1.0-i[j][k]);
			 //printf("    %d %d %e\n",j,k,ans);
			}
		}
		else{
			if(fabs(i[j][k])>ans){ans=fabs(i[j][k]);
			 //printf("   %d %d %e\n",j,k,ans);
			}
		}
	
	
	}}
	
	for(j=0;j<el;j++)delete [] i[j];
	delete [] i;
	
	return ans;

}

void invert_lapack(double **matin, double **min, int el,int verb){

//inverts the matrix matin[][] and stores the invers in min[][]
//uses lapack routines

	//double matrix[maxdata*maxdata],work[maxdata];
	int i,j;
	int info,lda,m,n,lwork;
	//int ipiv[maxdata];
	
	double *matrix, *work;
	int *ipiv;
	//printf("in invert_lapack\n");
	//hi
	/*matrix=new double [maxdata*maxdata];
	work=new double [maxdata];
	ipiv=new int [maxdata];*/
	
	if(verb==2)printf("in invert lapack el %d\n",el);
	
	lwork=4*el;
	matrix=new double[el*el];
	work=new double[lwork];
	ipiv=new int [el];
	
	if(verb==2)printf("about to assign to matrix\n");
	for(i=0;i<el;i++){
	for(j=0;j<el;j++){
		matrix[j*el+i]=matin[i][j]; //Fortran is backwards
					//look up "column major order"
					//on wikipedia
	}
	}
	if(verb==2)printf("assigned to matrix\n");
	
	
	m=el;
	n=el;
	lda=el;
	if(verb==2)printf("about to call dgetrf\n");
	dgetrf_(&m,&n,matrix,&lda,ipiv,&info);
	//cout<<"after dgetrf info is "<<info<<endl;
	//printf("after dgetrf info is %d\n",info);
	
	if(verb==2)printf("about to call dgetri\n");
	dgetri_(&n,matrix,&lda,ipiv,work,&lwork,&info);
	//cout<<"after dgetri info is "<<info<<endl;
 	//printf("after dgetri info is %d\n",info);

	for(i=0;i<el;i++){
	for(j=0;j<el;j++){
	
		min[i][j]=matrix[j*el+i];
	
	}
	}

	delete [] matrix;
	delete [] work;
	delete [] ipiv;


}

camb_wrapper_wmap.F900000744000437200043720000002723611732141370015200 0ustar  danielsfdanielsf!     Code for Anisotropies in the Microwave Background
!     by Antony Lewis (http://cosmologist.info/) and Anthony Challinor
!     See readme.html for documentation. This is a sample driver routine that reads
!     in one set of parameters and produdes the corresponding output. 

    subroutine camb_wrap(pparams,ell,cltt,clte,clee,clbb)
        use IniFile
        use CAMB
        use LambdaGeneral
        use Lensing
       !krig use RECFAST 
	use Reionization   !krig    
#ifdef NAGF95
        use F90_UNIX
#endif
        implicit none
	!int sfdel(3000)
	!real sfdcl(3000)
      
        Type(CAMBparams) P
        
        character(LEN=Ini_max_string_len) numstr, VectorFileName, &
            InputFile, ScalarFileName, TensorFileName, TotalFileName, LensedFileName
        integer i
        character(LEN=Ini_max_string_len) TransferFileNames(max_transfer_redshifts), &
               MatterPowerFileNames(max_transfer_redshifts), outroot
        real(dl) output_factor, Age
	real(dl), intent(in) :: pparams(6)
	real(dl), intent(out) :: cltt(3000),ell(3000)
	real(dl), intent(out) :: clte(3000),clee(3000),clbb(3000)

#ifdef WRITE_FITS
       character(LEN=Ini_max_string_len) FITSfilename
#endif

#ifndef NAGF95
#ifndef __INTEL_COMPILER_BUILD_DATE
        integer iargc
        !external iargc
#endif        
#endif
        logical bad

        InputFile = ''
   

       ! if (iargc() /= 0)  call getarg(1,InputFile)
       ! if (InputFile == '') stop 'No parameter input file'

       ! call Ini_Open(InputFile, 1, bad, .false.)
       ! if (bad) stop 'Error opening parameter file'

       ! Ini_fail_on_not_found = .false.
    
       ! outroot = Ini_Read_String('output_root')
       ! if (outroot /= '') outroot = trim(outroot) // '_'
        
        call CAMB_SetDefParams(P)

        P%WantScalars = .TRUE.
        P%WantVectors = .FALSE.
        P%WantTensors = .FALSE.
        
        P%OutputNormalization=outNone
        !if (Ini_Read_Logical('COBE_normalize',.false.))  P%OutputNormalization=outCOBE
        output_factor = 7.4311e12

        P%WantCls= P%WantScalars .or. P%WantTensors .or. P%WantVectors

        P%WantTransfer=.TRUE.
        
        P%NonLinear = 0
   
        P%DoLensing = .false.
        if (P%WantCls) then
          if (P%WantScalars  .or. P%WantVectors) then
           P%Max_l = 2000
           P%Max_eta_k = 4000
           if (P%WantScalars) then
             P%DoLensing = .false.
             if (P%DoLensing) lensing_method = 1
           end if
           if (P%WantVectors) then
            if (P%WantScalars .or. P%WantTensors) stop 'Must generate vector modes on their own'
            i = Ini_Read_Int('vector_mode')
            if (i==0) then 
              vec_sig0 = 1
              Magnetic = 0
            else if (i==1) then
              Magnetic = -1
              vec_sig0 = 0
            else
              stop 'vector_mode must be 0 (regular) or 1 (magnetic)'
            end if 
           end if
          end if

          if (P%WantTensors) then
           P%Max_l_tensor = 1500
           P%Max_eta_k_tensor =  3000
          end if
        endif

                
!  Read initial parameters.
       
       w_lam = -1.0
       cs2_lam = 1.0

       P%h0     = pparams(3)
 
       !if (Ini_Read_Logical('use_physical',.false.)) then 

        P%omegab = pparams(1)/(P%H0/100)**2
        P%omegac = pparams(2)/(P%H0/100)**2
        P%omegan = 0.0/(P%H0/100)**2
        !P%omegav = 1- Ini_Read_Double('omk') - P%omegab-P%omegac - P%omegan
	P%omegav = 1- 0.0 - P%omegab-P%omegac - P%omegan
        
	!write(*,*)'ombh2 ',pparams(1),' omch2 ',pparams(2),' h ',P%H0
  
      ! else
       
       ! P%omegab = Ini_Read_Double('omega_baryon')
       ! P%omegac = Ini_Read_Double('omega_cdm')
       ! P%omegav = Ini_Read_Double('omega_lambda')
       ! P%omegan = Ini_Read_Double('omega_neutrino')

       !end if

       P%tcmb   = 2.726
       P%yhe    = 0.24
       P%Num_Nu_massless  = 3.04
       P%Num_Nu_massive   = 0.0
   
       P%nu_mass_splittings = .true.
       P%Nu_mass_eigenstates = 1
       if (P%Nu_mass_eigenstates > max_nu) stop 'too many mass eigenstates'
       numstr = '0'
       if (numstr=='') then
         P%Nu_mass_degeneracies(1)= P%Num_nu_massive
       else
        read(numstr,*) P%Nu_mass_degeneracies(1:P%Nu_mass_eigenstates)
       end if
       numstr = '1'
       if (numstr=='') then
        P%Nu_mass_fractions(1)=1  
        if (P%Nu_mass_eigenstates >1) stop 'must give nu_mass_fractions for the eigenstates'
       else
        read(numstr,*) P%Nu_mass_fractions(1:P%Nu_mass_eigenstates)
       end if

       if (P%NonLinear==NonLinear_lens .and. P%DoLensing) then
          if (P%WantTransfer) &
             write (*,*) 'overriding transfer settings to get non-linear lensing'
          P%WantTransfer  = .true.
          call Transfer_SetForNonlinearLensing(P%Transfer)
          P%Transfer%high_precision=  .false.
       
       else if (P%WantTransfer)  then
        P%Transfer%high_precision=  .false.
        P%transfer%kmax          =  10
        P%transfer%k_per_logint  =  0
        P%transfer%num_redshifts =  1
       ! if (P%transfer%num_redshifts > max_transfer_redshifts) stop 'Too many redshifts'
       ! do i=1, P%transfer%num_redshifts
	
        !     write (numstr,*) i 
        !     numstr=adjustl(numstr)
        !     P%transfer%redshifts(i)  = Ini_Read_Double('transfer_redshift('//trim(numstr)//')',0._dl)
        !     TransferFileNames(i)     = Ini_Read_String('transfer_filename('//trim(numstr)//')')
        !     if (TransferFileNames(i)/= '') &
        !           TransferFileNames(i) = trim(outroot)//TransferFileNames(i)
        !     MatterPowerFilenames(i)  = Ini_Read_String('transfer_matterpower('//trim(numstr)//')')
        !     if (MatterPowerFilenames(i) /= '') &
        !         MatterPowerFilenames(i)=trim(outroot)//MatterPowerFilenames(i)
        
	!end do
        !P%transfer%kmax=P%transfer%kmax*(P%h0/100._dl)
                
       !else
       !  P%transfer%high_precision = .false.
       !endif
  
        !krig P%Reionization = .true.
        !krig P%use_optical_depth = P%Reionization .and. .true.
	P%Reion%use_optical_depth = .true.
  
       ! if ( P%use_optical_depth) then
        
	      P%Reion%optical_depth = pparams(4)
	      P%Reion%delta_redshift=0.5 !krig
        
	 !  else if (P%Reionization) then
         !     P%Reion%redshift = Ini_Read_Double('re_redshift')
         !     P%Reion%fraction = Ini_Read_Double('re_ionization_frac')
       ! end if 

           Ini_fail_on_not_found = .false. 
           
       ! RECFAST_fudge = 1.14

           i = 1
           if (i==2) then
           ! use_Dubrovich = .true.
           else if (i/=1) then
             stop 'Unknown recombination'
           end if

           P%InitPower%nn = 1
           if (P%InitPower%nn>nnmax) stop 'Too many initial power spectra - increase nnmax in InitialPower'
           P%InitPower%rat(:) = 1
           do i=1, P%InitPower%nn
              write (numstr,*) i 
              numstr=adjustl(numstr)
              P%InitPower%an(i) = pparams(5)!&
                   !Ini_Read_Double('scalar_spectral_index('//trim(numstr)//')')

              P%InitPower%n_run(i) = 0!&
                   !Ini_Read_Double('scalar_nrun('//trim(numstr)//')',0._dl)
    
            !  if (P%WantTensors) then
            !     P%InitPower%ant(i) = Ini_Read_Double('tensor_spectral_index('//trim(numstr)//')')
            !     P%InitPower%rat(i) = Ini_Read_Double('initial_ratio('//trim(numstr)//')')
            !  end if              

              P%InitPower%ScalarPowerAmp(i) = pparams(6)
              !Always need this as may want to set tensor amplitude even if scalars not computed
           end do
      
            if (P%WantScalars .or. P%WantTransfer) then
            P%Scalar_initial_condition = 1
            if (P%Scalar_initial_condition == initial_vector) then
                P%InitialConditionVector=0
              numstr = '-1 0 0 0 0'
              read (numstr,*) P%InitialConditionVector(1:initial_iso_neutrino_vel)
            end if
        end if

        
      ! if (P%WantScalars) then
      !    ScalarFileName = trim(outroot)//Ini_Read_String('scalar_output_file')
      !    LensedFileName =  trim(outroot) //Ini_Read_String('lensed_output_file')
      !  end if
      !  if (P%WantTensors) then
      !    TensorFileName =  trim(outroot) //Ini_Read_String('tensor_output_file')
      !   if (P%WantScalars) TotalFileName =  trim(outroot) //Ini_Read_String('total_output_file')
      !  end if
      !  if (P%WantVectors) then
      !    VectorFileName =  trim(outroot) //Ini_Read_String('vector_output_file')
      !  end if
         
#ifdef WRITE_FITS
       ! if (P%WantCls) then
       ! FITSfilename =  trim(outroot) //Ini_Read_String('FITS_filename',.true.)
       ! if (FITSfilename /='') then
       ! inquire(file=FITSfilename, exist=bad)
       ! if (bad) then
       !  open(unit=18,file=FITSfilename,status='old')
       !  close(18,status='delete')
       ! end if
       !end if
       ! end if
#endif        
       

       Ini_fail_on_not_found = .false. 

!optional parameters controlling the computation

       P%AccuratePolarization = .true.
       P%AccurateReionization = .false.
       P%AccurateBB = .false.
        
       !Mess here to fix typo with backwards compatibility
       if (Ini_Read_String('do_late_rad_trunction') /= '') then
         DoLateRadTruncation = .true.
        ! if (Ini_Read_String('do_late_rad_truncation')/='') stop 'check do_late_rad_xxxx'
       else
        DoLateRadTruncation = .true.
       end if
       DoTensorNeutrinos = .false.
       FeedbackLevel = 0
       
       P%MassiveNuMethod  = 3

       ThreadNum      = 0 !Ini_Read_Int('number_of_threads',0)
       AccuracyBoost  = 1 !Ini_Read_Double('accuracy_boost',1.d0)
       lAccuracyBoost = 1 !Ini_Read_Double('l_accuracy_boost',1.d0)
       lSampleBoost   = 1 !Ini_Read_Double('l_sample_boost',1.d0)

       !if (outroot /= '') then
       !  call Ini_SaveReadValues(trim(outroot) //'params.ini',1)
       !end if

       !call Ini_Close

       if (CAMB_ValidateParams(P)) then
       !write(*,*)'it is valid'
#ifdef RUNIDLE
       call SetIdle
#endif 

       if (FeedbackLevel > 0) then
         Age = CAMB_GetAge(P) 
         write (*,'("Age of universe/GYr  = ",f7.3)') Age  
       end if 

	!write(*,*)'about to call get results in inidriver '
       call CAMB_GetResults(P)
    
        if (P%WantTransfer .and. .not. (P%NonLinear==NonLinear_lens .and. P%DoLensing)) then
         !call Transfer_SaveToFiles(MT,TransferFileNames)
         !call Transfer_SaveMatterPower(MT,MatterPowerFileNames)
         if ((P%OutputNormalization /= outCOBE) .or. .not. P%WantCls)  call Transfer_output_sig8(MT)
        end if

        
  
         if (P%OutputNormalization == outCOBE) then

            if (P%WantTransfer) call Transfer_output_Sig8AndNorm(MT)
           
          end if
         
	! open(unit=55,file='sfdoutput.sav',form='formatted',status='unknown')
	! do i=lmin,CP%Max_l
	! write(55,'(1I5,3E15.5)')i ,output_factor*Cl_scalar(i,1,C_Temp:C_Cross)
	! end do
	! close(55)
      
        do i=lmin,CP%Max_l
	  ell(i)=i
	  cltt(i)=output_factor*Cl_scalar(i,1,C_Temp)
	  clee(i)=output_factor*Cl_scalar(i,1,C_E)
	  clte(i)=output_factor*Cl_scalar(i,1,C_Cross)
	  clbb(i)=0.0
	  !write(*,*)'writing ',i,output_factor*Cl_scalar(i,1,C_Temp)
	end do
      
     
      else
       write(*,*)'your parameters make no sense, but I do not care'
       !write(*,*)lmin,CP%Max_l
       do i=1,3000
         !write(*,*)'i ',i
	  ell(i)=i
	  cltt(i)=0.01
	  clee(i)=0.01
	  clte(i)=0.01
	  clbb(i)=0.01
	  !write(*,*)'writing ',i,output_factor*Cl_scalar(i,1,C_Temp)
	end do
	!write(*,*)'end do'
      end if
     ! write(*,*)'end if'
  end if
  !write(*,*)'leaving camb routine'
        end subroutine camb_wrap


wmap_wrapper.F900000644000437200043720000000235511732141370014210 0ustar  danielsfdanielsfmodule wmapwrap


use wmap_likelihood_7yr
use wmap_options
use wmap_util
logical::initswitch=.true.
contains

subroutine wmaplikeness_actual(cltt,clte,clee,clbb,likeout)

	integer i
	integer :: tt_npix, teeebb_npix
	real(8) cltt(3000),clte(3000),clee(3000),clbb(3000),ell(3000)
	real(8) like(num_WMAP)
	real(8) likeout
	
	!these settings will cause the likelihood function to consider only
	!the TT correlation spectrum, and to work directly in C_l space,
	!rather than pixel space.
	
	use_TT               = .true.
	use_TE               = .false. 
	use_lowl_TT          = .false.
	use_lowl_pol         = .false.
	use_TT_beam_ptsrc    = .true.
	
	use_gibbs=.false.
	!lowl_max=30
	
	if(initswitch) then
	call wmap_likelihood_init
	initswitch=.false.
	end if
	call wmap_likelihood_dof( tt_npix, teeebb_npix )
	call wmap_likelihood_compute(cltt,clte,clee,clbb,like)
	
	likeout=0.0
	
	likeout=2.0*like(1)+2.0*like(4)

	!do i=1,num_WMAP
	!write(*,*)'i ',i,' likepart ',like(i)
	!end do
	
	

end subroutine wmaplikeness_actual

end module

subroutine wmaplikeness(cltt,clte,clee,clbb,likeout)
use wmapwrap
real(8) cltt(3000),clte(3000),clee(3000),clbb(3000),ell(3000)
real(8) likeout
call wmaplikeness_actual(cltt,clte,clee,clbb,likeout)

end subroutine wmaplikeness
kd.h0000644000437200043720000000110111732141364011762 0ustar  danielsfdanielsf//this is the 27 February 2012 rewrite

class kd_tree{
 private:
  int **tree;
  int masterparent;
  int room,roomstep;
  double tol;
  
  void confirm(int,int,int,int);
  void organize(int*,int,int,int,int);
  int find_node(double*);
  void neigh_check(double*,int,int*,double*,int,int);
  
 public:
  int dim,pts,diagnostic,xplr;

  double **data,*maxs,*mins;

  kd_tree(int,int,double**,int*,double*,double*);
  ~kd_tree();
 
 void check_tree(int);
 double distance(double*,double*);

 void write_tree(char*);
 void add(double*);
 void nn_srch(double*,int,int*,double*);
};
kd.cpp0000644000437200043720000002431011732157551012331 0ustar  danielsfdanielsf//this is the 27 February 2012 rewrite

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "goto_tools.h"
#include "kd.h"

kd_tree::~kd_tree(){
  int i;
  delete [] mins;
  delete [] maxs;
  for(i=0;i<room;i++){
    delete [] tree[i];
    delete [] data[i];
  }
  delete [] tree;
  delete [] data;
  
}

kd_tree::kd_tree(int dd, int pp, double **mm, int *innin, double *nmin, \
double *nmax){
  
   int i,j,k,l,*inn,*use,*uinn,rct,lct,inp;
   double *tosort;
   
   tol=1.0e-7;
   
   room=pp;
   roomstep=10000;
   
   dim=dd;
   pts=pp;
   diagnostic=1;
   
   data=new double*[pts];
   for(i=0;i<pts;i++)data[i]=new double[dim];
   tree=new int*[pts];
   for(i=0;i<pts;i++){
     tree[i]=new int[4];
     //tree[i][0] will be the dimension being split
     //tree[i][1] will be left hand node (so, lt)
     //tree[i][2] will be right hand node (so, ge)
     //tree[i][3] will be parent
   }
   maxs=new double[dim];
   mins=new double[dim];
   
   for(i=0;i<dim;i++){
     mins[i]=nmin[i];
     maxs[i]=nmax[i];
   }
   
   for(i=0;i<pts;i++){
     for(j=0;j<dim;j++)data[i][j]=mm[i][j];
   }

   tosort=new double[pts];
   inn=new int[pts];

   for(i=0;i<pts;i++){
     tosort[i]=data[i][0];
     inn[i]=i;
   }
   
   sort(tosort,inn,pts);
   
   inp=pts/2;
   while(inp>0 && tosort[inp]-tosort[inp-1]<tol){
     //make sure that the division doesn't come in the middle of a bunch
     //of identical points
     inp--;
   }
   
   masterparent=inn[inp];
   //printf("inp %d masterparent %d\n",inp,masterparent);
   if(inp==0){
      rct=pts-1;
      lct=0; 
   }
   else{  
   rct=0;
   lct=0;
     for(j=0;j<pts;j++){
       if(masterparent!=j){
          if(data[j][0]<tosort[inp])lct++;
	  else rct++;
       }
     } 
   }//of masterparent!=0
   
   if(lct>0){
     use=new int[lct];
    
     
     for(j=0,i=0;i<pts;i++){
       if(i!=masterparent){
         if(data[i][0]<tosort[inp]){
	   use[j]=i;
	   j++;
	 }
       }
     }
     
     organize(use,masterparent,1,lct,1);
     
     delete [] use;
    
   }
   else tree[masterparent][1]=-1;
   
   if(rct>0){
     use=new int[rct];
     
     for(j=0,i=0;i<pts;i++){
       if(i!=masterparent){
         if(data[i][0]>=tosort[inp]){
	   use[j]=i;
	   j++;
	 } 
       }
     } 
     
     organize(use,masterparent,1,rct,2);
     
     delete [] use;
    
   }
   else tree[masterparent][2]=-1;
   
   tree[masterparent][3]=-1;
   tree[masterparent][0]=0;
   
   
   
   
  delete [] tosort;
  delete [] inn;

}

void kd_tree::organize(int *use, int parent, int idim, int ct, int dir){

    int i,j,k,l,*inn,newparent,ff,fb,inp;
    double *tosort,pivot,nn;
   
   if(idim>=dim)idim=0;
   
   /*printf("parent %d pdim %d dir %d\n",parent,idim-1,dir);
   for(i=0;i<ct;i++)printf("%d\n",use[i]);
   printf("\n");*/
   
  
   if(ct>2){
      inn=new int[ct];
      for(i=0;i<ct;i++)inn[i]=i;
      tosort=new double[ct];
      for(i=0;i<ct;i++)tosort[i]=data[use[i]][idim];
   
      sort(tosort,inn,ct);
   
      inp=ct/2;
      while(inp>0 && tosort[inp]-tosort[inp-1]<tol)inp--;
   
      newparent=use[inn[inp]];   
      pivot=data[newparent][idim];
   
      tree[parent][dir]=newparent;
      tree[newparent][3]=parent;
      tree[newparent][0]=idim;
      
      //now I need to re-order use[] so that I can pass it to another
      //call of ::organize and have the proper indices available
      
         k=use[inp];
	 use[inp]=newparent;
	 use[inn[inp]]=k;
	 
	 delete [] inn;
	 delete [] tosort;
      
      if(inp!=0){

	 ff=0;
	 fb=ct-1;
	 while(ff<inp && fb>inp){
	    if(data[use[ff]][idim]<pivot && \
	    data[use[fb]][idim]>=pivot){
	      ff++;
	      fb--;
	    }
	    else if(data[use[ff]][idim]>=pivot && \
	    data[use[fb]][idim]>=pivot){
	      fb--;
	    }
	    else if(data[use[ff]][idim]<pivot && \
	    data[use[fb]][idim]<pivot){
	      ff++;
	    }
	    else if(data[use[ff]][idim]>=pivot && \
	    data[use[fb]][idim]<pivot){
	      k=use[fb];
	      use[fb]=use[ff];
	      use[ff]=k;
	    }
	 }
	 
	 organize(use,newparent,idim+1,inp,1);
	 organize(&use[inp+1],newparent,idim+1,ct-inp-1,2);
	 
      
      }//if(inp!=0)
      else{
	
	tree[newparent][1]=-1;
	organize(&use[1],newparent,idim+1,ct-1,2);
	
      }//if(inp==0)
      
      
   }//if(ct>2)
   else if(ct==2){
     if(data[use[0]][idim]<data[use[1]][idim]){
       tree[parent][dir]=use[1];
       tree[use[1]][1]=use[0];
       tree[use[1]][2]=-1;
       tree[use[1]][3]=parent;
       tree[use[1]][0]=idim;
     }
     else{
       tree[parent][dir]=use[1];
       tree[use[1]][1]=-1;
       tree[use[1]][2]=use[0];
       tree[use[1]][3]=parent;
       tree[use[1]][0]=idim;
     }
     
     tree[use[0]][0]=idim+1;
     if(tree[use[0]][0]>=dim)tree[use[0]][0]=0;
     tree[use[0]][1]=-1;
     tree[use[0]][2]=-1;
     tree[use[0]][3]=use[1];
   
   }
   else if(ct==1){
      tree[parent][dir]=use[0];
      tree[use[0]][1]=-1;
      tree[use[0]][2]=-1;
      tree[use[0]][3]=parent;
      tree[use[0]][0]=idim;
   }
   else if(ct==0)printf("WARNING called organize with ct==0\n");
  
}

void kd_tree::write_tree(char *name){
   
   int i,j,k,l;
   FILE *output;

   output=fopen(name,"w");
   for(i=0;i<pts;i++){
     fprintf(output,"%d tree dim %d l %d r %d p %d ",\
     i,tree[i][0],tree[i][1],tree[i][2],tree[i][3]);
     fprintf(output,"    ");
     for(j=0;j<dim;j++)fprintf(output,"p%d %e ",j,data[i][j]);
     fprintf(output,"\n");
   
   
   }

   fclose(output);
}

void kd_tree::check_tree(int where){
   int i;
   

   
   if(where<0)where=masterparent;
   
      //printf("checking %d %d %d %d\n",where,tree[where][1],tree[where][2],tree[where][0]);
   
   if(tree[where][1]>-1)confirm(tree[where][0],where,1,tree[where][1]);
   
   //printf("confirmed 1\n");
   
   if(tree[where][2]>-1)confirm(tree[where][0],where,2,tree[where][2]);
   //printf("confirmed 2\n");
   
   if(tree[where][1]>-1)check_tree(tree[where][1]);
   if(tree[where][2]>-1)check_tree(tree[where][2]);
     
 
}

void kd_tree::confirm(int idim, int compareto, int dir, int where){
   
   //printf("confirm %d %d %d %d\n",idim,compareto,dir,where);
   
   if(dir==1){
     if(data[where][idim]>=data[compareto][idim])diagnostic=0;
     
     if(tree[where][1]>-1)confirm(idim,compareto,dir,tree[where][1]);
     if(tree[where][2]>-1)confirm(idim,compareto,dir,tree[where][2]);
     
   }
   else if(dir==2){
     if(data[where][idim]<data[compareto][idim])diagnostic=0;
     
     if(tree[where][1]>-1)confirm(idim,compareto,dir,tree[where][1]);
     if(tree[where][2]>-1)confirm(idim,compareto,dir,tree[where][2]);
   }
   
}

void kd_tree::add(double *v){
  
  int i,j,k,l,**tbuff,node,dir;
  double **dbuff;
  
  if(pts==room){
     dbuff=new double*[pts];
     for(i=0;i<pts;i++){
       dbuff[i]=new double[dim];
       for(j=0;j<dim;j++){
        dbuff[i][j]=data[i][j];
       }
       delete [] data[i];
     }
     delete [] data;
     
     tbuff=new int*[pts];
     for(i=0;i<pts;i++){
       tbuff[i]=new int[4];
       for(j=0;j<4;j++){
         tbuff[i][j]=tree[i][j];
       }
       delete [] tree[i];
     }
     delete [] tree;
     
     room+=roomstep;
     data=new double*[room];
     tree=new int*[room];
     for(i=0;i<room;i++){
       data[i]=new double[dim];
       tree[i]=new int[4];
     }
     
     for(i=0;i<pts;i++){
        for(j=0;j<dim;j++){
	  data[i][j]=dbuff[i][j];
	}
	delete [] dbuff[i];
	for(j=0;j<4;j++){
	  tree[i][j]=tbuff[i][j];
	}
	delete [] tbuff[i];
     }
     delete [] dbuff;
     delete [] tbuff;
     
  }//if(pts==room) and we have to make more room
  
  node=find_node(v);
  for(j=0;j<dim;j++)data[pts][j]=v[j];
  
  if(v[tree[node][0]]<data[node][tree[node][0]])dir=1;
  else dir=2;
  
  tree[node][dir]=pts;
  tree[pts][3]=node;
  tree[pts][0]=tree[node][0]+1;
  if(tree[pts][0]>=dim)tree[pts][0]=0;
  tree[pts][1]=-1;
  tree[pts][2]=-2;
  
  pts++;


}

int kd_tree::find_node(double *v){
   
    int i,j,k,l,nextstep,where;
    
    where=masterparent;
    
    if(v[tree[masterparent][0]]<data[masterparent][tree[masterparent][0]]){
      nextstep=tree[masterparent][1];
    }
    else nextstep=tree[masterparent][2];
    
    
    
    while(nextstep>-1){
      where=nextstep;   
      if(v[tree[where][0]]<data[where][tree[where][0]]){
        nextstep=tree[where][1];
      } 
      else nextstep=tree[where][2];
    }
    
    return where;
    
}

double kd_tree::distance(double *p1, double *p2){
  double dd;
  int i;
  
  dd=0.0;
  for(i=0;i<dim;i++)dd+=power((p1[i]-p2[i])/(maxs[i]-mins[i]),2);
  dd=sqrt(dd);
  return dd;
}

void kd_tree::neigh_check(double *v, int kk, int *neigh, double *dd,\
 int where, int wherefrom){

   int i,j,k,l,side;
   double dtry,dwhere;
   
   if(v[tree[where][0]]<data[where][tree[where][0]])side=1;
   else side=2;
   
   dtry=fabs((v[tree[where][0]]-data[where][tree[where][0]])/\
   (maxs[tree[where][0]]-mins[tree[where][0]]));
   
   if(dtry<=dd[kk-1]){
   
     dwhere=distance(data[where],v);
     if(dwhere<dd[kk-1]){
        for(i=kk-2;i>=0 && dd[i]>dwhere;i--){
           dd[i+1]=dd[i];
	   neigh[i+1]=neigh[i];
	}
	i++;
	
	
	dd[i]=dwhere;
	neigh[i]=where;
     }
     
     if(wherefrom==tree[where][3] || wherefrom==tree[where][side]){
      if(tree[where][3-side]>-1){
       neigh_check(v,kk,neigh,dd,tree[where][3-side],where);
       //check the other side of this node
      }
     }
     
   }
   
   if(wherefrom==tree[where][3]){
     if(tree[where][side]>-1){
       neigh_check(v,kk,neigh,dd,tree[where][side],where);
       //check the side of this node I am naturally on
     } 
   }
   else{
     if(tree[where][3]>-1){
       neigh_check(v,kk,neigh,dd,tree[where][3],where);
       //check the parent of this node, if that is not where I came from
     }
   }
    

}

void kd_tree::nn_srch(double *v, int kk, int *neigh, double *dd){

   int i,j,k,l,node,where,behind,*inn;
   double ddnode,ddtry;
  
   node=find_node(v);
   ddnode=distance(v,data[node]);
   dd[0]=ddnode;
   neigh[0]=node;
   
  // printf("node is %d dd %e\n",node,ddnode);
   
   for(i=1;i<kk;i++){
     dd[i]=dd[0]+double(i)*1.0e6;
     neigh[i]=-1;
   }
   
   /*j=1;
   for(i=0;j<kk;i++){
     
     l=1;
     for(k=0;k<j;k++){
       if(neigh[k]==i)l=0;
     }
     if(l==1){
       dd[j]=distance(data[i],v);
       neigh[j]=i;
       j++;
     }
     
   }*/
   
   
  // sort(dd,neigh,kk);
   
   
   
   neigh_check(v,kk,neigh,dd,tree[node][3],node);
  
}


likelihoodinterface.cpp0000644000437200043720000011214111744047501015733 0ustar  danielsfdanielsf#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#ifdef OPENMP
#include <omp.h>
#else
#define omp_get_num_threads() 1
#define omp_get_thread_num() 0
#endif

#include <time.h>
#include "likelihoodinterface.h"
#include "eigen_wrapper.h"

int ipower(int arg, int ee){
  int ans,i;
  
  ans=1;
  for(i=0;i<ee;i++)ans=ans*arg;
  
  return ans;
}

likelihood::~likelihood(){
  int i;
  double **dd,*d;

  //Call subroutines with delswit=-1 so that they delete
  //pointers that were allocated on their initial call  
  predict(dd,d,d,d,nparams,i,d,d,-1);
  sample_pts(-1);
 
  delete [] fpredicted;
  for(i=0;i<threads;i++)delete kptr[i];
  delete kptr;
  delete dice;
  delete [] mins;
  delete [] maxs;
  delete [] chisq;
  delete [] sigma;
  delete [] ddmin;
  delete [] ddavg;
  delete [] ddmax;
  delete [] suboptime;
  delete [] suboptpct;
  delete [] CLmax;
  delete [] CLmin;  
  delete [] kpa;
}

likelihood::likelihood(){
 
 printf("WARNING called the likelihood default constructor\n");
  printf("you shouldn't do that...");
  printf("I'm going to freeze the program now\n");
  scanf("%lf",&junk);
}

likelihood::likelihood(int nn,double *mns,double *mxs){

 //This routine doesn't actually do much.
 //It just builds some of the necessary matrices
 //and assigns the values you want to some of the
 //class member variables.  The real action of
 //beginning to explore parameter space is in
 //likelihoodinterface.initiailize()

 //nn is the number of parameters
 
 //*mns and *mxs are the min and max values allowed in parameter space
 
 int i,j,k,l;
 double rn;
 
 printf("starting\n");
 
 roomstep=100000; //When the arrays storing all of the known points
 		//need to be expanded, this is the amount they are expanded by

 lingerswitch=0; //default is that lingering will not happen
 lingermax=50;

 writevery=100; //call write_pts() every 100 sampled points
 
 nprinted=0; //keep track of how many points have been written with
 		//write_pts() so that you do not write the same point
		//twice.
 maxerr=-1.0;

 precision=0.001;//how close to the confidence limit is close enough

 spentlingering=0;
 lastsubop=0;

 printf("in likelihood constructor\n");
 #pragma omp parallel
 threads=omp_get_num_threads();
 printf("threads is %d\n",threads);

 
 i=int(time(NULL));
 
 printf("\n\nran seed is %d\n\n",i);

 dice=new Ran(i);
 npts=10000; //Number of random starting point sampled.
 	//This is a public parameter, so you can
	//set it to a different value before calling
	//likelihood.initialize()
 
 kk=15; //Number of nearest neighbors in the Gaussian Process 
 	//(also can be reset)

 nsamples=1000; //Number of samples to consider when using the Gaussian process

 target=1280.0; //Default value of chi squared for which to look
 chimin=1197.0;

 nparams=nn;
 kpmin=0.00001; //Value that kriging parameter gets set to in case of lingering

 CLmin=new double[nparams];//These will store the current upper and lower 
 CLmax=new double[nparams];//confidence limit bounds on parameters
 
 for(i=0;i<nparams;i++){
   CLmin[i]=1.0e10;
   CLmax[i]=-1.0e10;
 }
 

 kriging_parameter=1.0;
 rkp=1.0;

 mins=new double[nparams];
 maxs=new double[nparams];

 for(i=0;i<nparams;i++){
  mins[i]=mns[i];
  maxs[i]=mxs[i];
 }
 
 kptr=new kd_tree*[threads];
 //In order to allow parallelization without threads 
 //interfering with each other,
 //each thread needs its own kdtree storing the known points for the Gaussian
 //process.
 
}

void likelihood::initialize(double **guesses, int nguess){

//This routine will actually generate the initial set of random samples,
//evaluate their chi squared values,
//and store them in the kdtrees for the Gaussian process

 int i,j,k,l,*inn;
 double rn,tp,cc;
 double **base;
 
 //These are all parameters related to the timing statistics
 //recorded in timingfile.sav
 krigct=0;
 addct=0;
 krigtimewall=0.0;
 addtimewall=0.0;
 krigtimecpu=0.0;
 addtimecpu=0.0;
 
 
 ngood=0;//Number of points found within the confidence limit
 
 startpts=npts;//Number of points you are starting with
 
 room=npts;//This will be the initial size of the arrays keeping track
 	//of things like chi squared for all sampled points.
	//This is what gets incremented by roomstep when necessary.
	//It is not equal to npts to cut down on the number of times
	//we have to call delete [] and new []
 
 suboptimal=0;//Number of points chosen because of large 
 	     //sigma rather than small |target-mu|
 
 suboptct=0;//Number of total points chosen not counting lingering.
 	//suboptimal/suboptct is the percentage of points chosen for 
	//large sigma
		    
 fpredicted=new double[npts]; //This stores the value of chi squared
 			//predicted for each evaluated point, so
			//you can check if the Gaussian Process/Kriging
			//code is working properly
			
 sigma=new double[npts]; //This is the sqrt(variance) of each point predicted
 			//by the Gaussian Process/Kriging code

 chisq=new double[npts];
 
 ddmin=new double[npts];//These store the min, max, and average distances
 ddavg=new double[npts];//between evaluated candidate points and their kk
 ddmax=new double[npts];//nearest neighbors.  This information allows us to
 			//verify that the Gaussian process gets more accurate 
			//as ddavg gets smaller
 
 suboptime=new int[npts]; //this allows us to keep track of the number of 
 			//points chosen suboptimally
 			//as a function of number of points chosen 
			//(i.e. a time series)
			
 suboptpct=new double[npts]; //same as above for percentage of 
 				//suboptimal points

 kpa=new double[npts]; //this will store kriging parameter as a 
 			//function of time


 printf("initializing likelihood object\n");
 
 
 //below is the code to select the random base of points to start with
 
 base=new double*[npts];
 for(i=0;i<npts;i++){
  base[i]=new double[nparams];
 }
 
 printf("declared base npts %d params %d\n",npts,nparams);
 
 //first read in the supplied guesses
 for(j=0;j<nguess;j++){
   for(i=0;i<nparams;i++)base[j][i]=guesses[j][i];
 }
  
 for(;j<npts;j++){
  for(i=0;i<nparams;i++){
   rn=dice->doub();
   base[j][i]=(maxs[i]-mins[i])*rn+mins[i];
  }
  
 }
 
 
 printf("made base\n");

 
 inn=new int[npts];

   
 for(i=0;i<threads;i++){
     printf("about to make tree %d %d\n",nparams,npts);
     kptr[i]=new kd_tree(nparams,npts,base,inn,mins,maxs);
     kptr[i]->check_tree(-1);
     printf("diag %d\n",kptr[i]->diagnostic);

 }
 delete [] inn;
 
 //after this time, the list of points sampled from parameter space
 //will be stored in the arrays kptr[i]->data[][]
 
 printf("threads %d\n",threads);
 for(i=0;i<threads;i++){
  kptr[i]->check_tree(-1);
  printf("tree %d diag %d\n",i,kptr[i]->diagnostic);
 }
 

 printf("done with tree\n");
 //now calculate chi squared for each of those points
 for(j=0;j<npts;j++){
  cc=call_likelihood(kptr[0]->data[j]);
  printf("called on j %d\n",j);
  
  kpa[j]=0.0;
  fpredicted[j]=0.0;
  sigma[j]=0.0;
  chisq[j]=cc;
  ddmin[j]=0.0;
  ddmax[j]=0.0;
  ddavg[j]=0.0;
  suboptime[j]=0;
  suboptpct[j]=0.0;
 
 }
 
 printf("assigned f, fpred, chisq\n");
 
 for(i=0;i<npts;i++)delete [] base[i];
 delete [] base;
 
 printf("done with initializer\n");
 }

void likelihood::resume(char *inname){

//if your run gets interrupted, call this routine instead of initialize

//inname will be the name of the last output the code wrote
//the routine will read that file and then you will be set to pick up
//where you left off

 int i,j,k,l,*inn,th,*tsubopt,nn[2],oldsubopt=0;
 double rn,tp,cc,tkp;
 char injunk[letters],word[letters];
 double *tsig,*tfp,*tchi,*tdmin,*tdmax,*tdavg,tr,number;
 double *tsubpct,*tkpa;
 double **base,dd[2];
 FILE *input;

 krigct=0;
 addct=0;
 krigtimewall=0.0;
 addtimewall=0.0;
 krigtimecpu=0.0;
 addtimecpu=0.0;

 suboptimal=0;
 suboptct=0;
 
 printf("about to read %s\n",inname);
 npts=0;
 
 //open the file and read through it to see how many points it contains
 input=fopen(inname,"r");
 while(fscanf(input,"%s %le ",injunk,&junk)==2){
  for(i=1;i<nparams;i++)fscanf(input,"%s %le ",word,&junk);
  while(compare_char(word,"kp")==0){
   //for this reason 'kp' should always be the last entry in the 
   //output file
   
    fscanf(input,"%s %le",word,&junk);
  }
  
  npts++;
 }
 fclose(input);
 
 startpts=npts;
 room=npts;
 
 printf("resuming with %d points already sampled from %s\n",\
 npts,inname);

 //These are temporary storage for things like chi squared, mu, and sigma.
 tsig=new double[npts];
 tfp=new double[npts];
 tchi=new double[npts];
 tdmin=new double[npts];
 tdmax=new double[npts];
 tsubopt=new int[npts];
 tsubpct=new double[npts];
 tdavg=new double[npts];
 tkpa=new double[npts];
		    
 fpredicted=new double[npts]; 
 kpa=new double[npts];
			
 sigma=new double[npts];

 chisq=new double[npts];
 
 base=new double*[npts]; 

 ddmin=new double[npts];
 ddmax=new double[npts];
 ddavg=new double[npts];
 suboptime=new int[npts];
 suboptpct=new double[npts];
			
 for(i=0;i<npts;i++){
  base[i]=new double[nparams];
 }
 
 inn=new int[npts];
 
 //now actually read the file
 input=fopen(inname,"r");
 for(i=0;fscanf(input,"%s %lf",injunk,&base[i][0])==2;i++){
  for(j=1;j<nparams;j++)fscanf(input,"%s %lf",word,&base[i][j]);
  
  while(compare_char(word,"kp")==0){
    fscanf(input,"%s %lf",word,&number);
    
    if(compare_char(word,"sig")==1)tsig[i]=number;
    else if(compare_char(word,"chisq")==1){tchi[i]=number;
      if(tchi[i]<0.01)tchi[i]=10.0*target;
    }
    else if(compare_char(word,"predicted")==1)tfp[i]=number;
    else if(compare_char(word,"ddmin")==1)tdmin[i]=number;
    else if(compare_char(word,"ddmax")==1)tdmax[i]=number;
    else if(compare_char(word,"ddavg")==1)tdavg[i]=number;
    else if(compare_char(word,"suboptimal")==1){tsubopt[i]=int(number);
      
      suboptimal=tsubopt[i];
      
      /*if(tsubopt[i]>oldsubopt)suboptimal++;
      oldsubopt=tsubopt[i];
      tsubopt[i]=suboptimal;*/
      
    }
    else if(compare_char(word,"suboptpct")==1)tsubpct[i]=number;
  }
  tkp=number;
  if(tkp>0.0)tsig[i]=tsig[i]/sqrt(tkp);
  else tsig[i]=tsig[i];
      //this is because the root_out.sav file stores sigma accounting
      //for the Kriging parameter; we want the array sigma[] to store
      //a generic sigma that can be rescaled for different Kriging
      //parameters
  
  tkpa[i]=tkp;
  
  //inn[i]=i;
  junk=dice->doub(); //just in case you started with the same seed
 }
 fclose(input);
 
 kriging_parameter=number;//kriging parameter will get set to the last value
 rkp=number;		//stored in your old run
			//there is code at the end of this routine to deal with the
			//possibility that the old run was lingering (kriging parameter = very small)
			//when it got interrupted

 
 suboptct=int(double(suboptimal)/tsubpct[npts-1]);
 
 //we will now build a kd tree in which we will store the evaluated points
 //we do this so that we have access to a fast nearest neighbor
 //search each time Kriging asks for the 15 nearest neighbors of
 //the candidate test point
 

 
 
 for(th=0;th<threads;th++){

   //inn will allows you to associate the order of chi squared, mu, etc. 
   //as it is read in
   //with the order of the points after they are arranged into the kdtree
   
   //(this is relevant for an older version of kd_tree() that rearranged
   //its inputs; I am just too lazy to rewrite the code, even though kd_tree
   //no longer behaves that way)
   
   for(j=0;j<npts;j++)inn[j]=j;
   
   kptr[th]=new kd_tree(nparams,npts,base,inn,mins,maxs);
   

   if(th==0){
   //because we only need one copy of these arrays, only do it for the zeroth 
   //thread
   
    ngood=0;
    for(i=0;i<npts;i++){
     fpredicted[i]=tfp[inn[i]];
     sigma[i]=tsig[inn[i]];
     chisq[i]=tchi[inn[i]];
     ddmax[i]=tdmax[inn[i]];
     ddmin[i]=tdmin[inn[i]];
     ddavg[i]=tdavg[inn[i]];
     suboptime[i]=tsubopt[inn[i]];
     suboptpct[i]=tsubpct[inn[i]];
     kpa[i]=tkpa[inn[i]];
     
     if(chisq[i]<=target+precision){ngood++;
       for(j=0;j<nparams;j++){
         if(kptr[th]->data[i][j]<CLmin[j])CLmin[j]=kptr[th]->data[i][j];
	 if(kptr[th]->data[i][j]>CLmax[j])CLmax[j]=kptr[th]->data[i][j];
       }
     
     }
    }
   }
   printf("wrote\n");
    

 }
 
 
 delete [] tfp;
 delete [] tchi;
 delete [] tsig;
 delete [] inn;
 for(i=0;i<npts;i++)delete [] base[i];
 delete [] base;
 delete [] tdmin;
 delete [] tdmax;
 delete [] tdavg;
 delete [] tsubopt;
 delete [] tsubpct;
 delete [] tkpa;
 
 
 
 //this is the code to set the kriging parameter 
 //in the case where the old run
 //was lingering when it got interrupted
 if(kriging_parameter<10.0){
   if(npts>1100)i=npts-1000;
   else i=npts/2;
   set_kp(i);
 }
 
 if(compare_char(inname,masteroutname)==1){
 nprinted=npts;
 }
 else{
  nprinted=0;
  write_pts();
 }
 printf("done resuming\n");
 
}


//USER MODIFIED FUNCTION
double likelihood::call_likelihood(double *v){
  
  //this is the function that calls the likelihood function
  //to evaluate chi squared
  
  //*v is the point in parameter space you have chosen to evaluate
  
  //as written, it calls the CAMB likelihood function for the CMB anisotropy 
  //spectrum
  
  int i,start;
  double params[14],chisquared,omm;
  double d1,d2,sncc;
  double cltt[3000],clte[3000],clee[3000],clbb[3000],el[3000];
  FILE *output;
  
  //code below this comment is a cartoon likelihood function
  
  
   d1=v[0]+5.0;
   d1=d1*1.9*pi/20.0;
  
    chisquared=1197.0+500.0*log(1.0+sin(d1)*sin(d1));
    for(i=1;i<nparams;i++)chisquared+=300.0*log(fabs(v[i])+1.0);
  
  
  //code below this comment interfaces with CAMB and the WMAP7 likelihood
  //function
  
 
 /*
  
  while((v[0]+v[1])/(v[2]*v[2])>1.0){
    v[0]=0.9*v[0];
    v[1]=0.9*v[1];
    v[2]=1.1*v[2];
    //in the event that total omega_matter>1
    //printf("had to futz\n");
  }
  
  for(i=0;i<6;i++)params[i]=v[i];
 
  for(i=0;i<3000;i++)el[i]=0;
  params[2]=100.0*params[2];
  params[5]=exp(params[5])*1.0e-10;
  

  camb_wrap_(params,el,cltt,clte,clee,clbb); //a function to call CAMB
  
  for(start=0;el[start]<1.0;start++);

  wmaplikeness_(&cltt[start],&clte[start],&clee[start],\
  &clbb[start],&chisquared); //a function to call the WMAP likelihood code

 //printf("done with likelihood\n");

  if(chisquared<0.01)chisquared=10.0*target; 
  			//in case the model was so pathological that the
			//likelihood code crashed and returned
			//chisquared=0  (this has been known to happen)
  

  */
  
  return chisquared;
  
}

void likelihood::add_pt(double *v, double *dd, double fin, double sin, \
double chitrue){
//once the code has chosen the best point for evaluation and called
//call_likelihood on it, this routine will add that point to the 
//kd tree of evaluated points for future use in the Gaussian process

//*v is the point in parameter space that you are adding
//*dd stores the distance between it and its nearest neighbors
//fin stores the value of mu predicted by the Gaussian process
//sin stores the value of sigma returned by the Gaussian process
//chitrue is the value of chi squared returned by call_likelihood

 double *fpbuff,*sbuff,*cbuff,cc;
 double *ddnb,*ddxb,*ddab,*subpctbuff,*kpbuff;
 int *subbuff;
 int i,j,k,l,th,*inn;
 
 //printf("in adding\n");

if(npts==room){
//you have evaluated so many points that you need to make more room in
//the chi squared, mu, etc. arrays

 kpbuff=new double[npts];
 fpbuff=new double[npts];
 sbuff=new double[npts];
 cbuff=new double[npts];
 ddnb=new double[npts];
 ddxb=new double[npts];
 ddab=new double[npts];
 subbuff=new int[npts];
 subpctbuff=new double[npts];
 
 #pragma omp parallel for
 for(i=0;i<npts;i++){
 
  kpbuff[i]=kpa[i];
  fpbuff[i]=fpredicted[i];
  sbuff[i]=sigma[i];
  cbuff[i]=chisq[i];
  ddnb[i]=ddmin[i];
  ddxb[i]=ddmax[i];
  ddab[i]=ddavg[i];
  subbuff[i]=suboptime[i];
  subpctbuff[i]=suboptpct[i];
 }

 delete [] kpa;
 delete [] fpredicted;
 delete [] sigma;
 delete [] chisq;
 delete [] ddmin;
 delete [] ddmax;
 delete [] ddavg;
 delete [] suboptime;
 delete [] suboptpct;
 
 #pragma omp parallel for
 for(th=0;th<threads;th++){
  kptr[th]->add(v);
  //the kdtree code has its own infrastucture for dealing with
  //the size of the tree and when you need to make more room
 }

 npts=kptr[0]->pts;

 room+=roomstep;

 kpa=new double[room];
 fpredicted=new double[room];
 sigma=new double[room];
 chisq=new double[room];
 ddmin=new double[room];
 ddmax=new double[room];
 ddavg=new double[room];
 suboptime=new int[room];
 suboptpct=new double[room];

 #pragma omp parallel for
 for(i=0;i<npts-1;i++){
  
  kpa[i]=kpbuff[i];
  fpredicted[i]=fpbuff[i];
  sigma[i]=sbuff[i];
  chisq[i]=cbuff[i];
  ddmin[i]=ddnb[i];
  ddmax[i]=ddxb[i];
  ddavg[i]=ddab[i];
  suboptime[i]=subbuff[i];
  suboptpct[i]=subpctbuff[i];
 }
 
 kpa[npts-1]=kriging_parameter;
 fpredicted[npts-1]=fin;
 sigma[npts-1]=sin;
 
 suboptime[npts-1]=suboptimal;
 suboptpct[npts-1]=double(suboptimal)/double(suboptct);
 chisq[npts-1]=chitrue;
 
 if(chisq[npts-1]<=target+precision){
   ngood++;
   for(i=0;i<nparams;i++){
     if(v[i]<CLmin[i])CLmin[i]=v[i];
     if(v[i]>CLmax[i])CLmax[i]=v[i];
   }
 }
 
 ddmin[npts-1]=dd[0];
 ddmax[npts-1]=dd[0];
 ddavg[npts-1]=0.0;
 for(i=0;i<kk;i++){
  if(dd[i]<ddmin[npts-1])ddmin[npts-1]=dd[i];
  if(dd[i]>ddmax[npts-1])ddmax[npts-1]=dd[i];
  ddavg[npts-1]+=dd[i];
 }
 ddavg[npts-1]=ddavg[npts-1]/double(kk);
 
// printf("assigned new f\n");

 delete [] kpbuff;
 delete [] fpbuff;
 delete [] sbuff;
 delete [] cbuff;
 delete [] ddnb;
 delete [] ddxb;
 delete [] ddab;
 delete [] subbuff;
 delete [] subpctbuff;
 }//if pts==room
 else{
 
 //printf("plenty of room\n");
 
 //if you already have room in the relevant arrays
#pragma omp parallel for
   for(th=0;th<threads;th++){
    kptr[th]->add(v);

   }
 //printf("k added\n");
 
   npts=kptr[0]->pts;
 
   kpa[npts-1]=kriging_parameter;
   fpredicted[npts-1]=fin;
   sigma[npts-1]=sin;
   chisq[npts-1]=chitrue;
   suboptime[npts-1]=suboptimal;
   suboptpct[npts-1]=double(suboptimal)/double(suboptct);
 
 if(chisq[npts-1]<=target+precision){
   
   ngood++;
   
   for(i=0;i<nparams;i++){
     if(v[i]<CLmin[i])CLmin[i]=v[i];
     if(v[i]>CLmax[i])CLmax[i]=v[i];
   } 
 }
 
  
   ddmin[npts-1]=dd[0];
   ddmax[npts-1]=dd[0];
   ddavg[npts-1]=0.0;
   for(i=0;i<kk;i++){
    if(dd[i]<ddmin[npts-1])ddmin[npts-1]=dd[i];
    if(dd[i]>ddmax[npts-1])ddmax[npts-1]=dd[i];
    ddavg[npts-1]+=dd[i];
   }
   ddavg[npts-1]=ddavg[npts-1]/double(kk);

 
 }
 
 //printf("added %d lm %d\n",npts,lingermax);
 
 if(npts%(writevery)==0){
 //printf("\nwriting %d\n\n",npts);
   write_pts();
   if(suboptimal-lastsubop<10){
      rkp=10.0*rkp;
    }
   else if(suboptimal-lastsubop>9*writevery/10){
      rkp=0.5*rkp;
    }
   if(kriging_parameter>kpmin)kriging_parameter=rkp;
   lastsubop=suboptimal;
 }
 //printf("assessed writing\n");
 
}


void likelihood::sample_pts(int delswit){
//generates a random set of candidate points
//uses a Gaussian process to find which is the best for actual evaluation with
//call_likelihood() 

//adds the best scoring one to the sample of known data

//delswit>0 means that the routine will go ahead and sample points
//delswit<0 means that the routine will delete its allocated pointers
//and exit

 static int samplinit=0;
 int i,j,k,l,maxdex=0,lct;
 double *bestpt,mu,sig,strad,**vdummy;
 double beforewall,afterwall,chibar;
 double beforecpu,aftercpu,startchi;
 static int called=0;
 
 static int **neigh,*calledthread;
 static double **sample,**bestpts,*mufound,*sfound;
 static double **dd,***neighpts,**neighf,*stradmax;
 static double **gg,**ggin,*ggq,ldmufactor;
 
 static double **keptd,*bestkeptd,*closest,*dmu;
 double chitrue,dchi,dx,ldmu;
 
 int th,ddint;
 double stradfinal,sfinal,mufinal;

 FILE *output;
 
 vdummy=new double*[threads];
 for(i=0;i<threads;i++)vdummy[i]=new double[nparams];
 
 if(called==0){
   calledthread=new int[threads];
   gg=new double*[kk];
   ggin=new double*[kk];
   for(i=0;i<kk;i++){
     gg[i]=new double[kk];
     ggin[i]=new double[kk];
   }
   ggq=new double[kk];
   
   dmu=new double[nparams];
   closest=new double[threads];
   keptd=new double*[threads];
   bestkeptd=new double[kk];
   neigh=new int*[threads];
   dd=new double*[threads];
   neighpts=new double**[threads];
   neighf=new double*[threads];
   mufound=new double[threads];
   sfound=new double[threads];
   stradmax=new double[threads];
   sample=new double*[threads];
   bestpts=new double*[threads];

   for(th=0;th<threads;th++){
     keptd[th]=new double[kk];  
     neigh[th]=new int[kk];
     dd[th]=new double[kk];
     neighpts[th]=new double*[kk];
     for(i=0;i<kk;i++)neighpts[th][i]=new double[nparams];
     neighf[th]=new double[kk];
 
     bestpts[th]=new double[nparams];
     sample[th]=new double[nparams];
   }
  called=1;
 }//if called==0

 if(delswit>0){
 
 for(th=0;th<threads;th++){
   stradmax[th]=-1.0e20;
   closest[th]=1.0e20;
   calledthread[th]=0;
 }
 
 beforewall=double(time(NULL));
 beforecpu=double(clock());

 //printf("sampling again\n");

 #pragma omp parallel for private(j,k,l,strad,mu,sig,th,ddint)
 for(i=0;i<nsamples;i++){
 th=omp_get_thread_num();

  #pragma omp critical (flip_a_coin)
  for(j=0;j<nparams;j++){
   sample[th][j]=(maxs[j]-mins[j])*dice->doub()+mins[j];   
  }
  //the code above generates the random sample

  kptr[th]->nn_srch(sample[th],kk,neigh[th],dd[th]); 
                         //find the kk nearest neighbors
  			//the neighbors are stored in neigh[th] 
			//as a list of nodes on
			//the kdtree
  
  for(j=0;j<kk;j++){
   for(l=0;l<nparams;l++){
     neighpts[th][j][l]=kptr[th]->data[neigh[th][j]][l];  
     //here, there neighbors are stored as points in parameter space
   }
   neighf[th][j]=chisq[neigh[th][j]];
   //this stores the values of chi squared associated with the neighbors
  }

  //do the Kriging/Gaussian Process
  predict(neighpts[th],neighf[th],sample[th],dd[th],kk,nparams,&mu,&sig,1);

  if(fabs(mu-target)<closest[th])closest[th]=fabs(mu-target);
  //keeping track of the closest candidate point considered

  strad=1.96*sig*sqrt(kriging_parameter)-fabs((mu-target));
  //the Kriging parameter appears here because the subroutine predict()
  //yields the value of sigma if the kriging parameter were equal to unity
  
  if(!(sig>0.0)){
    //error message
    output=fopen("gp_error.sav","a");
    fprintf(output,"strad is %e mu %e sig %e npts %d\n",strad,mu,sig,npts);
    
    fprintf(output,"\n found nn\n");
   for(j=0;j<kk;j++){
    fprintf(output,"%d ",neigh[th][j]);
    for(l=0;l<nparams;l++)\
    fprintf(output,"(%.2e %.2e) ",\
    neighpts[th][j][l],kptr[th]->data[neigh[th][j]][l]);
    fprintf(output,"\n");
   }
   kptr[th]->check_tree(-1);
   fprintf(output,"tree check is %d\n",kptr[th]->diagnostic);
  
   fclose(output);

  }
  
  if(strad>stradmax[th] || calledthread[th]==0){
                         //each omp thread stores its own best point
   maxdex=i;		//they will compare notes after the calculation
   stradmax[th]=strad;  //is done
   mufound[th]=mu;
   sfound[th]=sig;
   calledthread[th]=1;
   for(l=0;l<nparams;l++)bestpts[th][l]=sample[th][l];
   for(l=0;l<kk;l++){
     keptd[th][l]=0.0;
     for(j=0;j<nparams;j++){
       keptd[th][l]+=power((sample[th][j]-neighpts[th][l][j])/(maxs[j]-mins[j]),2);
     }
     keptd[th][l]=sqrt(keptd[th][l]);
   } 
  }
  
 }


 //below is where the omp threads compare to see which one found the
 //absolute best straddle point
 bestpt=bestpts[0];
 mufinal=mufound[0];
 sfinal=sfound[0];
 stradfinal=stradmax[0];
 
 for(l=0;l<kk;l++){
  bestkeptd[l]=keptd[0][l];
 }
 
 
 for(th=1;th<threads;th++){
  if(stradmax[th]>stradfinal){
   stradfinal=stradmax[th];
   mufinal=mufound[th];
   sfinal=sfound[th];
   bestpt=bestpts[th];
  }
  for(l=0;l<kk;l++){
    bestkeptd[l]=keptd[th][l];
  }
 }
 
 //did a point get chosen based on sigma rather than
 //proximity to the target value?
 strad=closest[0];
 for(th=1;th<threads;th++){
  if(closest[th]<strad)strad=closest[th];
 }
 if(fabs(mufinal-target)>strad)suboptimal++; 
 
 afterwall=double(time(NULL));
 aftercpu=double(clock());
 if(samplinit>0){
   krigtimewall+=(afterwall-beforewall);
   krigtimecpu+=(aftercpu-beforecpu);
   krigct++;
 }

 beforewall=double(time(NULL));
 beforecpu=double(clock());
 chitrue=call_likelihood(bestpt);
 
 suboptct++; //how many points were selected while not lingering
 add_pt(bestpt,bestkeptd,mufinal,sfinal,chitrue);
 
// printf("out of adding\n");
 
 afterwall=double(time(NULL));
 aftercpu=double(clock());
 
 if(samplinit>0){
   addct++;
   addtimewall+=(afterwall-beforewall);
   addtimecpu+=(aftercpu-beforecpu);
 }
 else samplinit++;


 //below is the code to handle lingering around highly likely points
 
 if(lingerswitch==1 && fabs(chisq[npts-1]-target)<0.1*target \
 && chisq[npts-1]!=target){
    //printf("starting to linger %e %d %e\n",chisq[npts-1],lingermax,precision);
   
    kriging_parameter=kpmin;
    //this is a warning to output routines that these points were sampled
    //under lingering
       
   ldmufactor=0.001;
   
   startchi=chisq[npts-1];
   for(lct=0;lct<lingermax && fabs(startchi-chimin)>=0.4*precision\
   && ldmufactor>1.0e-7;lct++){
 
      for(j=0;j<nparams;j++)bestpts[0][j]=kptr[0]->data[npts-1][j];
     startchi=chisq[npts-1];
     
     kptr[0]->nn_srch(bestpts[0],kk,neigh[0],dd[0]); 
    
    chibar=0.0;                   
    for(j=0;j<kk;j++){
     for(l=0;l<nparams;l++){
       neighpts[0][j][l]=kptr[0]->data[neigh[0][j]][l];  
       //here, there neighbors are stored as points in parameter space
     }
     neighf[0][j]=chisq[neigh[0][j]];
     chibar+=neighf[0][j];
     //this stores the values of chi squared associated with the neighbors
    }
    chibar=chibar/double(kk);

   

    for(i=0;i<kk;i++){
      for(j=i;j<kk;j++){
        gg[i][j]=covariogram(neighpts[0][i],neighpts[0][j],nparams);
        if(i!=j)gg[j][i]=gg[i][j];
        else if(i==j)gg[i][j]=1.0001*gg[i][j];
      }
    }
  
    invert_lapack(gg,ggin,kk,1);
    for(i=0;i<kk;i++)ggq[i]=covariogram(neighpts[0][i],bestpts[0],nparams);
    
    //calculate the direction of the gradient of chi squared
    //using a Gaussian process
    for(k=0;k<nparams;k++){
      dmu[k]=0.0;
      for(i=0;i<kk;i++){
        for(j=0;j<kk;j++){
	 dmu[k]+=-1.0*(bestpts[0][k]-neighpts[0][j][k])\
	 *ggq[j]*ggin[j][i]*(neighf[0][i]-chibar)/(maxs[k]-mins[k]); 
	}
      }
      
    }
    
    chitrue=3.0*target;
    
     ldmu=0.0;
     for(j=0;j<nparams;j++)ldmu+=dmu[j]*dmu[j];
     ldmu=sqrt(ldmu);
    
    
    dchi=chimin-startchi;
   
    while(chitrue>2.0*target && ldmufactor>1.0e-7){
      //printf("chitrue %e ldmuf %e\n",chitrue,ldmufactor);
      dx=ldmufactor*dchi/ldmu;
    
      for(j=0;j<nparams;j++){
         sample[0][j]=bestpts[0][j]+dx*dmu[j];
      }
     
     j=1;
     for(k=0;k<nparams;k++){
      if(sample[0][k]<mins[k] || sample[0][k]>maxs[k]){j=0;
         //make sure the sample point is within the allowed bounds
	 //on parameter space
      }
     } 
     if(j==1){ 
      chitrue=call_likelihood(sample[0]);
      spentlingering++;
      
      for(k=0;k<kk;k++)dd[0][k]=0.0;//so that points found with
      				//gradient descent are easily identified
				//in output file
    
    add_pt(sample[0],dd[0],target,0.0,chitrue);
    
       if(chitrue>=2.0*target)ldmufactor=ldmufactor/2.0;
     
     }
     else{ ldmufactor=ldmufactor/2.0;}//this is if the sample point
     					//was outside the parameter
					//space bounds.
    }

  
  }
  //printf("done lingering %e %d %e\n\n",chisq[npts-1],lct,ldmufactor);
 

  kriging_parameter=rkp;
 
 
 }
 
 }//if delswit>0
 else if(called>0){
 
 delete [] calledthread;
 delete [] closest;
 delete [] mufound;
 delete [] sfound;
 delete [] stradmax;
 
 
 for(th=0;th<threads;th++){
   delete [] sample[th];
   delete [] bestpts[th];
   delete [] neigh[th];
   delete [] dd[th];
   for(i=0;i<kk;i++)delete [] neighpts[th][i];
   delete [] neighpts[th];
   delete [] neighf[th];
   
   delete [] keptd[th];
 
 }
 delete [] sample;
 delete [] bestpts;
 delete [] neigh;
 delete [] dd;
 delete [] neighpts;
 delete [] neighf;
 delete [] keptd;
 delete [] bestkeptd;
 delete [] dmu;
 for(i=0;i<kk;i++){
   delete [] gg[i];
   delete [] ggin[i];
 }
 delete [] gg;
 delete [] ggin;
 delete [] ggq;
 
   called=0;
 }

 for(i=0;i<threads;i++)delete [] vdummy[i];
 delete [] vdummy;

//ddfile fclose(output);
// printf("done sampling\n");
}

void likelihood::write_pts(){

//this routine will write the evaluated points along with
//*sigma, chi squared, *fpredicted, and the difference between
//chi squared and *fpredicted into the file *name

 FILE *output,*timefile,*goodfile;
 int i,j,th;
 int tot,trapped,good,tried,wayoff,dontcount;
 double trappedpct;
 
 
 for(i=0;i<threads;i++){
  //check to make the kdtrees are still properly constructed
  kptr[i]->check_tree(-1);
  //printf("tree %d\n",kptr[i]->diagnostic);
 }
 
 //printf("checked all trees\n");
 
 tot=0;
 trapped=0;

 
 
 good=0;
 tried=0;
 wayoff=0;
 dontcount=0;


  output=fopen(masteroutname,"a");

 tot=0;
 trapped=0;
 
 //goodfile=fopen("goodpts.sav","w");
 
 for(i=0;i<npts;i++){//printf("i %d np %d\n",i,nprinted);
  
  //this code determines what percentage of points (ignoring lingering)
  //were actually trapped within their 1-sigma kriging bounds 
  //note: this trapped percentage assumes that all points were sampled
  //with the current Kriging parameter
  //it is another way to judge whether or not our choice of Kriging parameter
  //was appropriate
   
     if(sigma[i]>0.0 && kpa[i]>kpmin){
      tot++;
      if(fabs(fpredicted[i]-chisq[i])\
      <=sqrt(kriging_parameter)*sigma[i])trapped++;
    }
   // printf("checked trapped and tot\n");
 
  if(tot>0){
   trappedpct=double(trapped)/double(tot);
  }
  else trappedpct=0.0;
  
  if(chisq[i]<target+precision){good++;

  }
  
  //printf("checked good\n");
  
  //below, ``tried'' keeps track of points that were chosen within 10% of target
  //``wayoff'' keeps track of points that were chosen that were more than 90%
  //different from target
  //This is another way to evaluate our choice of Kriging parameter
  
  if(fabs(fpredicted[i]-target)/target<0.1 && sigma[i]>0.0 \
  && kpa[i]>1.1*kpmin)\
  tried++;
  else if(fabs(fpredicted[i]-target)/target>0.9 && sigma[i]>0.0 && \
  kpa[i]>1.1*kpmin)wayoff++;
  else if(sigma[i]<=0.0 || kpa[i]<=1.1*kpmin)dontcount++;
 
  if(i>=nprinted){//no need to print points that have already been printed
  
  for(j=0;j<nparams;j++){
   fprintf(output,"%s %e ",pnames[j],kptr[0]->data[i][j]);
  }
  fprintf(output,"chisq %e ",chisq[i]);
  fprintf(output,"predicted %e sig %e ",\
  fpredicted[i],sqrt(kpa[i])*sigma[i]);
  fprintf(output,"diffpct %e ",fabs((chisq[i]-fpredicted[i])/chisq[i]));
  fprintf(output,"ddmin %e ddmax %e ddavg %e ",ddmin[i],ddmax[i],ddavg[i]);
  fprintf(output,\
  "goodpct %e triedpct %e wayoff %e suboptimal %d suboptpct %e trappedpct %e ",\
  double(good)/double(i+1),double(tried)/double(i+2-dontcount),\
  double(wayoff)/double(i+2-dontcount),\
  suboptime[i],suboptpct[i],trappedpct);
 
  fprintf(output,"kp %e\n",kpa[i]);
  }//if(i>=nprinted)
  
 }
 fclose(output);
 //fclose(goodfile);

 
 timefile=fopen("timingfile.sav","a");
 fprintf(timefile,"%s %d ling %d gd %d ",masteroutname\
 ,npts,spentlingering,ngood);
  fprintf(timefile,"avgkrigtimewall %e avgaddtimewall %e ns %d kk %d ",\
  krigtimewall/double(krigct),addtimewall/double(addct),nsamples,kk);
     fprintf(timefile,"avgkrigtimecpu %e avgaddtimecpu %e ",\
     krigtimecpu/double(krigct),addtimecpu/double(addct));
 fprintf(timefile,"suboptpct %e ct %d kp %e ",suboptpct[npts-1],\
 suboptct,kriging_parameter);
 for(th=0;th<threads;th++)fprintf(timefile,"kd%d %d ",th,kptr[th]->diagnostic);
 //for(th=0;th<nparams;th++)printf("%s %e %e ",pnames[th],CLmin[th],CLmax[th]);
 fprintf(timefile,"target %e\n",target);
 fclose(timefile);
 
 
 
  //reset variables related to timing performance
   krigct=0;
   addct=0;
   addtimecpu=0.0;
   addtimewall=0.0;
   krigtimewall=0.0;
   krigtimecpu=0.0;
   
   nprinted=npts;


}

double likelihood::covariogram(double *v1, double *v2, int coords){

 int i;
 double ans,d;
 
 d=0.0;
 for(i=0;i<coords;i++){
  d+=power((v1[i]-v2[i])/(maxs[i]-mins[i]),2);
 }
 ans=exp(-0.5*d);
 return ans;
 
}



void likelihood::predict(double **old, double *ffn, double *q, double *dd, \
int pts, int coords, double *mu, double *sig, int delswit){

//this subroutine actually does the Kriging/Gaussian Process

//**old stores the coordinates of the kk nearest neighbors in parameter space

//*ffn is the chisquared values associated with the neighbors

//*q is the test point being considered

//*dd stores the distances to the kk nearest neighbors

//pts is really just kk (the number of points in **old and *ffn)

//coords is the dimensionality of parameter space (i.e., nparams)

//*mu is where the predicted value of chi squared at *q gets stored

//*sig is the uncertainty on *mu

//delswit is an integer; delswit>0 means this routine will do Kriging
//delswit<0 means that it will delete its allocated pointers and exit

//delswit==2 causes the routine to print updates about where it is in the
//process (this is probably not very useful)

  double fbar,err,junk,d;
  int i,j,k,l;
  double static **gg,**ggin,*ggq;
  static int called=0;
  
  if(called==0){
    gg=new double*[pts];
  
    ggin=new double*[pts];
    for(i=0;i<pts;i++){
     gg[i]=new double[pts];
     ggin[i]=new double[pts];
  
   }
   ggq=new double[pts];
   called=1;
  }
  
  if(delswit<0 && called>0){
    delete [] ggq;
    for(i=0;i<pts;i++){
      delete [] gg[i];
      delete [] ggin[i];
    }
    
    delete [] gg;
    delete [] ggin;
    called=0;
  }
  else{
  
  for(i=0;i<pts;i++){
   for(j=i;j<pts;j++){
     gg[i][j]=covariogram(old[i],old[j],coords);
     if(i!=j)gg[j][i]=gg[i][j];
     else if(i==j)gg[i][j]=gg[i][j]*1.001;
   }
  }
  
  if(delswit==2){printf("made gg\n\n");
  
  }
  invert_lapack(gg,ggin,pts,delswit);
  
  if(delswit==2)printf("inverted gg\n");
  
  //the comments below are code to check that the matrix inversion worked;
  //you can generally trust lapack
  
  //err=check_inversion(gg,ggin,pts);
  //if(err>maxerr)maxerr=err;

  for(i=0;i<pts;i++)ggq[i]=covariogram(old[i],q,coords);

  fbar=0.0;
  for(i=0;i<pts;i++)fbar+=ffn[i];
  fbar=fbar/double(pts);
  
  *mu=fbar;
  
  for(i=0;i<pts;i++){
   for(j=0;j<pts;j++){
    *mu+=ggq[i]*ggin[i][j]*(ffn[j]-fbar);
   }
  }
    
  *sig=0.0;
  for(i=0;i<pts;i++){
   for(j=0;j<pts;j++){
     *sig+=ggq[i]*ggq[j]*ggin[i][j];

   }
  }
  
  *sig=covariogram(q,q,coords)-*sig;
      d=*sig;
     *sig=sqrt(fabs(*sig));

  }

}


double likelihood::test_kp(int pts, char *name){

//this is a subroutine to test how good a choice your Kriging parameter
//actually is.  It will generate 'pts' random points, run them through
//Kriging, and tell you what percentage of them actually fall within
//the 1-sigma range of chi squared.  The results are stored in the
//file specified by *name


   double err,cc;
   double **kx,*ky,*kmu,*ksig;
   double ***neighpts,**dd,**neighf;
   int i,j,k,l,th,**neigh;
   int trapped;
   FILE *output;


   trapped=0;
  

  neigh=new int*[threads];
  dd=new double*[pts];
  neighpts=new double**[threads];
  neighf=new double*[threads];

  for(th=0;th<threads;th++){
   neigh[th]=new int[kk];
   neighpts[th]=new double*[kk];
   for(i=0;i<kk;i++)neighpts[th][i]=new double[nparams];
   neighf[th]=new double[kk];
  }
  
  kx=new double*[pts];
  for(i=0;i<pts;i++){
     kx[i]=new double[nparams];
     dd[i]=new double[kk];
  }
  ky=new double[pts];
  kmu=new double[pts];
  ksig=new double[pts];

#pragma omp parallel for private(i,j,k,l,th)
  for(i=0;i<pts;i++){
    th=omp_get_thread_num();
    
#pragma omp critical (flip_a_coin_kp)
    for(j=0;j<nparams;j++){
     kx[i][j]=(maxs[j]-mins[j])*dice->doub()+mins[j];
    }
  

   kptr[th]->nn_srch(kx[i],kk,neigh[th],dd[i]);
   
    
    for(j=0;j<kk;j++){
      for(l=0;l<nparams;l++){
        neighpts[th][j][l]=kptr[th]->data[neigh[th][j]][l];  
      }
      neighf[th][j]=chisq[neigh[th][j]];
    }
    predict(neighpts[th],neighf[th],kx[i],dd[i],kk,nparams,&ky[i],&ksig[i],1);
  
  }
  
  printf("ready to call like\n");
  //cannot parallelize this part because of how CAMB works
  for(i=0;i<pts;i++){
    kmu[i]=call_likelihood(kx[i]);

    if(fabs(kmu[i]-ky[i])<=sqrt(rkp)*ksig[i]){
     trapped++;
    }
    
    
    add_pt(kx[i],dd[i],ky[i],ksig[i],kmu[i]);
    //add the point to the tree so you still learn something from having
    //called the likelihood function

  }

 output=fopen(name,"a");
  fprintf(output,\
  "%s kptest pts %d trapped %d pct %e npts %d kp %e subopt %d\n",\
  masteroutname,pts,trapped,double(trapped)/double(pts),npts,rkp,suboptimal);
 fclose(output);


  for(i=0;i<pts;i++){
   delete [] kx[i];
   delete [] dd[i];
  }
  
  for(th=0;th<threads;th++){
   for(i=0;i<kk;i++)delete [] neighpts[th][i];
   delete [] neighpts[th];
   delete [] neigh[th];
   delete [] neighf[th];
  
  }
  
  delete [] dd;
  delete [] neighf;
  delete [] neighpts;
  delete [] neigh;
  delete [] kx;
  delete [] ky;
  delete [] kmu;
  delete [] ksig;
  
  return double(trapped)/double(pts);

}

void likelihood::set_kp(int ignore){
   
//This routine will set the Kriging parameter heuristically.
//It looks through the points already sampled (ignoring
//the first set of points, the length of which is specified
//by the integer 'ignore') and sets the Kriging parameter
//to that value such that 68% of the points sampled
//would have been in the 1-sigma interval predicted by the
//Gaussian process
   
   double *rat;
   int *sind,i,ct;
 
   
   rat=new double[npts];
   sind=new int[npts];
   
   ct=0;
   for(i=0;i<npts;i++){
    if(sigma[i]>0.0 && i>ignore){
     
     rat[ct]=power((chisq[i]-fpredicted[i])/sigma[i],2);
     sind[ct]=ct;
     
     ct++;
     
    }
   }
   
   sort(rat,sind,ct);
   kriging_parameter=rat[ct*68/100];
   rkp=kriging_parameter;
   
   
   delete [] rat;
   delete [] sind;
   
}



likelihoodinterface.h0000644000437200043720000000455711741610741015412 0ustar  danielsfdanielsf#include "goto_tools.h"
#include "kd.h"
#define likeletters 500

//the two external function definitions below are only needed to
//interface with CAMB and the WMAP 7 likelihood function
//as provided in aps_cmb_module.cpp


extern "C" void \
camb_wrap_(double*,double*,double*,double*,double*,double*);
extern "C" void wmaplikeness_(double*,double*,double*,double*,double*);


class likelihood{
 private:
  
  double *fpredicted,*ddmin,*ddmax,*ddavg,*chisq,*sigma;
     //f is the distance from the confidence ball center
     //fpredicted is what Kriging predicted for f
     //sigma is the variance Kriging predicted for f
     
     

   
  //CAMB and WMAP necessary arrays

  int room,roomstep,nprinted;
  
 
  
  double *ldata,*lx;
  int lpts;
    
  double *mins,*maxs,rkp,*suboptpct,*kpa;
    //limits on the cosmologicalparameters
    
  int nok,nbad,datapts,startpts; 
     //npts is number of sample points accumulated 
     //nparams is the number of theoretical parameters
     
  int nparams;
  
  Ran *dice;
  kd_tree **kptr;
  
 
 public:
  char **pnames,masteroutname[100];
   //cmd will be a command that calls an external program (the likelihood
   //code) which takes a point in parameter space and outputs a curve
   //to be compared to the data
   
   //pstore will contain the parameter values to be input to the
   //likelihood code
   
   //dname is the name of the file that will contain the
   //data produced by the likelihood code
   
   double junk,kriging_parameter,maxerr;
   double *ans1,*ans2,db,kpmin,chimin;
     
   int npts,nsamples,threads,ngood,suboptimal,suboptct,*suboptime;
   int krigct,addct,lingerswitch,lingermax,kk,spentlingering;
   double krigtimewall,addtimewall,*CLmin,*CLmax;
   double krigtimecpu,addtimecpu;
   double target,precision;
   int lastsubop,writevery;
   
   
  
   
   likelihood();
   likelihood(int,double*,double*);
   ~likelihood();
   void initialize(double**,int);
   void resume(char*);
   
   double call_likelihood(double*);
   void add_pt(double*,double*,double,double,double);
     //now you feed true chi and f to add_pt
     //so that the convergence files can keep track of diffpct
     //on kriging
   double covariogram(double*,double*,int);
   void predict(double**,double*,double*,double*,int,int,double*,double*,int);
   void sample_pts(int);
   void write_pts();
   void set_kp(int);
   double test_kp(int,char*);
   

   

};
aps_interface.cpp0000664000437200043720000001273511741653445014553 0ustar  danielsfdanielsf#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include "likelihoodinterface.h"

main(int iargc, char *argv[]){
   
   int i,j,k,l,ii,dim,nstart=1000,nkp=1000,nend=10000,nc=200,llswit=0;
   int resumeswitch=0;
   char paramname[letters],keyword[letters],word[letters];
   char resumename[letters];
   double nn,mm,*min,*max;
   
   double **guess,**buff;
   int nguess;
   
   nguess=0;
   
   likelihood *likeness;
   
   FILE *input;
   
   for(j=0;argv[1][j]!=0;j++){
     paramname[j]=argv[1][j];
   }
   paramname[j]=0;
   
   input=fopen(paramname,"r");
   while(fscanf(input,"%s",keyword)==1){
     if(compare_char(keyword,"#dim")==1){
       
       fscanf(input,"%d",&dim);
       
       min=new double[dim];
       max=new double[dim];
       
     }
     else{
       if(compare_char(keyword,"#param")==1){
         fscanf(input,"%d %s %le %le",&ii,word,&nn,&mm);
	 
	 if(nn>mm){
	   max[ii]=nn;
	   min[ii]=mm;
	 }
	 else{
	   max[ii]=mm;
	   min[ii]=nn;
	 }
	 
	 //printf("set lims %d %e %e\n",ii,min[ii],max[ii]);
	 
       }
     }
   }
   fclose(input);
   likeness=new likelihood(dim,min,max);
   
   likeness->pnames=new char*[dim];
   for(i=0;i<dim;i++)likeness->pnames[i]=new char[letters];
   for(i=0;i<dim;i++){
     sprintf(likeness->pnames[i],"param%d",i);
   }
   
   input=fopen(paramname,"r");
   while(fscanf(input,"%s",keyword)==1){
    if(compare_char(keyword,"#param")==1){
      fscanf(input,"%d %s %le %le",&ii,word,&nn,&mm);
      
      for(i=0;word[i]!=0;i++){
        likeness->pnames[ii][i]=word[i];
      }
      likeness->pnames[ii][i]=0;
      
      //printf("named %d %s\n",ii,likeness->pnames[ii]);
      
    }
    else if(compare_char(keyword,"#Ns")==1){
      fscanf(input,"%d",&nstart);
      likeness->npts=nstart;
      
      //printf("set start to %d\n",likeness->npts);
      
    }
    else if(compare_char(keyword,"#Nk")==1){
      fscanf(input,"%d",&nkp);
      
      //printf("set nkp %d\n",nkp);
      
    }
    else if(compare_char(keyword,"#end")==1){
      fscanf(input,"%d",&nend);
      
      //printf("set end %d\n",nend);
      
    }
    else if(compare_char(keyword,"#Nc")==1){
      fscanf(input,"%d",&nc);
      
      //printf("set nc %d\n",nc);
      
    }
    else if(compare_char(keyword,"#linger")==1){
      llswit=1;
      fscanf(input,"%d",&ii);
      likeness->lingermax=ii;
      
      //printf("set lingermax %d\n",likeness->lingermax);
      
    }
    else if(compare_char(keyword,"#chitarget")==1){
      fscanf(input,"%le",&likeness->target);
      
      //printf("set target %e\n",likeness->target);
      
    }
    else if(compare_char(keyword,"#chimin")==1){
      fscanf(input,"%le",&likeness->chimin);
      
      //printf("set chimin %e\n",likeness->chimin);
      
    }
    else if(compare_char(keyword,"#Ng")==1){
      fscanf(input,"%d",&likeness->kk);
      
      //printf("set kk %d\n",likeness->kk);
      
    }
    else if(compare_char(keyword,"#resumefile")==1){
      resumeswitch=1;
      fscanf(input,"%s",word);
      for(i=0;word[i]!=0;i++){
        resumename[i]=word[i];
      }
      resumename[i]=0;
      
      //printf("resuming from %s\n",resumename);
      
    }
    else if(compare_char(keyword,"#write_every")==1){
      
      fscanf(input,"%d",&i);
      likeness->writevery=i;
     
    }
    else if(compare_char(keyword,"#outputfile")==1){
      fscanf(input,"%s",word);
      //printf("%s\n",word);
      for(i=0;word[i]!=0 && word[i]!=' ';i++){
        likeness->masteroutname[i]=word[i];
      }
      likeness->masteroutname[i]=0;
      
      //printf("set output %s\n",likeness->masteroutname);
    }
    else if(compare_char(keyword,"#guess")==1){
      
      nguess++;
      if(nguess==1){
        guess=new double*[1];
	guess[0]=new double[dim];
	for(i=0;i<dim;i++)fscanf(input,"%le",&guess[0][i]);
      
      }
      else{
        buff=new double*[nguess-1];
	for(i=0;i<nguess-1;i++){
	  buff[i]=new double[dim];
	  for(j=0;j<dim;j++){
	    buff[i][j]=guess[i][j];
	  }
	  delete [] guess[i];
	}
	delete [] guess;
	guess=new double*[nguess];
	for(i=0;i<nguess-1;i++){
	  guess[i]=new double[dim];
	  for(j=0;j<dim;j++)guess[i][j]=buff[i][j];
	  delete [] buff[i];
	}
	delete [] buff;
	guess[nguess-1]=new double[dim];
	for(i=0;i<dim;i++)fscanf(input,"%le",&guess[nguess-1][i]);
	
      }
      
    
    }
   }
   fclose(input);
   
   //printf("done with input\n");
   
   /*for(i=0;i<nguess;i++){
     printf("guess %d\n",i);
     for(j=0;j<dim;j++)printf("%.3e ",guess[i][j]);
     printf("\n");
   }*/
   
   if(resumeswitch==0){
     printf("about to initialize\n");
     likeness->initialize(guess,nguess);
     likeness->nsamples=1; //so that no Kriging is done when setting Kriging parameter
     for(i=0;i<nkp;i++)likeness->sample_pts(1);
     likeness->set_kp(nstart);
     printf("done with that\n");
   }
   else{
    likeness->resume(resumename);
    nstart=likeness->npts;
   }
   likeness->write_pts();

   likeness->lingerswitch=llswit;
   likeness->nsamples=nc;
   while(likeness->npts<nend){
     
     likeness->sample_pts(1);
   }
   
   likeness->write_pts();
   
   printf("after everything\n");
   printf("kk %d\n",likeness->kk);
   printf("ns %d\n",likeness->nsamples);
   printf("target %e\n",likeness->target);
   printf("chimin %e\n",likeness->chimin);
   printf("spent lingering %d\n",likeness->spentlingering);
   printf("lingermax %d\n",likeness->lingermax);
   for(i=0;i<dim;i++){
     printf("%s min %e max %e\n",likeness->pnames[i], \
     likeness->CLmin[i],likeness->CLmax[i]);
   }
   
}
params_CMB.sav0000664000437200043720000000244011744056424013710 0ustar  danielsfdanielsfA parameter file for the WMAP 7 likelihood function
Note: all keywords must be prefixed with a #

#dim 6
#linger 100 -- enable lingering; set the maximum number of iterations to 110

#Ns 1000 -- start with 1000 randomly sampled points
#Nk 1000 -- us an additional 1000 random points to set the Kriging parameter
#end 20000 -- #end after the algorithm has sampled 200000 points

#Nc 250 -- sample 250 candidates each time you use the Gaussian process
#Ng 15 -- use 15 neareast neighbors in the Gaussian process

#chitarget 1286.69 -- target chi squared
#chimin 1197

#write_every 300 -- write after sampling 300 points


#outputfile test_out_wmap_3.sav -- the name of the output file
!#resumefile test_out_wmap.sav -- if you wanted to resume old #output file 
			(this line is currently commented out, note the '!')

#param 0 ombh2 0.01 0.04 -- the density of baryons
#param 1 omch2 0.05 0.2 -- the density fo dark matter
#param 2 h 0.5 0.8 -- the hubble parameter in units of 100 km/s/Mpc
#param 3 tau 0.01 0.15 -- optical depth to last scattering
#param 4 0.8 1.1 -- scalar index controlling primordial power spectrum
#param 5 2.5 3.5 -- 10^10 ln(A) where A is normalizatin of primordial power spec.


below is a guess roughly approximating the known best fit point

#guess
0.02273 0.1099 0.72 0.087 0.963 3.1

params.sav0000664000437200043720000000202511741653261013225 0ustar  danielsfdanielsfA parameter file for the cartoon likelihood function
Note: all key words must be prefixed with a #


#dim 2
#linger 30 -- enable lingering; set the maximum number of iterations to 30

#Ns 100 -- start with 100 randomly sampled points
#Nk 100 -- us an additional 100 random points to set the Kriging parameter
#end 2000 -- end after the algorithm has sampled 2000 points

#Nc 45 -- sample 45 candidates each time you use the Gaussian process
#Ng 7 -- use 7 neareast neighbors in the Gaussian process

#chitarget 1300.0 -- target chi squared
#chimin 1000.0

#write_every 30 -- write after sampling 30 points

#outputfile test_out_2.sav -- the name of the output file
!#resumefile test_out.sav -- if you wanted to resume an old output file 
			(this line is currently commented out, note the '!')


#param 0 alice -10 10  -- explore two parameters named 'alice' and 'bob'
#param 1 bob -10 10       allow each to vary between -10 and 10


below are a series of guesses that you want to make sure the code samples

#guess
1 2

#guess
4 6

#guess 3 7