update_system_mat_fields

 UPDATE_SYSTEM_MAT_FIELDS: fields (elem to nodes) fraction of system mat
 Create the fields for a small change (fmdl1) from a base model (fmdl0)
 If the change is small, this computation is faster
 FC= system_mat_fields( fmdl0, fmdl1)
 input: 
   fmdl0 = base forward model
   fmdl1 = new (updated) forward model
 output:
   FC:        s_mat= C' * S * conduct * C = FC' * conduct * FC;

0001 function FC1 = update_system_mat_fields( fwd_model0, fwd_model1 )
0002 % UPDATE_SYSTEM_MAT_FIELDS: fields (elem to nodes) fraction of system mat
0003 % Create the fields for a small change (fmdl1) from a base model (fmdl0)
0004 % If the change is small, this computation is faster
0005 % FC= system_mat_fields( fmdl0, fmdl1)
0006 % input:
0007 %   fmdl0 = base forward model
0008 %   fmdl1 = new (updated) forward model
0009 % output:
0010 %   FC:        s_mat= C' * S * conduct * C = FC' * conduct * FC;
0011 
0012 % (C) 2008 Andy Adler. License: GPL version 2 or version 3
0013 % $Id: update_system_mat_fields.m 6483 2022-12-26 11:31:22Z aadler $
0014 
0015 if ischar(fwd_model0) && strcmp(fwd_model0,'UNIT_TEST'); do_unit_test; return; end
0016 
0017 FC0 = system_mat_fields(fwd_model0); % this should be a cache hit
0018 
0019 copt.cache_obj= {mk_cache_obj(fwd_model0), ...
0020                  mk_cache_obj(fwd_model1)};
0021 copt.fstr = 'update_system_mat_fields';
0022 dFC = eidors_cache(@calc_update_system_mat_fields,{fwd_model0, fwd_model1},copt );
0023 
0024 % updated system_mat_fields
0025 FC1 = FC0 + dFC;
0026 
0027 % now we fake out the caching system by telling it that we are doing
0028 % system_mat_fields(fwd_model1) with the new model, calls to system_mat_fields
0029 % will get the cached value from this function
0030 fstr = 'system_mat_fields';
0031 cache_obj = mk_cache_obj(fwd_model1);
0032 eidors_obj('set-cache', cache_obj, fstr, {FC1});
0033 
0034 
0035 % only cache stuff which is really relevant here
0036 function cache_obj = mk_cache_obj(fwd_model)
0037    cache_obj.elems       = fwd_model.elems;
0038    cache_obj.nodes       = fwd_model.nodes;
0039    try
0040    cache_obj.electrode   = fwd_model.electrode; % if we have it
0041    end
0042    cache_obj.type        = 'fwd_model';
0043    cache_obj.name        = ''; % it has to have one
0044 
0045 function dFC= calc_update_system_mat_fields( fwd_model0, fwd_model1 )
0046    p0= fwd_model_parameters( fwd_model0, 'skip_VOLUME' );
0047    p1= fwd_model_parameters( fwd_model1, 'skip_VOLUME' );
0048    d0= p0.n_dims+0;
0049    d1= p0.n_dims+1;
0050    e= p0.n_elem;
0051    n= p0.n_node;
0052 
0053    % find moved nodes, then the elements affected by those moved nodes
0054    [dn,~] = find(abs(fwd_model0.nodes - fwd_model1.nodes) > eps);
0055    dn = unique(dn);
0056    [de,~] = find(ismember(fwd_model1.elems,dn)); % our modified nodes touched these elements
0057    de = unique(de);
0058    assert(all(all(fwd_model0.elems == fwd_model1.elems)), 'expected fmdl0 and fmdl1 to have the same element structure');
0059    assert(all(de <= e), 'invalid elem# found for delta nodes');
0060 
0061    % TODO: use vectorized inverse from see system_mat_fields
0062 
0063 if 1 % Slow way
0064    FFjidx= floor([0:d0*e-1]'/d0)*d1*ones(1,d1) + ones(d0*e,1)*(1:d1);
0065    FFiidx= [1:d0*e]'*ones(1,d1);
0066    FFdata= zeros(d0*e,d1);
0067    dfact = (d0-1)*d0;
0068    for j=de(:)'
0069      a0=  inv([ ones(d1,1), p0.NODE( :, p0.ELEM(:,j) )' ]);
0070      a1=  inv([ ones(d1,1), p1.NODE( :, p0.ELEM(:,j) )' ]);
0071      idx= d0*(j-1)+1 : d0*j;
0072      FFdata(idx,1:d1)= a1(2:d1,:)/ sqrt(dfact*abs(det(a1))) - a0(2:d1,:)/ sqrt(dfact*abs(det(a0)));
0073    end %for j=1:ELEMs
0074    % Remove non-zero components
0075    idx = all(FFdata==0,2);
0076    FFdata(idx,:) = [];
0077    FFiidx(idx,:) = [];
0078    FFjidx(idx,:) = [];
0079    
0080 else
0081    ded0 = bsxfun(@plus,[de(:)-1]*d0,1:d0)'; 
0082    ded0 = ded0(:); lded0= length(ded0);
0083    FFjidx= floor([ded0-1]/d0)*d1*ones(1,d1) + ...
0084            ones(lded0,1)*(1:d1);
0085    FFiidx= ded0*ones(1,d1);
0086    FFdata= zeros(lded0,d1);
0087    dfact = (d0-1)*d0;
0088    for i=1:length(de); j=de(i);
0089      a0=  inv([ ones(d1,1), p0.NODE( :, p0.ELEM(:,j) )' ]);
0090      a1=  inv([ ones(d1,1), p1.NODE( :, p0.ELEM(:,j) )' ]);
0091      idx= d0*(i-1)+1 : d0*i;
0092      FFdata(idx,1:d1)= a1(2:d1,:)/ sqrt(dfact*abs(det(a1))) - a0(2:d1,:)/ sqrt(dfact*abs(det(a0)));
0093    end %for j=1:ELEMs
0094 end
0095 
0096    % TODO this could be better: we could look to see which, if any electrodes
0097    % have been modified and only update those ones... currently the
0098    % implementation here is pretty brain dead on the assumption this part is
0099    % fast
0100    % TODO currently can't handle electrode node number changes
0101    [F2data0,F2iidx0,F2jidx0, C2data0,C2iidx0,C2jidx0] = compl_elec_mdl(fwd_model0,p0,dn);
0102    [F2data1,F2iidx1,F2jidx1, C2data1,C2iidx1,C2jidx1] = compl_elec_mdl(fwd_model1,p1,dn);
0103 
0104 if 1
0105    FF= sparse([FFiidx(:);  F2iidx0(:);  F2iidx1(:)],...
0106               [FFjidx(:);  F2jidx0(:);  F2jidx1(:)],...
0107               [FFdata(:); -F2data0(:); +F2data1(:)]);
0108 else
0109 % New sparse assembly. Is it faster?
0110    F1=sparse(FFiidx, FFjidx, ...
0111              FFdata, d0*e, d1*e);
0112    F2=sparse(F2iidx0,F2jidx0, ...
0113            -F2data0, d0*e, d1*e);
0114    F3= sparse(F2iidx1,F2jidx1, ...
0115              F2data1,d0*e, d1*e);
0116    FF = F1 + F2 + F3;
0117 end
0118    
0119    CC= sparse([(1:d1*e)';     C2iidx0(:)], ...
0120               [p0.ELEM(:);    C2jidx0(:)], ...
0121               [ones(d1*e,1);  C2data0(:)]);
0122 
0123    dFC= FF*CC;
0124 
0125 
0126 % Add parts for complete electrode model
0127 function [FFdata,FFiidx,FFjidx, CCdata,CCiidx,CCjidx] = ...
0128              compl_elec_mdl(fwd_model,pp,dn)
0129    d0= pp.n_dims;
0130    FFdata= zeros(0,d0);
0131    FFd_block= sqrtm( ( ones(d0) + eye(d0) )/6/(d0-1) ); % 6 in 2D, 12 in 3D
0132    FFiidx= zeros(0,d0);
0133    FFjidx= zeros(0,d0);
0134    FFi_block= ones(d0,1)*(1:d0);
0135    CCdata= zeros(0,d0);
0136    CCiidx= zeros(0,d0);
0137    CCjidx= zeros(0,d0);
0138   
0139    sidx= d0*pp.n_elem;
0140    cidx= (d0+1)*pp.n_elem;
0141    for i= 1:pp.n_elec
0142       eleci = fwd_model.electrode(i);
0143       % contact impedance zc is in [Ohm.m] for 2D or [Ohm.m^2] for 3D
0144       zc=  eleci.z_contact;
0145       if any(find(ismember(eleci.nodes,dn)))
0146          [bdy_idx, bdy_area] = find_electrode_bdy( ...
0147              pp.boundary, fwd_model.nodes, eleci.nodes );
0148              % bdy_area is in [m] for 2D or [m^2] for 3D
0149    
0150          for j= 1:length(bdy_idx);
0151             bdy_nds= pp.boundary(bdy_idx(j),:);
0152    
0153             % 3D: [m^2]/[Ohm.m^2] = [S]
0154             % 2D: [m]  /[Ohm.m]   = [S]
0155             FFdata= [FFdata; FFd_block * sqrt(bdy_area(j)/zc)];
0156             FFiidx= [FFiidx; FFi_block' + sidx];
0157             FFjidx= [FFjidx; FFi_block  + cidx];
0158    
0159             CCiidx= [CCiidx; FFi_block(1:2,:) + cidx];
0160             CCjidx= [CCjidx; bdy_nds ; (pp.n_node+i)*ones(1,d0)];
0161             CCdata= [CCdata; [1;-1]*ones(1,d0)];
0162             sidx = sidx + d0;
0163             cidx = cidx + d0;
0164          end
0165       else
0166          [bdy_idx] = find_electrode_bdy( ...
0167              pp.boundary, fwd_model.nodes, eleci.nodes );
0168              % bdy_area is in [m] for 2D or [m^2] for 3D
0169    
0170          for j= 1:length(bdy_idx);
0171             bdy_nds= pp.boundary(bdy_idx(j),:);
0172    
0173             % 3D: [m^2]/[Ohm.m^2] = [S]
0174             % 2D: [m]  /[Ohm.m]   = [S]
0175             FFdata= [FFdata; FFd_block * 0];
0176             FFiidx= [FFiidx; FFi_block' + sidx];
0177             FFjidx= [FFjidx; FFi_block  + cidx];
0178    
0179             CCiidx= [CCiidx; FFi_block(1:2,:) + cidx];
0180             CCjidx= [CCjidx; bdy_nds ; (pp.n_node+i)*ones(1,d0)];
0181             CCdata= [CCdata; [1;-1]*ones(1,d0)];
0182             sidx = sidx + d0;
0183             cidx = cidx + d0;
0184          end
0185       end
0186       
0187    end
0188 
0189 function do_unit_test
0190    disp('running system_mat_fields UNIT_TEST');
0191    system_mat_fields('UNIT_TEST'); % we depend explicitly on system_mat_fields... check it works properly!
0192 
0193    eidors_cache('clear_name', 'system_mat_fields');
0194    eidors_cache('clear_name', 'update_system_mat_fields');
0195    disp('running update_system_mat_fields UNIT_TEST');
0196    imdl=  mk_geophysics_model('h2e',32);
0197    ndim = size(imdl.fwd_model.nodes,2);
0198    for i = 1:10
0199       fmdl0 = imdl.fwd_model;
0200       fmdl0.nodes(1,:) = fmdl0.nodes(1,:) + rand(1,ndim)*1e-10; % defeat cache
0201       t=tic; FC0 = system_mat_fields(fmdl0); t0(i)=toc(t);
0202    
0203       % perturb nodes
0204       fmdl1 = fmdl0;
0205       nn = fmdl1.electrode(1).nodes;
0206       nn = [nn(1); nn(end)];
0207       fmdl1.nodes(nn,:) = fmdl1.nodes(nn,:) + 1e-4+rand(2,ndim)*1e-10; % perturb
0208       FC0 = system_mat_fields(fmdl0);
0209       fwd_model_parameters(fmdl1,'skip_VOLUME'); % cache
0210       t=tic; FC1a = update_system_mat_fields(fmdl0, fmdl1); t1(i)=toc(t);
0211       t=tic; FC1b = system_mat_fields(fmdl1); t2(i)=toc(t);
0212       eidors_cache('clear_name', 'system_mat_fields');
0213       t=tic; FC1c = system_mat_fields(fmdl1); t3(i)=toc(t);
0214    end
0215    fprintf('system_mat_fields(fmdl0) = %0.2f sec\n',mean(t0));
0216    fprintf('update_system_mat_fields(fmdl0,fmdl1) = %0.3f sec [faster?]\n',mean(t1));
0217    fprintf('system_mat_fields(fmdl1) = %0.3f sec [cache hit?]\n',mean(t2));
0218    fprintf('system_mat_fields(fmdl1+delta) = %0.3f sec [recalculate]\n',mean(t3));
0219 %  disp(mean([t0;t1;t2;t3]'))
0220    unit_test_cmp('delta FC is fast                  ',mean(t1) < mean(t0)/2,1);
0221    unit_test_cmp('cache trick is really fast        ',mean(t2./t0) < 0.030,1);
0222    unit_test_cmp('cache was cleared before recalc   ',mean(t3) > mean(t2)*5,1); % did actually clear cache
0223    unit_test_cmp('cache trick returns correct result',FC1a,FC1b);
0224    unit_test_cmp('delta FC gives same result        ',FC1a,FC1c,3e2*eps); % Variation since random
0225    fprintf('speed-up: %0.2f\n',mean(t0./t1));
0226    if(sum(sum((FC1a - FC1c).^2)) > 100*eps)
0227       err = FC1a - FC1c
0228    end
0229 
0230    if 0
0231       disp('---- profiling ----');
0232       fmdl1.nodes(nn,:) = fmdl1.nodes(nn,:) + rand(2,ndim)*1e-3; % perturb
0233       t=tic; FC0 = system_mat_fields(fmdl0); t0=toc(t);
0234       profile clear;
0235       profile on;
0236       t = tic; FC1a = update_system_mat_fields(fmdl0, fmdl1); t1=toc(t);
0237       profile off;
0238       profview;
0239       profsave(profile('info'),'profile');
0240       fprintf('update_system_mat_fields(fmdl0,fmdl1) = %0.3f sec [profiled]\n',t1);
0241    end
0242 
0243 function time_test_code
0244     f0= mk_common_model('e3cr',4);
0245     f0 = f0.fwd_model; f1=f0;
0246     f1.nodes(1,:)= f1.nodes(1,:) + 1e-6*[1,1,1];
0247     Fa= system_mat_fields(f1);
0248     for sp=1:10;eidors_cache clear;switch sp;
0249        case 1;
0250     tic; F1= update_system_mat_fields(f0,f1); toc
0251        case 2;
0252     tic; F0= system_mat_fields(f0); toc
0253     tic; F1= update_system_mat_fields(f0,f1); toc
0254        case 3;
0255     sv = 'skip_VOLUME';
0256     tic; fwd_model_parameters(f0,sv); toc
0257     tic; fwd_model_parameters(f0,sv); toc
0258     tic; fwd_model_parameters(f1,sv); toc
0259     tic; fwd_model_parameters(f1,sv); toc
0260     disp ----
0261     tic; F0= system_mat_fields(f0); toc
0262     tic; F0= system_mat_fields(f0); toc
0263     disp ----
0264     tic; F1= update_system_mat_fields(f0,f1); toc
0265     tic; F1= update_system_mat_fields(f0,f1); toc
0266     tic; Fb= system_mat_fields(f1); toc
0267     unit_test_cmp('',Fa,F1,1e-10);
0268        otherwise; break;
0269     end; end
0270