function fitMatlabCellML1_5(display)
clear all;
clc;
tic;
%%%Set up some common variables needing to be changed between models
time = 'time'; %Some models use 't', some use 'time'
global lengthToRun
lengthToRun = 0:1:5000; %Length of time to run the model for
suffix = '_COR_CuminEdit'; %Mark at end of saved file to show it's being used in this function
global objF
objF = MyHELDTobjectiveFunction; %may write different objective functions and place here
boundSearch = false;
constantMultipliers = 20;
global modelname Y constantNames constantValues computedNames computedValues stateNames
if nargin == 0
display =1;
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Read in file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[filename path] = uigetfile('*.*');
t=textread([path,filename],'%s','delimiter','\r');
if(display)
disp(['Using model: ',filename(1:strfind(filename,'.')-1),' (time elapsed = ',num2str(toc),'s)']);
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Find variables & constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Set up to find where things are
textToFind = {'Initial conditions';'State variables';...
'Constants';'Computed variables';'State variables'};
textPosition = zeros(2,length(textToFind));
%%%Find where things are
%%%Fill in blank lines so when saves, gap remains
for i=1:length(t)
if strcmp(t{i},'')
t{i} = '%';
end
end
%%%Find start and end of sections
for i=1:length(t)
for t2f = 1:length(textToFind);
if (strfind(t{i}, textToFind{t2f}))
textPosition(1,t2f) = i+3;
for j=i+3:length(t)
if( strfind(t{j},'%-') )
textPosition(2,t2f) = j-2;
break;
end
end
end
end
end
%%%INITIAL CONDITIONS
%%%Get the initial conditions to run the model with
Y = (t{textPosition(1,1)});
Y = eval(Y(6:length(Y)));
%%%STATE VARIABLES
%%%Store the constants for later use
if( (textPosition(1,5) > 0) && (textPosition(2,5)>0) )
stateCount=1;
stateNames={};
for i=textPosition(1,5):textPosition(2,5)
lineText = t{i};
%%Save constant and 'original' value for later if needed.
stateNames{stateCount} = strtrim(lineText(strfind(lineText,':')+1:strfind(lineText,'(')-1));
stateCount = stateCount+1;
end
end
%%%CONSTANTS
%%%Store the constants for later use
if( (textPosition(1,3) > 0) && (textPosition(2,3)>0) )
constantCount=1;
constantNames={};
constantValues=[]; %%Need to initialise the matrix
for i=textPosition(1,3):textPosition(2,3)
lineText = t{i};
%%Save constant and 'original' value for later if needed.
constantNames{constantCount} = strtrim(lineText(1:strfind(lineText,'=')-1));
tempValue = str2num(lineText(strfind(lineText,'=')+1:...
strfind(lineText,';')-1));
constantValues(constantCount) = tempValue;
constantCount = constantCount+1;
end
end
%%%COMPUTED VARIABLES
%%%Store computed variable names
if( (textPosition(1,4) > 0) && (textPosition(2,4)>0) )
computedCount=1;
computedNames={};
for i=textPosition(1,4):textPosition(2,4)
lineText = t{i};
varName = lineText(2:strfind(lineText,'(')-2);
computedNames{computedCount} = strtrim(varName);
computedCount = computedCount+1;
end
end
if(display)
disp(['All varaibales read in. (time elapsed = ',num2str(toc),'s)']);
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Find constant(s) to use and edit line(s) in text
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Find the constants we want to change to fit data
constantsToUse = DCsGUI(constantNames);
drawnow; %%gets rid on the figure properly (otherwise computation prevents screen refresh)
count = 1;
for i=textPosition(1,3)+constantsToUse-1
%%edit the line of the constant to be used in the fminsearch
lineText = t{i};
lineText = [lineText(1:strfind(lineText,';')-1),...
' * constantMult(',num2str(count),');',lineText(strfind(lineText,'%'):length(lineText))];
t{i} = lineText;
count = count+1;
end
if(display)
disp(['Selected constants altered in temp file. (time elapsed = ',num2str(toc),'s)']);
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Edit the file to use later
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%COMPUTED VARS - so that we can easily plot them (or use later) if we want
%%%Working from the bottom up so that we don't 'cross over'
%%%Add variables to 'computed'
if( (textPosition(1,4) > 0) && (textPosition(2,4)>0) )
for i=1:length(computedNames)
lineText = strcat('computed(:,',num2str(i),')=',computedNames{i},';');
t=[t;lineText];
end
%%% Initialise 'computed'
t{14} = 'computed = [];';
%%% add 'computed' to the function line
t{13} = ['function [dY,computed] = ',filename(1:length(filename)-2),suffix,'(',time,',Y,constantMult)'];
else
%%% Don't initialise 'computed'
t{14} = '%';
%%% Don't add 'computed' to the function line
t{13} = ['function [dY] = ',filename(1:length(filename)-2),suffix,'(',time,',Y,constantMult)'];
end
%%%Write out new file
writeascii([path,filename(1:strfind(filename,'.')-1),suffix,'.m'],t,'%r\n')
%%Read in new file and use this one
filename=[filename(1:length(filename)-2),suffix,'.m'];
t=textread([path,filename],'%s','delimiter','\r');
modelname = filename(1:length(filename)-2);
if(display)
disp(['File ready to be used (time elapsed = ',num2str(toc),'s)']);
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Do the parameter fitting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
constantMultipliers = constantMultipliers.*ones(1,length(constantsToUse));
if(boundSearch)
%%%Use this one if you MUST bound the results - it is not as good as fminsearch
[xFits fval exit output] = fminsearchbnd(@model2Min,constantMultipliers,...
constantMultipliers,constantMultipliers...
optimset('MaxIter',1000));
else
[xFits fval exit output] = fminsearch(@model2Min,constantMultipliers*20,...
optimset('MaxIter',1000));
end
disp(['This took ', num2str(output.iterations), ' iterations.']);
%%Write out the results nicely
outString='The new values are:\n';
for i=1:length(constantsToUse)
outString = strcat(outString,'\n',constantNames{constantsToUse(i)},' = ');
value = constantValues(constantsToUse(i));
for j=1:length(constantsToUse)
if i==j
newValue(i) = value * xFits(i);
end
end
outString = strcat(outString, num2str(newValue(i)), ' (',num2str(xFits(i)),')');
end
outString = strcat(outString, '\n Fval = ',num2str(fval));
NewValues=sprintf(outString)
if(display)
disp(['Paramaters fitted (time elapsed = ',num2str(toc),'s)']);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Save the new values if wanted
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
button = questdlg('Would you like to save the new parameters?','Save fitted parameters?','Yes','No','No');
if strcmp(button,'Yes')
tNew=textread([path,filename],'%s','delimiter','\r');
for i=1:length(t)
%%Fill in blank lines so when saves, gap remains
if strcmp(tNew{i},'')
tNew{i} = '%';
end
end
count=1;
for i=constantsToUse-1
%edit the line of the constant to be used in the fminsearch
lineText = tNew{i+textPosition(1,3)};
lineText = [lineText(1:strfind(lineText,'=')+1), num2str(newValue(count)), lineText(strfind(lineText,';'):length(lineText))];
tNew{i+textPosition(1,3)} = lineText;
count=count+1;
end
%%Save new file
[Savename Savepath] = uiputfile('*.*');
%%% add 'computed' to the function line
tNew{13} = ['function [dY,computed] = ',Savename(1:strfind(Savename,'.')-1),'(t,Y,constantMult)'];
%%Write out new file
writeascii([Savepath,Savename],tNew,'%r\n');
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Delete the file once used
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
delete([path filename]);
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end %%%End of function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Objective function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function obj2Min = model2Min(guess,lengthToRun)
%%%get constants and computeds from above
global modelname Y constantNames constantValues computedNames computedValues stateNames
global lengthToRun objF
%%%Compute the model and variables
time=lengthToRun;
[tData, YData] = runModel(modelname, time, guess);
%%%Save computed variables as their names
for i=1:length(computedNames)
eval([computedNames{i},'=[',num2str(computedValues(:,i)'),'];']);
end
%%%Save computed variables as their names
for i=1:length(stateNames)
eval([stateNames{i},'=[',num2str(YData(:,i)'),'];']);
end
%%%Save the constants as their names
for i=1:length(constantValues)
eval([constantNames{i},'=',num2str(constantValues(i)),';']);
end
obj2Min = evalc(objF); %%%Calculate the objective function. Will store the characters 'ans = <answer>' in obj2Min
obj2Min = str2num(obj2Min(find(obj2Min=='=')+1:length(obj2Min))); %%%Retrieve the answer as a number not part of char
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Run the model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tData, YData] = runModel(modelname, time, constantsToUse)
global Y constantNames constantValues computedNames computedValues stateNames
constantMultipliers = ones(1,length(constantsToUse));
fhandle = str2func(modelname);
[tData, YData] = ode45(fhandle, time, Y, [], constantsToUse);
if(length(computedNames)>0)
for i=1:length(tData)
eval(['[dY,computedValues(i,:)] = ',modelname,'(tData(i),YData(i,:), constantsToUse);']);
end
%%%Save computed variables as their names
for i=1:length(computedNames)
eval([computedNames{i},'=[',num2str(computedValues(:,i)'),'];']);
end
%%%Save computed variables as their names
for i=1:length(stateNames)
eval([stateNames{i},'=[',num2str(YData(:,i)'),'];']);
end
else
for i=1:length(tData)
eval(['[dY] = ',modelname,'(tData(i),YData(i,:), constantsToUse);']);
end
end
%%%Save the constants as their names
for i=1:length(constantValues)
eval([constantNames{i},'=',num2str(constantValues(i)),';']);
end
end