test_c2f_jacobian

 Test calc of jacobian given coarse to fine mapping

 (C) Andy Adler 2008. Licenced under GPL v2 or v3
 $Id: test_c2f_jacobian.m 3213 2012-06-29 20:31:29Z aadler $

0001 function test_c2f_jacobian
0002 % Test calc of jacobian given coarse to fine mapping
0003 %
0004 % (C) Andy Adler 2008. Licenced under GPL v2 or v3
0005 % $Id: test_c2f_jacobian.m 3213 2012-06-29 20:31:29Z aadler $
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 % Fine model
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 % Create 2D FEM of all NODES with z=0
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 % Create f_mdl with coarse2fine mapping
0027 c2f= mk_coarse_fine_mapping( f1mdl, c_mdl );
0028 f2mdl= f1mdl;
0029 f2mdl.coarse2fine = c2f;
0030 
0031 % Just in case - define reconstruction
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]; % see iso_f_smooth: deg=1, w=1
0036 imdl.jacobian_bkgnd.value= .6;
0037 imdl.reconst_type= 'difference';
0038 imdl.fwd_model= f1mdl; % fine model
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 %%%%%%%%%% J2
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 %%%%%%%%%% J2p
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 %%%%%%%%%% J1p
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 %%%%%%%%%% J2
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 %%%%%%%%%% J2p
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 % Fine model
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 % Create 2D FEM of all NODES with z=0
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 % Create f_mdl with coarse2fine mapping
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