0001 function test_c2f_jacobian
0002
0003
0004
0005
0006 eidors_cache clear
0007
0008 test1_np_3d
0009 test2_aa_2d
0010 test3_aa_3d
0011 test4_npaa_3d
0012
0013 function test1_np_3d
0014
0015 imdl = mk_common_model('n3r2');
0016 f1mdl = imdl.fwd_model;
0017 f1mdl.solve = @np_fwd_solve;
0018 f1mdl.system_mat = @np_calc_system_mat;
0019 f1mdl.jacobian = @np_calc_jacobian;
0020
0021
0022 n2d = f1mdl.nodes( (f1mdl.nodes(:,3) == 0), 1:2);
0023 e2d = delaunayn(n2d);
0024 c_mdl = eidors_obj('fwd_model','2d','elems',e2d,'nodes',n2d);
0025
0026
0027 c2f= mk_coarse_fine_mapping( f1mdl, c_mdl );
0028 f2mdl= f1mdl;
0029 f2mdl.coarse2fine = c2f;
0030
0031
0032 imdl.solve= @np_inv_solve;
0033 imdl.hyperparameter.value = 1e-3;
0034 imdl.R_prior= @np_calc_image_prior;
0035 imdl.np_calc_image_prior.parameters= [3 1];
0036 imdl.jacobian_bkgnd.value= .6;
0037 imdl.reconst_type= 'difference';
0038 imdl.fwd_model= f1mdl;
0039
0040 img = calc_jacobian_bkgnd( imdl);
0041
0042 t=cputime;
0043 J1= np_calc_jacobian( f1mdl, img)*c2f;
0044 fprintf('np_calc - full: t=%f\n', cputime-t );
0045 n_J1 = norm(J1,'fro');
0046
0047 t=cputime;
0048 J2= np_calc_jacobian( f2mdl, img);
0049 fprintf('np_calc - c2f: t=%f\n', cputime-t );
0050
0051 n_J1_J2 = norm(J1-J2,'fro')/ n_J1;
0052 if n_J1_J2 > 1e-4;
0053 warning('J1-J2 is too big');
0054 end
0055
0056
0057 t=cputime;
0058 J2p= jacobian_perturb( f2mdl, img);
0059 fprintf('perturb_j - c2f: t=%f\n', cputime-t );
0060
0061 n_J1_J2p = norm(J1-J2p,'fro')/ n_J1;
0062 if n_J1_J2p > 1e-4;
0063 warning('J1-J2p is too big');
0064 end
0065
0066
0067
0068
0069 function test2_aa_2d
0070
0071 imdl= mk_common_model('a2c2',16);
0072 cmdl = imdl.fwd_model;
0073 imdl= mk_common_model('c2c0',16);
0074 f1mdl = imdl.fwd_model;
0075
0076 c2f= mk_coarse_fine_mapping(f1mdl,cmdl);
0077 img= calc_jacobian_bkgnd(imdl);
0078
0079 t=cputime;
0080 J1= jacobian_adjoint( f1mdl, img)*c2f;
0081 fprintf('aa_calc - full: t=%f\n', cputime-t );
0082 n_J1 = norm(J1,'fro');
0083
0084 f2mdl= f1mdl;
0085 f2mdl.coarse2fine = c2f;
0086
0087
0088 t=cputime;
0089 J2= jacobian_adjoint( f2mdl, img);
0090 fprintf('aa_calc - c2f: t=%f\n', cputime-t );
0091
0092 n_J1_J2 = norm(J1-J2,'fro')/ n_J1;
0093 if n_J1_J2 > 1e-4;
0094 warning('J1-J2 is too big');
0095 keyboard
0096 end
0097
0098
0099
0100 t=cputime;
0101 J2p= jacobian_perturb( f2mdl, img);
0102 fprintf('perturb_j - c2f: t=%f\n', cputime-t );
0103
0104 n_J1_J2p = norm(J1-J2p,'fro')/ n_J1;
0105
0106 if n_J1_J2p > 1e-4;
0107 warning('J1-J2p is too big');
0108 end
0109
0110
0111 t=cputime;
0112 J1p= jacobian_perturb( f1mdl, img)*c2f;
0113 fprintf('perturb_j - full: t=%f\n', cputime-t );
0114
0115 n_J1_J1p = norm(J1-J1p,'fro')/ n_J1;
0116 if n_J1_J1p > 1e-4;
0117 warning('J1-J1p is too big');
0118 end
0119
0120
0121 function test3_aa_3d
0122
0123 imdl= mk_common_model('a2c0',16);
0124 cmdl = imdl.fwd_model;
0125 imdl= mk_common_model('a3cr',16);
0126 f1mdl = imdl.fwd_model;
0127
0128 c2f= mk_coarse_fine_mapping(f1mdl,cmdl);
0129 img= calc_jacobian_bkgnd(imdl);
0130
0131 t=cputime;
0132 J1= jacobian_adjoint( f1mdl, img)*c2f;
0133 fprintf('aa_calc - full: t=%f\n', cputime-t );
0134 n_J1 = norm(J1,'fro');
0135
0136 f2mdl= f1mdl;
0137 f2mdl.coarse2fine = c2f;
0138
0139
0140 t=cputime;
0141 J2= jacobian_adjoint( f2mdl, img);
0142 fprintf('aa_calc - c2f: t=%f\n', cputime-t );
0143
0144 n_J1_J2 = norm(J1-J2,'fro')/ n_J1;
0145 if n_J1_J2 > 1e-4;
0146 warning('J1-J2 is too big');
0147 keyboard
0148 end
0149
0150
0151
0152 t=cputime;
0153 J2p= jacobian_perturb( f2mdl, img);
0154 fprintf('perturb_j - c2f: t=%f\n', cputime-t );
0155
0156 n_J1_J2p = norm(J1-J2p,'fro')/ n_J1;
0157
0158 if n_J1_J2p > 1e-4;
0159 warning('J1-J2p is too big');
0160 end
0161
0162 function test4_npaa_3d
0163
0164 imdl = mk_common_model('n3r2');
0165 f1mdl = imdl.fwd_model;
0166 f1mdl.solve = @fwd_solve_1st_order;
0167 f1mdl.jacobian = @jacobian_adjoint;
0168 f1mdl.system_mat = @system_mat_1st_order;
0169
0170
0171 n2d = f1mdl.nodes( (f1mdl.nodes(:,3) == 0), 1:2);
0172 e2d = delaunayn(n2d);
0173 c_mdl = eidors_obj('fwd_model','2d','elems',e2d,'nodes',n2d);
0174
0175
0176 c2f= mk_coarse_fine_mapping( f1mdl, c_mdl );
0177 f2mdl= f1mdl;
0178 f2mdl.coarse2fine = c2f;
0179
0180 imdl.jacobian_bkgnd.value= .6;
0181 img = calc_jacobian_bkgnd( imdl);
0182
0183 t=cputime;
0184 J1= jacobian_adjoint( f1mdl, img)*c2f;
0185 fprintf('aa_calc - full: t=%f\n', cputime-t );
0186 n_J1 = norm(J1,'fro');
0187
0188 t=cputime;
0189 J2= jacobian_adjoint( f2mdl, img);
0190 fprintf('aa_calc - c2f: t=%f\n', cputime-t );
0191
0192 n_J1_J2 = norm(J1-J2,'fro')/ n_J1;
0193 if n_J1_J2 > 1e-4;
0194 warning('J1-J2 is too big');
0195 end
0196
0197
0198 t=cputime;
0199 J2p= jacobian_perturb( f2mdl, img);
0200 fprintf('perturb_j - c2f: t=%f\n', cputime-t );
0201
0202 n_J1_J2p = norm(J1-J2p,'fro')/ n_J1;
0203 if n_J1_J2p > 1e-4;
0204 warning('J1-J2p is too big');
0205 end
0206