Module OptDesign

Set up the quantities necessary to run the optimal design
algorithms. NOTE: This function needs to be called first 
before any of the optimal design functions can be called.

paramfile: the name of a pickled file containing the 
best fit parameters in KeyedList format
calcobject: the calculation object for which we are doing the 
optimal design. (Note that in general, may be searching a 
design over many different calculations, but here we only 
consider one. Thus, we set design_sentraj equal to senstraj)
senstrajfile: the name of the file containing the pickled 
sensitivity trajectory for the calculation, calcobject, 
for the set of parameters in paramfile.
jtjfile: the name of the file containing the pickled Fisher
Information Matrix (J^t J) for the current set of data and
for the parameters in paramfile. 
NOTE: The derivatives computed for J^tJ need to be with respect
to the *log* of the parameters

Make the sensitivity trajectory for the calculation
calcoject (same as in setup(...) above). 
params: parameters as a KeyedList, sensitivity traj is 
calculated at these parameters (should be same as in paramfile 
in setup(...) above)
times: the timepoints in the sensitivity trajectory (1-d array)
senstrajfilename: the file to save the sensitivity trajectory to

Note that if times is very finely spaced, the 
sensitivity trajectory will need a lot of storage space 

chemnames = list of unmeasurable chemicals
designchemnames = list of measurable chemicals
logprior = prior on params, e.g. logprior = log(1000.0) means
parameter standard deviation will be less than a factor of 1000.0

Out of the list chemnames, find the best chemical and 
best time point, that most reduces the integrated variance
over designchemnames 

sensvect : a sensitivity vector (length = # of params) of 
unmeasurable quantity of interest
designchemnames : list of measurable chemicals

sensvect could be the sensitivity of a single chemical at a
single timepoint; then can use method get_sens_vect (see elsewhere
in this file) to compute this sensitivity vector. In that
case we are designing over the species variance at that single point

chemnames : list of chemical names for which the 
variance at all timepoints will be computed
logprior : prior on parameters. logprior = log(1000.0)
means params allowed to vary by about a factor of 1000.0
return values : 
times: times of the trajectory
bestfit: a dictionary of best fit trajectories (keys are entries in chemnames)
var: a dictionary of variances (keys are entries in chemnames)

Get the variance for a single function of parameters 
that has a sensitivity vector sensvect. Useful for looking at
variances in parameter combinations, or simple functions of
parameters. Note that if we are concerned with ratios and 
products of parameters, it's often best to consider sensvect
as a sensitivity w.r.t. log parameters 

chemnames : list of chemical names at which we will look 
at variance
sensvect_design : the sensitivity vector (one by no. params array) at
the new design point.
returns : (times, varchange) 
the times and the change in variances at those times (should
be negative) for each of the chemicals in chemnames, after the 
addition of the new timepoint. varchange is a dictionary 
indexed by entries in chemnames.

sensvect : given a single function f(p) of parameters, this is the
derivative w.r.t. each of the parameters (in log parameters). For
ratios or products of rate constants, f(p) is a linear function
sensvect_design : the sensitivity vector of the new point in the 
design you wish to add
returns: the variance change of the quantity f(p), given the 
addition of the new data point, with sensitivity vector sensvect_design.

get a sensitivity vector for a chemical "chemname" at a
time, time 

get an array of sens_vects for all the times the chemical is defined
and convert to log sensitivities 

This is similar to var_change except now we pass in a sensarray
instead of sensvect --- this is a matrix of sensvects aligned rowwise.
Row i will be multiplied by sqrt(weights[i]) where sum(weights)=1 and 
each weight is a number between zero and one. We will return the 
change in variance for all the chemicals in chemnames 

For this cost function we're going to assume unconstrained 
variables are being passed in, so we need to convert them to 
a range between 0 and 1. The sum of the weights should also = 1 

weights : a vector of positive numbers with length the same as the number of 
rows of sensarray_design. The weights should sum to 1 
chemnames: a list of unmeasurable chemical names over which we wish
to design experiments 
sensarray_design: an array of sensitivities of measurable chemicals
or just an array of sensitivity vectors, each row a different 
sensitivity vector 
logprior : prior on parameters. logprior = log(1000.0) allows parameters
to fluctuate by a factor of 1000 

chemnames: list of chemical names 
logprior: prior on params. logprior = log(1000.0) means parameters
allowed to fluctuate by a factor of 1000 

chemnames: list of chemical names
logprior: prior on params
Plots the standard deviation of the chemicals when the variance
is computed using logs of the chemical trajectories. This 
makes sure the final plots do not have best_fit+-stddev that
do not become negative 

chemnames: list of chemical names
sensvect_design: a sensivity vector of a quantity that is 
measurable
This will plot the old and new variances of the chemicals in 
chemnames, given a new measurement that has sensitivity vector
sensvect_design

weights : a proposed set of weights for each of the row vectors in 
sensarray_design
chemnames : a list of chemicals for which we will plot the variance
logprior : as before 
This will plot the old and new variances on chemnames, similar to
above.
NOTE: the weights that are passed in do not necessarily have to sum to
one. e.g. if the weights are normalized such that max(weights) = 1, then
by scaling all the weights by 1/sigma, you are then assuming that
the most accurate measurement has an error of size sigma. sigma for 
example could be 20% of the maximum value of a trajectory.