calc_meas_icov

 meas_icov = calc_meas_icov( inv_model )
 CALC_MEAS_ICOV: calculate inverse covariance of measurements
   The meas_icov is a matrix n_meas x n_meas of the
     inverse of measurement covariances. Normally measurements
     are assumed to be independant, so the meas_icov is 
     a diagonal matrix of 1/var(meas)

 calc_meas_icov can be called as
    meas_icov= calc_meas_icov( inv_model )

 meas_icov   is the calculated data prior
 inv_model    is an inv_model structure

 if:
    inv_model.meas_icov    is a function
          -> call it to calculate meas_icov
    inv_model.meas_icov    is a matrix
          -> return it as meas_icov
    inv_model.meas_icov    does not exist
          -> use I, or 1./homg (for normalized difference)

0001 function meas_icov = calc_meas_icov( inv_model )
0002 % meas_icov = calc_meas_icov( inv_model )
0003 % CALC_MEAS_ICOV: calculate inverse covariance of measurements
0004 %   The meas_icov is a matrix n_meas x n_meas of the
0005 %     inverse of measurement covariances. Normally measurements
0006 %     are assumed to be independant, so the meas_icov is
0007 %     a diagonal matrix of 1/var(meas)
0008 %
0009 % calc_meas_icov can be called as
0010 %    meas_icov= calc_meas_icov( inv_model )
0011 %
0012 % meas_icov   is the calculated data prior
0013 % inv_model    is an inv_model structure
0014 %
0015 % if:
0016 %    inv_model.meas_icov    is a function
0017 %          -> call it to calculate meas_icov
0018 %    inv_model.meas_icov    is a matrix
0019 %          -> return it as meas_icov
0020 %    inv_model.meas_icov    does not exist
0021 %          -> use I, or 1./homg (for normalized difference)
0022 
0023 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0024 % $Id: calc_meas_icov.m 6718 2024-04-02 18:15:20Z aadler $
0025 
0026 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0027 
0028 if isfield(inv_model,'meas_icov')
0029    if isnumeric(inv_model.meas_icov);
0030       meas_icov = inv_model.meas_icov;
0031    else
0032       try inv_model.meas_icov = str2func(inv_model.meas_icov); end
0033       meas_icov= eidors_cache( inv_model.meas_icov, {inv_model});
0034    end
0035 else
0036    meas_icov= eidors_cache(@default_meas_icov,{inv_model} );
0037 end
0038 
0039  
0040 % Calculate a data prior for an assumption of uniform noise
0041 % on each channel
0042 %
0043 function meas_icov = default_meas_icov( inv_model )
0044 
0045    fwd_model= inv_model.fwd_model;
0046    fwd_model = mdl_normalize(fwd_model,mdl_normalize(fwd_model));
0047 
0048    n =  calc_n_meas( fwd_model );
0049 
0050    if ~mdl_normalize(fwd_model);
0051       meas_icov= speye( n );
0052    else
0053       homg_data=  solve_homg_image( fwd_model );
0054 % if we normalize, then small data get increased
0055 %   this means that noise on small data gets increased,
0056 %    so the covariance is large when data are small
0057 %   so the icov is small when data are small
0058 % sig = k/h -> std = k/h -> 1/std = kh
0059       meas_icov = sparse(1:n, 1:n, ( homg_data.meas ).^2 );
0060    end
0061 
0062 function n_meas = calc_n_meas( fwd_model )
0063 
0064    n_meas = 0;
0065    for i= 1:length(fwd_model.stimulation );
0066        n_meas = n_meas + size(fwd_model.stimulation(i).meas_pattern,1);
0067    end
0068 
0069 % create homogeneous image + simulate data
0070 function homg_data = solve_homg_image( fwd_mdl )
0071     n_elems= size( fwd_mdl.elems , 1);
0072     mat= ones( n_elems, 1);
0073     homg_img= eidors_obj('image', 'homogeneous image', ...
0074                          'elem_data', mat, 'fwd_model', fwd_mdl );
0075     homg_data=fwd_solve( homg_img);
0076 
0077 function do_unit_test
0078    mdl= getfield( mk_common_model('a2c2',16), 'fwd_model');;
0079    run_dataprior_test( mdl );
0080 
0081    mdl = mdl_normalize(mdl,1); mdl.name = [mdl.name,':normalize'];
0082    run_dataprior_test( mdl );
0083 
0084    mdl= getfield( mk_common_model('n3r2',[16,2]), 'fwd_model');;
0085    run_dataprior_test( mdl );
0086    mdl = mdl_normalize(mdl,1); mdl.name = [mdl.name,':normalize'];
0087    run_dataprior_test( mdl );
0088 
0089 % test dataprior
0090 %     Difference dataprior should be 1
0091 %     normalized difference dataprior should be homg_data.^2
0092 function ok= run_dataprior_test( mdl )
0093     img= mk_image(mdl,1);
0094     homg_data=fwd_solve( img);
0095 
0096     DP= calc_meas_icov( img );
0097 
0098     % difference dataprior
0099     testvec= diag(DP);
0100     if mdl_normalize(mdl)
0101 %  from calc_meas_icov, we have the following
0102 %     meas_icov = sparse(1:n, 1:n, ( homg_data.meas ).^2 );
0103         testvec = homg_data.meas.^2 ./ diag(DP);
0104     end
0105 
0106     unit_test_cmp(['dataprior ',mdl.name],testvec,1,1e-10);
0107