% 
% This is an example on how to run an parameter fitting task.
% The example creates a simple model and runs a time course simulation on it.
% The timecourse data is written to file with some noise added to it.
% This data is used to fit the original parameters.
% 

COPASI

MODEL_STRING="<?xml version=\"1.0\" encoding=\"UTF-8\"?>\
<!-- Created by COPASI version 4.5.30 (Debug) on 2009-03-30 08:01 with libSBML version 3.3.2. -->\
<sbml xmlns=\"http://www.sbml.org/sbml/level2\" level=\"2\" version=\"1\">\
  <model metaid=\"COPASI1\" id=\"Model_1\" name=\"Model\">\
    <listOfUnitDefinitions>\
      <unitDefinition id=\"volume\">\
        <listOfUnits>\
          <unit kind=\"litre\" scale=\"-3\"/>\
        </listOfUnits>\
      </unitDefinition>\
      <unitDefinition id=\"substance\">\
        <listOfUnits>\
          <unit kind=\"mole\" scale=\"-3\"/>\
        </listOfUnits>\
      </unitDefinition>\
    </listOfUnitDefinitions>\
    <listOfCompartments>\
      <compartment id=\"compartment_1\" name=\"compartment\" size=\"1\"/>\
    </listOfCompartments>\
    <listOfSpecies>\
      <species id=\"species_1\" name=\"A\" compartment=\"compartment_1\" initialConcentration=\"5\"/>\
      <species id=\"species_2\" name=\"B\" compartment=\"compartment_1\" initialConcentration=\"0\"/>\
      <species id=\"species_3\" name=\"C\" compartment=\"compartment_1\" initialConcentration=\"0\"/>\
    </listOfSpecies>\
    <listOfReactions>\
      <reaction id=\"reaction_1\" name=\"reaction\" reversible=\"false\">\
        <listOfReactants>\
          <speciesReference species=\"species_1\"/>\
        </listOfReactants>\
        <listOfProducts>\
          <speciesReference species=\"species_2\"/>\
        </listOfProducts>\
        <kineticLaw>\
          <math xmlns=\"http://www.w3.org/1998/Math/MathML\">\
            <apply>\
              <times/>\
              <ci> compartment_1 </ci>\
              <ci> k1 </ci>\
              <ci> species_1 </ci>\
            </apply>\
          </math>\
          <listOfParameters>\
            <parameter id=\"k1\" name=\"k1\" value=\"0.03\"/>\
          </listOfParameters>\
        </kineticLaw>\
      </reaction>\
      <reaction id=\"reaction_2\" name=\"reaction_1\" reversible=\"false\">\
        <listOfReactants>\
          <speciesReference species=\"species_2\"/>\
        </listOfReactants>\
        <listOfProducts>\
          <speciesReference species=\"species_3\"/>\
        </listOfProducts>\
        <kineticLaw>\
          <math xmlns=\"http://www.w3.org/1998/Math/MathML\">\
            <apply>\
              <times/>\
              <ci> compartment_1 </ci>\
              <ci> k1 </ci>\
              <ci> species_2 </ci>\
            </apply>\
          </math>\
          <listOfParameters>\
            <parameter id=\"k1\" name=\"k1\" value=\"0.004\"/>\
          </listOfParameters>\
        </kineticLaw>\
      </reaction>\
    </listOfReactions>\
  </model>\
</sbml>";

%!assert( swig_this(CRootContainer.getRoot()) != 0);
% create a datamodel
dataModel = CRootContainer.addDatamodel();
%!assert( CRootContainer.getDatamodelList().size() == 1);
% first we load a simple model
% load the model 
if (dataModel.importSBMLFromString(MODEL_STRING) == false)
  error("Error while importing the model.\n");
  quit(1)
endif
% now we need to run some time course simulation to get data to fit
% against

% get the trajectory task object
trajectoryTask = dataModel.getTask("Time-Course");
%!assert( swig_this(trajectoryTask) != 0);
% if there isn't one
if (swig_this(trajectoryTask) == 0)
  % create a one
  trajectoryTask = CTrajectoryTask();

  % add the time course task to the task list
  % this method makes sure that the object is now owned 
  % by the list and that it does not get deleted by SWIG
  dataModel.getTaskList().addAndOwn(trajectoryTask);
endif

% run a deterministic time course
trajectoryTask.setMethodType(CTaskEnum_Method_deterministic);

% pass a pointer of the model to the problem
trajectoryTask.getProblem().setModel(dataModel.getModel());

% activate the task so that it will be run when the model is saved
% and passed to CopasiSE
trajectoryTask.setScheduled(true);

% get the problem for the task to set some parameters
problem = trajectoryTask.getProblem();

% simulate 4000 steps
problem.setStepNumber(4000);
% start at time 0
dataModel.getModel().setInitialTime(0.0);
% simulate a duration of 400 time units
problem.setDuration(400);
% tell the problem to actually generate time series data
problem.setTimeSeriesRequested(true);

