resistor_model

 DEMO to show really simple application of EIDORS framework

 This code models a resistor (resitance R), with one electrode at each end. 

 Stimulation patterns apply current I1, I2, I3, and measure voltages
 Forward Model: V= IR    => [V1;V2;V3] = [I1;I2*I3]*(R+2*Zc)
 Jacobian:      J= dV/dR =I = [I1; I2; I3]
 Inverse Model: R= inv(J'*J)*J'*V
    This corresponds to the least squares solution

0001 % DEMO to show really simple application of EIDORS framework
0002 %
0003 % This code models a resistor (resitance R), with one electrode at each end.
0004 %
0005 % Stimulation patterns apply current I1, I2, I3, and measure voltages
0006 % Forward Model: V= IR    => [V1;V2;V3] = [I1;I2*I3]*(R+2*Zc)
0007 % Jacobian:      J= dV/dR =I = [I1; I2; I3]
0008 % Inverse Model: R= inv(J'*J)*J'*V
0009 %    This corresponds to the least squares solution
0010 
0011 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0012 % $Id: resistor_model.m 3630 2012-11-18 18:29:52Z aadler $
0013 
0014 function resistor_model;
0015 
0016 % Set log level to show all messages
0017 eidors_msg('set_level',4);
0018 
0019 %
0020 % Step 1: create FEM model structure
0021 %
0022 % Fwd model:
0023 %  Two nodes are in space at [1,1,1] and [2,2,2]
0024 %  The resistor is connected between them
0025 
0026 r_mdl.name = 'demo resistor model';
0027 r_mdl.nodes= [1,1,1;  2,2,2];
0028 r_mdl.elems= [1;2];
0029 r_mdl.solve=      @f_solve;
0030 r_mdl.jacobian=   @c_jacobian;
0031 
0032 %
0033 % Step 2: create FEM model electrode definitions
0034 %
0035 
0036 r_mdl.electrode(1).z_contact= 10; % ohms
0037 r_mdl.electrode(1).nodes=     1;
0038 r_mdl.gnd_node= 2;
0039 r_mdl.normalize_measurements = 0;
0040 
0041 %
0042 % Step 3: create stimulation and measurement patterns
0043 % patterns are 0.010,0.020,0.030 Amp
0044 
0045 for i=1:3
0046     r_mdl.stimulation(i).stimulation= 'Amp';
0047     r_mdl.stimulation(i).stim_pattern= ( 0.010*i );
0048     r_mdl.stimulation(i).meas_pattern= 1; % measure electrode 1
0049 end
0050 
0051 r_mdl= eidors_obj('fwd_model', r_mdl);
0052 
0053 %
0054 % Step 4: simulate data for medium with R=1 kohms
0055 % This medium is called an 'image'
0056 %
0057 
0058 img_1k = eidors_obj('image', 'homogeneous image', ...
0059                      'elem_data', 1e3, ...
0060                      'fwd_model', r_mdl );
0061 
0062 data_1k =fwd_solve( r_mdl, img_1k );
0063 
0064 %
0065 % Step 5: add noise to simulated data
0066 %
0067 
0068 data_noise= eidors_obj('data', 'noisy data', ...
0069                        'meas', data_1k.meas + 1e-3*randn(3,1));
0070 
0071 %
0072 % Step 7: create inverse model
0073 %
0074 
0075 % create an inv_model structure of name 'demo_inv'
0076 r_inv.name=  'Resistor Model inverse';
0077 r_inv.solve= @i_solve;
0078 r_inv.reconst_type= 'static';
0079 r_inv.fwd_model= r_mdl;
0080 r_inv= eidors_obj('inv_model', r_inv);
0081 
0082 %
0083 % solve inverse model');
0084 %
0085 
0086 R= inv_solve( r_inv, data_1k );
0087 fprintf('R calculated with clean data= %5.3f kOhm\n', R.elem_data / 1000 );
0088 
0089 R= inv_solve( r_inv, data_noise );
0090 fprintf('R calculated with noisy data= %5.3f kOhm\n', R.elem_data / 1000 );
0091 
0092 
0093 % Forward Model:
0094 % For each stimulation there is I1 into Node1
0095 %  Node2 is connected to gnd with Zcontact
0096 %
0097 % each stim has one measurement pattern where
0098 %  Vmeas= Meas_pat * Node1
0099 %       = Meas_pat * I1 * ( Zcontact*2 + R )
0100 %
0101 % Thus
0102 %  V= IR    => [V1;V2;V3] = [I1;I2*I3]*(R + 2*Zcontact)
0103 function data =f_solve( f_mdl, img )
0104   R= img.elem_data;
0105 
0106   n_stim= length( f_mdl.stimulation );
0107   V= zeros(n_stim, 1);
0108 
0109   for i=1:n_stim
0110     if ~strcmp( f_mdl.stimulation(i).stimulation, 'Amp' )
0111        error('f_solve expects current in Amp');
0112     end
0113 
0114     I        = f_mdl.stimulation(i).stim_pattern;
0115     meas_pat = f_mdl.stimulation(i).meas_pattern;
0116 
0117     stim_elec= find( I );
0118     zc       = f_mdl.electrode( stim_elec ).z_contact;
0119 
0120     V(i)= meas_pat * I * ( R + 2*zc);
0121   end
0122 
0123   data.name='resistor model data';
0124   data.meas= V;
0125   
0126 % Jacobian:      J= dV/dR =I = [I1; I2; I3]
0127 function J= c_jacobian( f_mdl, img)
0128   n_stim= length( f_mdl.stimulation );
0129   J= zeros(n_stim, 1);
0130   for i=1:n_stim
0131     J(i)     = f_mdl.stimulation(i).stim_pattern; % Amp
0132   end
0133 
0134 % Inverse Model: R= inv(J'*J)*J'*V
0135 %    This corresponds to the least squares solution
0136 function img= i_solve( i_mdl, data )
0137   % Normally the Jacobian depends on an image. Create a dummy one here
0138   i_img= mk_image(i_mdl,1);
0139   J = calc_jacobian( i_img); 
0140 
0141   img.name= 'solved by i_solve';
0142   img.elem_data= (J'*J)\J'* data(:);
0143   img.inv_model= i_mdl;
0144   img.fwd_model= i_mdl.fwd_model;
0145 
0146