-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathIris_active.m
169 lines (136 loc) · 5.14 KB
/
Iris_active.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
function Iris_active()
%% Initialize
clc
close all
warning off
parpool('local',24)
if ~isfolder('temp')
mkdir 'temp'
end
diary on
fprintf('Active Iris Curve Fitting\nWritten by Babak N. Safa (Ethier Lab)\n')
fprintf('\nMatlab version \t %s \n\n',version)
fprintf('Starting iris_active.m ...\n\n')
fprintf('%s\n',datetime)
N = 25;
DataFit_Output.N = N;
fprintf('Iris_active.m run, with grid size of N = %d \n',N)
DataFit_Output.Start_date = datetime('now');
%% Read experimental data
%% Seg
raw_seg = readtable('Iris_seg_summary.csv');
time = raw_seg.Time_sec_;
e_r_p = raw_seg.e_pupil;
%linear interpolation
index = (time>17)&(time<20);
e_r_p_max_0 = mean(e_r_p(index));
time_resampled = [0;.1; 1];% the .1 step is added to make sure the reference state of elements is logged in the log file
T_s_lc = [0;0;1];
T_d_lc = [0;0;0];
load = [time_resampled,T_s_lc,T_d_lc];
step_size = .1;
figure
hold on
plot(time,e_r_p,'o')
plot(time(index),e_r_p(index),'-*r')
%% Define cost function
%parameters E (MPa), v, T_sphincter(MPa), a_sphincter(mm)
lb = [0, 0, 0, .1 ];
ub = [1, .49, 1, 2];
if sum(ub<lb)>0
error('The bounds do not match')
end
%set the sphincter width (see
%/home/asixbabak/Dropbox (GaTech)/project-IrisBiomechanics/Iris sphincter
%dimension estimation from literature)
% test cases are a_s = [.4, .7, 1, 1.3]
fun_par_normal = @(x) [x(1), x(2), x(3), (1.3-lb(4))/(ub(4)-lb(4))]; %set the function to assign sphincter width
%in the variable names I have kept the usual FEBio units when using mm as dimension, which are based
%on mm, MPa and sec
DataFit_Output.par_name = {'E (MPa)','\nu','T_sphincter(MPa)', 'a_sphincter (mm)'};
DataFit_Output.lb = lb;
DataFit_Output.ub = ub;
%% Define the objective function
%% parameters for testing
E_test = .01;% in MPa, which is equivalent of 10 kPa
v_test = 0.4;
T_s_test = .01;% in MPa, which is equivalent of 10 kPa
a_s_test = 1;% in mm
r_p = 3.4; %pupil radius
x_test = [E_test, v_test, T_s_test, a_s_test];
temp = (x_test-lb)./(ub-lb);
x_test_normal = fun_par_normal(temp([1,2,3]));
cleanup = true;
[e_r_test1, ~, detail] = FEBio_run_Iris_Active(x_test_normal,lb,ub,r_p,load,step_size,cleanup);
rho = detail.rho;
%% Definitions go in here
cleanup = true;
%Fit only maximum constriction pupil strain
fun = @(x) sqrt((min(FEBio_run_Iris_Active(fun_par_normal(x),lb,ub,r_p,...
load,step_size,cleanup))-e_r_p_max_0).^2);
%% run tests
test1 = sqrt(sum((min(e_r_test1)-e_r_p_max_0).^2)/length(e_r_p_max_0));
test2 = fun(x_test_normal([1,2,3]));
if isempty(test1)||isempty(test2)||abs(test1-test2)>1e-12
error('FEBio test failed!')
end
%% Curvefitting
s=rng('shuffle');%for reproducing the random number
DataFit_Output.rand_state = s;
M_par_normal_master = lhsdesign(N,3); %the list of random guesses the values are picked from
M_par_normal = M_par_normal_master(1:N,:);
lb_normal = zeros(1,size(M_par_normal,2));
ub_normal = ones(1,size(M_par_normal,2));
options = optimoptions(@fmincon,'Algorithm','interior-point',...
'Display','Iter','DiffMinChange',.01,...
'ObjectiveLimit',1e-5,'StepTolerance',1e-20);
tmp_swap=struct([]);
fprintf('Optimization started ...\n')
parfor i=1:N
sprintf('Starting cycle number %d/%d',[i,N])
tic
x0_temp = M_par_normal(i,:);
[x_fit,fval,exitflag,~,~,grad,hessian] = fmincon(fun,x0_temp,[],[],[],[],...
lb_normal,ub_normal,[],options);
t_elapsed = toc;
fprintf('Optimization for cycle %d ended in %.3f hours\n',[i, t_elapsed/3600])
%% Evaluate the fit response
cleanup = false;
[e_r_p_max_fit, ~, detail] = FEBio_run_Iris_Active(fun_par_normal(x_fit),lb,ub,r_p,...
load,step_size,cleanup);
RMSE_e_r = fval;
elapsed_time = toc
pause(2)
if isfile('Iris_active.log')
movefile('Iris_active.log',sprintf('./temp/temp_%d_iris_active.log',i),'f');
end
if isfile('Iris_active.feb')
movefile('Iris_active.feb',sprintf('./temp/temp_%d_iris_active.feb',i),'f');
end
if isfile('Iris_active.xplt')
movefile('Iris_active.xplt',sprintf('./temp/temp_%d_iris_active.xplt',i),'f');
end
tmp_swap(i).elapsed_time = elapsed_time;
tmp_swap(i).x_par0 = lb + fun_par_normal(x0_temp).*(ub-lb);
tmp_swap(i).x_fit = lb + fun_par_normal(x_fit).*(ub-lb);
tmp_swap(i).fval = fval;
tmp_swap(i).exitflag = exitflag;
tmp_swap(i).grad = grad;
tmp_swap(i).hessian = hessian;
tmp_swap(i).time = time_resampled;
tmp_swap(i).rho = rho;
tmp_swap(i).e_r_p_max_fit = e_r_p_max_fit;
tmp_swap(i).e_r_p_max_0 = e_r_p_max_0;
tmp_swap(i).RMSE_e_r = RMSE_e_r;
tmp_swap(i).detail = detail;
pause(1)
end
%% Save Outputs
DataFit_Output.active_iris = tmp_swap;
fprintf('%s \n',datetime)
fprintf('Analyses concluded! \n')
formatOut = 'mm_dd_yyyy_HH_MM';
date_marker = datestr(now,formatOut);
save([sprintf('temp/DataFit_Output_N%d_',N),date_marker,'.mat'],'DataFit_Output')
diary off
end