Figure A.10 (page 658):

Fit of model to measurements using estimated parameters.

Code for Figure A.10

Text of the GNU GPL.

main.m

%%
%% Estimate parameters for the A->B->C model
%% jbr,  11/2008

clear all
close all

% ensure sundials/cvode is in the path
if (~ exist ('CVode'))
  if (exist ('startup_STB'))
    %%startup_STB;
    startup_STB ("/usr/share/octave/packages/3.2/sundialsTB");
  else
    error ('%s: requires the Sundials Toolbox, please see the README file for installation instructions');
  end
end

global model 

k1 = 2;
k2 = 1;
ca0 = 1;
cb0 = 0;
cc0 = 0;
thetaac = [k1; k2];


tfinal = 6;
nplot = 100;
model.odefcn  = @massbal_ode;
model.tplot   = linspace(0, tfinal, nplot)';
model.param   = thetaac;
model.ic = [ca0; cb0; cc0];

objective.estflag = [1, 2];
objective.paric   = [0.5; 3];
objective.parlb   = [1e-4; 1e-4;];
objective.parub   = [10; 10];

measure.states = [1,2,3];
%% load measurements from a file
table = load ('ABC_data.dat')
measure.time = table(:,1);
measure.data = table(:,1+measure.states);

%% estimate the parameters
estimates = parest(model, measure, objective);


disp('Estimated Parameters and Bounding Box')
[estimates.parest estimates.bbox]


%%plot the model fit to the noisy measurements
figure(1);
plot(model.tplot, estimates.x, measure.time, measure.data, 'o');

%% if 2-d, plot confidence ellipse
if (length(objective.estflag) == 2)
  [xe, ye] = ellipse(estimates.Pinv, estimates.b, 100, estimates.parest);
  figure(2); 
  plot(xe, ye)
end

