reciprocity_idx

 RECIPROCITY_IDX: find indices of stim, meas pairs that are recirocal
     ie. stimulation/measurement is same as measurement/stimulation on other
 usage:  idx = reciprocity_idx( fwd_model);
  fmdl: an eidors fwd_model structure
  idx:  index of corresponding reciprocal pairs
   ie measurement #3 is reciprocal to idx(3)

  if a measurement has no reciprocal pair. idx(m) = NaN;

 example 
     imdl= mk_common_model('a2c2',8);
     idx = reciprocity_idx(imdl.fwd_model);

 With an options vector
     idxr = reciprocity_idx( fmdl, 'reduce');
   calculates a 'reduced' idx which can be used to keep 
   only the unique parts of a measurement. For example
   to average pair-drive measurements
     vh_reduce = sparse(idxr,1:length(idxr),0.5)*vh

0001 function idxr = reciprocity_idx( fmdl, opt);
0002 % RECIPROCITY_IDX: find indices of stim, meas pairs that are recirocal
0003 %     ie. stimulation/measurement is same as measurement/stimulation on other
0004 % usage:  idx = reciprocity_idx( fwd_model);
0005 %  fmdl: an eidors fwd_model structure
0006 %  idx:  index of corresponding reciprocal pairs
0007 %   ie measurement #3 is reciprocal to idx(3)
0008 %
0009 %  if a measurement has no reciprocal pair. idx(m) = NaN;
0010 %
0011 % example
0012 %     imdl= mk_common_model('a2c2',8);
0013 %     idx = reciprocity_idx(imdl.fwd_model);
0014 %
0015 % With an options vector
0016 %     idxr = reciprocity_idx( fmdl, 'reduce');
0017 %   calculates a 'reduced' idx which can be used to keep
0018 %   only the unique parts of a measurement. For example
0019 %   to average pair-drive measurements
0020 %     vh_reduce = sparse(idxr,1:length(idxr),0.5)*vh
0021 
0022 % (C) 2010-2018 Andy Adler. License: GPL version 2 or version 3
0023 % $Id: reciprocity_idx.m 5823 2018-09-12 06:26:13Z aadler $
0024 
0025 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0026 
0027 sm_0= []; sm_r = [];
0028 for stim = fmdl.stimulation(:)';
0029   sp = stim.stim_pattern';
0030   nsp = sum(abs(sp))/2; %norm so max unidirectional current is 1
0031   mp = stim.meas_pattern;
0032   sp = ones(size(mp,1),1)*sp;
0033   sm_0= [sm_0; [sp/nsp, mp*nsp]];
0034   sm_r= [sm_r; [mp, sp]];
0035 end
0036 idxr = ones(size(sm_0,1),1);
0037 for i=1:size(idxr)
0038   sm0i= ones(size(sm_0,1),1)*sm_0(i,:);
0039   mm= all( abs(sm0i - sm_r) < 1e-10, 2);
0040   mp= all( abs(sm0i + sm_r) < 1e-10, 2);
0041   m = [find(mm), find(mp)];
0042   if length(m)>1;
0043       meas= sprintf('%d,',m); 
0044       error('More than one reciprocal measure %d=>(%s). Giving up', i,meas);
0045   elseif length(m)==0;
0046       idxr(i) = NaN;
0047   else
0048       idxr(i) = m;
0049   end
0050 end
0051 
0052 if nargin>=2 && strcmp(lower(opt),'reduce');
0053 % Reduce idx to the minimal vector output
0054    idxk= min([idxr,idxr(idxr)],[],2);
0055    [~,~,idxr] = unique(idxk);
0056 end
0057 
0058 function do_unit_test
0059 
0060 %     [01] [12] [23] [34] [45] [56 [67] [70]
0061 % [01] X    X    15   19   23   27  31   X
0062 % [12] X    X    X    20   24   28  32   36
0063 % [23] 1    X    X    X    25   29  33   37
0064 % [34] 2    6    X    X    X    30  34   38
0065 % [45] 3    7    11   X    X    X   35   39
0066 % [56] 4    8    12   16   X    X   X    40
0067 % [67] 5    9    13   17   21   X   X    X
0068 % [70] X    10   14   18   22   26  X    X
0069 tst.stimulation = mk_stim_patterns(8,1,[0,1],[0,1],{'rotate_meas'},1);
0070 idx = reciprocity_idx( tst ); idx = reshape(idx,5,8);
0071 unit_test_cmp('8-[01]-[01]-rotate',idx(:,[1,5]), ...
0072     [15,19,23,27,31;35,39,3,7,11]');
0073 
0074 tst.stimulation = mk_stim_patterns(8,1,[0,1],[0,1],{'rotate_meas'},10);
0075 idx = reciprocity_idx( tst ); idx = reshape(idx,5,8);
0076 unit_test_cmp('8(10)-[01]-[01]-rotate',idx(:,[1,5]), ...
0077     [15,19,23,27,31;35,39,3,7,11]');
0078 
0079 %     [01]   [12]   [23]   [34]   [45]   [50]
0080 % [01] X      X      7  9   10 11  13 13  X
0081 % [12] X      X      X      11 12  14 14  16 16
0082 % [23] 1  1   X      X      X      15 15  17 17
0083 % [34] 2  2   4  4   X      X      X      18 18
0084 % [45] 3  3   5  5   8  7   X      X      X
0085 % [50] X      6  6   9  8   12 10  X      X
0086 tst.stimulation = mk_stim_patterns(6,1,[0,1],[0,1],{'rotate_meas'},1);
0087 idx = reciprocity_idx( tst ); idx = reshape(idx,3,6); 
0088 unit_test_cmp('6-[01]-[01]-rotate',idx(:,[1,4]), [9,11,13;18 2 4]');
0089 
0090 
0091 idxr= reciprocity_idx( tst,'reduce' );
0092 unit_test_cmp('6-[01]-[01]-reduce-rotate',idxr,[1;2;3;4;5;6;7;8;1;9;2;4;3;5;7;6;8;9]);
0093 
0094 
0095 tst.stimulation = mk_stim_patterns(6,1,[0,1],[0,1],{'no_rotate_meas'},1);
0096 idx = reciprocity_idx( tst ); idx = reshape(idx,3,6); 
0097 unit_test_cmp('6-[01]-[01]-no_rotate',idx(:,[1,4]), [7,10,13;2 4 18]');
0098     
0099 %     [02]   [13]   [24]   [35]   [40]   [51]
0100 % [02] X      4   6  X      10 11  X      16 16
0101 % [13] 1  1   X      7  9   X      13 14  X
0102 % [24] X      5   4  X      11 12  X      17 17
0103 % [35] 2  2   X      8  7   X      14 15  X
0104 % [40] X      6   5  X      12 10  X      18 18
0105 % [51] 3  3   X      9  8   X      15 13  X
0106 tst.stimulation = mk_stim_patterns(6,1,[0,2],[0,2],{'rotate_meas'},1);
0107 idx = reciprocity_idx( tst ); idx = reshape(idx,3,6);
0108 unit_test_cmp('6-[02]-[02]-rotate',idx(:,[1,4]), [6,11,16;15,2,7]');
0109 
0110 tst.stimulation = mk_stim_patterns(6,1,[0,2],[0,2],{'no_rotate_meas'},1);
0111 idx = reciprocity_idx( tst ); idx = reshape(idx,3,6); 
0112 unit_test_cmp('6-[02]-[02]-no_rotate',idx(:,[1,4]), [4,10,16;2,8,14]');
0113 
0114 idxr= reciprocity_idx( tst,'reduce' );
0115 unit_test_cmp('6-[02]-[02]-reduce-no_rotate',idxr,[1;2;3;1;4;5;4;6;7;2;6;8;5;8;9;3;7;9]);
0116 
0117 
0118 %     [02]   [13]   [24]   [35]   [40]   [51]
0119 % [02] 1  1   7  12  13 17  19 22  25 27  31 32
0120 % [13] 2  2   8   7  14 18  20 23  26 28  32 33
0121 % [24] 3  3   9   8  15 13  21 24  27 29  33 34
0122 % [35] 4  4   10  9  16 14  22 19  28 30  34 35
0123 % [40] 5  5   11 10  17 15  23 20  29 25  35 36
0124 % [51] 6  6   12 11  18 16  24 21  30 26  36 31
0125 tst.stimulation = mk_stim_patterns(6,1,[0,2],[0,2],{'meas_current','rotate_meas'},1);
0126 idx = reciprocity_idx( tst ); idx = reshape(idx,6,6);
0127 unit_test_cmp('6-[02]-[02]-mc-rotate',idx(:,[1,4]),  ...
0128             [ 1,12,17,22,27,32;19,30,35,4,9,14]');
0129 
0130 tst.stimulation = mk_stim_patterns(6,1,[0,2],[0,2],{'meas_current','no_rotate_meas'},1);
0131 idx = reciprocity_idx( tst ); idx = reshape(idx,6,6); 
0132 unit_test_cmp('6-[02]-[02]-mc-no_rotate',idx(:,[1,4]),  ...
0133             [ 1,7,13,19,25,31; 4,10,16,22,28,34]');
0134 
0135 tst.stimulation = mk_stim_patterns(6,1,[0,4],[0,4],{'meas_current','no_rotate_meas'},1);
0136 idx = reciprocity_idx( tst ); idx = reshape(idx,6,6); 
0137 unit_test_cmp('6-[04]-[04]-mc-no_rotate',idx(:,[1,4]),  ...
0138             [ 1,7,13,19,25,31; 4,10,16,22,28,34]');