% set some parameters for the LSODA method through the method
method = trajectoryTask.getMethod();

result=true;
try
    % now we run the actual trajectory
    result=trajectoryTask.process(true);
catch
    % check if there are additional error messages
    printf("%s\n",trajectoryTask.getProcessWarning());
    printf("%s\n",trajectoryTask.getProcessError());
    error("Error. Running the time course simulation failed.\n");
    quit(1)
end_try_catch

if(result == false)
    % check if there are additional error messages
    dataModel.saveModel('example6.cps', true);
    printf("%s\n",trajectoryTask.getProcessWarning());
    printf("%s\n",trajectoryTask.getProcessError());
    error("An error occured while running the time course simulation.\n");
    quit(1)
endif

% we write the data to a file and add some noise to it
% This is necessary since COPASI can only read experimental data from
% file.
timeSeries = trajectoryTask.getTimeSeries();
% we simulated 100 steps, including the initial state, this should be
% 101 step in the timeseries
%!assert( timeSeries.getRecordedSteps() == 4001);
iMax = timeSeries.getNumVariables()-1;
% there should be four variables, the three metabolites and time
%!assert( iMax == 5);
lastIndex = timeSeries.getRecordedSteps() - 1;
% open the file
% we need to remember in which order the variables are written to file
% since we need to specify this later in the parameter fitting task
indexSet={};
metabVector={};

% write the header
% the first variable in a time series is a always time, for the rest
% of the variables, we use the SBML id in the header
r=0.0;
os=fopen("fakedata_example6.txt","w");
fprintf(os,"# time ");
keyFactory=CRootContainer.getKeyFactory();
%!assert( swig_this(keyFactory) != 0);
for i = 1:iMax
  key=timeSeries.getKey(i);
  object=keyFactory.get(key);
  %!assert( swig_this(object) != 0);
  % only write header data or metabolites
  if ( strcmp(swig_type(object),"CMetab") == 1 )
    fprintf(os,", ");
    fprintf(os,"%s",timeSeries.getSBMLId(i,dataModel));
    indexSet{length(indexSet)+1}=i;
    metabVector{length(metabVector)+1}=object;
  endif
endfor
fprintf(os,"\n");
data=0.0;
for i = 0:lastIndex-1
  s="";
  for j = 0:iMax-1
    % we only want to  write the data for metabolites
    % the compartment does not interest us here

    % Could not find an easy way to check for the existence 
    % of an element in a cell array, so we do it manually
    cMax=length(indexSet);
    found=false;
    for e = 1:cMax
      if (indexSet{e} == j)
        found=true;
        break;
      endif
    endfor
    if (j==0 || (found == true))
      % write the data with some noise (+-5% max)
      r=rand();
      data=timeSeries.getConcentrationData(i, j);
      % don't add noise to the time
      if j != 0
        data=data+data*(r*0.1-0.05);
      endif
      s=strcat(s,num2str(data),", ");
    endif
  endfor  
  % remove the last two characters again
  fprintf(os,"%s",substr(s,1,length(s)-2));
  fprintf(os,"\n");
endfor
fclose(os);


% now we change the parameter values to see if the parameter fitting
% can really find the original values
r=rand()*10;
reaction=dataModel.getModel().getReaction(0);
% we know that it is an irreversible mass action, so there is one
% parameter
%!assert( reaction.getParameters().size() == 1);
%!assert( reaction.isLocalParameter(0) == true);
% the parameter of a irreversible mass action is called k1
reaction.setParameterValue("k1",r);

reaction=dataModel.getModel().getReaction(1);
% we know that it is an irreversible mass action, so there is one
% parameter
%!assert( reaction.getParameters().size() == 1);
%!assert( reaction.isLocalParameter(0) == true);
reaction.setParameterValue("k1",r);

fitTask=dataModel.addTask(CTaskEnum_Task_parameterFitting);
%!assert( swig_this(fitTask) != 0);
% the method in a fit task is an instance of COptMethod or a subclass of
% it.
fitMethod=fitTask.getMethod();
%!assert( swig_this(fitMethod) != 0);
% the object must be an instance of COptMethod or a subclass thereof
% (CFitMethod)
fitProblem=fitTask.getProblem();
%!assert( swig_this(fitProblem) != 0);

experimentSet=fitProblem.getParameter("Experiment Set");
%!assert( swig_this(experimentSet) != 0);

% first experiment (we only have one here)
experiment=CExperiment(dataModel);
%!assert( swig_this(experiment) != 0);
% tell COPASI where to find the data
% reading data from string is not possible with the current C++ API
experiment.setFileName("fakedata_example6.txt");
% we have to tell COPASI that the data for the experiment is a komma
% separated list (the default is TAB separated)
experiment.setSeparator(",");
% the data start in row 1 and goes to row 4001
experiment.setFirstRow(1);
%!assert( experiment.getFirstRow()==1);
experiment.setLastRow(4001);
%!assert( experiment.getLastRow()==4001);
experiment.setHeaderRow(1);
%!assert( experiment.getHeaderRow()==1);
experiment.setExperimentType(CTaskEnum_Task_timeCourse);
%!assert( experiment.getExperimentType() == CCopasiTask_timeCourse);
experiment.setNumColumns(4);
%!assert( experiment.getNumColumns()==4);
objectMap=experiment.getObjectMap();
%!assert( swig_this(objectMap) != 0);
result=objectMap.setNumCols(4);
%!assert( result == true);
result=objectMap.setRole(0,CExperiment_time);
%!assert( result == true);
%!assert( objectMap.getRole(0) == CExperiment_time);

