calc_jacobian_test

 Verify Jacobian Calculation by small derivative from forward problem
 Also calculate dataprior
     Difference dataprior should be 1
     normalized difference dataprior should be 1./ homg_data

0001 function ok= calc_jacobian_test
0002 % Verify Jacobian Calculation by small derivative from forward problem
0003 % Also calculate dataprior
0004 %     Difference dataprior should be 1
0005 %     normalized difference dataprior should be 1./ homg_data
0006 
0007 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0008 % $Id: calc_jacobian_test.m 3733 2013-01-15 20:59:22Z aadler $
0009 
0010 ok= 1;
0011 delta = 1e-5;
0012 testvec= [5,20,40,130];
0013 
0014 mdl= make_aa_mdl2;
0015 %
0016 disp('test jacobian_adjoint (2D) for difference data')
0017 ok=ok & run_jacobian_test( mdl, delta, testvec );
0018 ok=ok & run_dataprior_test( mdl );
0019 
0020 %
0021 disp('test jacobian_adjoint (2D) for normalized difference data');
0022 mdl = mdl_normalize(mdl,1);
0023 ok=ok & run_jacobian_test( mdl, delta, testvec );
0024 ok=ok & run_dataprior_test( mdl );
0025 
0026 mdl= make_aa_mdl3;
0027 %
0028 disp('test jacobian_adjoint (3D) for difference data')
0029 ok=ok & run_jacobian_test( mdl, delta, testvec );
0030 ok=ok & run_dataprior_test( mdl );
0031 
0032 mdl = mdl_normalize(mdl,1);
0033 disp('test jacobian_adjoint (3D) for normalized difference data')
0034 ok=ok & run_jacobian_test( mdl, delta, testvec );
0035 ok=ok & run_dataprior_test( mdl );
0036 
0037 mdl= make_np_mdl;
0038 %
0039 disp('test np_calc_jacobian for difference data')
0040 ok=ok & run_jacobian_test( mdl, delta, testvec );
0041 ok=ok & run_dataprior_test( mdl );
0042 
0043 mdl = mdl_normalize(mdl,1);
0044 disp('test np_calc_jacobian for normalized difference data')
0045 ok=ok & run_jacobian_test( mdl, delta, testvec );
0046 ok=ok & run_dataprior_test( mdl );
0047 
0048 
0049 % run the jacobian test
0050 function ok= run_jacobian_test( mdl, delta, testvec ); 
0051     calc_norm = 0;
0052     if mdl_normalize(mdl)
0053         calc_norm = 1;    
0054     end
0055 
0056     img= eidors_obj('image', 'homg image');
0057     img.fwd_model= mdl;
0058 
0059     img.elem_data= ones( size(mdl.elems,1) ,1);
0060     homg_data=fwd_solve( img);
0061 
0062     J= calc_jacobian( img );
0063 
0064     % J = dF/dx = [F(x+d)  - F(x)]/d
0065     sumdiff= 0;
0066     bkgnd_elem_data= img.elem_data;
0067     for testelem = testvec
0068        img.elem_data= bkgnd_elem_data;
0069        img.elem_data(testelem)= bkgnd_elem_data(testelem)+delta;
0070        inh_data=fwd_solve( img);
0071 
0072        if calc_norm
0073           simJ= 1/delta* (inh_data.meas ./ homg_data.meas - 1);
0074        else
0075           simJ= 1/delta* (inh_data.meas-homg_data.meas);
0076        end
0077        
0078        plot([J(:,testelem) simJ]);
0079        sumdiff = sumdiff + std( J(:,testelem) - simJ );
0080     end
0081 
0082     tol= 1e-4*std(J(:));
0083     dif= sumdiff/length(testvec);
0084     if sumdiff/length(testvec) > tol
0085        eidors_msg('Jacobian calculation error (%s) tol(%3.2g)>diff(%3.2g)', ...
0086                 mdl.name,tol,dif, 1);
0087        ok=0;
0088     else
0089        ok=1;
0090     end
0091 
0092 %
0093 % test dataprior
0094 %     Difference dataprior should be 1
0095 %     normalized difference dataprior should be homg_data.^2
0096 function ok= run_dataprior_test( mdl )
0097     img= eidors_obj('image', 'homg image');
0098     img.fwd_model= mdl;
0099 
0100     img.elem_data= ones( size(mdl.elems,1) ,1);
0101     homg_data=fwd_solve( img);
0102 
0103     DP= calc_meas_icov( img );
0104 
0105     % difference dataprior
0106     testvec= diag(DP);
0107     if mdl_normalize(mdl)
0108 %  from calc_meas_icov, we have the following
0109 %     meas_icov = sparse(1:n, 1:n, ( homg_data.meas ).^2 );
0110         testvec = homg_data.meas.^2 ./ diag(DP);
0111     end
0112 
0113     mdiff = full(max(abs(diff( testvec ))));
0114     if mdiff > 1e-12 
0115        ok=0;
0116        eidors_msg('Dataprior calculation error (%s) = %f', mdl.name, mdiff, 1);
0117 keyboard
0118     else
0119        ok=1;
0120     end
0121 
0122 
0123 
0124 % 2D model with point electrodes
0125 %
0126 function mdl= make_aa_mdl2;
0127     n_elec= 16;
0128     n_rings= 1;
0129     options = {'no_meas_current','no_rotate_meas'};
0130     params= mk_circ_tank(8, [], n_elec);
0131 
0132     params.stimulation= mk_stim_patterns(n_elec, n_rings, '{ad}','{ad}', ...
0133                                 options, 10);
0134     params.solve=      'fwd_solve_1st_order';
0135     params.system_mat= 'system_mat_1st_order';
0136     params.jacobian=   'jacobian_adjoint';
0137     params.normalize_measurements = 0;
0138     mdl = eidors_obj('fwd_model', params);
0139     mdl.name= 'AA_1996 mdl';
0140 
0141 function mdl= make_aa_mdl3;
0142     i_mdl = mk_common_model('b3cz2',16);
0143     mdl= i_mdl.fwd_model;
0144     mdl.name= 'AA_1996 mdl';
0145     mdl.solve=      'fwd_solve_1st_order';
0146     mdl.system_mat= 'system_mat_1st_order';
0147     mdl.jacobian=   'jacobian_adjoint';
0148     
0149     
0150 
0151 function mdl= make_np_mdl;
0152     i_mdl = mk_common_model('n3r2',16);
0153     mdl= i_mdl.fwd_model;
0154     mdl.name=     'NP_2003 mdl';
0155     mdl.solve=      'np_fwd_solve';
0156     mdl.system_mat= 'np_calc_system_mat';
0157     mdl.jacobian=   'np_calc_jacobian';
0158 
0159 function mdl= make_ms_mdl;
0160     i_mdl = mk_common_model('n3r2',16);
0161     mdl= i_mdl.fwd_model;
0162     mdl.name=       'MS_2005 mdl';
0163     mdl.solve=      'np_fwd_solve';
0164     mdl.system_mat= 'np_calc_system_mat';
0165     mdl.jacobian=   'ms_calc_jacobian';