stim_meas_list

STIM_MEAS_LIST: mk stimulation pattern from list of electrodes
 [stim, meas_sel]= stim_meas_list( sp_mp , Nelec, amplitude, gain);
 [sp_mp]         = stim_meas_list( stim);

 stim =    EIDORS stimulation structure
     use: fwd_model.stimulation = stim;
 meas_sel =EIDORS meas_select (select all values specified)
                to form part of a fwd_model object

 sp_mp = matrix [stim+, stim-, meas+, meas-] x N_patterns

 Nelec  = Num of electrodes         DEFAULT = max sp_mp
 current= drive current levels,     DEFAULT = .010 Amp
 gain   = gain on voltage channels  DEFAULT = 1

 example: stim_meas_list([1,2,3,4; 1,2,4,5])
    create stim pattern in to elecs 1,2 with differential
    measurements on electrodes 3,4 and 4,5;
    then convert the stim struct back to a list of stim/meas pairs
  s=stim_meas_list([1,2,3,4; 1,2,4,5]);
  stim_meas_list(s)
  ans=[1,2,3,4; 1,2,4,5]

0001 function [stim, meas_sel]= stim_meas_list( sp_mp , Nelec, current, gain);
0002 %STIM_MEAS_LIST: mk stimulation pattern from list of electrodes
0003 % [stim, meas_sel]= stim_meas_list( sp_mp , Nelec, amplitude, gain);
0004 % [sp_mp]         = stim_meas_list( stim);
0005 %
0006 % stim =    EIDORS stimulation structure
0007 %     use: fwd_model.stimulation = stim;
0008 % meas_sel =EIDORS meas_select (select all values specified)
0009 %                to form part of a fwd_model object
0010 %
0011 % sp_mp = matrix [stim+, stim-, meas+, meas-] x N_patterns
0012 %
0013 % Nelec  = Num of electrodes         DEFAULT = max sp_mp
0014 % current= drive current levels,     DEFAULT = .010 Amp
0015 % gain   = gain on voltage channels  DEFAULT = 1
0016 %
0017 % example: stim_meas_list([1,2,3,4; 1,2,4,5])
0018 %    create stim pattern in to elecs 1,2 with differential
0019 %    measurements on electrodes 3,4 and 4,5;
0020 %    then convert the stim struct back to a list of stim/meas pairs
0021 %  s=stim_meas_list([1,2,3,4; 1,2,4,5]);
0022 %  stim_meas_list(s)
0023 %  ans=[1,2,3,4; 1,2,4,5]
0024 
0025 % (C) 2010,2015 Andy Adler, Alistair Boyle. License: GPL version 2 or version 3
0026 % $Id: stim_meas_list.m 5838 2018-10-23 03:33:21Z aadler $
0027 
0028 if ischar(sp_mp) && strcmp(sp_mp,'UNIT_TEST'); do_unit_test; return; end
0029 
0030 if isstruct(sp_mp)
0031   stim = sp_mp;
0032   meas_sel = [];
0033   nst = length(stim);
0034   nvt = 0;
0035   for i = 1:nst;
0036     nvt = nvt + size(stim(i).stim_pattern, 2) * size(stim(i).meas_pattern, 1);
0037   end
0038   stim = flatten_stim(stim, nst, nvt);
0039   stim = stim(:,[2 1 4 3]);
0040   return
0041 end
0042 
0043 if nargin <2; Nelec = max(sp_mp(:)); end
0044 if nargin <3; current = 1;           end
0045 if nargin <4; gain    = 1;           end
0046 
0047 if any(sp_mp(:) > Nelec);
0048     error('Electrode patterns require more electrodes than Nelec');
0049 end
0050 stim = struct([]);
0051 Npat = size(sp_mp,1);
0052 is = 0;
0053 current = current.*ones(Npat,1); % Extend if required
0054 for i=1:Npat
0055    cur = sp_mp(i,1:2);
0056    new_stim = sparse( cur, 1, current(i)*[-1,1], Nelec,1);
0057    % create a new stim if it isn't the same as the last one
0058    if (is < 1) || any(any(stim(is).stim_pattern ~= new_stim))
0059      is = is + 1;
0060    end
0061    stim(is).stimulation = 'Amp';
0062    stim(is).stim_pattern = new_stim;
0063    mes = sp_mp(i,3:4);
0064    if isfield(stim(is),'meas_pattern') % append pattern if required
0065      stim(is).meas_pattern = [ stim(is).meas_pattern; sparse( 1, mes, gain *  [-1,1], 1, Nelec)];
0066    else
0067      stim(is).meas_pattern = sparse( 1, mes, gain *  [-1,1], 1, Nelec);
0068    end
0069 end
0070 
0071 % take in a stim/meas struct
0072 % return a matrix of stim/meas pairs, per row with drive current 'i' and measurement gain 'g' calculated
0073 % [ +s -s +m -m i g ]
0074 function stim_flat = flatten_stim(stim, nst, nvt)
0075   stim_flat = zeros(nvt, 6);
0076   idx = 1;
0077   % TODO calculate 'gain' when it matched (m+ == - m-)
0078   % TODO calculate 'gain' when it is unmatched (m+ ~= m-)
0079   % TODO calculate 'current' when it matched (s+ == - s-)
0080   % TODO calculate 'current' when it is unmatched (s+ ~= s-)
0081   for i = 1:nst
0082       nmp= size(stim(i).meas_pattern, 1); % number of measurement patterns for this stim pair
0083       [sp, jnk, spv]= find(stim(i).stim_pattern>0);
0084       [sn, jnk, snv]= find(stim(i).stim_pattern<0);      
0085       [order, mp, mpv]= find(stim(i).meas_pattern>0);  
0086       [~, idx2sort] = sort(order); mp = mp(idx2sort); mpv = mpv(idx2sort);
0087       [order, mn, mnv]= find(stim(i).meas_pattern<0);    
0088       [~, idx2sort] = sort(order); mn = mn(idx2sort); mnv = mnv(idx2sort);
0089       % expand s+/s- to match the size of m+/m-
0090       sp  = zeros(nmp,1)+sp;
0091       sn  = zeros(nmp,1)+sn;
0092       spv = zeros(nmp,1)+spv;
0093       snv = zeros(nmp,1)+snv;
0094       stim_flat(idx:idx+nmp-1,:) = ...
0095         [ sp sn ... % stim pairs
0096             mp mn spv mpv];  % meas pairs
0097       idx = idx + nmp;
0098   end
0099 
0100 function  do_unit_test
0101    imdl = mk_common_model('a2c0',16);
0102    img = mk_image(imdl);
0103    list_in = [1,2,3,4;1,2,4,5];
0104    img.fwd_model.stimulation = stim_meas_list(list_in,16);
0105    list_out = stim_meas_list(img.fwd_model.stimulation);
0106    unit_test_cmp('pattern#1', list_in, list_out);
0107 
0108 
0109    vh = fwd_solve(img);
0110    list_in = [6,7,3,4;1,2,4,5];
0111    img.fwd_model.stimulation = stim_meas_list(list_in,16);
0112    list_out = stim_meas_list(img.fwd_model.stimulation);
0113    unit_test_cmp('pattern#2', list_in, list_out);
0114 
0115    vh = fwd_solve(img);
0116    
0117    
0118    nElecs = 16;
0119    stim = mk_stim_patterns( nElecs, 1, '{ad}', '{ad}');
0120    sp_mp = stim_meas_list(stim);
0121    inj_diff = mod(sp_mp(:,2) - sp_mp(:,1), nElecs);
0122    meas_diff = mod(sp_mp(:,3) - sp_mp(:,4), nElecs);
0123    unit_test_cmp('pattern#3', true, all(inj_diff == 1) && all(meas_diff == 1));
0124    
0125    stim = mk_stim_patterns( nElecs, 1, '{ad}', '{ad}', {'no_meas_current', 'no_rotate_meas'});
0126    sp_mp = stim_meas_list(stim);
0127    inj_diff = mod(sp_mp(:,2) - sp_mp(:,1), nElecs);
0128    meas_diff = mod(sp_mp(:,3) - sp_mp(:,4), nElecs);
0129    unit_test_cmp('pattern#4', true, all(inj_diff == 1) && all(meas_diff == 1));
0130    
0131    stim = mk_stim_patterns( nElecs, 1, '{ad}', '{ad}', {'no_meas_current'});
0132    sp_mp = stim_meas_list(stim);
0133    inj_diff = mod(sp_mp(:,2) - sp_mp(:,1), nElecs);
0134    meas_diff = mod(sp_mp(:,3) - sp_mp(:,4), nElecs);
0135    unit_test_cmp('pattern#5', true, all(inj_diff == 1) && all(meas_diff == 1));
0136    
0137    Skip = 3;
0138    stim = mk_stim_patterns( nElecs, 1, [0 1+Skip], [0 1+Skip], {'no_meas_current', 'no_rotate_meas'});
0139    sp_mp = stim_meas_list(stim);
0140    inj_diff = mod(sp_mp(:,2) - sp_mp(:,1), nElecs);
0141    meas_diff = mod(sp_mp(:,3) - sp_mp(:,4), nElecs);
0142    unit_test_cmp('pattern#6', true, all(inj_diff == Skip+1) && all(meas_diff == Skip+1));
0143    
0144    nElecs = 32;
0145    Skip = 7;
0146    stim = mk_stim_patterns( nElecs, 1, [0 1+Skip], [0 1+Skip], {'meas_current', 'rotate_meas'});
0147    sp_mp = stim_meas_list(stim);
0148    inj_diff = mod(sp_mp(:,2) - sp_mp(:,1), nElecs);
0149    meas_diff = mod(sp_mp(:,3) - sp_mp(:,4), nElecs);
0150 
0151 
0152    unit_test_cmp('pattern#7', true, all(inj_diff == Skip+1) && all(meas_diff == Skip+1));
0153    
0154    
0155    
0156