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 FC = 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 6747 2024-04-05 18:29:27Z aadler $
0014 
0015 if ischar(fwd_model0) && strcmp(fwd_model0,'UNIT_TEST'); do_unit_test; return; end
0016 
0017 copt.cache_obj= {mk_cache_obj(fwd_model0), ...
0018                  mk_cache_obj(fwd_model1)};
0019 copt.fstr = 'update_system_mat_fields';
0020 dFC = eidors_cache(@calc_update_system_mat_fields,{fwd_model0, fwd_model1},copt );
0021 
0022 % update system_mat_fields
0023 FC = system_mat_fields(fwd_model0); % this should be a cache hit
0024 FC = FC + dFC;
0025 
0026 % now we fake out the caching system by telling it that we are doing
0027 % system_mat_fields(fwd_model1) with the new model, calls to system_mat_fields
0028 % will get the cached value from this function
0029 fstr = 'system_mat_fields';
0030 cache_obj = mk_cache_obj(fwd_model1);
0031 eidors_obj('set-cache', cache_obj, fstr, {FC});
0032 
0033 
0034 % only cache stuff which is really relevant here
0035 function cache_obj = mk_cache_obj(fwd_model)
0036    cache_obj.elems       = fwd_model.elems;
0037    cache_obj.nodes       = fwd_model.nodes;
0038    try
0039    cache_obj.electrode   = fwd_model.electrode; % if we have it
0040    end
0041    cache_obj.type        = 'fwd_model';
0042    cache_obj.name        = ''; % it has to have one
0043 
0044 function dFC= calc_update_system_mat_fields( fwd_model0, fwd_model1 )
0045    p0= fwd_model_parameters( fwd_model0, 'skip_VOLUME' );
0046    p1= fwd_model_parameters( fwd_model1, 'skip_VOLUME' );
0047    e= p0.n_elem;
0048    n= p0.n_node;
0049    d0= p0.n_dims+0;
0050    d1= p0.n_dims+1;
0051 
0052    % find moved nodes, then the elements affected by those moved nodes
0053    [dn,~] = find(abs(fwd_model0.nodes - fwd_model1.nodes) > eps);
0054    dn = unique(dn);
0055    [de,~] = find(ismember(fwd_model1.elems,dn)); % our modified nodes touched these elements
0056    de = unique(de);
0057    assert(all(all(fwd_model0.elems == fwd_model1.elems)), 'expected fmdl0 and fmdl1 to have the same element structure');
0058    assert(all(de <= e), 'invalid elem# found for delta nodes');
0059    [FFiidx,FFjidx,FFdata] = assemble_elements(p1,p0,de(:),3); % Use vectorized mode=3
0060 
0061    % TODO this could be better: we could look to see which, if any electrodes
0062    % have been modified and only update those ones... currently the
0063    % implementation here is pretty brain dead on the assumption this part is
0064    % fast
0065    % TODO currently can't handle electrode node number changes
0066    [F2data0,F2iidx0,F2jidx0, C2data0,C2iidx0,C2jidx0] = compl_elec_mdl(fwd_model0,p0,dn);
0067    [F2data1,F2iidx1,F2jidx1, C2data1,C2iidx1,C2jidx1] = compl_elec_mdl(fwd_model1,p1,dn);
0068 
0069 Cmax = max([d0*e; max(FFiidx(:)); max(F2iidx0(:)); F2iidx1(:)]);
0070 Rmax = max([d1*e; max(FFjidx(:)); max(F2jidx0(:)); F2jidx1(:)]);
0071 if 1
0072    FF= sparse([FFiidx(:);  F2iidx0(:);  F2iidx1(:)],...
0073               [FFjidx(:);  F2jidx0(:);  F2jidx1(:)],...
0074               [FFdata(:); -F2data0(:); +F2data1(:)],Cmax,Rmax);
0075 else
0076 % New sparse assembly. Is it faster?
0077    F1=sparse(FFiidx, FFjidx,  FFdata,  Cmax, Rmax);
0078    F2=sparse(F2iidx0,F2jidx0,-F2data0, Cmax, Rmax);
0079    F3= sparse(F2iidx1,F2jidx1,F2data1, Cmax, Rmax);
0080    FF = F1 + F2 + F3;
0081 end
0082    
0083 if 0 % don't assemble whole CC matrix
0084    CC= sparse([(1:d1*e)';     C2iidx0(:)], ...
0085               [p0.ELEM(:);    C2jidx0(:)], ...
0086               [ones(d1*e,1);  C2data0(:)]);
0087 else % only assemble the parts we need
0088    ded0 = reshape(((de-1)*d1 + [1:d1])', [],1);
0089    deels= reshape(p0.ELEM(:,de),[],1);
0090    CC= sparse([ded0;     C2iidx0(:)], ...
0091               [deels;    C2jidx0(:)], ...
0092               [ded0*0+1; C2data0(:)],Rmax,max([n;max(C2jidx0(:))]));
0093 end
0094    dFC= FF*CC;
0095 
0096 function [FFiidx,FFjidx,FFdata] = assemble_elements(p1,p0,de,approach)
0097    if nargin<4; approach=1; end
0098 
0099    d0= p0.n_dims+0;
0100    d1= p0.n_dims+1;
0101    e= p0.n_elem;
0102    n= p0.n_node;
0103    switch approach
0104    case 3; % Use vectorize
0105       ded0 = reshape(((de-1)*d0 + [1:d0])', [],1);
0106       FFjidx= floor((ded0-1)/d0)*d1 + (1:d1);
0107       FFiidx= ded0*ones(1,d1);
0108       FFdata = assemble_vectorize(p1,de) - assemble_vectorize(p0,de);
0109 
0110    case 1; %slow way
0111       FFjidx= floor([0:d0*e-1]'/d0)*d1*ones(1,d1) + ones(d0*e,1)*(1:d1);
0112       FFiidx= [1:d0*e]'*ones(1,d1);
0113       dfact = (d0-1)*d0;
0114       FFdata= zeros(d0*e,d1);
0115       for j=de(:)'
0116         a0=  inv([ ones(d1,1), p0.NODE( :, p0.ELEM(:,j) )' ]);
0117         a1=  inv([ ones(d1,1), p1.NODE( :, p0.ELEM(:,j) )' ]);
0118         idx= d0*(j-1)+1 : d0*j;
0119         FFdata(idx,1:d1)= a1(2:d1,:)/ sqrt(dfact*abs(det(a1))) -  ...
0120                           a0(2:d1,:)/ sqrt(dfact*abs(det(a0)));
0121       end %for j=1:ELEMs
0122    
0123       % Remove non-zero components
0124       idx = all(FFdata==0,2);
0125       FFdata(idx,:) = [];
0126       FFiidx(idx,:) = [];
0127       FFjidx(idx,:) = [];
0128    case 2; %improved
0129       ded0 = bsxfun(@plus,[de(:)-1]*d0,1:d0)'; 
0130       ded0 = ded0(:); lded0= length(ded0);
0131       FFjidx= floor([ded0-1]/d0)*d1*ones(1,d1) + ...
0132               ones(lded0,1)*(1:d1);
0133       FFiidx= ded0*ones(1,d1);
0134       FFdata= zeros(lded0,d1);
0135       dfact = (d0-1)*d0;
0136       for i=1:length(de); j=de(i);
0137         a0=  inv([ ones(d1,1), p0.NODE( :, p0.ELEM(:,j) )' ]);
0138         a1=  inv([ ones(d1,1), p1.NODE( :, p0.ELEM(:,j) )' ]);
0139         idx= d0*(i-1)+1 : d0*i;
0140         FFdata(idx,1:d1)= a1(2:d1,:)/ sqrt(dfact*abs(det(a1))) - a0(2:d1,:)/ sqrt(dfact*abs(det(a0)));
0141       end %for j=1:ELEMs
0142    otherwise error('no such path');
0143    end
0144 
0145 function FFdata = assemble_vectorize(p,de);
0146    d0= p.n_dims+0;
0147    d1= p.n_dims+1;
0148    M3 = ones(d1,d1,length(de));
0149    M3(2:end,:,:)= reshape(p.NODE(:,p.ELEM(:,de)),d0,d1,[]);
0150 % To stabilize inverse, remove average
0151 % This is a multiple of row1 = 1
0152 % This affects row one only
0153    M3(2:end,:,:) = bsxfun(@minus, M3(2:end,:,:), ...
0154               mean(M3(2:end,:,:),2));
0155    [I,D] = vectorize_inverse(M3);
0156    dfact = (d0-1)*d0; % d factorial for 2 or 3d
0157    Is = bsxfun(@times, I, sqrt(abs(D)/dfact));
0158    Is = permute(Is(:,2:d1,:),[2,3,1]);
0159    FFdata= reshape(Is,[],d1);
0160 
0161 % Add parts for complete electrode model
0162 function [FFdata,FFiidx,FFjidx, CCdata,CCiidx,CCjidx] = ...
0163              compl_elec_mdl(fwd_model,pp,dn)
0164    d0= pp.n_dims;
0165    FFdata= zeros(0,d0);
0166    FFd_block= sqrtm( ( ones(d0) + eye(d0) )/6/(d0-1) ); % 6 in 2D, 12 in 3D
0167    FFiidx= zeros(0,d0);
0168    FFjidx= zeros(0,d0);
0169    FFi_block= ones(d0,1)*(1:d0);
0170    CCdata= zeros(0,d0);
0171    CCiidx= zeros(0,d0);
0172    CCjidx= zeros(0,d0);
0173   
0174    sidx= d0*pp.n_elem;
0175    cidx= (d0+1)*pp.n_elem;
0176    for i= 1:pp.n_elec
0177       eleci = fwd_model.electrode(i);
0178       % contact impedance zc is in [Ohm.m] for 2D or [Ohm.m^2] for 3D
0179       zc=  eleci.z_contact;
0180       if any(find(ismember(eleci.nodes,dn)))
0181          [bdy_idx, bdy_area] = find_electrode_bdy( ...
0182              pp.boundary, fwd_model.nodes, eleci.nodes );
0183              % bdy_area is in [m] for 2D or [m^2] for 3D
0184    
0185          for j= 1:length(bdy_idx);
0186             bdy_nds= pp.boundary(bdy_idx(j),:);
0187    
0188             % 3D: [m^2]/[Ohm.m^2] = [S]
0189             % 2D: [m]  /[Ohm.m]   = [S]
0190             FFdata= [FFdata; FFd_block * sqrt(bdy_area(j)/zc)];
0191             FFiidx= [FFiidx; FFi_block' + sidx];
0192             FFjidx= [FFjidx; FFi_block  + cidx];
0193    
0194             CCiidx= [CCiidx; FFi_block(1:2,:) + cidx];
0195             CCjidx= [CCjidx; bdy_nds ; (pp.n_node+i)*ones(1,d0)];
0196             CCdata= [CCdata; [1;-1]*ones(1,d0)];
0197             sidx = sidx + d0;
0198             cidx = cidx + d0;
0199          end
0200       else
0201          [bdy_idx] = find_electrode_bdy( ...
0202              pp.boundary, fwd_model.nodes, eleci.nodes );
0203              % bdy_area is in [m] for 2D or [m^2] for 3D
0204    
0205          for j= 1:length(bdy_idx);
0206             bdy_nds= pp.boundary(bdy_idx(j),:);
0207    
0208             % 3D: [m^2]/[Ohm.m^2] = [S]
0209             % 2D: [m]  /[Ohm.m]   = [S]
0210             FFdata= [FFdata; FFd_block * 0];
0211             FFiidx= [FFiidx; FFi_block' + sidx];
0212             FFjidx= [FFjidx; FFi_block  + cidx];
0213    
0214             CCiidx= [CCiidx; FFi_block(1:2,:) + cidx];
0215             CCjidx= [CCjidx; bdy_nds ; (pp.n_node+i)*ones(1,d0)];
0216             CCdata= [CCdata; [1;-1]*ones(1,d0)];
0217             sidx = sidx + d0;
0218             cidx = cidx + d0;
0219          end
0220       end
0221       
0222    end
0223 
0224 function do_unit_test
0225    disp('running system_mat_fields UNIT_TEST');
0226    system_mat_fields('UNIT_TEST'); % we depend explicitly on system_mat_fields... check it works properly!
0227 
0228    eidors_cache('clear_name', 'system_mat_fields');
0229    eidors_cache('clear_name', 'update_system_mat_fields');
0230    disp('running update_system_mat_fields UNIT_TEST');
0231 
0232    fmdl0 = getfield(mk_common_model('a2c2',4),'fwd_model');
0233    fmdl1 = fmdl0; fmdl1.nodes([5,6],2)= fmdl1.nodes([5,6],2) + .01;
0234    SM0= system_mat_fields(fmdl0);
0235    SM1= system_mat_fields(fmdl1);
0236    SM1a = update_system_mat_fields(fmdl0,fmdl1);
0237    unit_test_cmp('simple test',SM1a,SM1,1e-14)
0238 
0239    imdl=  mk_geophysics_model('h2e',32);
0240    ndim = size(imdl.fwd_model.nodes,2);
0241    for i = 1:10
0242       fmdl0 = imdl.fwd_model;
0243       fmdl0.nodes(1,:) = fmdl0.nodes(1,:) + rand(1,ndim)*1e-10; % defeat cache
0244       t=tic; FC0 = system_mat_fields(fmdl0); t0(i)=toc(t);
0245    
0246       % perturb nodes
0247       fmdl1 = fmdl0;
0248       nn = fmdl1.electrode(1).nodes;
0249       nn = [nn(1); nn(end)];
0250       fmdl1.nodes(nn,:) = fmdl1.nodes(nn,:) + 1e-4+rand(2,ndim)*1e-10; % perturb
0251       FC0 = system_mat_fields(fmdl0);
0252       fwd_model_parameters(fmdl1,'skip_VOLUME'); % cache
0253       t=tic; FC1a = update_system_mat_fields(fmdl0, fmdl1); t1(i)=toc(t);
0254       t=tic; FC1b = system_mat_fields(fmdl1); t2(i)=toc(t);
0255       eidors_cache('clear_name', 'system_mat_fields');
0256       t=tic; FC1c = system_mat_fields(fmdl1); t3(i)=toc(t);
0257    end
0258    fprintf('system_mat_fields(fmdl0) = %0.2f sec\n',mean(t0));
0259    fprintf('update_system_mat_fields(fmdl0,fmdl1) = %0.3f sec [faster?]\n',mean(t1));
0260    fprintf('system_mat_fields(fmdl1) = %0.3f sec [cache hit?]\n',mean(t2));
0261    fprintf('system_mat_fields(fmdl1+delta) = %0.3f sec [recalculate]\n',mean(t3));
0262 %  disp(mean([t0;t1;t2;t3]'))
0263    fprintf('Time ratio: %5.3f / %5.3f\n', mean(t1), mean(t0));
0264    unit_test_cmp('delta FC is faster                ',mean(t1) < mean(t0),1);
0265 
0266    unit_test_cmp('cache trick is really fast        ',mean(t2./t0) < 0.030,1);
0267    unit_test_cmp('cache was cleared before recalc   ',mean(t3) > mean(t2)*5,1); % did actually clear cache
0268    unit_test_cmp('cache trick returns correct result',FC1a,FC1b);
0269    unit_test_cmp('delta FC gives same result        ',FC1a,FC1c,3e2*eps); % Variation since random
0270    fprintf('speed-up: %0.2f\n',mean(t0./t1));
0271    if(sum(sum((FC1a - FC1c).^2)) > 100*eps)
0272       err = FC1a - FC1c
0273    end
0274 
0275    if 0
0276       disp('---- profiling ----');
0277       fmdl1.nodes(nn,:) = fmdl1.nodes(nn,:) + rand(2,ndim)*1e-3; % perturb
0278       t=tic; FC0 = system_mat_fields(fmdl0); t0=toc(t);
0279       profile clear;
0280       profile on;
0281       t = tic; FC1a = update_system_mat_fields(fmdl0, fmdl1); t1=toc(t);
0282       profile off;
0283       profview;
0284       profsave(profile('info'),'profile');
0285       fprintf('update_system_mat_fields(fmdl0,fmdl1) = %0.3f sec [profiled]\n',t1);
0286    end
0287 
0288 function time_test_code
0289     f0= mk_common_model('e3cr',4);
0290     f0 = f0.fwd_model; f1=f0;
0291     f1.nodes(1,:)= f1.nodes(1,:) + 1e-6*[1,1,1];
0292     Fa= system_mat_fields(f1);
0293     for sp=1:10;eidors_cache clear;switch sp;
0294        case 1;
0295     tic; F1= update_system_mat_fields(f0,f1); toc
0296        case 2;
0297     tic; F0= system_mat_fields(f0); toc
0298     tic; F1= update_system_mat_fields(f0,f1); toc
0299        case 3;
0300     sv = 'skip_VOLUME';
0301     tic; fwd_model_parameters(f0,sv); toc
0302     tic; fwd_model_parameters(f0,sv); toc
0303     tic; fwd_model_parameters(f1,sv); toc
0304     tic; fwd_model_parameters(f1,sv); toc
0305     disp ----
0306     tic; F0= system_mat_fields(f0); toc
0307     tic; F0= system_mat_fields(f0); toc
0308     disp ----
0309     tic; F1= update_system_mat_fields(f0,f1); toc
0310     tic; F1= update_system_mat_fields(f0,f1); toc
0311     tic; Fb= system_mat_fields(f1); toc
0312     unit_test_cmp('',Fa,F1,1e-10);
0313        otherwise; break;
0314     end; end
0315