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 = 'mA'
     (#).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.5mA
   if all electrodes are used for measurement
 stim(1).stim_pattern= [0.5;-0.5;0;0]
 stim(1).meas_pattern= [1, 0, 0,-1 
                       -1, 1, 0, 0 
                        0,-1, 1, 0 
                        0, 0,-1, 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
      '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 = 1mA

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