Eidors在编程中使用了cache,主要是用来重用缓存的值对于从而提高数值效率。但是这样就对新手阅读代码造成了一定的困扰。
关于如何将数据写进缓存,我们大可不必关心。
以fwd_model_parameters.m的一段程序为例。
还有calc_jacobian_bkgnd.m
给定输入参数(逆问题求解结构体,正问题结构体,电压数据),varargout = eidors_obj('get-cache', cache_obj, fstr );
会先查看缓存里是否计算过重构电导率,倘若计算过,就不会再计算。直接将varargout输出。
倘若没计算过,varargout应该为空。
倘若计算过,varargout为非空,可以看到elem_data(这就是待求的电导率)是有值的,所以不用再进行逆问题求解,从而提高成像速度。
一、逆问题求解
这里以一个官网的最简单的单电阻模型例程为例。
考虑最简单的EIT系统。我们有一个电阻,我们想知道它的值。因此,我们将一个电极连接到每个端子上,并施加几个不同的测试电流(I1, I2, I3)和测量电压:(V1, V2, V3)从这些测量值计算出电阻的最小二乘估计。
r_mdl= eidors_obj('fwd_model','demo resistor model');
r_mdl = mdl_normalize( r_mdl, 0);
% Geometry
r_mdl.nodes= [1,1,1; 2,2,2];
r_mdl.elems= [1,2];
r_mdl.boundary= [1,2];
% Define Electrodes (there is only one)
r_mdl.electrode(1).z_contact= 10; % ohms
r_mdl.electrode(1).nodes= 1;
r_mdl.gnd_node= 2;
show_fem(r_mdl); view(-12,24);
% Define stimulation patterns
for i=1:3
r_mdl.stimulation(i).stimulation= 'Amps';
r_mdl.stimulation(i).stim_pattern= ( 0.001*i );
r_mdl.stimulation(i).meas_pattern= 1; % measure electrode 1
end
r_mdl.solve= @tutorial020_f_solve;
% Define an 'image'
resistor = eidors_obj('image', 'resistor');
resistor.elem_data= 1000;
resistor.fwd_model= r_mdl;
% Calculate data for 1k resistor
data_1k0 =fwd_solve( resistor );
% Now change resistor to be 1.2k
resistor.elem_data= 1800;
data_1k2 =fwd_solve( resistor );
% Solve resistor model
% $Id: tutorial020c.m 3127 2012-06-08 16:19:25Z bgrychtol $
% Now we complete the fwd_model
r_mdl.jacobian= @jacobian_perturb;
% Now create an inverse model
i_mdl= eidors_obj('inv_model','resistor inverse');
i_mdl.fwd_model= r_mdl;
i_mdl.jacobian_bkgnd.value= 1000;
% regulatization not needed for this problem
i_mdl.RtR_prior= @prior_tikhonov;
i_mdl.hyperparameter.value= 0;
i_mdl.reconst_type= 'difference';
i_mdl.solve= @inv_solve_diff_GN_one_step;
% Reconstruct resistor change
reconst= inv_solve(i_mdl, data_1k0, data_1k2)
这里对如何实现逆问题重构,将eidors主要计算过程展开。
从这里开始断点调试。
1、inv_solve
2、 imgc= eidors_cache( inv_model.solve, {inv_model, fdata1, fdata2},'inv_solve');
输入:
inv_model.solve:@inv_solve_diff_GN_one_step
{inv_model, fdata1, fdata2}:
‘inv_solve’:表示类型为逆问题求解
输出:imgc(该图像被进一步 处理为eidors_obj,作为最终的输出img)
1)eval(sprintf(‘%s = %s’, output, ‘cache_shorthand(command, varargin{:});’));
- sprintf(‘%s = %s’, output, ‘cache_shorthand(command, varargin{:});’)
ans =
'[varargout{1} ] = cache_shorthand(command, varargin{:});'
输入:
command:调用的函数或者说求解的方法
varargin{:}:包含了逆模型结构体
、用于差分成像的电压数据
(vh,vi)、求解类型
。
输出:
varargout{1}:计算结果
- cache_shorthand(command, varargin{:})进一步探究之eval(sprintf(‘%s = %s’, output, ‘feval(fhandle,args{:});’));
输入:
sprintf(‘%s = %s’, output, ‘feval(fhandle,args{:});’)
ans =
'[varargout{1} ] = feval(fhandle,args{:});'
%fhandle:@inv_solve_diff_GN_one_step
%args{:}:
args{:}:
-
eval(sprintf(‘%s = %s’, output, ‘feval(fhandle,args{:});’));
- 调用inv_solve_diff_GN_one_step.m
- sol = eidors_cache(@get_RM, inv_model,‘inv_solve_diff_GN_one_step’ ) * dv;
eidors_cache(@get_RM, inv_model,‘inv_solve_diff_GN_one_step’ ) 计算出重构矩阵RM
所以sol = RM*dv-
eval(sprintf(‘%s = %s’, output, ‘cache_shorthand(command, varargin{:});’));
- eval(sprintf(‘%s = %s’, output, ‘feval(fhandle,args{:});’));
两个输入
执行完的结果
- eval(sprintf(‘%s = %s’, output, ‘feval(fhandle,args{:});’));
-
- 调用inv_solve_diff_GN_one_step.m
输出:
J= jacobian_adjoint( fwd_model, img)
注:eval的用法