Figure 9.34 (page 553):

Concentration of species B versus time.

Code for Figure 9.34

Text of the GNU GPL.

main.m

%% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%%
%% 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.

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 vrdiclo k1 k2 vro teaf flowrate tflow xsmall idx fcn sx0 ...
    timeout tcrit x xdot sx sxdot sens lcmeas resid grad Hgn ...
    lcpred theta
%%
%% global variables                     
%%
%% vrdiclo                      initial moles of dichloro
%% k1                           first rate constant cu dm/kgmol min
%% k2                           second rate constant cu dm/kgmol min
%% vro	                        initial reactor volume  cu dm
%% teef	                        TEA feed concentration kgmol/cu dm
%% teaden                       TEA density kg/cu dm
%% tflow, flowrate              measured flowrate at time tflow
%%                              read from "flow.dat" containing time
%%                              (min) and flowrate (kg/min) converted to
%%                              (cu dm/min)
%% tlc, lcmeas                  lc measurement at time tlc (min)
%%

if (~ exist ('CVode'))
  if (exist ('startup_STB'))
    startup_STB;
  else
    error ('%s: requires the Sundials Toolbox, please see the README file for installation instructions');
  end
end

xsmall    = 1e-6;
%% reaonably optimal values
%% phi is 0.0011810 for these
% vrdiclo     =  2.3497
% k1          = 2200.
% k2          = 1120.
%%
%% other initial guesses tried
% vrdiclo   = 2.38;
% k1        = 2000;
% k2        = k1/3;
% vrdiclo   = 2.33;
% k1        = 2117;
% k2        = 690;
% vrdiclo   = 2.33;
% k1        = 2040;
% k2        = 682;
%%
%% Sankash's optimal parameters, 5/29/00
%% these would not solve for me.
% vrdiclo =      2.352674877
% k1      = 8360353635
% k2      = 4025874934
%Obj.    0.000622978
%% Sanakash's best 2 parameter fit
vrdiclo = 2.350634;
k2      = 1250;
k1      = 1.997264*k2;
%Obj.    0.00113
vro       = 2370;
teaf      = 0.00721;
teaden    = 0.728;
flow = load ('flow.dat')
lc = load ('lc.dat')
tflow  = [0. ; flow(:,1)];
flowrate = [0. ; flow(:,2)./teaden];
tlc    = lc(:,1);
lcmeas = lc(:,2);

theta = [vrdiclo; k1; k2];

%% this version depends on the following functions defined in .cc files
%%
% function res(x, time, theta)
% function jac = j_fcn(x, time, theta)
% function bcol = b_fcn (x, t, theta, col)

%tcrit = tflow(end);
tcrit = flow(:,1);

ntimes    = 200;
tlin = linspace (0, tflow(end), ntimes)';
timeout = sort (unique ([tlin; tflow; tlc]));

%% I can't think of a way to do this without a loop.
idx = zeros (size (tlc));
for i = 1:length (tlc)
  %% This assignment will fail if find returns more than one element,
  %% but that's OK, because each element of tlc should only appear once
  %% in timeout.
  idx(i) = find (timeout == tlc(i));
end

%% sanity check
if (any (timeout (idx) ~= tlc))
  error ('extracting index vector for tlc');
end



model.cvodesfcn = @res_revised;
model.tplot = timeout;
model.ic = [vro; xsmall; xsmall; xsmall; theta(1)];
model.stateic_par = zeros(size(model.ic));
model.stateic_par(5) = 1;
model.param = theta;

measure.data = lcmeas';
measure.statefcn = @lcfcn;
measure.dstatedx = @dlcdx;
measure.time = tlc;

objective.estflag = [1;2;3];
objective.parlb   = zeros(size(objective.estflag));
objective.parub   = 1e20*ones(size(objective.estflag));
objective.paric   = theta(objective.estflag);
objective.printlevel = 0;

estimates = parest(model, measure, objective);

estimates.info
estimates.obj
estimates.parest
estimates.bbox

table = [model.tplot, estimates.x', estimates.y'];

figure(1);
%plot measurement and fit
plot (measure.time, measure.data, '*', model.tplot, estimates.y)
axis ([0,model.tplot(end),0,0.35]);
figure(2);
%plot c_A, c_C, c_D
plot (model.tplot, estimates.x([2,4,5],:))
figure(3);
%plot c_B
plot (model.tplot, estimates.x(3,:))
title ('Figure 9.34')

