calc_noise_figure

 CALC_NOISE_FIGURE: calculate the noise amplification (NF) of an algorithm
 [NF,SE] = calc_noise_figure( inv_model, hp, iterations)
    inv_model  => inverse model object
    hp         => value of hyperparameter to use (if not spec
         then use the value of inv_model.hyperparameter.value)
    iterations => number of iterations (default 10)
       (for calculation of noise figure using random noise)
    if hp is [] => use the hp on the inv_model
  NF = calculated NF. SE = standard error on NF

 hp is specified, it will be used for the hyperparameter.
    Otherwise the inv_model.hyperparameter will be used.

 Noise Figure must be defined for a specific measurement
 In order to specify data, use
     inv_model.hyperparameter.tgt_data.meas_t1
     inv_model.hyperparameter.tgt_data.meas_t2
   to use a temporal solver (or the Kalman filter), the
   measurement to perform the NF calc must also be specified,
   using:
     inv_model.hyperparameter.tgt_data.meas_select
   otherwise, the middle measurement will be used

 In order to automatically simulate data, specify tgt_elems,
   containing a vector of elements to use
 
     inv_model.hyperparameter.tgt_elems

 [NF,SE] = calc_noise_figure( inv_model, vh, vi)
    interface provided for compatibility with the deprecated
    calc_noise_params

