nodal_solve

 NODAL_SOLVE inverse solver using approach of Adler&Guardo 1996
 img= nodal_solve( inv_model, data1, data2)
 img        => output image (or vector of images)
 inv_model  => inverse model struct
 data1      => differential data at earlier time
 data2      => differential data at later time

 both data1 and data2 may be matrices (MxT) each of
  M measurements at T times
 if either data1 or data2 is a vector, then it is expanded
  to be the same size matrix

0001 function img= nodal_solve( inv_model, data1, data2)
0002 % NODAL_SOLVE inverse solver using approach of Adler&Guardo 1996
0003 % img= nodal_solve( inv_model, data1, data2)
0004 % img        => output image (or vector of images)
0005 % inv_model  => inverse model struct
0006 % data1      => differential data at earlier time
0007 % data2      => differential data at later time
0008 %
0009 % both data1 and data2 may be matrices (MxT) each of
0010 %  M measurements at T times
0011 % if either data1 or data2 is a vector, then it is expanded
0012 %  to be the same size matrix
0013 
0014 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0015 % $Id: nodal_solve.m 4833 2015-03-29 21:32:08Z bgrychtol-ipa $
0016 
0017 
0018 dv = calc_difference_data( data1, data2, inv_model.fwd_model);
0019 
0020 sol = eidors_cache(@get_RM, inv_model, 'nodal_solve' ) * dv;
0021 
0022 % create a data structure to return
0023 img.name= 'solved by nodal_solve';
0024 img.node_data = sol;
0025 img.fwd_model= inv_model.fwd_model;
0026 
0027 function one_step_inv = get_RM( inv_model );
0028    fwd_model= inv_model.fwd_model;
0029 
0030    img_bkgnd= calc_jacobian_bkgnd( inv_model );
0031    J = calc_jacobian(img_bkgnd);
0032 
0033    RtR = calc_RtR_prior( inv_model );
0034    W   = calc_meas_icov( inv_model );
0035    hp  = calc_hyperparameter( inv_model );
0036 
0037    [e2n,Ne, Nn] = elem2node( fwd_model.elems );
0038    if size(J,2) == Ne;
0039       J= J*e2n;
0040    end
0041    if size(RtR,2) == Ne;
0042       RtR = e2n'*RtR*e2n;
0043       RtR = RtR +  1e-4*spdiags(diag(RtR),0,Nn,Nn); %Need just a little regularization
0044    end
0045 
0046    one_step_inv= (J'*W*J +  hp^2*RtR)\J'*W;
0047 
0048 function [e2n, Ne, Nn] = elem2node( elems )
0049    [Ne,d] = size(elems);
0050    elemo = (1:Ne)'*ones(1,d);
0051    e2n = sparse(elemo,elems,1/d); 
0052    Nn = size(e2n,2);