From f0675ed53718cff0cdfe5a2b34640ca88bcb701a Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Sun, 25 Dec 2022 09:00:58 -0600 Subject: [PATCH] Delete calibrator.m --- Legacy Code/calibrator.m | 348 --------------------------------------- 1 file changed, 348 deletions(-) delete mode 100644 Legacy Code/calibrator.m diff --git a/Legacy Code/calibrator.m b/Legacy Code/calibrator.m deleted file mode 100644 index 2149631..0000000 --- a/Legacy Code/calibrator.m +++ /dev/null @@ -1,348 +0,0 @@ -% -% /-------------------------------------------------------\ -% | Date: 7/16/2018 | -% | Author: Jonathan Cody | -% | Institution: Purdue University | -% | Weldon School of Biomedical Engineering | -% | Pienaar Lab | -% \-------------------------------------------------------/ -% -%% To Do List ============================================================= - -% -> switch to fullname variables -% -> update formatting -% -> rewrite introduction -% -> write calling script template - -% -> add parallel processing -% -> add nan error catch to costfun? -% -> collect model solutions? - -%% FUNCTION =============================================================== - -function calibrator(x,Y,mod,R,suf,runs,... - p,l,cf,W,tout,mix,r,Np,Rp,Ns,xu,Yu,Yr) - -% -> add function outline here - -%% Defining some global variables ----------------------------------------- - -Ny = length(Y) ;% number of outputs -xx = unique(cat(2,x{:})) ;% time points for evaluation -lR = log10(R) ;% log-scale parameter range - -%% Setting default optional inputs ---------------------------------------- - -if isempty(p) - p = find(R(1,:)~=R(2,:) & ... - (R(2,:)<10*R(1,:) | R(1,:)<=0)) ;% sample indices -end -if isempty(l) - l = find(R(2,:)>=10*R(1,:) & R(1,:)>0) ;% log-sample indices -end -if isempty(cf) - cf = repmat({@(y,Y) sqrt(mean2((y - Y).^2))},size(Y)) ;% cost fun -end -if isempty(W) ; W = ones(size(Y))' ; end % cost weights of outputs -if isempty(Np) ; Np = 1 ; end % number of plots -if isempty(Rp) ; Rp = [1 0 1] ; end % plot parameter range -if isempty(xu) ; xu = ('Input') ; end % x units for plots -if isempty(Yu) - for m = 1:Ny ; Yu{m} = ['Output ' num2str(m)] ; end % Y units -end -if isempty(Yr) ; Yr = repmat({[-inf inf]},size(Y)) ; end % Y plot range - -%% Defining inputs for 'nlmefit' (if requested) --------------------------- - -if strcmp(mix,'yes') - XX = [] ;% predictor matrix - YY = [] ;% response vector - GG = [] ;% group vector - for n = 1:Ny - Xout = zeros(numel(Y{n}),Ny-1) ; - if n > 1 - Xout(:,n-1) = ones(numel(Y{n}),1) ; - end - XX = [XX ; [repmat(x{n}',size(Y{n},1),1) Xout] ] ; - for m = 1:size(Y{n},1) - YY = [YY ; Y{n}(m,:)'] ; - GG = [GG ; m*ones(size(x{n}))'] ; - end - end -end - -%% Creating results data structure or loading previous structure ---------- - -if exist(fullfile(cd,['C' suf '.mat']),'file') == 2 - load(['C' suf],'C') ;% load output structure - disp('Resuming runs.') ; disp(' ') -else - LHS = lhsdesign(runs,length([p l])) ;% latin hypercube sample - for n = 1:runs - C(n).L = LHS(n,:) ;% create output structure - end - C(1).in = struct('x',{x},'Y',{Y},'mod',{mod},'R',R,...% required inputs - 'p',p,'l',l,'cf',{cf},'W',W,'tout',tout,... % optional inputs - 'n',0) ;% last run completed - disp('Starting runs.') ; disp(' ') -end - -%% Calibrating from each initial parameter guess -------------------------- - -for n = C(1).in.n+1:runs - if n > runs ; break ; end % skipping to printing if done - tic % starting algorithm timer - - try - if strcmp(mix,'yes') - %% Calling 'nlmefit' algorithm (if requested) ----------------- - - % Generating fixed-effect guess with fmincon - [C(n).O] = fmincon(@(L)costfun(L,0),C(n).L,[],[],[],[],... - zeros(1,size(R,2)),ones(1,size(R,2)),[],... - optimoptions(@fmincon,'Display','none')) ; - - % Building mixed-effect model with nlmefit - C(n).P = getpars(C(n).O) ;% guess for fixed effects - [C(n).P(r),C(n).PSI,stats,C(n).B] = nlmefit(XX,YY,GG,[],... - @(phi,X)modelfun(C(n).P,phi,X),C(n).P(r)') ; - - %% Saving 'nlmefit' results (if requested) -------------------- - - C(n).B = C(n).B' ;% random effects - C(n).M = repmat(C(n).P,GG(end),1) ; - C(n).M(:,r) = C(n).M(:,r) + C(n).B ;% mixed effects - C(n).gof = stats.logl ;% log-likelihood - C(n).stats = stats ;% statistical outputs - - else - %% Calling 'fmincon' algorithm -------------------------------- - - [C(n).O,C(n).gof] = ... - fmincon(@(L)costfun(L,0),C(n).L,[],[],[],[],... - zeros(1,size(R,2)),ones(1,size(R,2)),[],... - optimoptions(@fmincon,'Display','none')) ; - - %% Saving 'fmincon' results ----------------------------------- - - C(n).P = getpars(C(n).O) ;% optimized parameters - C(n).cost = costfun(C(n).O,1) ;% cost for each output - - end - disp(['Run ' num2str(n) ' succeeded at ' datestr(now,'HH:mm') '.']) - - catch ME - %% Displaying error (if one occurs) ------------------------------ - - C(n).gof = nan ;% setting goodness-of-fit to NaN - disp(['Error: ' ME.message]) - for m = 1:length(ME.stack) - disp(['Function "' ME.stack(m).name ... - '" line ' num2str(ME.stack(m).line)]) - end - disp(['Run ' num2str(n) ' failed at ' datestr(now,'HH:mm') '.']) - end - - %% Saving results ----------------------------------------------------- - C(1).in.n = n ;% run number - save(['C' suf '.mat'],'C') % saving results -end - -%% Retrieving parameters from normalized values --------------------------- - -function P = getpars(L) - P = R(1,:) ;% parameter vector - if ~isempty(p) - P(p) = R(1,p) + L(1:length(p)).*(R(2,p)-R(1,p)) ;% adding sample - end - if ~isempty(l) - P(l) = 10.^(lR(1,l) + L(length(p)+1:end).*(lR(2,l)-lR(1,l))) ; - end -end - -%% Defining target function for nlmefit ----------------------------------- - -function yfit = modelfun(P,phi,X) - if ~isempty(tout) && (toc >= 60*tout) - error('Optimization algorithm timed out.') - end - P(r) = phi ;% adding mixed effects parameters - Yout = mod(xx,P) ;% solving (see model) - - for m = 1:Ny - yofX(:,m) = interp1(xx,Yout(:,m),X(:,1)) ; - end - for m = 1:size(X,1) - yfit(m,:) = yofX(m,1) ; - for o = 2:Ny - if X(m,o) == 1 - yfit(m,:) = yofX(m,o) ;% YY output - end - end - end -end - -%% Defining target function for fmincon ----------------------------------- - -function cost = costfun(L,all) - if ~isempty(tout) && (toc >= 60*tout) - error('Optimization algorithm timed out.') - end - P = getpars(L) ;% parameter vector - Yout = mod(xx,P) ;% solving (see model) - - for m = 1:Ny % per output: - y = interp1(xx,Yout(:,m),x{m}) ; - y = repmat(y,size(Y{m}, 1),1) ;% model value - cost(m) = cf{m}(y,Y{m}) ;% output cost - end - if all == 1 ; return ; end % returning output costs - cost = sum(cost.*W)/sum(W) ;% weighted mean cost -end - -%% Sorting results by goodness-of-fit ------------------------------------- - -if strcmp(mix,'yes') - [~,o] = sort(-cat(2,C(:).gof)) ;% order by maximum loglikelihood -else - [~,o] = sort(cat(2,C(:).gof)) ;% order by minimum cost -end -[C(:).in] = deal(C(1).in) ;% copying input structure to all elements -C = C(o) ;% ordering output structure elements - -fits = find(~isnan(cat(2,C(:).gof))) ; -if isempty(fits) - disp('No fits were found.') ; disp(' ') ; return -end -fits = fits(end) ;% number of fits found - -save(['C' suf '.mat'],'C') % saving results - -%% Limiting plotted results per Rp (if requested) ------------------------- - -O = cat(1,C(:).O) ; -b = 1:size(O,1) ; -for n = 1:size(Rp,1) - a = find([p l]==Rp(n,1)) ; - b = intersect(b,find(Rp(n,2)<=O(:,a) & O(:,a)<=Rp(n,3))) ; -end - -if isempty(b) - disp('No results meet parameter restrictions.') ; disp(' ') ; return -end - -%% Defining plot variables ------------------------------------------------ - -disp('Plotting figures.') - -if isempty(Ns) ; Ns = fits ; end % default fitted parameter sets -if fits < Ns ; Ns = fits ; end % summarizing all fitted sets -if fits < Np ; Np = fits ; end % plotting all cost minimums - -switch Ny + 1 + (Rp(1,3)-Rp(1,2)~=1) - case 2 ; plotgrid = [2 1] ;% size of subplot grid - case {3,4} ; plotgrid = [2 2] ;% size of subplot grid - case {5,6} ; plotgrid = [3 2] ;% size of subplot grid - case {7,8,9} ; plotgrid = [3 3] ;% size of subplot grid - case {10,11,12} ; plotgrid = [4 3] ;% size of subplot grid - otherwise - Ny = 10 ; plotgrid = [4 3] ;% size of subplot grid - disp('Too many ouputs. Plotting first 10.') -end - -mark = { ':o' ':^' ':s' ':d' ':p' ':h' ':+' ':x' ':*' ... - '-.o' '-.^' '-.s' '-.d' '-.p' '-.h' '-.+' '-.x' '-.*' ... - '--o' '--^' '--s' '--d' '--p' '--h' '--+' '--x' '--*' }; -for m = 1:Ny - if size(Y{m},1) > length(mark) - Y{m} = mean(Y{m}) ; - disp(['Output ' num2str(m) ... - ' has too many series to plot. Plotting mean.']) - end -end - -color = 'brgmcy' ;% colors for output plots - -xx = xx(1):getprec(xx(end)-xx(1)):xx(end) ;% denser vector for plotting - -function p = getprec(x) % returns precision of number - f = 14-9*isa(x,'single') ;% Don't know how... - s = sprintf('%.*e',f,x) ; - v = [f+2:-1:3,1] ; - s(v) = '0'+diff([0,cumsum(s(v)~='0')>0]); - p = str2double(s) ; -end - -%% Generating Plots ------------------------------------------------------- - -for n = 1:Np - n = b(n) ; - name = ['Run' num2str(n) suf] ;% figure filename - F = figure('Name',name,'Color',[1 1 1],... - 'defaultAxesColorOrder',[0 0 0;0 0 0]); - - %% Generating model solutions ----------------------------------------- - - if strcmp(mix,'yes') - for o = 1:GG(end) - Ymix{o} = mod(xx,C(n).M(o,:)) ;% model solutions - end - end - Yout = mod(xx,C(n).P) ;% model solution - - %% Plotting best fits against data ------------------------------------ - - for m = 1:Ny - subplot(plotgrid(1),plotgrid(2),m) - hold on - for o = 1:size(Y{m},1) - plot(x{m},Y{m}(o,:),mark{o},'Color',color(o),'MarkerSize',5) - end - if strcmp(mix,'yes') - for o = 1:GG(end) - plot(xx,Ymix{o}(:,m),'Color',color(o),'LineWidth',2) - end - end - plot(xx,Yout(:,m),'Color','k','LineWidth',2) - axis([xx(1) xx(end) Yr{m}(1) Yr{m}(2)]) - ylabel(Yu{m}) - xlabel(xu) - set(gca,'FontSize',11,'FontName','Arial','LineWidth',1.5) - grid on - end - - %% Summarizing parameter GoF and parameters --------------------------- - - for o = 1:2 - subplot(plotgrid(1),plotgrid(2),Ny+o) - hold on - for m = Ns:-1:1 - if ~any(b==m) && o==2 ; continue ; end - bluescale = [0 0 1]+[1 1 0]*(m-1)/(Ns-1) ; - grayscale = [1 1 1]*(m-1)/(Ns-1) ; - yyaxis left - plot(0,C(m).gof,'o','MarkerEdgeColor',bluescale,... - 'MarkerFaceColor',bluescale) - yyaxis right - plot([p l],C(m).O,'*','MarkerEdgeColor',grayscale) - end - yyaxis left - plot(0,C(n).gof,'o','MarkerEdgeColor','r','MarkerSize',10) - ylabel('Model GoF (blue)') - yyaxis right - plot([p l],C(n).O,'o','MarkerEdgeColor','r','MarkerSize',10) - ylabel('Normalized Parameters') - axis([0 size(R,2)+1 0 1]) - xticks(1:2:size(R,2)) - if o == 2 ; xlabel('(Restricted Range)') ; end - set(gca,'FontSize',11,'FontName','Arial','LineWidth',1.5) - if (Rp(1,3)-Rp(1,2)==1) ; break ; end - end - saveas(F,name,'fig') -end - -disp(' ') ; disp(['Done! Found ' num2str(fits) ' fits.']) -%#ok<*AGROW> -%#ok<*FXUP> -%#ok<*FXSET> -end