lcfcn.m

function retval = lcfcn(x)
  d = x(4) + 2*x(2);
  if (d == 0)
    retval = 1;
  else
    retval = x(4) / d;
  end

dlcdx.m

function retval = dlcdx(x)
  dlcdx = zeros(size(x'));
  d = (x(4)+2*x(2))^2;
  if (d != 0)
    dlcdx(2) = -2*x(4)/d;
    dlcdx(4) = 2*x(2)/d;
  end
  retval = dlcdx;

../../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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Octave dynamically-linked functions (C++ language)

res_revised.cc

/*

Copyright (C) 2001, James B. Rawlings

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.

*/

#include 
#include 
#include 
#include 

DEFUN_DLD (res_revised, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {[@var{xdot}, @var{flag}, @var{new_data}] =} res_revised (@var{t}, @var{x}, @var{data})\n\
ODE RHS function for CVODES for BVSM model.\n\
@end deftypefn")

//	Passed variable                  Local equivalent
//
//	x(1)                             reactor contents volume
//	x(2)                             total BVSM moles
//	x(3)                             total TEA  moles
//	x(4)                             total monochloro moles
//	x(5)                             total dichloro moles

{
  octave_value_list retval;
  retval(2) = Matrix ();
  retval(1) = 0.0;
  retval(0) = 0.0;
  if (args.length () == 3)
    {
      double time  = args(0).double_value ();
      Array x = args(1).vector_value ();
      octave_scalar_map data = args(2).scalar_map_value ();
      if (! error_state)
	{
	if (x.numel () == 5)
	  {
  	    double vr = x(0);
	    if (vr <= 0.0)
	    {
	      error ("(res_revised).... vr is zero, time= %f", time);
	      return retval;
	    }
	  double teaf = get_global_value ("teaf") . double_value ();
	  ColumnVector tflow    (get_global_value ("tflow") . vector_value ());
	  ColumnVector flowrate (get_global_value ("flowrate") . vector_value ());
	  if (error_state)
	    {	  
	      error ("res_revised: unable to locate global variables");
	      return retval;
	    }
	  double tea    = x(2) / vr;
	  double monocl = x(3) / vr;
	  double dicl   = x(4) / vr;

          Cell theta_cell = data.contents ("theta");
          if (theta_cell.numel () == 1)
            {
              Array theta = theta_cell(0).vector_value ();
              if (theta.numel () >= 3)
                {
                  double k1 = theta(1);
                  double k2 = theta(2);
                  // Set up local parameter names.
                  double k1bdvr = k1 * tea * dicl * vr;
                  double k2bmvr = k2 * tea * monocl * vr;
                  // Look up flowrate at current time.  Do it by calling the M-file
                  // (or .oct file) find_time_index.
                  octave_value_list tmp_args;
                  tmp_args(1) = time;
                  tmp_args(0) = tflow;
                  int tmp_nargout = 2;
                  octave_value_list tmp_result = feval ("find_time_index", tmp_args,
                                                        tmp_nargout);
                  if (! error_state)
                    {
                      int index = tmp_result(0) . int_value ();
                      double fteaf = flowrate(index-1) * teaf;
                      ColumnVector xdot (5);
                      xdot(0) = flowrate(index-1);
                      xdot(1) = k2bmvr;
                      xdot(2) = fteaf - (k1bdvr + k2bmvr);
                      xdot(3) = k1bdvr - k2bmvr;
                      xdot(4) = -k1bdvr;
                      retval(0) = xdot;
                    }
                }
              else
                error ("res_revised: expecting theta to be a 3-element vector");
            }
          else
            error ("res_revised: data.theta is invalid");
          }
	else
	  error ("res_revised: expecting x to be a 5-element vector");
	}
      else
	error ("res_revised: invalid argument values");
    }
  else
    print_usage ();

  if (error_state)
    retval(1) = -1.0;

  return retval;
}

Data files

lc.dat

    414        0.17050
    424        0.16040
    434        0.13250
    444        0.10840
    493        0.10140
    503        0.10420
    513        0.10450
    523        0.09700
    533        0.08240
    582        0.06780
    592        0.06970
    602        0.07020
    612        0.06920
    621        0.06840
    631        0.06750
    641        0.06730
    651        0.06730
    661        0.06540
    671        0.05770
    681        0.05710
    690        0.05630
    700        0.04660
    710        0.04550
    720        0.04490
    730        0.04560
    740        0.04500
    750        0.04340
    759        0.04280
    769        0.04260
    779        0.04090
    789        0.03990
    799        0.03970
    848        0.03940
    858        0.03740
    868        0.03710

flow.dat

   9.000000      9.7817838E-02   4.786000    
   19.00000       1.290601       14.21800    
   29.00000       1.262779       13.99800    
   39.00000       1.240016       13.81800    
   49.00000       1.215736       13.62600    
   59.00000       1.187156       13.40000    
   69.00000       1.160599       13.19000    
   79.00000      0.5545961       8.398001    
   89.00000       1.109508       12.78600    
   99.00000       1.361925       14.78200    
   109.0000       1.332333       14.54800    
   119.0000       1.294648       14.25000    
   129.0000       1.271885       14.07000    
   139.0000       1.244569       13.85400    
   149.0000       1.220794       13.66600    
   159.0000       1.197778       13.48400    
   169.0000       1.170463       13.26800    
   179.0000       1.150988       13.11400    
   189.0000       1.128983       12.94000    
   199.0000       1.105968       12.75800    
   209.0000       1.092563       12.65200    
   219.0000       1.360155       14.76800    
   229.0000       1.344473       14.64400    
   239.0000       1.329045       14.52200    
   249.0000       1.305524       14.33600    
   259.0000       1.275426       14.09800    
   269.0000       1.259239       13.97000    
   279.0000       1.244063       13.85000    
   289.0000       1.219783       13.65800    
   299.0000       1.191961       13.43800    
   309.0000       1.173751       13.29400    
   319.0000       1.158322       13.17200    
   329.0000      0.3163429       6.514000    
   339.0000      0.5495375       8.358001    
   349.0000       1.224841       13.69800    
   359.0000      0.3228684       6.565600    
   369.0000      9.4845297E-04   4.020000    
   379.0000      4.6794413E-04   4.016201    
   389.0000      5.4383278E-04   4.016800    
   399.0000      4.9326418E-04   4.016400    
   409.0000      3.6680698E-04   4.015400    
   419.0000      0.2701340       6.148601    
   429.0000       1.185891       13.39000    
   439.0000      0.9299334       11.36600    
   449.0000      0.5259147       8.171202    
   459.0000      8.9793204E-04   4.019600    
   469.0000      8.7256433E-04   4.019400    
   479.0000      7.9672335E-04   4.018800    
   489.0000      4.9319270E-04   4.016400    
   499.0000      4.6799183E-04   4.016200    
   509.0000      3.6678315E-04   4.015400    
   519.0000      0.2330809       5.855600    
   529.0000      0.7784327       10.16800    
   539.0000      7.5244429E-03   4.072000    
   549.0000      7.5244186E-03   4.072000    
   559.0000      5.7539940E-03   4.058000    
   569.0000      8.0302954E-03   4.076001    
   579.0000      0.5019882       7.982001    
   589.0000     -3.0983209E-03   3.988000    
   599.0000      9.4850064E-04   4.020000    
   609.0000     -6.3252446E-05   4.012000    
   619.0000     -1.3278484E-03   4.002000    
   629.0000     -4.6158792E-03   3.976000    
   639.0000     -6.3204767E-05   4.012000    
   649.0000     -2.0865917E-03   3.996000    
   659.0000      0.1752122       5.398000    
   669.0000      0.3024322       6.404000    
   679.0000      4.9952744E-03   4.052001    
   689.0000      3.9835931E-03   4.044001    
   699.0000      0.5114223       8.056601    
   709.0000      3.4151078E-04   4.015200    
   719.0000      3.1619071E-04   4.015000    
   729.0000      1.1386871E-04   4.013400    
   739.0000      6.3300133E-05   4.013000    
   749.0000      3.4151078E-04   4.015200    
   759.0000      0.1059619       4.850400    
   769.0000      6.7019463E-04   4.017800    
   779.0000      6.8099439E-02   4.551000    
   789.0000      8.7261200E-04   4.019400    
   799.0000      9.7372534E-04   4.020200    
   809.0000      8.7256433E-04   4.019400    
   819.0000      7.9669955E-04   4.018800    
   829.0000      7.2083471E-04   4.018200    
   839.0000      5.6912901E-04   4.017000    
   849.0000      4.6801567E-04   4.016200    
   859.0000      5.1856041E-04   4.016600    
   869.0000      4.4267176E-04   4.016000