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);

0001 function idxr = reciprocity_idx( fmdl);
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 % (C) 2010 Andy Adler. License: GPL version 2 or version 3
0016 % $Id: reciprocity_idx.html 2819 2011-09-07 16:43:11Z aadler $
0017 
0018 if isstr(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0019 
0020 sm_0= []; sm_r = [];
0021 for stim = fmdl.stimulation(:)';
0022   sp = stim.stim_pattern';
0023   nsp = sum(abs(sp))/2; %norm so max unidirectional current is 1
0024   mp = stim.meas_pattern;
0025   sp = ones(size(mp,1),1)*sp;
0026   sm_0= [sm_0; [sp/nsp, mp*nsp]];
0027   sm_r= [sm_r; [mp, sp]];
0028 end
0029 idxr = ones(size(sm_0,1),1);
0030 for i=1:size(idxr)
0031   sm0i= ones(size(sm_0,1),1)*sm_0(i,:);
0032   mm= all( abs(sm0i - sm_r) < 1e-10, 2);
0033   mp= all( abs(sm0i + sm_r) < 1e-10, 2);
0034   m = [find(mm), find(mp)];
0035   if length(m)>1;
0036       meas= sprintf('%d,',m); 
0037       error('More than one reciprocal measure %d=>(%s). Giving up', i,meas);
0038   elseif length(m)==0;
0039       idxr(i) = NaN;
0040   else
0041       idxr(i) = m;
0042   end
0043 end
0044 
0045 function do_unit_test
0046 
0047 %     [01] [12] [23] [34] [45] [56 [67] [70]
0048 % [01] X    X    15   19   23   27  31   X
0049 % [12] X    X    X    20   24   28  32   36
0050 % [23] 1    X    X    X    25   29  33   37
0051 % [34] 2    6    X    X    X    30  34   38
0052 % [45] 3    7    11   X    X    X   35   39
0053 % [56] 4    8    12   16   X    X   X    40
0054 % [67] 5    9    13   17   21   X   X    X
0055 % [70] X    10   14   18   22   26  X    X
0056 tst.stimulation = mk_stim_patterns(8,1,[0,1],[0,1],{'rotate_meas'},1);
0057 idx = reciprocity_idx( tst ); idx = reshape(idx,5,8);
0058 unit_test_cmp('8-[01]-[01]-rotate',idx(:,[1,5]), ...
0059     [15,19,23,27,31;35,39,3,7,11]');
0060 
0061 tst.stimulation = mk_stim_patterns(8,1,[0,1],[0,1],{'rotate_meas'},10);
0062 idx = reciprocity_idx( tst ); idx = reshape(idx,5,8);
0063 unit_test_cmp('8(10)-[01]-[01]-rotate',idx(:,[1,5]), ...
0064     [15,19,23,27,31;35,39,3,7,11]');
0065 
0066 %     [01]   [12]   [23]   [34]   [45]   [50]
0067 % [01] X      X      7  9   10 11  13 13  X
0068 % [12] X      X      X      11 12  14 14  16 16
0069 % [23] 1  1   X      X      X      15 15  17 17
0070 % [34] 2  2   4  4   X      X      X      18 18
0071 % [45] 3  3   5  5   8  7   X      X      X
0072 % [50] X      6  6   9  8   12 10  X      X
0073 tst.stimulation = mk_stim_patterns(6,1,[0,1],[0,1],{'rotate_meas'},1);
0074 idx = reciprocity_idx( tst ); idx = reshape(idx,3,6); 
0075 unit_test_cmp('6-[01]-[01]-rotate',idx(:,[1,4]), [9,11,13;18 2 4]');
0076 
0077 tst.stimulation = mk_stim_patterns(6,1,[0,1],[0,1],{'no_rotate_meas'},1);
0078 idx = reciprocity_idx( tst ); idx = reshape(idx,3,6); 
0079 unit_test_cmp('6-[01]-[01]-no_rotate',idx(:,[1,4]), [7,10,13;2 4 18]');
0080     
0081 %     [02]   [13]   [24]   [35]   [40]   [51]
0082 % [02] X      4   6  X      10 11  X      16 16
0083 % [13] 1  1   X      7  9   X      13 14  X
0084 % [24] X      5   4  X      11 12  X      17 17
0085 % [35] 2  2   X      8  7   X      14 15  X
0086 % [40] X      6   5  X      12 10  X      18 18
0087 % [51] 3  3   X      9  8   X      15 13  X
0088 tst.stimulation = mk_stim_patterns(6,1,[0,2],[0,2],{'rotate_meas'},1);
0089 idx = reciprocity_idx( tst ); idx = reshape(idx,3,6);
0090 unit_test_cmp('6-[02]-[02]-rotate',idx(:,[1,4]), [6,11,16;15,2,7]');
0091 
0092 tst.stimulation = mk_stim_patterns(6,1,[0,2],[0,2],{'no_rotate_meas'},1);
0093 idx = reciprocity_idx( tst ); idx = reshape(idx,3,6); 
0094 unit_test_cmp('6-[02]-[02]-no_rotate',idx(:,[1,4]), [4,10,16;2,8,14]');
0095 
0096 %     [02]   [13]   [24]   [35]   [40]   [51]
0097 % [02] 1  1   7  12  13 17  19 22  25 27  31 32
0098 % [13] 2  2   8   7  14 18  20 23  26 28  32 33
0099 % [24] 3  3   9   8  15 13  21 24  27 29  33 34
0100 % [35] 4  4   10  9  16 14  22 19  28 30  34 35
0101 % [40] 5  5   11 10  17 15  23 20  29 25  35 36
0102 % [51] 6  6   12 11  18 16  24 21  30 26  36 31
0103 tst.stimulation = mk_stim_patterns(6,1,[0,2],[0,2],{'meas_current','rotate_meas'},1);
0104 idx = reciprocity_idx( tst ); idx = reshape(idx,6,6);
0105 unit_test_cmp('6-[02]-[02]-mc-rotate',idx(:,[1,4]),  ...
0106             [ 1,12,17,22,27,32;19,30,35,4,9,14]');
0107 
0108 tst.stimulation = mk_stim_patterns(6,1,[0,2],[0,2],{'meas_current','no_rotate_meas'},1);
0109 idx = reciprocity_idx( tst ); idx = reshape(idx,6,6); 
0110 unit_test_cmp('6-[02]-[02]-mc-no_rotate',idx(:,[1,4]),  ...
0111             [ 1,7,13,19,25,31; 4,10,16,22,28,34]');
0112 
0113 tst.stimulation = mk_stim_patterns(6,1,[0,4],[0,4],{'meas_current','no_rotate_meas'},1);
0114 idx = reciprocity_idx( tst ); idx = reshape(idx,6,6); 
0115 unit_test_cmp('6-[04]-[04]-mc-no_rotate',idx(:,[1,4]),  ...
0116             [ 1,7,13,19,25,31; 4,10,16,22,28,34]');