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