0001 function [NF,SE] = calc_noise_figure( inv_model, hp, iterations)
0002 % CALC_NOISE_FIGURE: calculate the noise amplification (NF) of an algorithm
0003 % [NF,SE] = calc_noise_figure( inv_model, hp, iterations)
0004 %    inv_model  => inverse model object
0005 %    hp         => value of hyperparameter to use (if not spec
0006 %         then use the value of inv_model.hyperparameter.value)
0007 %    iterations => number of iterations (default 10)
0008 %       (for calculation of noise figure using random noise)
0009 %    if hp is [] => use the hp on the inv_model
0010 %  NF = calculated NF. SE = standard error on NF
0011 %
0012 % hp is specified, it will be used for the hyperparameter.
0013 %    Otherwise the inv_model.hyperparameter will be used.
0014 %
0015 % Noise Figure must be defined for a specific measurement
0016 % In order to specify data, use
0017 %     inv_model.hyperparameter.tgt_data.meas_t1
0018 %     inv_model.hyperparameter.tgt_data.meas_t2
0019 %   to use a temporal solver (or the Kalman filter), the
0020 %   measurement to perform the NF calc must also be specified,
0021 %   using:
0022 %     inv_model.hyperparameter.tgt_data.meas_select
0023 %   otherwise, the middle measurement will be used
0024 %
0025 % In order to automatically simulate data, specify tgt_elems,
0026 %   containing a vector of elements to use
0027 %
0028 %     inv_model.hyperparameter.tgt_elems
0029 %
0030 % [NF,SE] = calc_noise_figure( inv_model, vh, vi)
0031 %    interface provided for compatibility with the deprecated
0032 %    calc_noise_params
0033 
0034 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0035 % $Id: calc_noise_figure.m 5867 2018-12-09 18:36:39Z aadler $
0036 
0037 % A normal definition of noise power is based on power:
0038 %      NF = SNR_z / SNR_x
0039 %    SNR_z = sumsq(z0) / var(z) = sum(z0.^2) / trace(Rn)
0040 %    SNR_x = sumsq(x0) / var(x) = sum(x0.^2) / trace(ABRnB'A')
0041 % but his doesn't work, because the signal spreads like
0042 % the amplitude, not like the power, thus
0043 %      NF = SNR_z / SNR_x
0044 %    SNR_z = sum(|z0|/len_z) / std(z/len_z)
0045 %    SNR_x = sum(|x0|/len_x) / std(x/len_x)
0046 
0047 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'), do_unit_test, return, end
0048 
0049 if nargin>=2 && numel(hp) == 1
0050    inv_model.hyperparameter.value= hp;
0051 % Remove function parameter because it will recurse
0052    try; inv_model.hyperparameter = rmfield(inv_model.hyperparameter,'func'); end
0053 end
0054 if nargin == 3 && numel(hp) > 1
0055     h_data = hp;
0056     c_data = iterations;
0057     if isstruct(c_data) && strcmp(c_data.type,'data'); 
0058        c_data = c_data.meas;
0059     end
0060 elseif nargin == 3 && isstruct(hp);
0061     if strcmp(hp.type,'data'); 
0062        h_data = hp.meas;
0063     else
0064        error('expecting object of type "data"')
0065     end
0066     c_data = iterations;
0067     if isstruct(c_data) && strcmp(c_data.type,'data'); 
0068        c_data = c_data.meas;
0069     end
0070 else
0071     [inv_model, h_data, c_data] = process_parameters( inv_model );
0072 end
0073 
0074 % disable solution checking (slow)
0075 inv_model.inv_solve.calc_solution_error = false;
0076 
0077 %NF= nf_calc_use_matrix( inv_model, h_data, c_data);
0078 %NF= nf_calc_iterate( inv_model, h_data, c_data);
0079 if nargin<3; iterations= 10; end
0080 solver = inv_model.solve;
0081 if ischar(solver) && strcmp(solver, 'eidors_default')
0082     solver = eidors_default('get','inv_solve');
0083 end
0084 if isa(solver,'function_handle')
0085     solver = func2str(solver);
0086 end
0087 switch solver
0088     case {'inv_solve_backproj'
0089             'inv_solve_conj_grad'
0090             'inv_solve_diff_GN_one_step'
0091             'inv_solve_trunc_iterative'
0092             'inv_solve_TSVD'
0093             'solve_use_matrix'}
0094         NF = nf_calc_linear( inv_model, h_data, c_data);
0095         SE = 0;
0096     otherwise
0097         [NF,SE]= nf_calc_random( inv_model, h_data, c_data, iterations);
0098 end
0099 eidors_msg('@@ NF= %f', NF, 2);
0100 
0101 function [inv_model, h_data, c_data] = process_parameters( inv_model );
0102 
0103    if     isfield(inv_model.hyperparameter,'tgt_elems')
0104       [h_data, c_data]= simulate_targets( inv_model.fwd_model, ...
0105            inv_model.hyperparameter.tgt_elems);
0106    elseif isfield(inv_model.hyperparameter,'tgt_data')
0107       tgt_data= inv_model.hyperparameter.tgt_data;
0108       h_data= tgt_data.meas_t1;
0109       c_data= tgt_data.meas_t2;
0110    else
0111       error('unsure how to get data to measure signal');
0112    end
0113 
0114    % if hp is specified, then use that value
0115 
0116    if isfield(inv_model.hyperparameter,'func')
0117       funcname= inv_model.hyperparameter.func;
0118       if strcmp( class(funcname), 'function_handle')
0119          funcname= func2str(funcname);
0120       end
0121 
0122       if strcmp(funcname, 'choose_noise_figure')
0123          error('specifying inv_model.hp.func = choose_noise_figure will recurse');
0124       end
0125    end
0126 
0127 function params = nf_calc_linear(imdl, vh, vi )
0128 % params = GREIT_noise_params(imdl, homg_voltage, sig_voltage)
0129 %  params(1,:) = Noise Figure = SNR(image) / SNR(data)
0130 %
0131 %  see also: eval_GREIT_fig_merit or using test_performance
0132 
0133 % (C) 2008 Andy Adler. Licensed under GPL v2 or v3
0134 % $Id: calc_noise_figure.m 5867 2018-12-09 18:36:39Z aadler $
0135 
0136 % NOTE THAT WE ASSUME A LINEAR ALGORITHM FOR THIS MEASURE
0137 
0138 if 0 % old code with random noise
0139      % Keep this to validate while we test it
0140    Nnoise = 1000;
0141    noise = 0.01*std(vh)*randn(size(vh,1),Nnoise);
0142    vhn= vh*ones(1,Nnoise) + noise;
0143 else % use independent noise model on each channel
0144    noise = 0.01*std(vh)*speye(size(vh,1));
0145    vhn= vh*ones(1,size(vh,1)) + noise;
0146 end
0147 
0148 signal_y = calc_difference_data( vh, vi,  imdl.fwd_model);
0149 noise_y  = calc_difference_data( vh, vhn, imdl.fwd_model);
0150 
0151 signal_x = inv_solve(imdl, vh, vi);  
0152 signal_x = data_mapper(signal_x); signal_x = signal_x.elem_data;
0153 noise_x  = inv_solve(imdl, vh, vhn); 
0154 noise_x  = data_mapper(noise_x);  noise_x  = noise_x.elem_data;
0155 
0156 use_rec = 1;
0157 try 
0158    use_rec = ~imdl.prior_use_fwd_not_rec;
0159 end
0160 if use_rec
0161    try
0162       VOL = get_elem_volume(imdl.rec_model);
0163    catch
0164       VOL = get_elem_volume(imdl.fwd_model);
0165    end
0166 else
0167    VOL = get_elem_volume(imdl.fwd_model);
0168 end
0169 VOL = spdiags(VOL,0, length(VOL), length(VOL));
0170 
0171 signal_x = VOL*signal_x;
0172 noise_x = VOL*noise_x;
0173 
0174 signal_x = mean(abs(signal_x),1);
0175 noise_x  = mean(std(noise_x));
0176 snr_x = signal_x / noise_x;
0177 
0178 signal_y = mean(abs(signal_y),1);
0179 noise_y  = mean(std(noise_y)); 
0180 snr_y = signal_y / noise_y;
0181 
0182 params= [snr_y(:)./snr_x(:)]';
0183 
0184 try
0185 eidors_msg('@@ NF= %f (hp=%e)', NF, imdl.hyperparameter.value, 2);
0186 end
0187 
0188 
0189 % NOTES on the calculations: AA - Feb 20, 2012
0190 % SNR = mean(abs(x)); VAR =
0191 % ym= E[y]
0192 % Sy= E[(y-ym)*(y-ym)'] = E[y*y'] - ym*ym'
0193 % ny = sqrt(trace(Sy))
0194 % xm= E[x]  = E[R*y] = R*E[y] = R*ym
0195 % Sx= E[(x-xm)*(x-xm)'] = E[x*x'] - xm*xm'
0196 %   = E[R*ym*ym'*R'] = R*E[ym*ym']*R' = R*Sy*R'
0197 % nx = sqrt(trace(Sx))
0198 %
0199 % signal = mean(abs(x));
0200 %
0201 % In this case, these are exactly the same:
0202 %    noise_x  = mean(std(noise_x));
0203 %    noise_x  = sqrt(mean(noise_x.^2,2));
0204    
0205    
0206 function NF= nf_calc_use_matrix( inv_model, h_data, c_data)
0207 % To model std(z) we use z=z0+n
0208 % so that std(z) = sqrt(var(z)) = sqrt(1/L * E[n'*n])
0209 % we know a priori that the mean noise is zero, thus
0210 % we do not need to divide by L-1 in the variance
0211 % E[n'*n] = E[ trace(n*n') ] = trace( cov_N )
0212 % Thus, std(z) = sqrt( trace( cov_N )/L )
0213 %              = sqrt( mean( diag( cov_N )))
0214 % And,  std(x) = sqrt( mean( diag( cov_X )))
0215 %
0216 % To model cov_N, we consider independant noise
0217 %  on each channel, cov_N= N*N', N=diag( sigma_i )
0218 % And,              cov_X= X*X', X=reconst(N)
0219 % where X= reconst(mdl, z0,z0+N)
0220 %
0221 % To run efficiently mean(diag(cov_N))=mean(sum(N.^2,2))
0222 % The sum over X is actually weighted by the volume of
0223 %  each element, so VOL.^2*(sum(X.^2,2)
0224    try 
0225        VOL = get_elem_volume(inv_model.rec_model)';
0226    catch
0227        VOL = get_elem_volume(inv_model.fwd_model)';
0228    end
0229 
0230    % calculate signal
0231    d_len   = size(h_data,1);
0232    delta   = 1e-2* mean(h_data);
0233    c_noise = c_data*ones(1,d_len) + eye(d_len);
0234    h_full  = h_data*ones(1,d_len);
0235 
0236    sig_data = mean(abs( ...
0237          calc_difference_data( h_data, c_data , inv_model.fwd_model ) ...
0238                        ));
0239    var_data = mean(sum( ...
0240          calc_difference_data( h_full, c_noise, inv_model.fwd_model ) ...
0241                        .^2, 2)); 
0242                       
0243 
0244    % calculate image
0245    % Note, this won't work if the algorithm output is not zero biased
0246    [img0, img0n] = get_images( inv_model, h_data, c_data, ...
0247                                h_full, c_noise);
0248 
0249    i_len = length(img0);
0250    sig_img= VOL*abs(img0) / i_len;;
0251    var_img= VOL.^2*sum(img0n.^2 ,2) / i_len;
0252    
0253    NF = ( sig_data/ sqrt(var_data) ) / ( sig_img / sqrt(var_img)  );
0254 
0255    % For the record, the expression for var_img is derived as:
0256    % Equiv expresssions for var_img % given: A= diag(pp.VOLUME);
0257    % var_img= trace(A*RM*diag(Rn.^2)*RM'*A');
0258    % vv=A*RM*diag(Rn);var_img=trace(vv*vv'); var_img= sum(sum(vv.^2));
0259    % var_img= VOL2* (RM*diag(Rn)).^2
0260    % var_img= VOL2* RM.^2 * Rn.^2
0261 
0262    
0263 % simulate homg data and a small target in centre
0264 function [h_data, c_data]= simulate_targets( fwd_model, ctr_elems)
0265 
0266    homg= 1; % homogeneous conductivity level is 1
0267 
0268    %Step 1: homogeneous image
0269    sigma= homg*ones( size(fwd_model.elems,1) ,1);
0270 
0271    img= eidors_obj('image', 'homogeneous image', ...
0272                    'elem_data', sigma, ...
0273                    'fwd_model', fwd_model );
0274    h_data=fwd_solve( img );
0275    h_data= h_data.meas;
0276 
0277    %Step 1: inhomogeneous image with contrast in centre
0278    delta = 1e-2;
0279    sigma(ctr_elems) = homg*(1 + delta);
0280    img.elem_data = sigma;
0281    c_data=fwd_solve( img );
0282    c_data= c_data.meas;
0283 
0284 function [img0, img0n] = get_images( inv_model, h_data, c_data, ...
0285                                h_full, c_noise);
0286    if isa(inv_model.solve,'function_handle')
0287       solve= func2str(inv_model.solve);
0288    else
0289       solve= inv_model.solve;
0290    end
0291 
0292 % Test for special functions and solve them specially
0293    switch solve
0294    case 'ab_tv_diff_solve'
0295       error('Dont know how to calculate TV noise figure')
0296 
0297    case 'inv_solve_diff_kalman'
0298       inv_model.inv_solve_diff_kalman.keep_K_k1= 1;
0299       stablize = 6;
0300       img0 = inv_solve( inv_model, h_data, ...
0301                                    c_data*ones(1,stablize) );
0302       K= img0.inv_solve_diff_kalman.K_k1;
0303       img0.elem_data = K*calc_difference_data( h_data , c_data , inv_model.fwd_model);
0304       img0n.elem_data= K*calc_difference_data( h_full , c_noise, inv_model.fwd_model);
0305 
0306    otherwise
0307       img0 = inv_solve( inv_model, h_data, c_data);
0308       if nargin>4
0309       img0n= inv_solve( inv_model, h_full, c_noise);
0310       end
0311    end
0312 
0313    % Need elem or nodal data
0314    if isfield(img0,'node_data');
0315       img0 = img0.node_data;
0316    else
0317       img0 = img0.elem_data;
0318    end
0319 
0320    if isfield(img0n,'node_data');
0321       img0n = img0n.node_data;
0322    else
0323       img0n = img0n.elem_data;
0324    end
0325 
0326 % OLD CODE - iterate
0327 function NF= nf_calc_iterate( inv_model, h_data, c_data);
0328    try 
0329        VOL = get_elem_volume(inv_model.rec_model)';
0330    catch
0331        VOL = get_elem_volume(inv_model.fwd_model)';
0332    end
0333    % calculate signal
0334    d_len   = size(h_data,1);
0335    delta   = 1e-2* mean(h_data);
0336    sig_data = mean(abs( ...
0337          calc_difference_data( h_data, c_data , inv_model.fwd_model ) ...
0338                        ));
0339    % calculate image
0340    % Note, this won't work if the algorithm output is not zero biased
0341 
0342    [img0] = get_images( inv_model, h_data, c_data);
0343 %  sig_img= mean(VOL'.*abs(img0.elem_data));
0344    sig_img= VOL*abs(img0) / length(img0);
0345 
0346    % Now do noise
0347    var_data= 0;
0348    var_img = 0;
0349    for i=1:d_len
0350       this_noise = -ones(d_len, size(c_data,2))/(d_len-1);
0351       this_noise(i,:) = 1;
0352       c_noise = c_data + this_noise;
0353       [imgn] = get_images( inv_model, h_data, c_noise);
0354       if 1
0355          var_data = var_data + mean(sum( ...
0356             calc_difference_data( h_data, c_noise, inv_model.fwd_model ) ...
0357                           .^2, 2)); 
0358 %        var_img= var_img +  mean( (VOL'.*imgn.elem_data).^2 );
0359          var_img= var_img +  (VOL.^2)*sum(imgn.elem_data.^2,2 ) / length(imgn.elem_data); 
0360       else
0361          % OLD APPROACH BASED ON variance, rather than matrix calcs
0362          var_data = var_data + var( ...
0363             calc_difference_data( h_data, c_noise, inv_model.fwd_model ) ...
0364                                  ); 
0365          var_img= var_img + var( VOL'.*imgn.elem_data ); 
0366       end
0367    end
0368    var_data = var_data / d_len;
0369    var_img  = var_img  / d_len;
0370    NF = ( sig_data/ sqrt(var_data) ) / ( sig_img / sqrt(var_img)  );
0371 
0372 function [NF,SE]= nf_calc_random( rec, vh, vi, N_RUNS);
0373    eidors_cache('boost_priority',-2); % low priority values
0374 
0375    imgr= inv_solve(rec, vh, vi);
0376 
0377    if isfield(imgr,'node_data');
0378       img0 = imgr.node_data;
0379       try
0380           VOL = get_elem_volume(rec.rec_model, 1);
0381       catch
0382           VOL = get_elem_volume(rec.fwd_model, 1);
0383       end
0384    else
0385       n_els = length(rec.fwd_model);
0386       img0 = imgr.elem_data(1:n_els); % elem_data can also contain movement
0387       try
0388           VOL = get_elem_volume(rec.rec_model, 0);
0389       catch
0390           VOL = get_elem_volume(rec.fwd_model, 0);
0391       end
0392    end
0393 
0394    sig_ampl = mean( abs( VOL .* img0 )) / ...
0395               mean( abs( calc_difference_data( vh, vi, rec.fwd_model )));
0396 
0397 % Estimate Signal Amplitude
0398    for i=1:N_RUNS
0399       vn= addnoise(vh, vi, 1.0);
0400 
0401       imgr= inv_solve(rec, vh, vn);
0402 
0403       if isfield(imgr,'node_data'); img0 = imgr.node_data;
0404       else;                         img0 = imgr.elem_data(1:n_els);
0405       end
0406 
0407       noi_imag(i) = std( VOL .* img0 );
0408       noi_sgnl(i) = std( calc_difference_data( vh, vn, rec.fwd_model ));
0409    end
0410 
0411    noi_ampl = noi_imag./noi_sgnl;
0412    NF =  mean(noi_ampl/sig_ampl);
0413    SE =  std(noi_ampl/sig_ampl)/sqrt(N_RUNS);
0414    eidors_msg('NF= %f+/-%f', NF, SE, 3);
0415 
0416    eidors_cache('boost_priority',2);
0417 
0418 function noise= addnoise( vh, vi, SNR);
0419       if isstruct(vh); vh= vh.meas; end
0420       if isstruct(vi); vi= vi.meas; end
0421       noise = randn(size(vh));
0422       noise = noise*std(vh-vi)/std(noise);
0423       noise = vh + SNR*noise;
0424 
0425 function do_unit_test
0426     ll = eidors_msg('log_level',1);
0427     test1; % can we deal with c2f ?
0428     imdl = mk_common_model('a2t2',16); test2(imdl);
0429     imdl = mk_common_model('d2t2',16); test2(imdl);
0430     imdl = mk_common_model('a2c2',16); test2(imdl);
0431     imdl = mk_common_model('d2c2',16); test2(imdl);
0432     ll = eidors_msg('log_level',ll);
0433 
0434 function test1
0435     % big model with c2f and supposedly an NF of 0.5
0436     fmdl = mk_library_model('pig_23kg_16el');
0437     [fmdl.stimulation fmdl.meas_select] = mk_stim_patterns(16,1,'{ad}','{ad}');
0438     fmdl = mdl_normalize(fmdl, 1);  % Use normalized difference imaging
0439     opt.noise_figure = 0.5; opt.imgsz = [64 64];
0440     imdl = mk_GREIT_model(fmdl, 0.25, [], opt);
0441     % homogeneous measurement
0442     img = mk_image(fmdl,1);
0443     vh = fwd_solve(img);
0444     % inhomogeneous measurement
0445     select_fcn = inline('(x-0).^2+(y-0).^2+(z-0.5).^2<0.1^2','x','y','z');
0446     mfrac = elem_select(fmdl, select_fcn);
0447     img.elem_data = img.elem_data + mfrac*0.1;
0448     vi = fwd_solve(img);
0449     
0450     nf1 = calc_noise_params(imdl, vh.meas, vi.meas);
0451     
0452     % We use different measurements so naturally we don't get 0.5 here
0453     imdl.hyperparameter.tgt_data.meas_t1 = vh.meas;
0454     imdl.hyperparameter.tgt_data.meas_t2 = vi.meas;
0455     try
0456         % calc_noise_figure doens't support dual models
0457         nf2 = calc_noise_figure(imdl);
0458     catch
0459         nf2 = 0;
0460     end
0461     unit_test_cmp('Noise fig implementations',nf1, nf2, 1e-2);
0462 
0463 function test2(imdl)
0464     fmdl = imdl.fwd_model;
0465     % homogeneous measurement
0466     img = mk_image(fmdl,1);
0467     vh = fwd_solve(img);
0468     % inhomogeneous measurement
0469     select_fcn = inline('(x-0).^2+(y-0).^2.^2<15^2','x','y','z');
0470     mfrac = elem_select(fmdl, select_fcn);
0471     img.elem_data = img.elem_data + mfrac*0.1;
0472     vi = fwd_solve(img);
0473 
0474     nf1 = calc_noise_params(imdl, vh.meas, vi.meas);
0475     eidors_msg(nf1,0);
0476 
0477 imdl.hyperparameter.tgt_data.meas_t1 = vh.meas;
0478 imdl.hyperparameter.tgt_data.meas_t2 = vi.meas;
0479 % calc_noise_figure doens't support dual models
0480 nf2 = calc_noise_figure(imdl,[],1000);
0481 unit_test_cmp('Noise fig implementations',nf1, nf2, 1e-2);
0482