table = [model.tplot, estimates.x'];
save ABC.dat table

massbal_ode.m

function xdot = massbal_ode(x, t, theta)
  ca  = x(1);
  cb  = x(2);
  cc  = x(3);
  k1  = theta(1);
  k2  = theta(2);
  r1 = k1*ca;
  r2 = k2*cb;
  xdot = [-r1; r1-r2; r2];

../../util/common/ellipse.m

%% Copyright (C) 2001, James B. Rawlings and John W. Eaton
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License as
%% published by the Free Software Foundation; either version 2, or (at
%% your option) any later version.
%%
%% This program is distributed in the hope that it will be useful, but
%% WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%% General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; see the file COPYING.  If not, write to
%% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
%% MA 02111-1307, USA.

%% [x, y, major, minor, bbox] = ellipse (amat, level, n, shift)
%%
%% Given a 2x2 matrix, generate ellipse data for plotting.  The
%% arguments N and SHIFT are optional.  If N is an empty matrix, a
%% default value of 100 is used.

function [x, y, major, minor, bbox] = ellipse (amat, level, n, shift)

  if (nargin < 3)
    n = 100;
  end

  if (isempty (n))
    n = 100;
  end

  if (nargin < 4)
    shift = [0, 0];
  end

  ss = size (shift);

  if (any (ss ~= [1, 2]))
    if (ss == [2, 1])
      shift = shift';
    else
      error ('shift must be a 2-element row vector');
    end
  end

  if (nargin > 1)

    [v, l] = eig (amat / level);

    dl = diag(l);
    if (any (imag (dl)) || any (dl <= 0))
      error ('ellipse: amat must be positive definite');
    end

    %% Generate contour data.

    a = 1 / sqrt (l(1,1));
    b = 1 / sqrt (l(2,2));

    t = linspace (0, 2*pi, n)';

    xt = a * cos (t);
    yt = b * sin (t);

    %% Rotate the contours.

    ra = atan2 (v(2,1), v(1,1));

    cos_ra = cos (ra);
    sin_ra = sin (ra);

    x = xt * cos_ra - yt * sin_ra + shift(1);
    y = xt * sin_ra + yt * cos_ra + shift(2);

    %% Endpoints of the major and minor axes.

    minor = (v * diag ([a, b]))';
    major = minor;

    major(2,:) = -major(1,:);
    minor(1,:) = -minor(2,:);

    t = [1; 1] * shift;

    major = major + t;
    minor = minor + t;

    %% Bounding box for the ellipse using magic formula.

    ainv = inv (amat);
    xbox = sqrt (level * ainv(1,1));
    ybox = sqrt (level * ainv(2,2));

    bbox = [xbox ybox; xbox -ybox; -xbox -ybox; -xbox ybox; xbox ybox];

    t = [1; 1; 1; 1; 1] * shift;
    bbox = bbox + t;

  else
    error ('usage: ellipse (amat, level, n, shift)');
  end

%% endfunction

../../util/common/parest.m

function retval = parest(model, measure, objective)
    
%      This program is free software: you can redistribute it and/or modify
%      it under the terms of the GNU General Public License as published by
%      the Free Software Foundation, either version 3 of the License, or
%      (at your option) any later version.
%  
%      This program is distributed in the hope that it will be useful,
%      but WITHOUT ANY WARRANTY; without even the implied warranty of
%      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%      GNU General Public License for more details.
%  
%      You should have received a copy of the GNU General Public License
%      along with this program.  If not, see. 
%
%
%  Function File:	estimates = @parest(model, measure, objective)
%  
%  Estimates parameters p for differential equation models
%		dx/dt = f(x,t,p)
%		x(0) = x0(p)
%  for given masurements y_i = h(x(t_i)) at time points t_i.
% 
%  The estimates are obtained by solving a weighted least squares problem
% 
%		min Sum ( h(x(t_i)) - y_i )' W ( h(x(t_i)) - y_i )
%		 p   i
%
%  with weight matrix W.
%
%  Input Arguments
%  ===============
%  The three input arguments model, measure, objective 
%  are fields and contain:
%
%  1. model
%  --------
%  model.odefcn (required)
%		Provides the right hand side of the 
%		differential equation in the form 
%		xdot = f(x, t, parameters).
%
%  model.cvodesfcn (optional)
%		Provides the right hand side of the 
%		differential equation in the form 
%		expected by CVODES,
%		[xdot, flag,new_data] = cvodesrhs(t, x, data).
%		If you provide this function and model.odefcn,
%		then the CVODES RHS function is used and 
%		model.odefcn is not called.
%
%  model.dodedp, model.dodedx (recommended)
%  		Partial derivatives of f with respect to the 
%		model parameters and the state respectively.
%  		Necessary for exact calculation of gradients 
%		of the objective function.
%  
%  model.ic (required)
%  		Column vector with initial condition of the 
%		system states.
%		If several data sets for different experiments 
%		corresponding to different initial conditions 
%		are given, this is a matrix in which each the
%		i-th column is the state initial condition 
%		corresponding to the i-th experiment.
%
%  model.param (required)
%  		Values of model parameters or initial guesses for 
%		the parameters which are estimated.
%  
%  model.stateic_par (optional)
%  		Vector or matrix of the same size as model.ic
%		which describes which initial conditions of
%		the state are parameters (and possibly have 
%		to be estimated).
%		E.g. stateic_par = [0; 1; 0; 2]
%		states that the second initial condition is
%		parameter 1 and the fourth initial condition
%		is parameter 2.
%  
%  model.tplot (optional)
%  		Time points for which estimated states are 
%		returned, default: measure.time (cf. below).
%
%  model.atol
%  model.rtol (optional)
%		Absolute and relative tolerances for the
%		integration of the differential equations
%		with CVodes. Default: sqrt(eps)
%  model.odesteps (optional)
%		Maximum number of steps for solving the 
%		differential equations with CVodes.
%		Default: 100,000
%
%  2. measure
%  ----------
%  measure.states
%  or
%  measure.statefcn, measure.dstatedx
%  or
%  measure.measurematrix (recommended)
%		The user can either provide the measured states
%		as a vector or a function y = h(x) along with
%		its derivative dh/dx or the measurematrix C
%		for y = C*x.
%		If neither is provided, parest assumes all 
%		states to be measured.
%		If two or more are provided, an error is given.
%
%  measure.time (required)
%		Column vector of the measurement times.
%		If several data sets for different initial 
%		conditions are given, this is a matrix with
%		the second dimension corresponding to the 
%		index of the experiment.
%  
%  measure.data (required)
%  		Matrix of data, the columns correspond to the 
%		measured states or output functions h 
%		respectively and the rows correspond to 
%		the times provided in measure.time.
%		If several data sets for different initial 
%		conditions are given, this is a 3-dimensional 
%		structure with the third dimension corresponding 
%		to the index of the experiment.
%
%  measure.weight (optional)
%		Matrix W for weighted least squares or covariance
%		matrix of measurement noise, default: eye(nmeas).
%		Matrix W should be symmetric, of size nmeas x nmeas
%		and is same for all timepoints and experiments.
%
%  measure.genwt (optional)
%		Generalized weights. Matrix G of size nmeas x nmeas x ntimepoints x nexpts.
%		Parest still needs all experiments to have the same number of time points.
%		If measure.genwt is provided, measure.weight is ignored.
%  
%  measure.alpha (optional)
% 		Alpha for confidence intervals, default: 0.95
%  
%  3. objective
%  ------------
%  objective.estflag (recommended)
%  		Vector of indices of parameters in model.param 
%		that have to be estimated.
%		If estflag is not given, all parameters are 
%		estimated. If estflag is an empty vector [], 
%		the output values are calculated for the given 
%		initial guess of parameters (cf. objective.paric).
%  
%  objective.paric, objective.parlb, objective.parub (recommended)
%  		Vectors of length of estflag with initial guess, 
%		lower and upper bounds for the parameters which are 
%		estimated. When no estflag is provided, these vectors 
%		have to be of the same length as model.param.
%		For empty / missing vectors the following defaults 
%		are used:
%		Default for objective.paric: values in model.param
%		Default for objective.parlb: realmin
%		Default for objective.parub: realmax
%		If only scalars are given for objective.parlb and
%		objective.parub, those bounds are used for all 
%		parameters.
%
%  objective.logtr (optional)
%  		Vector of logicals of length of estflag which 
%		indicates for which parameters a logarithmic tranformation 
%		should be performed.
%  
%  objective.sensitivityAdjointMethod (optional)
%  		Logical, if true use adjoint method for gradient
%		calculations, if false forward sensitivity 
%		analysis is used.
%
%  objective.printlevel (optional)
%		Decides which informations about the current
%		optimization process are displayed on the screen.
%		The options are currently available:
%		0: (default) no output on screen is given
%		1: parameter values, value of objective function and
%			gradient are displayed on screen during
%			optimization
%
%  Return Values
%  =============
%  The return value estimates is a field and contains:
%
%  estimates
%  ---------
%  estimates.parest, estimates.parestlogtr
%  		Estimated parameter values or its logarithm for
%  		those specified in objective.logtr.
%  
%  estimates.obj
%  		Value of the objective function.
%  
%  estimates.grad, estimates.gradlogtr
%  estimates.Hgn, estimates.Hgnlogtr
%  estimates.HFD
%		Gradient and Hessian of objective function with 
%		respect to the parameters or their logarithms for
%		those specified in objective.logtr.
%		The Hessian in Hgn and Hgnlogtr is calculated by 
%		using a Gauss-Newton Approximation, whereas HFD 
%		is calculated with a two sided finite difference 
%		approach applied on the gradients.
%
%  estimates.info  
%		Information flag returned by the optimization program.
%  
%  estimates.mu
%		Matrix of adjoint variable in form mu(:,measure.time)
%  		if adjoint method is used for gradient calculation
%		(cf. objective.sensitivityAdjointMethod).
%		If several data sets for different initial conditions 
%		are given, this is a 3-dimensional structure with the
%		third dimension corresponding to the index of the 
%		experiment.
%  
%  estimates.x
%  estimates.sx
%  		State and sensitivities for the optimal parameter
%		values at time tplot or by default the measurement
%		times.
%		x is a matrix with the ith column being the state at
%		the ith time point, i.e.
%		x(:,i) = x(t_i).
%		sx is a three-dimensional structure, with the last
%		argument corresponding to time, i.e.
%		sx(:,:,i) = dx/dp(t_i).
%		If several data sets for different initial conditions
%		are given, x (sx) is a 3 (4)-dimensional structure
%		with the third (fourth) dimension corresponding to
%		the index of the experiment.
%
%  estimates.y
%		Output of system for the optimal parameter
%		values at time tplot or by default the measurement
%		times (cf. estimates.x).
%
%  estimates.measpred
%  		Matrix of predicted measurements with columns
%		corresponding to the measurement times.
%  estimates.resid
%  		Matrix of residuals in the same format as 
%		estimates.measpred.
%		If several data sets for different initial 
%		conditions are given, resid and measpred are 
%		3-dimensional structures with the third dimension 
%		corresponding to the index of the experiment.
%
%  estimates.measvar
%  		Estimate of measurement error covariance matrix.
%
%
%  estimates.bbox, estimates.bboxFD
%		Bounding box for confidence intervals for given 
%		measure.alpha confidence level using the two methods of 
%		finding the Hessian matrix.
%
%  estimates.Pinv, estimates.b
%               Confidence ellipse:
%  (theta - estimates.parest)' * (estimates.Pinv) * (theta - estimates.parest)  \leq b
%
%  See also: -

  global mod_ meas_ obj_
  global solveroption

  %% create local copies of these objects to pass to likelihood
  meas_  = measure; obj_ = objective;  mod_ = model;

  %% 1. Check model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if(~isfield(mod_,'odefcn') && ~isfield(mod_,'cvodesfcn'))
	error('parest: model.odefcn and model.cvodesfcn missing: no rhs of ode provided');
  end

  if(~isfield(mod_,'param'))
	error('parest: model.param missing');
  end

  if(~isfield(mod_,'ic'))
	error('parest: model.ic missing');
  end

  %% if user provides no stateic_par or par_stateic
  if(~isfield(mod_,'stateic_par') || isempty(mod_.stateic_par))
	mod_.stateic_par = zeros(size(mod_.ic));
  end

  %% if entry in model.stateic_par is too large
  if(max(mod_.stateic_par) > numel(mod_.param))
	error('parest: entry in model.stateic_par exceeds number of parameters')
  end

  %% if sizes don't match
  %% e.g. if only model.stateic_par(3) = 2 was used
  %% for a system with 4 states
  %% -> fill up with zeros
  if (~isequal( size(mod_.ic), size(mod_.stateic_par) ) )
	warning('parest: sizes of model.ic and model.stateic_par do not match')
	stateic_par = zeros(size(mod_.ic));
	for i = 1:size(stateic_par,1)
		for j = 1:size(stateic_par,2)
			stateic_par(i,j) = mod_.stateic_par(i,j);
		end
	end
	mod_.stateic_par = stateic_par;
  end


  if(isfield(mod_,'tplot'))
  if(size(mod_.tplot,1) == 1 && size(mod_.tplot,2) ~= 1 )
	% if tplot is provided as row vector -> transpose it
	mod_.tplot = mod_.tplot(:);
	warning('parest: model.tplot is row vector')
  end
  end

  if(~isfield(mod_,'atol'))
	% absolute tolerance for CVodes
	mod_.atol = sqrt(eps);
  elseif(~isscalar(mod_.atol))
	% absolute tolerance for CVodes
	error('parest: model.atol has to be scalar')
  end

  if(~isfield(mod_,'rtol'))
	% relative tolerance for CVodes
	mod_.rtol = sqrt(eps);
  elseif(~isscalar(mod_.rtol))
	% relative tolerance for CVodes
	error('parest: model.rtol has to be scalar')
  end

  if(~isfield(mod_,'odesteps'))
	% maximum number of steps for CVodes
	mod_.odesteps = 100000;
  elseif(~isscalar(mod_.odesteps))
	% max number of steps tolerance for CVodes
	error('parest: model.atol has to be scalar')
  end

  %% 2. Check objective %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  %% if user provides no estflag, estimate all parameters
  if(~isfield(obj_,'estflag'))
	obj_.estflag = 1:numel(mod_.param);
	warning('parest: objective.estflag is missing. All parameters are estimated')
  end
  
  %% if estflag is empty, don't estimate any parameters 
  %% but calculate objective function etc
  if(isempty(obj_.estflag))
	obj_.estflag = 1:numel(mod_.param);
	estflagempty = true;
	warning('parest: objective.estflag is empty. No optimization is performed')
  else
	estflagempty = false;	
  end

  %% make column vector
	if ( size(obj_.estflag,2) > 1 )
		obj_.estflag = obj_.estflag';
	end

  %% number of parameters to estimate
  nparest = numel(obj_.estflag);

  %% if user provides no special initial guess take values out of model.param
  if(~isfield(obj_,'paric') || isempty(obj_.paric) )
	obj_.paric = mod_.param(obj_.estflag);
	warning('parest: objective.paric is missing. Model.param used as initial guess')
  end

  %% if user doesn't provide logtr
  if(~isfield(obj_,'logtr') || isempty(obj_.logtr))
	%as default don't use log transform
	obj_.logtr = logical(zeros(nparest,1));
  end
		
  %% check if sizes in objective are correct
  if (numel(obj_.logtr) ~= nparest)
	error('parest: Number of elements in objective.logtr must be equal to length of objective.estflag') 
  end
  if (~isfield(obj_,'parlb') || isempty(obj_.parlb))
	% if no lower bound is given
	obj_.parlb = realmin*ones(nparest,1);%1e-5*obj_.paric;
	warning('parest: objective.parlb is missing')
	%% error('parest: objective.parlb missing')
  elseif (numel(obj_.parlb) ~= nparest)
	if (isscalar(obj_.parlb))
		obj_.parlb = obj_.parlb*ones(1,nparest);
	else
	error('parest: Number of elements in objective.parlb must be equal to length of objective.estflag')
	end
  end
  if (~isfield(obj_,'parub') || isempty(obj_.parub))
	% if no upper bound is given
	obj_.parub = realmax*ones(nparest,1);%1e5*obj_.paric;
	warning('parest: objective.parub is missing')
	%% error('parest: objective.parub missing')
  elseif (numel(obj_.parub) ~= nparest)
	if (isscalar(obj_.parub))
		obj_.parub = obj_.parub*ones(1,nparest);
	else
	error('parest: Number of elements in objective.parub must be equal to length of objective.estflag')
	end
  end
  % ensure that parub, parlb, and paric are column vectors
  obj_.parub = obj_.parub(:);
  obj_.parlb = obj_.parlb(:);
  obj_.paric = obj_.paric(:);
  % ensure ub > lb
  if (any( obj_.parub - obj_.parlb <= 0))
	error('parest: Upper bound has to be greater than lower bound')
  end
  if (numel(obj_.paric) ~= nparest)
	error('parest: Number of elements in objective.paric must be equal to length of objective.estflag')
  end

  if(~isfield(obj_,'printlevel'))
	obj_.printlevel = 0;
  end

  %% 3. Check measure %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if(isfield(meas_,'states') + isfield(meas_,'statefcn') + isfield(meas_,'measurematrix') > 1)
	error('parest: more than one of measure.statefcn, states and measurematrix provided')
  end

  if(~isfield(meas_,'states') && ~isfield(meas_,'statefcn') && ~isfield(meas_,'measurematrix'))
	meas_.states = 1:size(mod_.ic,1);
	warning('parest: none of measure.statefcn, states and measurematrix provided')
  end

  if(isfield(meas_,'statefcn') && ~isfield(meas_,'dstatedx'))
	error('parest: measure.dstatedx missing');
  end

  if(~isfield(meas_,'data'))
	error('parest: measure.data missing');
  end

  if(~isfield(meas_,'time'))
	error('parest: measure.time missing');
  end

  if(size(meas_.time,1) == size(mod_.ic,2) && size(meas_.time,2) ~= size(mod_.ic,2))
	% if measurement time is provided as column vector -> transpose it
	meas_.time = meas_.time(:)';
	warning('parest: measure.time seems to be transposed')
  end

  %% if transpose of measurement.data is given
  %% data is supposed to be in format (states,time)
  if (size(meas_.data,2) ~= size(meas_.time,1) && size(meas_.data,1) == size(meas_.time,1)) 
	for k = 1:size(meas_.data,3)
	measdata(:,:,k) = (meas_.data(:,:,k))';
	end
	meas_.data = measdata;
%	warning('parest: measure.data seems to be transposed')
  end

  if (size(meas_.time,2) ~= size(meas_.data,3) || size(meas_.time,2) ~= size(mod_.ic,2) || size(meas_.time,1) ~= size(meas_.data,2))
	error('parest: sizes of measure.data, measure.time and model.ic are not compatible')
  end

  %% number of measurements at each time point, e.g. =2 if x1 and x2 are measured
  nmeas = size(meas_.data, 1);

  if(isfield(meas_,'states'))
	if( numel(meas_.states) ~= nmeas)
		error('parest: sizes of data and measurement.states do not match')
	end
  elseif(isfield(meas_,'statefcn'))
	y0 = meas_.statefcn(mod_.ic);
	dy0dx = meas_.dstatedx(mod_.ic);
	if( size(y0,1) ~= nmeas)
		error('parest: sizes of data and measurement.statefcn do not match')
	end
	if( size(dy0dx,1) ~= nmeas || size(dy0dx,2) ~= numel(mod_.ic) )
		error('parest: size of measurement.dstatedx seems wrong')
	end
  elseif(isfield(meas_,'measurematrix'))
	if( size(meas_.measurematrix,1) ~= nmeas || size(meas_.measurematrix,2) ~= numel(mod_.ic) )
		error('parest: size of measurematrix seems wrong')
	end
  end

  %% ensure that measure weight matrix is symmetric
  %% and has appropriate size
  if(isfield(meas_,'weight'))
	if (size(meas_.weight) == nmeas*ones(1,2))
  		meas_.weight = 0.5*(meas_.weight+meas_.weight');
	else
		error('parest: size of measurement.weight does not fit to measurement.data')
	end
  else
	meas_.weight = eye(nmeas);
	%warning('parest: weight matrix missing. Identity matrix used')
  end

  %% ensure that measure general weight matrix is of appropriate size
  %% and G(:,:,i,k) is symmetric
  if(isfield(meas_,'genwt'))
      if isequal(size(meas_.genwt),[nmeas nmeas size(meas_.data,2) size(mod_.ic,2)])
	for i=1:size(meas_.data,2)
	  for k=1:size(mod_.ic,2)
	    meas_.genwt(:,:,i,k) = 0.5*(meas_.genwt(:,:,i,k)+meas_.genwt(:,:,i,k)');
	  end
	end
      else
	error('parest: size of measure.genwt does not fit to measure.data or model.ic)')
      end
  end

	
  if(size(mod_.ic,2) ~= size(meas_.data,3))
	error('parest: number of data sets with different ICs does not fit with meaurement.data')
  end

  if(size(mod_.ic,2) ~= size(meas_.time,2))
	error('parest: number of initial condition sets does not fit with meaurement.time')
  end

  %% 4. prepare optimization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  %% log transform the specified parameters
	% ensure objective.logtr to be logical
	obj_.logtr = logical(obj_.logtr);
	trind =  obj_.logtr;
	obj_.paric(trind) = log(obj_.paric(trind));
	obj_.parlb(trind) = log(obj_.parlb(trind));
	obj_.parub(trind) = log(obj_.parub(trind));

  %% initialize the par sensitivities; note the pars that are ICs
  %% sx0 = 0 for parameters that are not state IC
  %% sx0 = 1 for parameters that are state IC
  %% third dimension corresponds to different experiments
	mod_.sx0 = zeros(size(mod_.ic,1), numel(mod_.param), size(mod_.ic,2));
	
	for i = 1:size(mod_.ic,1) %for all states
		for k = 1:size(mod_.ic,2) % for all experiments
			if (mod_.stateic_par(i,k) > 0)
				mod_.sx0(i,mod_.stateic_par(i,k) ,k) = 1;
			end
		end %for k = ...
	end % for i = ...
        %% chop out the part of sx0 (2nd index) for parameters not estimated
        mod_.sx0 = mod_.sx0(:,obj_.estflag,:);
  %% 5. call optimizer %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if(estflagempty)
	%% if estflag is empty -> no parameters to estimate, no optimization
	theta = obj_.paric;
	info = 'estflag empty, no optimization performed';
  else
	%% optimize in order to
	%% find best fit of parameters
	%% solveroption determines if/how gradient is provided for optimization
	solveroption = 2;
	if (isfield(obj_,'sensitivityAdjointMethod') && obj_.sensitivityAdjointMethod)
	        if(isfield(mod_,'dodedx') && isfield(mod_,'dodedp'))
			solveroption = 3;
        	else
        		warning('parest: adjoint method cannot be used because model.dodedp or dodedx missing')
        	end
	end

	%% minimize negative logarithm of likelihood function
	switch solveroption
	case 1 % solve without provided gradients - currently not used
	  [theta, phi, info] = fmincon (@likelihood, obj_.paric, ...
        		[], [], [], [], obj_.parlb, obj_.parub);
	case 2 % solve using forward sensitivities
	  options = optimset('GradObj','on');
	  [theta, phi, info] = fmincon ({@likelihood,@gradi}, obj_.paric, ...
        		[], [], [], [], obj_.parlb, obj_.parub,[],options);
	case 3	% solve using adjoint method for gradients
	  options = optimset('GradObj','on');
	  [theta, phi, info] = fmincon ({@likelihoodA,@gradi}, obj_.paric, ...
        		[], [], [], [], obj_.parlb, obj_.parub,[],options);
	end %end switch solveroption
  end %if(estflagempty)

  %% 6. work on return values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  %% (Insurance) call to likelihood with optimal parameters to compute solution,
  %% sensitivities, gradient and Hessian. Not necessary if optimizer's
  %% last call was with optimal values.
	phi = likelihood(theta);

  % if tplot is given, solve ODE for all time points in tplot
  if(isfield(mod_,'tplot'))
	data.theta = mod_.param;
	%% set the model parameters; note the log transformation
	data.theta(obj_.estflag) = theta;
	trind =  obj_.estflag(obj_.logtr);
	data.theta(trind) = exp(theta(obj_.logtr));
	%% solve the model and sensitivites at tplot
		if (~isfield(mod_,'cvodesfcn'))
			if (~isfield(mod_,'dodedx') || ~isfield(mod_,'dodedp') )
				[x, sx, status] = cvode_sens(@oderhs, [], mod_.ic, ...
					obj_.estflag, data, mod_.sx0, mod_.tplot);
			else
				[x, sx, status] = cvode_sens(@oderhs, @sensrhs, mod_.ic, ...
					obj_.estflag, data, mod_.sx0, mod_.tplot);  
			end
		else
			if (~isfield(mod_,'dodedx') || ~isfield(mod_,'dodedp') )
				[x, sx, status] = cvode_sens(mod_.cvodesfcn, [], mod_.ic, ...
				obj_.estflag, data, mod_.sx0, mod_.tplot);
			else
				[x, sx, status] = cvode_sens(mod_.cvodesfcn, @sensrhs, mod_.ic, ...
					obj_.estflag, data, mod_.sx0, mod_.tplot);  
			end
		end		
	% if cvode_sens failed
	if (status ~= 0)
		error('parest: CVode failed with status = %d', status);
	end
	retval.x = x;
	retval.sx = sx;
  else
	%% if no tplot is provided by user return estimates and 
	%% sensitivities at available measurement time point
	retval.x = obj_.x;
	retval.sx = obj_.sx;
  end %if(isfield(mod_,'tplot'))

  nics = size(retval.x, 3);
  ntimes = size(retval.x,2);
  for k = 1:nics %for each IC = dataset
	for i = 1: ntimes
	if (isfield(meas_,'states'))
		% if user provides vector of measured states
		output(:,i,k) = retval.x(meas_.states,i,k);
	elseif (isfield(meas_,'statefcn'))
		% if user provides a measurement state function
		output(:,i,k) = meas_.statefcn(retval.x(:,i,k));
	elseif (isfield(meas_,'measurematrix'))
		% if output matrix y = measurematrix*x is provided
		output(:,i,k) = meas_.measurematrix*(retval.x(:,i,k));
	end
	end % for i = ...
  end % for k = ...

  %% store output for return
	retval.y = output;
	retval.obj  = phi;
	retval.info = info;
  %% log transformed parameters
	retval.parestlogtr = theta;
  %% gradient and Hessian; pull these out of common; computed in likelihood
	retval.resid = obj_.resid;
	retval.measpred = obj_.measpred;
	retval.grad  = obj_.grad;
	retval.gradlogtr = obj_.gradlogtr;
  %% undo log transform the specified parameters
	thetawolog = theta;
	thetawolog(obj_.logtr) = exp(theta(obj_.logtr));
	retval.parest = thetawolog;

  %% if user provides no alpha for confidence intervals
  if(~isfield(meas_,'alpha'))
	meas_.alpha = 0.95;
  end

  if(isfield(obj_,'Hgn'))
	retval.Hgn  = obj_.Hgn;
	retval.Hgnlogtr = obj_.Hgnlogtr;
	%%compute confidence intervals
	p     = length(theta);
	ndata = length(meas_.time)*size(meas_.data,3);
	Fstat = finv(meas_.alpha, p, ndata-p);
	Hgn   = retval.Hgn;
	if (isscalar(Hgn))
	  invHgn = 1./Hgn;
	else
	  invHgn = inv(Hgn);
	end
        retval.bbox = sqrt(2*p/(ndata-p)*Fstat*phi*diag(invHgn));
%  	keyboard
        samvar = retval.obj/(ndata-p);
        retval.Pinv = retval.Hgn/(2*samvar);
        retval.b = p*Fstat;
  end

  %% adjoint variable mu if computed in likelihoodA
  %% matrix of form mu(:,time)
  if(isfield(obj_,'mu'))
	retval.mu = obj_.mu;
  end

  %% Estimate of Measurement Variance Matrix
	residv = retval.resid;
	%% first find and delete bad data
	residv(find(~isfinite(residv))) = 0;
	%% formula for Variance
	retval.measvar = 0;
	for k = 1:size(mod_.ic,2)
	retval.measvar = retval.measvar + residv(:,:,k)*residv(:,:,k)'/size(residv,2);
	end

   %% approximate Hessian via Finite Difference
	HFD = zeros(nparest,nparest);
	for j =1:nparest
		FDperturbation = 1*sqrt(1e-3)*thetawolog(j);
		thetaFD = thetawolog;
		thetaFD(j) = thetawolog(j)+FDperturbation; %perturb one parameter
		thetaFD2 = thetawolog;
		thetaFD2(j) = thetawolog(j)-FDperturbation; %perturb one parameter
		thetaFD2(obj_.logtr) = log(thetaFD2(obj_.logtr)); %log trans
		likelihood(thetaFD2);
		grad2 = obj_.grad;
		thetaFD(obj_.logtr) = log(thetaFD(obj_.logtr)); %log trans
		likelihood(thetaFD);
		%HFD(:,j) = (obj_.grad - retval.grad)/FDperturbation;
		HFD(:,j) = (obj_.grad - grad2)/(2*FDperturbation);
    	end
	HFD = 0.5*(HFD+HFD'); % make HFD symmetric
	retval.HFD = HFD;
	  %%compute confidence intervals with FD
	p     = length(theta);
	ndata = length(meas_.time)*size(meas_.data,3);
	Fstat = finv(meas_.alpha, p, ndata-p);
	if (isscalar(HFD))
	  invHFD = 1./HFD;
	else
	  invHFD = inv(HFD);
    	end
        retval.bboxFD = sqrt(2*p/(ndata-p)*Fstat*retval.obj*diag(invHFD));

end %parest.m
%% END OF FUNCTION PAREST.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function phi = likelihood(theta)
%% Calculates negative logarithm of likelihood function,
%% this function is minimized by parest.m.
%% Gradients and sensitivities are calculated by forward calculations
%% using cvode_sens.
  global mod_ meas_ obj_ likelihood_last_theta

  % print on screen if printlevel is set accordingly
  if(obj_.printlevel == 1)
	disp('----------------------');
	disp('examined parameter values:');
	disp(theta);
  end

  %% data.theta contains all parameters, not only those which are estimated!
  data.theta = mod_.param;
  %% set the model parameters; note the log transformation
  data.theta(obj_.estflag) = theta;
  trind =  obj_.estflag(obj_.logtr);
  data.theta(trind) = exp(theta(obj_.logtr));

  %% set the state ICs that are parameters
  for i = 1:numel(mod_.stateic_par)
	if (mod_.stateic_par(i) > 0)
		mod_.ic(i) =  data.theta(mod_.stateic_par(i));
	end
  end

  %% solve the model and sensitivites at the measurement times;
  %% predict measurements
	if (~isfield(mod_,'cvodesfcn'))
		if (~isfield(mod_,'dodedx') || ~isfield(mod_,'dodedp') )
			[x, sx, status] = cvode_sens(@oderhs, [], mod_.ic, ...
			obj_.estflag, data, mod_.sx0,  meas_.time);
		else
			[x, sx, status] = cvode_sens(@oderhs, @sensrhs, mod_.ic, ...
			obj_.estflag, data, mod_.sx0,  meas_.time);  
		end
	else
		if (~isfield(mod_,'dodedx') || ~isfield(mod_,'dodedp') )
			[x, sx, status] = cvode_sens(mod_.cvodesfcn, [], mod_.ic, ...
			obj_.estflag, data, mod_.sx0,  meas_.time);
		else
			[x, sx, status] = cvode_sens(mod_.cvodesfcn, @sensrhs, mod_.ic, ...
			obj_.estflag, data, mod_.sx0,  meas_.time);  
		end
	end
  % if cvode_sens failed
	if (status ~= 0)
		phi = realmax;
		return;
	end

  nmeas = size(meas_.data, 1);
  ntimes = size(meas_.data, 2);
  nics = size(mod_.ic, 2);
  npar = numel(theta);
  resid = zeros(nmeas,ntimes,nics);
  measpred = zeros(nmeas,ntimes,nics);
  nx = size (mod_.ic,1);

  %% is a loop best here?
  phi = 0;
  grad = zeros(npar, 1);
  Hgn  = zeros(npar, npar);
  for k = 1:nics %for each IC = dataset
	for i = 1: ntimes
		% for current time point i and data set k:
	
		if (isfield(meas_,'states'))
			% if user provides vector of measured states
			measpred(:,i,k) = x(meas_.states,i,k);
			resid(:,i,k) = meas_.data(:,i,k) - x(meas_.states,i,k);
			IdentityMatrix = eye(nx);
			% Jacobian at of h
			sqdhdx = IdentityMatrix(meas_.states,:);
		elseif (isfield(meas_,'statefcn'))
			% if user provides a measurement state function
			measpred(:,i,k) = meas_.statefcn(x(:,i,k));
			resid(:,i,k) = meas_.data(:,i,k) - meas_.statefcn(x(:,i,k));
			% Jacobian of h at time point i if statefcn provided by user
			sqdhdx = meas_.dstatedx(x(:,i,k));
		elseif (isfield(meas_,'measurematrix'))
			% if output matrix y = measurematrix*x is provided
			measpred(:,i,k) = meas_.measurematrix*(x(:,i,k));
			resid(:,i,k) = meas_.data(:,i,k) - meas_.measurematrix*(x(:,i,k));
			% Jacobian of h at time point i if statefcn provided by user
			sqdhdx = meas_.measurematrix;
		end	

		% squeeze matrix / vector out of struct
		sqsx = sx(:,:,i,k);
		r = resid(:,i,k);
		%% look for bad data (inf, NaN,...) and do not use it for calculation
		badData = find(~isfinite(r));
		r(badData) = 0;
		sqsx(badData,:) = 0;
		% tmp = dh/dp via chain rule
		tmp = sqdhdx*sqsx;
		% update for current time point i and data set k
		if isfield(meas_,'genwt')
		  phi = phi + r'*meas_.genwt(:,:,i,k)*r;
		  grad = grad + tmp'*meas_.genwt(:,:,i,k)*r;
		  Hgn = Hgn + tmp'*meas_.genwt(:,:,i,k)*tmp;
		else
		  phi = phi + r'*meas_.weight*r;
		  grad = grad + tmp'*meas_.weight*r;
		  Hgn = Hgn + tmp'*meas_.weight*tmp;
		end

	end %i = 1: ntimes
  end %k = 1:nics

  %% gradient with respect to untransformed parameters
	grad = -2*grad';
	obj_.grad = grad;
  %% gradient with respect to log transformed parameters
  %% this is the optimizer's gradient (see gradi function)
	scale = ones(npar,1);
	scale(obj_.logtr) = data.theta(trind);
	obj_.gradlogtr = obj_.grad*diag(scale);
  %% Gauss Newtown approximation of Hessian with untransformed parameters
	obj_.Hgn =  2*Hgn;
	obj_.Hgnlogtr =  2*diag(scale)*Hgn*diag(scale);
	obj_.x = x;
	obj_.sx = sx;
	obj_.resid = resid;
	obj_.measpred = measpred;
  
  % print on screen if printlevel is set accordingly
  if(obj_.printlevel == 1)
	disp('objective function:');
	disp(phi);
	disp('gradient:');
	disp(grad);
  end

%% check on where the optimizer is looking
%theta
%phi
end %Likelihood
%% END OF FUNCTION LIKELIHOOD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x, sx, status] = cvode_sens(f, fp, x0, estflag, data, sx0, tmeas)
%% Call the cvode integrator and calculate states, x, and 
%% sensitivities, sx, at the requested times.
%% f is rhs of ode, fp is rhs of sensitivity ode
%% status = 0 -> no errors
global mod_

  t0 = 0;
  nics = size(x0,2);

  %% output size
  nx = size(x0,1);
  nt = size(tmeas,1);
  np = size(sx0,2);
  %% output initialization
  x = zeros (nx,nt,nics);
  sx = zeros (nx, np,nt,nics);

  for k = 1:nics %for each IC =  dataset
	
	%% Set options for integrator.
	options = CVodeSetOptions ('UserData', data, ...
				   'RelTol', mod_.rtol, ...
				   'AbsTol', mod_.atol, ...
				   'MaxNumSteps',mod_.odesteps);
	
	%% Allocate storage and set initial conditions for integrator.
%-	CVodeMalloc (f, t0, x0(:,k), options, data);
        CVodeInit(f,'BDF','Newton',t0,x0(:,k),options);
	%% Number of paramters
	np = numel(estflag);

	fsa_options = CVodeSensSetOptions('method','Simultaneous', ...
					 'ErrControl', true,...
					 'ParamField', 'theta',...
					 'ParamList', estflag,...
					 'ParamScales', ones(1, np));

	%% Set options for forward sensitivity problem.
	if (isempty(fp))
	   CVodeSensInit (np, [], sx0(:,:,k), fsa_options);		      
	else
	   CVodeSensInit (np, fp, sx0(:,:,k), fsa_options);		      
	end %if (isempty(fp))

	%% Allocate storage for forward sensitivity problem.
%-	CVodeSensMalloc (np, 'Simultaneous', sx0(:,:,k), fsa_options);

	if (tmeas(1) == t0)  
		% initial condition for this data set k
		x(:,1,k) = x0(:,k);
		sx(:,:,1,k) = sx0(:,:,k);
		iout = 1; 
	else
		iout = 0;
	end

	tout = tmeas(1+iout);
	nsteps = nt;

	while (true)
		[status, t, x_step, sx_step] = CVode (tout, 'Normal');
		if (status == 0)
			iout = iout + 1;
			x(:,iout,k) = x_step;
			sx(:,:,iout,k) = sx_step;
		elseif (status < 0)
			warning ('parest: CVode failed with status = %d', status);
			break;
		end %if (status == 0)
		if (iout == nsteps)
			break;
		end
		tout = tmeas(1 + iout);
	end % while(true)
	%% Free integrator storage.
	CVodeFree ();

  end % for k = 1:nics

end %cvode_sens
%% END OF FUNCTION cvode_sens %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function phi = likelihoodA(theta)
%% Calculates negative logarithm of likelihood function,
%% this function is minimized by parest.m.
%% Gradients are calculated by adjoint method.
  global mod_ meas_ obj_ data

  if(obj_.printlevel == 1)
	disp('----------------------');	
	disp('examined parameter values:');
	disp(theta);
  end

  %% data.theta contains all parameters, not only those which are estimated!
  data.theta = mod_.param;
  %% set the model parameters; note the log transformation
  data.theta(obj_.estflag) = theta;
  trind =  obj_.estflag(obj_.logtr);
  data.theta(trind) = exp(theta(obj_.logtr));

  %% set the state ICs that are parameters
  for i = 1:numel(mod_.stateic_par)
	if (mod_.stateic_par(i) > 0)
		mod_.ic(i) =  data.theta(mod_.stateic_par(i));
	end
  end

  t0 = 0;
  nmeas = size(meas_.data, 1);
  ntimes = size(meas_.data, 2);
  nx = size(mod_.ic,1);
  npar = numel(obj_.estflag);

  %% initialize objective and gradient
  phi = 0;
  grad = zeros (1, npar);

  %% Set options for integrator.
  	options = CVodeSetOptions ('UserData', data, ...
				   'RelTol', mod_.rtol, ...
				   'AbsTol', mod_.atol, ...
				   'MaxNumSteps', mod_.odesteps);

  %% number of IC's = data sets
  nics = size(mod_.ic, 2);

  x = zeros (nx,ntimes,nics);
  resid = zeros (nx,ntimes,nics);
  measpred = zeros (nx,ntimes,nics);
  mu = zeros (nx,ntimes,nics);

for k = 1:nics

  %% Allocate storage and set initial conditions for integrator.
  if (~isfield(mod_,'cvodesfcn'))
  CVodeInit (@oderhs, 'BDF', 'Newton', 0, mod_.ic(:,k), options);
  else
  CVodeInit (mod_.cvodesfcn, 'BDF', 'Newton', 0, mod_.ic(:,k), options);
  end

  %% Allocate memory for checkpointing of forward solution for
  %% backward integration of adjoint problem.
  CVodeAdjInit (150, 'Hermite');

  tmeas = meas_.time(:,k);
  t0 = 0;
  t1 = tmeas(end);

  x(:,:,k) = zeros (nx,ntimes);
  mu(:,:,k) = zeros (nx,ntimes);

  %% Integrate forward, saving intermediate points for calculation of
  %% the objective function.
	if (tmeas(1) == t0)  
		% initial condition for this data set k
		x(:,1,k) = mod_.ic(:,k);
		sx(:,:,1,k) = sx0(:,:,k);
		iout = 1; 
	else
		iout = 0;
	end

	tout = tmeas(1+iout);
	nsteps = nt;

	while (true)
		[status, t, x_step, sx_step] = CVode (tout, 'Normal');
		if (status == 0)
			iout = iout + 1;
			x(:,iout,k) = x_step;
		elseif (status < 0)
			warning ('parest: CVode failed with status = %d', status);
			CVodeFree ();
			phi = realmax;
			return;
		end %if (status == 0)
		if (iout == nsteps)
			break;
		end
		tout = tmeas(1 + iout);
	end % while(true)

  if(status < 0)
	CVodeFree ();
	phi = realmax;
	return;
  end

  %%compute residual and objective function for these param values
  %%compute gradient, here NO Hessian approximation; have to store in common
  if (isfield(meas_,'states'))
		% if user provides a vector of measured states
		statefcnExists = false;
		measpred(:,:,k) = x(meas_.states,:,k);
		resid(:,:,k) = meas_.data(:,:,k) - x(meas_.states,:,k);
		IdentityMatrix = eye(nx);
		dhdx = IdentityMatrix(meas_.states,:);
	elseif (isfield(meas_,'measurematrix'))
		% if user provides measurement matrix y = measurematrix*x
		statefcnExists = false;
		measpred(:,:,k) = meas_.measurematrix*x(:,:,k);
		resid(:,:,k) = meas_.data(:,:,k) - meas_.measurematrix*x(:,:,k);
		dhdx = meas_.measurematrix;
	else
		% if user provides a measurement state function
		statefcnExists = true;
		for t_ind = 1:ntimes	%loop over all time points
			resid(:,t_ind,k) = meas_.data(:,t_ind,k) - meas_.statefcn(x(:,t_ind,k));
			measpred(:,t_ind,k) = meas_.statefcn(x(:,t_ind,k));
		end
  end

  obj_.resid = resid;
  obj_.measpred = measpred;
  %% find and delete bad data
  resid(find(~isfinite(resid))) = 0;
	
  %% update objective function
  if isfield(meas_,'genwt')
    for i=1:size(resid,2)
      phi = phi + resid(:,i,k)'*meas_.genwt(:,:,i,k)*resid(:,i,k)
    end
  else
    phi = phi + trace(resid(:,:,k)'*meas_.weight*resid(:,:,k));
  end

  %% Free storage for the integration.
  %CVodeFree ();

  %% Set up reverse integration.  Must be done after the forward
  %% integration is complete.

  %% Initial condition for the backward quadrature.
  w = zeros (1, npar);

  %% Set options for the backward integration of the adjoint equation and
  %% quadrature.
  optionsb = CVodeSetOptions ('RelTol', mod_.rtol, ...
			      'AbsTol', mod_.atol);

  optionsQB = CVodeQuadSetOptions('ErrControl', true, ...
				  'RelTol', sqrt(eps), ...
				  'AbsTol', sqrt(eps));



  %% Initial condition for the backward integration of the adjoint
  %% equations
  if (statefcnExists)
	% Jacobian at last time point if statefcn provided by user
	sqdhdx = meas_.dstatedx(x(:,end,k));
	%% find and delete bad data
	sqdhdx(find(~isfinite(meas_.data(:,end,k))),:) = 0;
  else
	% Jacobian at last time point if only measurement provided by user
	sqdhdx = dhdx;
	%% find and delete bad data
	sqdhdx(find(~isfinite(meas_.data(:,end,k))),:) = 0;
  end
  
  if isfield(meas_,'genwt')
    mu_init = -2*( resid(:,end,k)'*meas_.genwt(:,:,end,k)*sqdhdx )';
  else
    mu_init = -2*( resid(:,end,k)'*meas_.weight*sqdhdx )';
  end
  mu(:,end,k) = mu_init;
  %% Allocate storage and set initial conditions for integrator.
  idxB = CVodeInitB (@adjrhs, 'BDF', 'Newton', t1, mu_init, optionsb);

  CVodeQuadInitB(idxB, @sensrhsA, w, optionsQB);
  %% Integrate backward, accumulating in w the gradient of phi, which is
  %% the least-squares objective \sum resid*meas_.weight*resid.
  for (i = ntimes-1:-1:1)
	tinit = tmeas(i+1);
	tout = tmeas(i);
	CVodeReInitB (idxB, tinit, mu_init, optionsb);

   	while (true)
		[status, t_step, mu_step, w_step] = CVodeB (tout, 'Normal');

		if (status < 0)
			warning ('parest: CVodeB failed with status = %d', status);
			CVodeFree ();
			phi = realmax;
			return;
		end

		%% We are integrating backward, and we are storing the results in
		%% reverse order too.
    		if (status == 0 || (status == 1 && t_step == 0))
		if (statefcnExists)
			% Jacobian at time point i if statefcn provided by user
			sqdhdx = meas_.dstatedx(x(:,i,k));
			%% find and delete bad data
			sqdhdx(find(~isfinite(meas_.data(:,i,k))),:) = 0;
		else
			% Jacobian at time point i if no measurement fcn provided by user
			sqdhdx = dhdx;
			%% find and delete bad data
			sqdhdx(find(~isfinite(meas_.data(:,i,k))),:) = 0;
		end
			if isfield(meas_,'genwt')
			  mu_init = mu_step - 2*(resid(:,i,k)'*meas_.genwt(:,:,i,k)*sqdhdx )';
			else
			  mu_init = mu_step - 2*(resid(:,i,k)'*meas_.weight*sqdhdx )';
			end
			mu(:,i,k) = mu_init;
			w = w + w_step';
			break;
		end % if (status == ...
	end % while(true)
  end % for i = ntimes:...

  if( tmeas(1) > t0) %% integrate till t = 0
	CVodeReInitB (idxB, tmeas(1), mu_init, optionsb);
	[status, t_s, mu_0, w_0] = CVodeB (t0, 'Normal');
	if (status ~= 0)
		warning ('parest: CVodeB failed with status = %d', status);
		return;
	else
		w = w + w_step';
   	end
  else
  	mu_0 = mu_init;
  end

  %% Free storage for the integration.
	CVodeFree ();
  %CVadjFree ();

    %% gradient with respect to untransformed parameters
    %% Eqn (3.19) of CVodeS user guide
    %% sx0 = 0 for parameters that are not state IC
    %% sx0 = 1 for parameters that are state IC
    grad = grad + w + mu_0'*mod_.sx0(:,:,k);

end % for k = 1:nics

    obj_.grad = grad;
    %% gradient with respect to log transformed parameters
    %% this is the optimizer's gradient (see gradi function)
    scale = ones(npar,1);
    scale(trind) = data.theta(trind);
    obj_.gradlogtr = obj_.grad*diag(scale);
    %% values of x at tmeas
    obj_.x = x;
    obj_.mu = mu;

  % print on screen if printlevel is set accordingly
  if(obj_.printlevel == 1)
	disp('objective function:');
	disp(phi);
	disp('gradient:');
	disp(grad);
  end

end %LikelihoodA
%% END OF FUNCTION LIKELIHOODA.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function retval = hineq(theta)
%% provides h_ineq(x) for inequality constraint h_ineq(x) >= 0
%% to ensure that parameters are between lower and upper bound
%% parlb <= theta <= parub
  global obj_
  retval = [theta - obj_.parlb; -theta + obj_.parub];
end %hineq
%% END OF FUNCTION hineq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function retval = gradi(theta)
%% optimizer should have just called the function likelihood with the same
%% theta value, and likelihood has computed obj_.gradlogtr.
%% So pull the gradient out of the global variables and return
  global obj_ likelihood_last_theta mod_ solveroption

	if (~ isequal (theta, likelihood_last_theta))
		if(solveroption == 3)
		%% resolve model to update cache
		obj = likelihoodA (theta);
		else
		obj = likelihood (theta);
		end
	end
  retval = obj_.gradlogtr';
end %gradi
%% END OF FUNCTION gradi %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xdot, flag, new_data] = oderhs(t, x, data)
%% Right hand side of differential equation
%% dx/dt = f(x,t,theta)
  global mod_

  [xdot] = mod_.odefcn(x, t, data.theta);
  flag = 0;
  new_data = [];
end %oderhs
%% END OF FUNCTION oderhs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xsd, flag, new_data] = sensrhs(t, x, xdot, xs, data)
%% Right hand side of forward sensitivity equation
%% sx = dx/dp
%% dsx / dt = df/dx * sx + df/dp
  global mod_ obj_

  dfdp = mod_.dodedp(x, t, data.theta);
  dfdp = dfdp(:,obj_.estflag);
  xsd = mod_.dodedx(x, t, data.theta)*xs + dfdp;
  flag = 0;
  new_data = [];
end %sensrhs
%% END OF FUNCTION sensrhs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [qbd, flag, new_data] = sensrhsA (t, x, mu, data)
%% Backward quadrature function use to compute gradient of objective
%% function.  In this problem, it computes Equation 3.19 from the CVodes
%% user guide:
%%   dg/dp(t1) = mu(t0)s(t0) + g_p(t1) + \int_t0^t1 mu'*f_p
%% 
%% compute the integral by integrating
%%    dw/dt = - mu' f_p
%%    w(t1) = 0
%% 
%%  backward over the time horizon.
  global mod_ obj_

  dfdp = mod_.dodedp(x, t, data.theta);
  dfdp = dfdp(:,obj_.estflag);
  qbd = -(mu'*dfdp)';
  flag = 0;
  new_data = [];
end %sensrhsA
%% END OF FUNCTION sensrhsA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [mudot, flag, new_data] = adjrhs (t, x, mu, data)
%% RHS of adjoint system.  In this problem, it computes Equation 3.20
%% from the CVodes user guide:
%%    mudot = - f_x'*mu
%%    mu(t1) = g_x' at t=t1
  global mod_

  mudot = -mod_.dodedx(x, t, data.theta)'*mu;
  flag = 0;
  new_data = [];
end %adjrhs
%% END OF FUNCTION adjrhs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Data files

ABC.dat

# Created by Octave 3.6.2, Tue May 14 00:33:52 2013 CDT 
# name: table
# type: matrix
# rows: 100
# columns: 4
 0 1 0 0
 0.06060606060606061 0.8848233223536865 0.1117109106715674 0.00346576674982394
 0.1212121212121212 0.7829123067621313 0.2040308391257652 0.01305685333214935
 0.1818181818181818 0.6927390515446727 0.2795742520067726 0.02768669699743825
 0.2424242424242424 0.6129516533762045 0.3406321955366671 0.04641615342224718
 0.303030303030303 0.5423539003189545 0.3892109152113909 0.06843518707418547
 0.3636363636363636 0.4798873783706807 0.4270658117792184 0.09304681178898236
 0.4242424242424243 0.4246155536122287 0.4557315732483217 0.1196528761795098
 0.4848484848484849 0.3757097607235237 0.4765486222849156 0.147741618139771
 0.5454545454545454 0.3324367605108264 0.4906865609223258 0.1768766759613666
 0.6060606060606061 0.2941478026148869 0.4991647030456889 0.2066874912256227
 0.6666666666666667 0.2602688382145835 0.5028703841751919 0.2368607725058967
 0.7272727272727273 0.2302919403461645 0.5025749852091307 0.2671330684236419
 0.7878787878787878 0.2037676819150678 0.4989480923810672 0.2972842185720354
 0.8484848484848485 0.1802983990466281 0.4925699809054059 0.3271316119784785
 0.9090909090909092 0.1595322300904848 0.4839426098832218 0.3565251511700092
 0.9696969696969697 0.1411578392275535 0.4734993153967521 0.3853428358349933
 1.03030303030303 0.1248997494932536 0.4616133450121498 0.4134868959757724
 1.090909090909091 0.1105142123319374 0.4486053718024498 0.440880406332665
 1.151515151515152 0.09778555347556032 0.4347501109952634 0.4674643259897164
 1.212121212121212 0.08652293917381268 0.4202821400363677 0.4931949112595115
 1.272727272727273 0.07655751526901554 0.4054010175850353 0.5180414576100107
 1.333333333333333 0.06773987569066836 0.3902757841788888 0.5419843305974378
 1.393939393939394 0.05993782245920518 0.3750489171865766 0.5650132508198895
 1.454545454545455 0.05303438371934838 0.3598398046140462 0.5871258021327808
 1.515151515151515 0.04692606004373614 0.344747794544232 0.6083261358780804
 1.575757575757576 0.04152127529637509 0.3298548658708658 0.6286238492982726
 1.636363636363636 0.03673899902808295 0.3152279826728423 0.6480330087641719
 1.696969696969697 0.03250752766784495 0.3009211569053662 0.666571305891807
 1.757575757575758 0.0287634215447189 0.2869772260723309 0.6842593428479988
 1.818181818181818 0.0254505487295796 0.273429437021086 0.7011200047143696
 1.878787878787879 0.0225192420879477 0.2603028248742175 0.7171779235028554
 1.939393939393939 0.01992555274887522 0.2476154182449457 0.7324590194712091
 2 0.01763059565213198 0.2353792797107512 0.7469901151021496
 2.060606060606061 0.01559996409694179 0.2236014239799888 0.760798602388098
 2.121212121212121 0.01380321351448919 0.2122846137534227 0.773912163197117
 2.181818181818182 0.01221340668503915 0.2014280478927336 0.7863585358872579
 2.242424242424243 0.01080670818611693 0.1910279641481316 0.7981653181307816
 2.303030303030303 0.00956202859312319 0.1810781561130601 0.8093598057588465
 2.363636363636364 0.008460706864600162 0.1715704250203074 0.8199688585801224
 2.424242424242424 0.007486231561312076 0.1624949658209353 0.8300187930827828
 2.484848484848485 0.006623993071107175 0.1538407014702757 0.8395352959236472
 2.545454545454545 0.005861064157564161 0.1455955703433666 0.8485433559640994
 2.606060606060606 0.005186007334546378 0.1377467698179541 0.8570672133125298
 2.666666666666667 0.004588703355060324 0.1302809658908966 0.8651303212190734
 2.727272727272728 0.004060197123240804 0.123184475085305 0.8727553182564846
 2.787878787878788 0.003592562735229833 0.1164434162746161 0.8799640114551844
 2.848484848484849 0.003178787295816458 0.110043833410022 0.8867773697591919
 2.909090909090909 0.002812667507365943 0.1039717991115374 0.893215523846127
 2.96969696969697 0.002488715704403829 0.09821350461439358 0.8992977701462329
 3.03030303030303 0.002202075998794117 0.09275533510864685 0.9050425793575894
 3.090909090909091 0.0019484505319788 0.08758393047485338 0.9104676094581982
 3.151515151515151 0.00172403636848565 0.08268623093772329 0.9155897231588214
 3.212121212121212 0.001525469017260709 0.07804951394239697 0.9204250075053726
 3.272727272727273 0.001349771782547499 0.07366142369277344 0.9249887949897093
 3.333333333333333 0.001194310723615486 0.06950999366184318 0.9292956860795717
 3.393939393939394 0.001056755081575004 0.06558366229547975 0.9333595730879757
 3.454545454545455 0.000935042495305569 0.06187128300402878 0.9371936649656961
 3.515151515151515 0.0008273482036626089 0.05836212972212509 0.9408105125392426
 3.575757575757576 0.0007320576735484587 0.05504589890307102 0.9442220338884109
 3.636363636363637 0.0006477423251742669 0.05191270842770073 0.9474395397121553
 3.696969696969697 0.000573138083491335 0.04895309381353757 0.9504737585680015
 3.757575757575758 0.0005071264462787704 0.04615800217331581 0.9533348618454358
 3.818181818181818 0.0004487177434942088 0.0435187843949971 0.9560324883265391
 3.878787878787879 0.0003970363054917097 0.0410271859408679 0.9585757682186709
 3.939393939393939 0.0003513073221811402 0.03867533655544366 0.9609733465874056
 4 0.0003108452166740596 0.03645573910429449 0.9632334061440618
 4.060606060606061 0.000275043374996885 0.03436125774581806 0.9653636893442155
 4.121212121212121 0.0002433650449941336 0.03238510564746584 0.9673715197725705
 4.181818181818182 0.0002153352864188676 0.03052083239941569 0.969263822779196
 4.242424242424242 0.0001905338753260695 0.02876231126695965 0.9710471453227448
 4.303030303030303 0.0001685889856266183 0.02710372640871676 0.9727276750706871
 4.363636363636363 0.0001491716150808266 0.02553956012480547 0.9743112587251442
 4.424242424242425 0.0001319908614232179 0.02406458028617328 0.9758034193174341
 4.484848484848485 0.0001167890976768277 0.02267382778599436 0.9772093735813595
 4.545454545454546 0.0001033389882104947 0.02136260462685521 0.978534046849965
 4.606060606060606 9.14378805278004e-05 0.0201264612500457 0.9797820913344573
 4.666666666666667 8.090757478514186e-05 0.01896118576718694 0.9809578971230587
 4.727272727272728 7.158994492727467e-05 0.01786279226414261 0.982065608255961
 4.787878787878788 6.334501357293608e-05 0.01682750962108882 0.9831091358303692
 4.848484848484849 5.604959178138912e-05 0.01585177199981472 0.9840921688734349
 4.909090909090909 4.959427931086259e-05 0.01493220748355515 0.985018188702165
 4.96969696969697 4.388242543173056e-05 0.01406562913079985 0.9858904789087994
 5.03030303030303 3.882861272312608e-05 0.01324902565188155 0.9867121362004263
 5.090909090909091 3.435679684879674e-05 0.01247955135604449 0.9874860823121377
 5.151515151515151 3.040002373996354e-05 0.01175451886969726 0.9882150715715938
 5.212121212121212 2.689895752361209e-05 0.01107139019918357 0.9889017013083239
 5.272727272727272 2.38010032332678e-05 0.01042776851226528 0.9895484209495325
 5.333333333333334 2.105985210702013e-05 0.009821391758337066 0.9901575388545871
 5.393939393939394 1.863438081141093e-05 0.009250124277957012 0.9907312318062628
 5.454545454545455 1.64882442193286e-05 0.008711950747443038 0.9912715514733689
 5.515151515151516 1.458934124994323e-05 0.008204969797741716 0.9917804313260395
 5.575757575757576 1.290910709951106e-05 0.007727386892669207 0.9922596944652625
 5.636363636363637 1.14223862348502e-05 0.007277509614456029 0.9927110584643404
 5.696969696969697 1.010690303213876e-05 0.006853741559190599 0.9931361420028085
 5.757575757575758 8.942889794258136e-06 0.006454576649864608 0.9935364709253723
 5.818181818181818 7.912945309928674e-06 0.006078595112405932 0.9939134824073154
 5.878787878787879 7.001615058067173e-06 0.005724457755399133 0.994268531094574
 5.939393939393939 6.195236896463615e-06 0.005390901975323785 0.9946028932528109
 6 5.481749089223554e-06 0.005076737636557397 0.9949177710793846


ABC_data.dat

   0.000   0.957  -0.031  -0.015
   0.263   0.557   0.330   0.044
   0.526   0.342   0.512   0.156
   0.789   0.224   0.499   0.310
   1.053   0.123   0.428   0.454
   1.316   0.079   0.396   0.556
   1.579   0.035   0.303   0.651
   1.842   0.029   0.287   0.658
   2.105   0.025   0.221   0.750
   2.368   0.017   0.148   0.854
   2.632  -0.002   0.182   0.845
   2.895   0.009   0.116   0.893
   3.158  -0.023   0.079   0.942
   3.421   0.006   0.078   0.899
   3.684   0.016   0.059   0.942
   3.947   0.014   0.036   0.991
   4.211  -0.009   0.014   0.988
   4.474  -0.030   0.036   0.941
   4.737   0.004   0.036   0.971
   5.000  -0.024   0.028   0.985