mk_stim_patterns

MK_STIM_PATTERNS: create an EIDORS stimulation pattern structure
                to form part of a fwd_model object
 [stim, meas_sel] = mk_stim_patterns( n_elec, n_rings, ...
                                      inj, meas, options, amplitude)

 where
 stim(#).stimulation = 'Amp'
     (#).stim_pattern= [vector n_elec*n_rings x 1 ]
     (#).meas_pattern= [matrix n_elec*n_rings x n_meas_patterns]

 for example, for an adjacent pattern for 4 electrodes, with 0.5 Amp
   if all electrodes are used for measurement
 stim(1).stim_pattern= [0.5;-0.5;0;0]
 stim(1).meas_pattern= [1,-1, 0, 0 
                        0, 1,-1, 0 
                        0, 0, 1,-1 
                       -1, 0, 0, 1]

 meas_sel: when not using data from current injection electrodes,
           it is common to be given a full measurement set.
           For example 16 electrodes gives 208 measures, but 256
           measure sets are common. 'meas_sel' indicates which
           electrodes are used

 PARAMETERS:
   n_elec:   number of electrodes per ring
   n_rings:  number of electrode rings (1 for 2D)

   inj: injection pattern
      '{ad}'        -> adjacent drive: equivalent to [0 1]
      '{op}'        -> opposite drive: equivalent to [0, n_elec/2]
      '{trig}'      -> trigonometric drive [cos,sin,cos,sin ...]
                       '{trig}' implies the 'meas_current' option.
      '{trigcscs}'  -> trigonometric drive [cos,sin,cos,sin ...]
      '{trigccss}'  -> trigonometric drive [cos,cos, ... sin,sin, ...]
      '{mono}'      -> Drive via each elec, current leaves by ground
      Bi-polar injection patterns:
        [x y]: First pattern is [x,y] next is [x+1,y+1] 
      Mono-polar injection patterns:
        [x]:   First pattern is [x]   next is [x+1] 

   meas: measurement pattern
      '{ad}'        -> adjacent measurement
      '{op}'        -> opposite drive: equivalent to [0, n_elec/2]
      '{trig}'      -> trigonometric drive [sin,cos,sin,cos ...]
      '{mono}'      -> Meas at each elec (Reminder, you may want to set meas_current);
      Bi-polar measurement patterns:
        [x y]: First pattern is [x,y] next is [x+1,y+1] 
      Mono-polar measurement patterns:
        [x]:   First pattern is [x]   next is [x+1] 

   options: cell array of options, eg {'no_meas_current'}
     if contradictory options are specified, only the last applies
      'no_meas_current' / 'meas_current'
         -> do / don't make mesurements on current carrying electrodes 
      'no_meas_current_next0' is same as 'no_meas_current'
         -> to not make measurements on the electrodes next to the current carrying,
            use 'no_meas_current_next1'
         -> to not make measurements on the two electrodes next to the current carrying,
            use 'no_meas_current_next2'
         -> to not make measurements on the three electrodes next to the current carrying,
            use 'no_meas_current_next3'
      'rotate_meas' / 'no_rotate_meas'
         -> do / don't rotate measurements with stimulation pattern
      'do_redundant' / 'no_redundant'
         -> do / don't make reciprocally redundant measures
      'balance_inj' / 'no_balance_inj'
         -> do / don't draw current from all electrodes so total
            injection is zero (useful for mono patterns)
      'balance_meas' / 'no_balance_meas'
         -> do / don't subtrant measurement from all electrodes so total
            average measurement is zero (useful for mono patterns)

   amplitude: drive current levels, DEFAULT = 0.010 Amp

0001 function [stim, meas_sel]= mk_stim_patterns( ...
0002             n_elec, n_rings, inj, meas, options, amplitude)
0003 %MK_STIM_PATTERNS: create an EIDORS stimulation pattern structure
0004 %                to form part of a fwd_model object
0005 % [stim, meas_sel] = mk_stim_patterns( n_elec, n_rings, ...
0006 %                                      inj, meas, options, amplitude)
0007 %
0008 % where
0009 % stim(#).stimulation = 'Amp'
0010 %     (#).stim_pattern= [vector n_elec*n_rings x 1 ]
0011 %     (#).meas_pattern= [matrix n_elec*n_rings x n_meas_patterns]
0012 %
0013 % for example, for an adjacent pattern for 4 electrodes, with 0.5 Amp
0014 %   if all electrodes are used for measurement
0015 % stim(1).stim_pattern= [0.5;-0.5;0;0]
0016 % stim(1).meas_pattern= [1,-1, 0, 0
0017 %                        0, 1,-1, 0
0018 %                        0, 0, 1,-1
0019 %                       -1, 0, 0, 1]
0020 %
0021 % meas_sel: when not using data from current injection electrodes,
0022 %           it is common to be given a full measurement set.
0023 %           For example 16 electrodes gives 208 measures, but 256
0024 %           measure sets are common. 'meas_sel' indicates which
0025 %           electrodes are used
0026 %
0027 % PARAMETERS:
0028 %   n_elec:   number of electrodes per ring
0029 %   n_rings:  number of electrode rings (1 for 2D)
0030 %
0031 %   inj: injection pattern
0032 %      '{ad}'        -> adjacent drive: equivalent to [0 1]
0033 %      '{op}'        -> opposite drive: equivalent to [0, n_elec/2]
0034 %      '{trig}'      -> trigonometric drive [cos,sin,cos,sin ...]
0035 %                       '{trig}' implies the 'meas_current' option.
0036 %      '{trigcscs}'  -> trigonometric drive [cos,sin,cos,sin ...]
0037 %      '{trigccss}'  -> trigonometric drive [cos,cos, ... sin,sin, ...]
0038 %      '{mono}'      -> Drive via each elec, current leaves by ground
0039 %      Bi-polar injection patterns:
0040 %        [x y]: First pattern is [x,y] next is [x+1,y+1]
0041 %      Mono-polar injection patterns:
0042 %        [x]:   First pattern is [x]   next is [x+1]
0043 %
0044 %   meas: measurement pattern
0045 %      '{ad}'        -> adjacent measurement
0046 %      '{op}'        -> opposite drive: equivalent to [0, n_elec/2]
0047 %      '{trig}'      -> trigonometric drive [sin,cos,sin,cos ...]
0048 %      '{mono}'      -> Meas at each elec (Reminder, you may want to set meas_current);
0049 %      Bi-polar measurement patterns:
0050 %        [x y]: First pattern is [x,y] next is [x+1,y+1]
0051 %      Mono-polar measurement patterns:
0052 %        [x]:   First pattern is [x]   next is [x+1]
0053 %
0054 %   options: cell array of options, eg {'no_meas_current'}
0055 %     if contradictory options are specified, only the last applies
0056 %      'no_meas_current' / 'meas_current'
0057 %         -> do / don't make mesurements on current carrying electrodes
0058 %      'no_meas_current_next0' is same as 'no_meas_current'
0059 %         -> to not make measurements on the electrodes next to the current carrying,
0060 %            use 'no_meas_current_next1'
0061 %         -> to not make measurements on the two electrodes next to the current carrying,
0062 %            use 'no_meas_current_next2'
0063 %         -> to not make measurements on the three electrodes next to the current carrying,
0064 %            use 'no_meas_current_next3'
0065 %      'rotate_meas' / 'no_rotate_meas'
0066 %         -> do / don't rotate measurements with stimulation pattern
0067 %      'do_redundant' / 'no_redundant'
0068 %         -> do / don't make reciprocally redundant measures
0069 %      'balance_inj' / 'no_balance_inj'
0070 %         -> do / don't draw current from all electrodes so total
0071 %            injection is zero (useful for mono patterns)
0072 %      'balance_meas' / 'no_balance_meas'
0073 %         -> do / don't subtrant measurement from all electrodes so total
0074 %            average measurement is zero (useful for mono patterns)
0075 %
0076 %   amplitude: drive current levels, DEFAULT = 0.010 Amp
0077 
0078 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0079 % $Id: mk_stim_patterns.m 6636 2024-01-11 15:22:43Z aadler $
0080 
0081 if ischar(n_elec) && strcmp(n_elec,'UNIT_TEST'); do_unit_test; return; end
0082 
0083 if nargin<6; amplitude= .01; end
0084 if nargin<5; options= {};  end
0085 v = process_args(n_elec, n_rings, inj, meas, options, amplitude );
0086 
0087 [stim,mpat ] = calc_stim(v, n_elec, n_rings);
0088 
0089 % Meas_sel selects meas which would exist if we measured all of them
0090 v.do_redundant = 1; v.use_meas_current = 1; 
0091 [jnk ,mpat0] = calc_stim(v, n_elec, n_rings);
0092 
0093 %meas_sel= meas_select_old( n_elec, v.inj, v);
0094  meas_sel= meas_select( mpat, mpat0);
0095 
0096 function [stim,mpat] = calc_stim(v, n_elec, n_rings)
0097    curr_pat = v.inj(:) * ones(1,n_elec);
0098    meas_pat = v.meas(:) * ones(1,n_elec);
0099    offset   = [1;1]*(0:n_elec-1);
0100 
0101    stim= struct([]);
0102    mpat= struct([]);
0103    i=1; j=1;
0104    for ring = 0:v.n_rings-1
0105       seen_patterns= struct;
0106       for elec= 0:v.n_elec-1
0107           if v.trig_inj && elec == v.n_elec-1 ; continue; end % No indep patterns
0108           s_pat= mk_stim_pat(v, elec, ring );
0109           m_pat= mk_meas_pat(v, elec, ring );
0110 
0111           if v.do_redundant == 0 % elim redudant
0112              [m_pat, seen_patterns] = elim_redundant(m_pat, s_pat, seen_patterns);
0113           end
0114 
0115           if ~isempty(m_pat) 
0116               stim(i).stimulation = 'Amp';
0117               stim(i).stim_pattern= sparse(s_pat);
0118               stim(i).meas_pattern= sparse(m_pat);
0119               i=i+1;
0120           end
0121               mpat(j).meas_pattern= sparse(m_pat);
0122               j=j+1;
0123       end
0124    end
0125 
0126 % when not using data from current injection electrodes,
0127 % it is common to be given a full measurement set.
0128 % For example 16 electrodes gives 208 measures, but 256
0129 % measure sets are common.
0130 % This function calculates a selector matrix to remove the extra
0131 %
0132 % reshape(meas_select( 6, [0,1], v),6,6) 0 0 1 1 1 0
0133 %                                        0 0 0 1 1 1
0134 %                                        1 0 0 0 1 1
0135 %                                        1 1 0 0 0 1
0136 %                                        1 1 1 0 0 0
0137 %                                        0 1 1 1 0 0
0138 %
0139 % reshape(meas_select( 6, [0,2], v),6,6) 0 0 1 1 0 0
0140 %                                        0 0 0 1 1 0
0141 %                                        0 0 0 0 1 1
0142 %                                        1 0 0 0 0 1
0143 %                                        1 1 0 0 0 0
0144 %                                        0 1 1 0 0 0
0145 %
0146 % reshape(meas_select( 6, [0,3], v),6,6) 0 0 1 0 0 1
0147 %                                        1 0 0 1 0 0
0148 %                                        0 1 0 0 1 0
0149 %                                        0 0 1 0 0 1
0150 %                                        1 0 0 1 0 0
0151 %                                        0 1 0 0 1 0
0152 function meas_sel= meas_select( mpat, mpat0);
0153    meas_sel = [];
0154    for i=1:length(mpat);
0155       [mset,  err ] = mk_meas_set( mpat(i).meas_pattern );
0156       [mset0, err0] = mk_meas_set( mpat0(i).meas_pattern );
0157       if err || err0; msel_i = 1 + 0*mset; % Set to 1 for error
0158       else            msel_i = ismember(mset0,mset);
0159       end
0160       meas_sel = [meas_sel; msel_i];
0161    end
0162 
0163    meas_sel = logical(meas_sel(:));
0164 
0165 function [mset,err] = mk_meas_set( meas_pat )
0166    mpats= size(~meas_pat,1);
0167    mset = zeros(mpats,1);
0168    err=0;
0169    for i=1:size(meas_pat,1);
0170       mpat = meas_pat(i,:);
0171       fp = find(mpat>0); if length(fp)==0; fp = 0; end
0172       fn = find(mpat<0); if length(fn)==0; fn = 0; end
0173       if length(fp)>1 || length(fn)>1 ; err=1; return; end
0174       mset(i) =  fp*1e7 + fn;
0175    end
0176 
0177 function stim_pat = mk_stim_pat(v, elec, ring );
0178    stim_pat = sparse(v.tn_elec, 1);
0179    if v.balance_inj
0180       stim_pat= stim_pat - sum(v.i_factor)/ (v.tn_elec-1);
0181    elseif v.trig_inj
0182       stim_pat = trig_pat( elec, v.tn_elec, v.trig_inj) * v.amplitude;
0183       return;
0184    end
0185 
0186    stim_idx = rem( v.inj + elec, v.n_elec) + 1 + v.n_elec*ring;
0187    stim_pat( stim_idx ) = v.amplitude*v.i_factor;
0188 
0189 % Measurement config can stay static, or can rotate with
0190 % the stim pattern. This code keeps measurements static
0191 function meas = mk_meas_pat(v, elec, ring );
0192    meas= sparse(v.tn_elec, v.tn_elec);
0193    if v.balance_meas
0194       meas= meas - sum(v.m_factor)/ (v.tn_elec-1);
0195    elseif v.trig_meas
0196       meas= trig_pat( 0:v.tn_elec-2, v.tn_elec)';
0197       return;
0198    end
0199 
0200    if v.rotate_meas
0201       ofs = elec;
0202    else
0203       ofs = 0;
0204    end 
0205 
0206    mseq= 0:v.tn_elec-1;
0207    within_ring = rem(v.n_elec +  mseq , v.n_elec);
0208    ouside_ring = floor( mseq / v.n_elec) * v.n_elec;
0209    meas_seq    = mseq *v.tn_elec + 1;
0210 
0211    for i=1:length(v.meas)
0212       meas_pat = rem( v.meas(i) + within_ring + ofs, v.n_elec ) + ...
0213                        ouside_ring + meas_seq;
0214       meas(meas_pat) = v.m_factor(i);
0215    end
0216 
0217    if v.use_meas_current == 0
0218    % each column of meas is a measurement pattern
0219    % Test whether each col has contribution from stim
0220        stim_idx = rem( v.inj + elec, v.n_elec) + 1 + v.n_elec*ring;
0221 
0222        if v.use_meas_current_next
0223    % If we want to eliminate the ones next to the stimulation electrodes
0224    %  ie. Swisstom device, then use this
0225      
0226           for ni= -v.use_meas_current_next:v.use_meas_current_next;
0227             stim_idx = [stim_idx, ...
0228                         rem( v.inj + elec + ni, v.n_elec) + 1 + v.n_elec*ring];
0229           end
0230           stim_idx = unique(stim_idx);
0231           stim_idx(stim_idx<=0) = stim_idx(stim_idx<=0) + v.n_elec;
0232        end
0233        elim= any(meas(stim_idx,:),1);
0234        meas(:,elim) = [];
0235    end
0236 
0237    meas= meas';
0238 
0239 
0240 function v = process_args(n_elec, n_rings, inj, meas, options, amplitude )
0241 
0242 % SET DEFAULTS
0243    v.trig_meas= 0;
0244    v.trig_inj = 0;
0245 
0246    v.use_meas_current = 0;
0247    v.use_meas_current_next = 0;
0248    v.rotate_meas = 0;
0249    v.do_redundant = 1;
0250    v.balance_inj = 0;
0251    v.balance_meas= 0;
0252 
0253 % Stimulation (injection) pattern
0254 % This currently does not handle complicated injection patterns
0255 if ischar(inj)
0256    if      strcmp(inj,'{ad}')
0257       inj= [0, 1];
0258       rel_ampl= [-1;1];
0259    elseif  strcmp(inj,'{op}')
0260       inj= [0, floor(n_elec/2)];
0261       rel_ampl= [-1;1];
0262    elseif  strcmp(inj,'{trig}')
0263       v.trig_inj = 1;
0264       v.use_meas_current = 1; % We need to measure on the electrodes
0265       rel_ampl= [];
0266    elseif  strcmp(inj,'{trigcscs}')
0267       v.trig_inj = 1;
0268       v.use_meas_current = 1; % We need to measure on the electrodes
0269       rel_ampl= [];
0270    elseif  strcmp(inj,'{trigccss}')
0271       v.trig_inj = 2;
0272       v.use_meas_current = 1; % We need to measure on the electrodes
0273       rel_ampl= [];
0274    elseif  strcmp(inj,'{mono}')
0275       inj= [0];
0276       rel_ampl= [1];
0277    else
0278       error(['parameter inj=',inj,' not understood']);
0279    end
0280 elseif prod(size(inj))==1
0281       rel_ampl= [1];
0282 elseif prod(size(inj))==2
0283       rel_ampl= [-1;1];
0284 else
0285       error(['parameter inj not understood']);
0286 end
0287 
0288 v.inj= inj;
0289 v.i_factor=      rel_ampl;
0290 
0291 % Measurement configuration.
0292 % All systems I know of use adjacent measurement,
0293 % are there actually any others?
0294 if ischar(meas)
0295    if      strcmp(meas,'{ad}')
0296       meas=     [0, 1];
0297       rel_ampl= [ 1;-1];
0298    elseif  strcmp(meas,'{op}')
0299       meas= [0, floor(n_elec/2)];
0300       rel_ampl= [ 1;-1];
0301    elseif  strcmp(meas,'{trig}')
0302       v.trig_meas= 1;
0303       rel_ampl= [ 1;-1];
0304    elseif  strcmp(meas,'{mono}')
0305       meas= [0];
0306       rel_ampl= [1];
0307    else
0308       error(['parameter meas=',meas,' not understood']);
0309    end
0310 elseif prod(size(meas))==1
0311       rel_ampl= [1];
0312 elseif prod(size(meas))==2
0313       rel_ampl= [ 1;-1];
0314 else
0315       error(['parameter meas not understood']);
0316 end
0317 
0318 v.meas=          meas;
0319 v.m_factor=      rel_ampl;
0320 
0321 v.n_elec = n_elec;
0322 v.n_rings= n_rings;
0323 v.tn_elec= n_rings * n_elec;
0324 v.amplitude = amplitude;
0325 
0326 v= parse_options(v, options);
0327 
0328 function v= parse_options(v, options);
0329 % iterate through the options cell array
0330 for opt = options
0331    switch opt{1};
0332    case {'no_meas_current', ...
0333          'no_meas_current_next0'};
0334       v.use_meas_current = 0;
0335       v.use_meas_current_next = 0;
0336    case 'no_meas_current_next1';
0337       v.use_meas_current = 0;
0338       v.use_meas_current_next = 1;
0339    case 'no_meas_current_next2';
0340       v.use_meas_current = 0;
0341       v.use_meas_current_next = 2;
0342    case 'no_meas_current_next3';
0343       v.use_meas_current = 0;
0344       v.use_meas_current_next = 3;
0345    case 'meas_current';
0346       v.use_meas_current = 1;
0347    case 'rotate_meas';
0348       v.rotate_meas = 1;
0349    case 'no_rotate_meas';
0350       v.rotate_meas = 0;
0351    case 'do_redundant';
0352       v.do_redundant = 1;
0353    case 'no_redundant';
0354       v.do_redundant = 0;
0355    case 'balance_inj';
0356       v.balance_inj = 1;
0357    case 'no_balance_inj';
0358       v.balance_inj = 0;
0359    case 'balance_meas';
0360       v.balance_meas= 1;
0361    case 'no_balance_meas';
0362       v.balance_meas= 0;
0363    otherwise
0364       error(['option parameter opt=',opt,' not understood']);
0365    end
0366 end
0367 
0368 function [m_pat, seen_patterns] = elim_redundant(m_pat, s_pat, seen_patterns);
0369    m_pat_new= sparse([]);
0370    s_pat_str= ['s',sprintf('%d_', find(s_pat) ),'m'];
0371    for j=1:size(m_pat,1);
0372       this_m_pat= m_pat(j,:);
0373       pat_str= [s_pat_str, sprintf('%d_', find(this_m_pat))];
0374       if ~isfield(seen_patterns,pat_str);
0375          m_pat_new= [m_pat_new;this_m_pat];
0376          % we've seen this pattern
0377          seen_patterns.(pat_str)= 1;
0378          % and it's dual by reciprocity
0379          pat_str= ['s',sprintf('%d_', find(this_m_pat) ), ...
0380                    'm',sprintf('%d_', find(s_pat))];
0381          seen_patterns.(pat_str)= 1;
0382       end
0383    end
0384    m_pat= m_pat_new;
0385 
0386 % create trig patterns.
0387 %
0388 % n_elecs is total number of electrodes
0389 % elec    is the electrodes selected (can be multiple)
0390 %         (elec goes from 0 to n_elecs-1
0391 % if sel = 1 -> cos|sin|cos|sin (default)
0392 % if sel = 2 -> cos|cos|sin|sin
0393 %
0394 function pat= trig_pat( elec_sel, n_elecs, sel);
0395     if nargin<3; sel=1; end
0396     idx= linspace(0,2*pi,n_elecs+1)'; idx(end)= [];
0397     omega= idx*[1:n_elecs/2];
0398     meas_pat= [cos(omega), sin(omega) ];
0399     if sel==1;
0400        % reorder so we get cos|sin|cos|sin
0401        order = reshape(1:n_elecs,[],2)';
0402        meas_pat= meas_pat(:,order(:));
0403     end
0404     meas_pat= meas_pat(:,1:end-1); % only n_elecs-1 independent patterns
0405     pat  = meas_pat(:, elec_sel+1);
0406 
0407 function trig_tests;
0408    stim = mk_stim_patterns(8,1,'{trig}',[0,1],{},2);
0409    t= linspace(0,2*pi,8+1)'; t(end)= [];
0410    unit_test_cmp('trig: t1',[stim(1:4).stim_pattern], ...
0411          2*[cos(t),sin(t),cos(2*t),sin(2*t)], 1e-10);
0412 
0413    stim = mk_stim_patterns(8,1,'{trigcscs}',[0,1],{},2);
0414    unit_test_cmp('trig: t2',[stim(1:4).stim_pattern], ...
0415          2*[cos(t),sin(t),cos(2*t),sin(2*t)], 1e-10);
0416 
0417    stim = mk_stim_patterns(8,1,'{trigccss}',[0,1],{},2);
0418    test = 2*[cos(t),cos(2*t),cos(3*t), cos(4*t), sin(t),sin(2*t),sin(3*t)];
0419    unit_test_cmp('trig: t2',[stim(:).stim_pattern], test, 1e-10);
0420 
0421    stim = mk_stim_patterns(8,1,'{trigccss}','{trig}',{},2);
0422    mp1 = stim(1).meas_pattern;
0423    mp2 = stim(2).meas_pattern;
0424    unit_test_cmp('trig: m1',mp1,mp2,1e-15);
0425    test = [cos(t),sin(t),cos(2*t),sin(2*t), ...
0426        cos(3*t),sin(3*t),cos(4*t)];
0427    unit_test_cmp('trig: m2',mp1',test,1e-15);
0428 
0429 function do_unit_test
0430    trig_tests;
0431    stim = mk_stim_patterns(4,1,[0,1],[0,1],{},1);
0432    unit_test_cmp('t1',stim(1).stim_pattern, [-1;1;0;0]);
0433    unit_test_cmp('t2',stim(4).stim_pattern, [1;0;0;-1]);
0434    unit_test_cmp('t3',stim(1).meas_pattern, [0,0,1,-1]);
0435    unit_test_cmp('t4',stim(4).meas_pattern, [0,1,-1,0]);
0436 
0437 %      'no_meas_current' / 'meas_current'
0438 %         -> do / don't make mesurements on current carrying electrodes
0439    stim = mk_stim_patterns(4,1,[0,1],[0,1],{'meas_current'},1);
0440    unit_test_cmp('meas_current: t1',stim(1).meas_pattern,  ...
0441          [1,-1,0,0; 0,1,-1,0; 0,0,1,-1; -1,0,0,1]);
0442    unit_test_cmp('meas_current: t2',stim(4).stim_pattern, [1;0;0;-1]);
0443    stim = mk_stim_patterns(4,1,[0,1],[0,1],{'no_meas_current'},1);
0444    unit_test_cmp('meas_current: t3',stim(1).meas_pattern,  [0,0,1,-1]);
0445    unit_test_cmp('meas_current: t2',stim(4).stim_pattern, [1;0;0;-1]);
0446 
0447 %      'rotate_meas' / 'no_rotate_meas'
0448 %         -> do / don't rotate measurements with stimulation pattern
0449 
0450    stim = mk_stim_patterns(6,1,[0,1],[0,1],{'no_rotate_meas'},1);
0451    unit_test_cmp('no_rotate_meas: t1',stim(2).stim_pattern, [0;-1;1;0;0;0]);
0452    unit_test_cmp('no_rotate_meas: t2',stim(2).meas_pattern, ...
0453          [0,0,0,1,-1,0; 0,0,0,0,1,-1; -1,0,0,0,0,1]);
0454    unit_test_cmp('no_rotate_meas: t3',stim(3).stim_pattern, [0;0;-1;1;0;0]);
0455    unit_test_cmp('no_rotate_meas: t4',stim(3).meas_pattern, ...
0456          [1,-1,0,0,0,0;0,0,0,0,1,-1; -1,0,0,0,0,1]);
0457 
0458    stim = mk_stim_patterns(6,1,[0,1],[0,1],{'rotate_meas'},1);
0459    unit_test_cmp('rotate_meas: t1',stim(2).stim_pattern, [0;-1;1;0;0;0]);
0460    unit_test_cmp('rotate_meas: t2',stim(2).meas_pattern, ...
0461          [0,0,0,1,-1,0; 0,0,0,0,1,-1; -1,0,0,0,0,1]);
0462    unit_test_cmp('rotate_meas: t3',stim(3).stim_pattern, [0;0;-1;1;0;0]);
0463    unit_test_cmp('rotate_meas: t4',stim(3).meas_pattern, ...
0464          [0,0,0,0,1,-1; -1,0,0,0,0,1; 1,-1,0,0,0,0]);
0465 
0466 %      'do_redundant' / 'no_redundant'
0467 %         -> do / don't make reciprocally redundant measures
0468    stim = mk_stim_patterns(6,1,[0,1],[0,1],{'do_redundant'},1);
0469 
0470    unit_test_cmp('do_redundant: t0',length(stim), 6);
0471    unit_test_cmp('do_redundant: t1',stim(2).stim_pattern, [0;-1;1;0;0;0]);
0472    unit_test_cmp('do_redundant: t2',stim(2).meas_pattern, ...
0473          [0,0,0,1,-1,0; 0,0,0,0,1,-1; -1,0,0,0,0,1]);
0474    unit_test_cmp('do_redundant: t3',stim(3).stim_pattern, [0;0;-1;1;0;0]);
0475    unit_test_cmp('do_redundant: t4',stim(3).meas_pattern, ...
0476          [1,-1,0,0,0,0;0,0,0,0,1,-1; -1,0,0,0,0,1]);
0477 
0478    stim = mk_stim_patterns(6,1,[0,1],[0,1],{'no_redundant'},1);
0479    unit_test_cmp('no_redundant: t0',length(stim), 4);
0480    unit_test_cmp('no_redundant: t1',stim(2).stim_pattern, [0;-1;1;0;0;0]);
0481    unit_test_cmp('no_redundant: t2',stim(2).meas_pattern, ...
0482          [0,0,0,1,-1,0; 0,0,0,0,1,-1; -1,0,0,0,0,1]);
0483    unit_test_cmp('no_redundant: t3',stim(3).stim_pattern, [0;0;-1;1;0;0]);
0484    unit_test_cmp('no_redundant: t4',stim(3).meas_pattern, ...
0485          [0,0,0,0,1,-1;-1,0,0,0,0,1]);
0486    unit_test_cmp('no_redundant: t5',stim(4).meas_pattern, ...
0487          [-1,0,0,0,0,1]);
0488 
0489 %      'balance_inj' / 'no_balance_inj'
0490 %         -> do / don't draw current from all electrodes so total
0491 %            injection is zero (useful for mono patterns)
0492    stim = mk_stim_patterns(4,1,'{mono}',[0,1],{'balance_inj','meas_current'},1);
0493    unit_test_cmp('balance_inj: t0',length(stim), 4,1);
0494    unit_test_cmp('balance_inj: t1',stim(2).stim_pattern, -[1;-3;1;1]/3);
0495    unit_test_cmp('balance_inj: t2',stim(2).meas_pattern, ...
0496          [1,-1,0,0;0,1,-1,0;0,0,1,-1;-1,0,0,1]);
0497 
0498    stim = mk_stim_patterns(4,1,'{mono}',[0,1],{'no_balance_inj','no_meas_current'},1);
0499    unit_test_cmp('no_balance_inj: t0',length(stim), 4,1);
0500    unit_test_cmp('no_balance_inj: t1',stim(2).stim_pattern, [0;1;0;0]);
0501    unit_test_cmp('no_balance_inj: t2',stim(2).meas_pattern, ...
0502          [0,0,1,-1;-1,0,0,1]);
0503 
0504    stim = mk_stim_patterns(4,1,'{mono}',[0,1],{},1);
0505    unit_test_cmp('no_balance_inj: t0',length(stim), 4,1);
0506    unit_test_cmp('no_balance_inj: t1',stim(2).stim_pattern, [0;1;0;0]);
0507 
0508 %      'balance_meas' / 'no_balance_meas'
0509 %         -> do / don't subtrant measurement from all electrodes so total
0510 %            average measurement is zero (useful for mono patterns)
0511    stim = mk_stim_patterns(4,1,[0,1],'{mono}',{'no_balance_meas','meas_current'},1);
0512    unit_test_cmp('no_balance_meas: t0',length(stim), 4,1);
0513    unit_test_cmp('no_balance_meas: t1',stim(2).stim_pattern, [0;-1;1;0]);
0514    unit_test_cmp('no_balance_meas: t1',stim(2).meas_pattern, eye(4));
0515 
0516    stim = mk_stim_patterns(4,1,[0,1],'{mono}',{'meas_current'},1);
0517    unit_test_cmp('no_balance_meas: t0',length(stim), 4,1);
0518    unit_test_cmp('no_balance_meas: t1',stim(2).stim_pattern, [0;-1;1;0]);
0519    unit_test_cmp('no_balance_meas: t1',stim(2).meas_pattern, eye(4));
0520 
0521    stim = mk_stim_patterns(4,1,[0,1],'{mono}',{'no_meas_current'},1);
0522    unit_test_cmp('no_balance_meas: t0',length(stim), 4,1);
0523    unit_test_cmp('no_balance_meas: t1',stim(2).stim_pattern, [0;-1;1;0]);
0524    unit_test_cmp('no_balance_meas: t1',stim(2).meas_pattern, [1,0,0,0;0,0,0,1]);
0525 
0526    stim = mk_stim_patterns(4,1,[0,1],'{mono}',{},1); % DO WE WANT THIS AS DEFAULT??
0527    unit_test_cmp('no_balance_meas: t0',length(stim), 4,1);
0528    unit_test_cmp('no_balance_meas: t1',stim(2).stim_pattern, [0;-1;1;0]);
0529    unit_test_cmp('no_balance_meas: t1',stim(2).meas_pattern, [1,0,0,0;0,0,0,1]);
0530 
0531    stim = mk_stim_patterns(4,1,[0,1],'{mono}',{'balance_meas','meas_current'},1);
0532    unit_test_cmp('balance_meas: t0',length(stim), 4);
0533    unit_test_cmp('balance_meas: t1',stim(2).stim_pattern, [0;-1;1;0]);
0534    unit_test_cmp('balance_meas: t1',stim(2).meas_pattern, (4*eye(4)-ones(4))/3);
0535 
0536    stim = mk_stim_patterns(4,1,[0,1],[0,1],{},2);
0537    unit_test_cmp('amplitude: t1',stim(2).stim_pattern, [0;-2;2;0]);
0538    stim = mk_stim_patterns(4,1,'{ad}',[0,1],{},2);
0539    unit_test_cmp('amplitude: t2',stim(2).stim_pattern, [0;-2;2;0]);
0540    stim = mk_stim_patterns(4,1,'{mono}',[0,1],{'no_balance_inj'},2);
0541    unit_test_cmp('amplitude: t3',stim(2).stim_pattern, [0;2;0;0]);
0542    stim = mk_stim_patterns(4,1,'{mono}',[0,1],{},2);
0543    unit_test_cmp('amplitude: t4',stim(2).stim_pattern, [0;2;0;0]);
0544   
0545    [stim,msel] = mk_stim_patterns(6,1,[0,2],[0,1],{'meas_current','no_rotate_meas'},2);
0546    msel = reshape(msel, 6, 6);
0547    unit_test_cmp('meas_sel: t1',msel(:,6), [1;1;1;1;1;1]);
0548 
0549    [stim,msel] = mk_stim_patterns(6,1,[0,2],[0,1],{'meas_current','rotate_meas'},2);
0550    msel = reshape(msel, 6, 6);
0551    unit_test_cmp('meas_sel: t1',msel(:,6), [1;1;1;1;1;1]);
0552 
0553    [stim,msel] = mk_stim_patterns(6,1,[0,1],[0,1],{'no_meas_current','no_rotate_meas'},2);
0554    msel = reshape(msel, 6, 6);
0555    unit_test_cmp('meas_sel: t1',msel(:,6), [0;1;1;1;0;0]);
0556 
0557    [stim,msel] = mk_stim_patterns(6,1,[0,2],[0,1],{'meas_current','no_rotate_meas'},2);
0558    msel = reshape(msel, 6, 6);
0559    unit_test_cmp('meas_sel: t1',msel(:,6), [1;1;1;1;1;1]);
0560 
0561    [stim,msel] = mk_stim_patterns(6,1,[0,1],[0,1],{'no_meas_current','no_rotate_meas'},2);
0562    msel = reshape(msel, 6, 6);
0563    unit_test_cmp('meas_sel: nnp01',msel(:,[4,5]), [1,1;1,1;0,1;0,0;0,0;1,0]);
0564 
0565    [stim,msel] = mk_stim_patterns(6,1,[0,2],[0,1],{'no_meas_current','no_rotate_meas'},2);
0566    msel = reshape(msel, 6, 6);
0567    unit_test_cmp('meas_sel: nnp02',msel(:,[4,5]), [1,0;1,1;0,1;0,0;0,0;0,0]);
0568 
0569    [stim,msel] = mk_stim_patterns(6,1,[0,3],[0,1],{'no_meas_current','no_rotate_meas'},2);
0570    msel = reshape(msel, 6, 6);
0571    unit_test_cmp('meas_sel: nnp03',msel(:,[4,5]), [0,0;1,0;0,1;0,0;1,0;0,1]);
0572 
0573    [stim,msel] = mk_stim_patterns(6,1,[1,2],[0,1],{'no_meas_current','no_rotate_meas'},2);
0574    msel = reshape(msel, 6, 6);
0575    unit_test_cmp('meas_sel: nnp12',msel(:,[4,5]), [1,0;1,1;1,1;0,1;0,0;0,0]);
0576 
0577    [stim,msel] = mk_stim_patterns(6,1,[2,4],[0,1],{'no_meas_current','no_rotate_meas'},2);
0578    msel = reshape(msel, 6, 6);
0579    unit_test_cmp('meas_sel: nnp24',msel(:,[4,5]), [0,0;0,0;1,0;1,1;0,1;0,0]);
0580 
0581    [stim,msel] = mk_stim_patterns(6,1,[2,4],[3,6],{'no_meas_current','no_rotate_meas'},2);
0582    msel = reshape(msel, 6, 6);
0583    unit_test_cmp('meas_sel: nnp2436',msel(:,[4,5]), [1,0;0,1;0,0;1,0;0,1;0,0]);
0584 
0585    [stim,msel] = mk_stim_patterns(6,1,[0,1],[0,1],{'no_meas_current','rotate_meas'},2);
0586    msel = reshape(msel, 6, 6);
0587    unit_test_cmp('meas_sel: nrp01',msel(:,6), [0;0;1;1;1;0]);
0588 
0589    [stim,msel] = mk_stim_patterns(6,1,[0,2],[0,1],{'no_meas_current','rotate_meas'},2);
0590    msel = reshape(msel, 6, 6);
0591    unit_test_cmp('meas_sel: nrp02',msel(:,6), [0;0;0;1;1;0]);
0592 
0593    [stim,msel] = mk_stim_patterns(6,1,[0,2],[3,4],{'no_meas_current','rotate_meas'},2);
0594    msel = reshape(msel, 6, 6);
0595    unit_test_cmp('meas_sel: nrp0234',msel(:,6), [1;1;0;0;0;0]);
0596 
0597    [stim,msel] = mk_stim_patterns(6,1,[0,1],[0,1],{'no_meas_current','no_rotate_meas','no_redundant'},2);
0598    msel = reshape(msel, 6, 6);
0599    unit_test_cmp('meas_sel: nnp01',msel(:,[2,5]), [0,0;0,0;0,0;1,0;1,0;1,0]);
0600 
0601    [stim,msel] = mk_stim_patterns(6,1,[0,2],'{mono}',{'no_meas_current','no_rotate_meas'},2);
0602    msel = reshape(msel, 6, 6);
0603    unit_test_cmp('meas_sel: nnp2436',msel(:,[4,5]), [1,0;1,1;1,1;0,1;1,0;0,1]);
0604 
0605    [stim,msel] = mk_stim_patterns(6,1,[0,2],'{mono}',{'no_meas_current_next0'},2);
0606    msel = reshape(msel, 6, 6);
0607    unit_test_cmp('meas_sel: nnp2436',msel(:,[4,5]), [1,0;1,1;1,1;0,1;1,0;0,1]);
0608 
0609    [stim,msel] = mk_stim_patterns(16,1,[0,1],[0,1],{'no_meas_current_next1','no_rotate_meas'},1);
0610    unit_test_cmp('meas_sel: next1a',msel(48+(1:16)), [1;0;0;0;0;0;1;1;1;1;1;1;1;1;1;1]);
0611    [stim,msel] = mk_stim_patterns(16,1,[0,5],[0,5],{'no_meas_current_next1','no_rotate_meas'},1);
0612    unit_test_cmp('meas_sel: next1b',msel(48+(1:16)), [1;1;0;0;0;1;1;0;0;0;1;1;1;0;0;0]);
0613    [stim,msel] = mk_stim_patterns(16,1,[0,5],[0,5],{'no_meas_current_next1','rotate_meas'},1);
0614    unit_test_cmp('meas_sel: next1c',msel(48+(1:16)), [0;0;1;1;0;0;0;1;1;1;0;0;0;1;1;0]);
0615 
0616    [stim,msel] = mk_stim_patterns(16,1,[0,1],[0,1],{'no_meas_current_next2','no_rotate_meas'},1);
0617    unit_test_cmp('meas_sel: next2a',msel(48+(1:16)), [0;0;0;0;0;0;0;1;1;1;1;1;1;1;1;1]);
0618    [stim,msel] = mk_stim_patterns(16,1,[0,5],[0,5],{'no_meas_current_next2','no_rotate_meas'},1);
0619    unit_test_cmp('meas_sel: next2b',msel(48+(1:16)), [0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0]);
0620    [stim,msel] = mk_stim_patterns(16,1,[0,5],[0,5],{'no_meas_current_next2','rotate_meas'},1);
0621    unit_test_cmp('meas_sel: next2c',msel(48+(1:16)), [0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0]);
0622 
0623    [stim,msel] = mk_stim_patterns(16,1,[0,1],[0,1],{'no_meas_current_next3','no_rotate_meas'},1);
0624    unit_test_cmp('meas_sel: next3a',msel(48+(1:16)), [0;0;0;0;0;0;0;0;1;1;1;1;1;1;1;0]);
0625    [stim,msel] = mk_stim_patterns(16,1,[0,5],[0,5],{'no_meas_current_next3','no_rotate_meas'},1);
0626    unit_test_cmp('meas_sel: next3b',msel(48+(1:16)), [0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0]);
0627    [stim,msel] = mk_stim_patterns(16,1,[0,5],[0,5],{'no_meas_current_next3','rotate_meas'},1);
0628    unit_test_cmp('meas_sel: next3c',msel(48+(1:16)), [0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0]);
0629 
0630 
0631    mselr=[];mseln=[]; for i=1:4; switch i;
0632       case 1;nmc= 'no_meas_current';
0633       case 2;nmc= 'no_meas_current_next1';
0634       case 3;nmc= 'no_meas_current_next2';
0635       case 4;nmc= 'no_meas_current_next3';
0636       otherwise; error 'huh?';
0637       end
0638       [~,mselr(:,i)] = mk_stim_patterns(32,1,[0,5],[0,5],{nmc,'rotate_meas'},1);
0639       [~,mseln(:,i)] = mk_stim_patterns(32,1,[0,5],[0,5],{nmc,'no_rotate_meas'},1);
0640    end
0641    ccv = zeros(32,1); ccv([1,2,end])=1/3;
0642    rsp = reshape(mseln,32,32,4); nxl= 0*rsp;
0643    for i=2:4; for j=1:32
0644       nxl(:,j,i) = cconv(rsp(:,j,i-1),ccv,32)>1-.01;
0645    end; end
0646    unit_test_cmp('meas_sel: next*c',nxl(:,:,2:4),rsp(:,:,2:4));
0647    rsp = reshape(mselr,32,32,4); nxl= 0*rsp;
0648    for i=2:4; for j=1:32
0649       nxl(:,j,i) = cconv(rsp(:,j,i-1),ccv,32)>1-.01;
0650    end; end
0651    unit_test_cmp('meas_sel: next*c',nxl(:,:,2:4),rsp(:,:,2:4));
0652       
0653 
0654 
0655 % TESTS FROM OLD mk_stim_patterns_test CODE
0656    pat= mk_stim_patterns(16,1,'{ad}','{ad}',{}, 1);
0657    test_adj(pat);
0658 
0659    options= {'no_rotate_meas'};
0660    pat= mk_stim_patterns(16,1,'{ad}','{ad}', options,1);
0661    test_adj(pat);
0662 
0663    options= {'no_rotate_meas', 'no_meas_current'};
0664    pat= mk_stim_patterns(16,1,'{ad}','{ad}', options,1);
0665    test_adj(pat);
0666 
0667    options= {'no_rotate_meas', 'meas_current'};
0668    pat= mk_stim_patterns(16,1,'{ad}','{ad}', options,1);
0669    test_adj_full(pat);
0670 
0671    options= {'meas_current'};
0672    pat= mk_stim_patterns(16,1,'{ad}','{ad}', options,1);
0673    test_adj_full(pat);
0674 
0675    options= {'rotate_meas'};
0676    pat= mk_stim_patterns(16,1,'{ad}','{ad}', options,1);
0677    test_adj_rotate(pat);
0678 
0679    options= {'rotate_meas', 'no_meas_current'};
0680    pat= mk_stim_patterns(16,1,'{ad}','{ad}', options,1);
0681    test_adj_rotate(pat);
0682 
0683    options= {'rotate_meas','no_redundant', 'no_meas_current'};
0684    pat= mk_stim_patterns(16,1,'{ad}','{ad}', options,1);
0685    test_adj_no_redund(pat);
0686 
0687 function ok= test_adj(pat)
0688    %%%  test adjacent current pattern
0689 
0690    unit_test_cmp('pt#01', length(pat), 16);
0691    unit_test_cmp('pt#02', pat(1).stimulation, 'Amp');
0692    unit_test_cmp('pt#03', pat(1).stim_pattern, [-1;1;zeros(14,1)]); % Stim pattern # 1
0693 
0694    meas= pat(1).meas_pattern;
0695    unit_test_cmp('pt#04', size(meas), [13 16] );
0696    unit_test_cmp('pt#05', meas(1,:), [0,0,1,-1,zeros(1,12)] );
0697    unit_test_cmp('pt#06', meas(13,:), [zeros(1,14),1,-1] );
0698    unit_test_cmp('pt#07', pat(10).stim_pattern , [zeros(9,1);-1;1;zeros(5,1)] );
0699 
0700    meas= pat(10).meas_pattern;
0701    unit_test_cmp('pt#08', size(meas), [13 16] );
0702    unit_test_cmp('pt#09', meas(1,:), [1,-1,zeros(1,14)] );
0703    unit_test_cmp('pt#10', meas(13,:), [-1,zeros(1,14),1] );
0704 
0705 function ok= test_adj_full(pat)
0706    %%% test adjacent current pattern (full)
0707 
0708    unit_test_cmp('pt#11', length(pat), 16);
0709    unit_test_cmp('pt#12', pat(1).stimulation, 'Amp');
0710    unit_test_cmp('pt#13', pat(1).stim_pattern, [-1;1;zeros(14,1)]);
0711 
0712    meas= pat(1).meas_pattern;
0713    unit_test_cmp('pt#14', size(meas), [16 16] );
0714    unit_test_cmp('pt#15', meas(1,:), [1,-1,zeros(1,14)] );
0715    unit_test_cmp('pt#16', meas(13,:), [zeros(1,12),1,-1,0,0] );
0716    unit_test_cmp('pt#17', pat(10).stim_pattern, [zeros(9,1);-1;1;zeros(5,1)] );
0717 
0718    meas= pat(10).meas_pattern;
0719    unit_test_cmp('pt#18', size(meas), [16 16] );
0720    unit_test_cmp('pt#19', meas(1,:), [1,-1,zeros(1,14)] );
0721    unit_test_cmp('pt#20', meas(13,:), [zeros(1,12),1,-1,0,0] );
0722 
0723 
0724 function ok= test_adj_rotate(pat)
0725    %%%% test adjacent current pattern (rotate)
0726 
0727    unit_test_cmp('pt#21', length(pat), 16);
0728    unit_test_cmp('pt#22', pat(1).stimulation, 'Amp');
0729    unit_test_cmp('pt#23', pat(1).stim_pattern, [-1;1;zeros(14,1)]);
0730 
0731    meas= pat(1).meas_pattern;
0732    unit_test_cmp('pt#24', size(meas), [13 16] );
0733    unit_test_cmp('pt#25', meas(1,:), [0,0,1,-1,zeros(1,12)] );
0734    unit_test_cmp('pt#26', meas(13,:), [zeros(1,14),1,-1] );
0735    unit_test_cmp('pt#27', pat(10).stim_pattern, [zeros(9,1);-1;1;zeros(5,1)] );
0736 
0737    meas= pat(10).meas_pattern;
0738    unit_test_cmp('pt#28', size(meas), [13 16] );
0739    unit_test_cmp('pt#29', meas(1,:), [zeros(1,11),1,-1,zeros(1,3)] );
0740    unit_test_cmp('pt#30', meas(13,:), [zeros(1,7),1,-1,zeros(1,7)] );
0741 
0742 function ok= test_adj_no_redund(pat)
0743    %%% test adjacent current pattern (rotate)
0744 
0745    unit_test_cmp('pt#31', length(pat), 14);
0746    unit_test_cmp('pt#32', pat(1).stimulation, 'Amp');
0747    unit_test_cmp('pt#33', pat(1).stim_pattern, [-1;1;zeros(14,1)]);
0748 
0749    meas= pat(1).meas_pattern;
0750    unit_test_cmp('pt#34', size(meas), [13 16] );
0751    unit_test_cmp('pt#35', meas(1,:), [0,0,1,-1,zeros(1,12)] );
0752    unit_test_cmp('pt#36', meas(13,:), [zeros(1,14),1,-1] );
0753    unit_test_cmp('pt#37', pat(10).stim_pattern, [zeros(9,1);-1;1;zeros(5,1)] );
0754 
0755    meas= pat(10).meas_pattern;
0756    unit_test_cmp('pt#38', size(meas), [5 16] );
0757    unit_test_cmp('pt#39', meas(1,:), [zeros(1,11),1,-1,zeros(1,3)] );
0758    unit_test_cmp('pt#40', meas(5,:), [-1,zeros(1,14),1] );
0759