model=dataModel.getModel();
%!assert( swig_this(model) != 0);
timeReference=model.getValueReference();
%!assert( swig_this(timeReference) != 0);
objectMap.setObjectCN(0,timeReference.getCN().getString());

% now we tell COPASI which column contain the concentrations of
% metabolites and belong to dependent variables
objectMap.setRole(1,CExperiment_dependent);
metab=metabVector{1};
%!assert( swig_this(metab) != 0);
particleReference=metab.getConcentrationReference();
%!assert( swig_this(particleReference) != 0);
objectMap.setObjectCN(1,particleReference.getCN().getString());

objectMap.setRole(2,CExperiment_dependent);
metab=metabVector{2};
%!assert( swig_this(metab) != 0);
particleReference=metab.getConcentrationReference();
%!assert( swig_this(particleReference) != 0);
objectMap.setObjectCN(2,particleReference.getCN().getString());

objectMap.setRole(3,CExperiment_dependent);
metab=metabVector{3};
%!assert( swig_this(metab) != 0);
particleReference=metab.getConcentrationReference();
%!assert( swig_this(particleReference) != 0);
objectMap.setObjectCN(3,particleReference.getCN().getString());

experimentSet.addExperiment(experiment);
%!assert( experimentSet.getExperimentCount()==1);
% addExperiment makes a copy, so we need to get the added experiment
% again
experiment=experimentSet.getExperiment(0);
%!assert( swig_this(experiment) != 0);

% now we have to define the two fit items for the two local parameters
% of the two reactions
reaction=model.getReaction(0);
%!assert( swig_this(reaction) != 0);
%!assert( reaction.isLocalParameter(0)==true);
parameter=reaction.getParameters().getParameter(0);
%!assert( swig_this(parameter) != 0);

% define a CFitItem
parameterReference=parameter.getValueReference();
%!assert( swig_this(parameterReference) != 0);
fitItem1=CFitItem(dataModel);
%!assert( swig_this(fitItem1) !=0);
fitItem1.setObjectCN(parameterReference.getCN());
fitItem1.setStartValue(4.0);
fitItem1.setLowerBound(CCommonName("0.00001"));
fitItem1.setUpperBound(CCommonName("10"));
% add the fit item to the correct parameter group
optimizationItemGroup=fitProblem.getParameter("OptimizationItemList");
%!assert( swig_this(optimizationItemGroup) != 0);
optimizationItemGroup.addParameter(fitItem1);

reaction=model.getReaction(1);
%!assert( swig_this(reaction) != 0);
%!assert( reaction.isLocalParameter(0) == true);
parameter=reaction.getParameters().getParameter(0);
%!assert( swig_this(parameter) != 0);

% define a CFitItem
parameterReference=parameter.getValueReference();
%!assert( swig_this(parameterReference) != 0);
fitItem2=CFitItem(dataModel);
%!assert( swig_this(fitItem2) != 0);
fitItem2.setObjectCN(parameterReference.getCN());
fitItem2.setStartValue(4.0);
fitItem2.setLowerBound(CCommonName("0.00001"));
fitItem2.setUpperBound(CCommonName("10"));
% add the fit item to the correct parameter group
optimizationItemGroup.addParameter(fitItem2);

result=true;

try
  % running the task for this example will probably take some time
  printf("This can take some time...\n");
  result=fitTask.process(true);
catch
  error("Error. Parameter fitting failed.\n");
  quit(1);
end_try_catch

%!assert( result == true);
% assert that there are two optimization items
%!assert( length(fitProblem.getOptItemList()) == 2);
% the order should be the order in whih we added the items above
optItem1 = fitProblem.getOptItemList(){1};
optItem2 = fitProblem.getOptItemList(){2};
% the actual results are stored in the fit problem
%!assert( fitProblem.getSolutionVariables().size() == 2);
printf("value for %s: %f\n" , optItem1.getObject().getCN().getString() , fitProblem.getSolutionVariables().get(0));
printf("value for %s: %f\n" , optItem2.getObject().getCN().getString(), fitProblem.getSolutionVariables().get(1));
% depending on the noise, the fit can be quite bad, so we are a litle
% relaxed here (we should be within 3% of the original values)
%!assert( (abs(fitProblem.getSolutionVariables().get(0) - 0.03) / 0.03) < 3e-2);
%!assert( (abs(fitProblem.getSolutionVariables().get(1) - 0.004) / 0.004) < 3e-2);
%%%   

