function change_thresh(path1,proc,walkspeed,muscles,chngHA)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% change_thresh.m allows one to re-do EMG stimulation thresholding
% for a given set of muscles and trials.  Note that chngHA is not
% currently implemented.
%
% Inputs:
% path1     - data directory (subject name and date).
% proc      - processing method (e.g. 'SangZaj' or 'lowpass').
% walkspeed - string of walking speed to process. 
% muscles   - cell of muscles to have thresholds re-done.
% chngHA    - 0 to keep HAM, HAB muscle selections as is; 1 to reasses.
%
% Outputs:
% Re-saves avg_act_thresh file and thresh_values file.
%
% Example Usage:
% change_thresh('EM072611','low100_mvc_lo','walk125',{'TA'})
%
% Jared Markowitz
% Spring 2013
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

addpath('~/Dropbox/harvard/processMOT');

% Load Perry muscles for reference:
load ~/Documents/MATLAB/data/Perry_muscles/stim_act.mat
p_sVAS  = p_sVLAT;   clear p_sVLAT;
p_sADDL = p_sAL;     clear p_sAL;
p_sADDM = p_sAM;     clear p_sAM;

eval(['load ''~/Documents/MATLAB/data/' path1 '/processed/EMG/avg_act/' proc '/' walkspeed '.mat'';']);

numpts = length(m_aTA);
Fs = 1000;
dt = 1/Fs;

for i = 1:length(muscles)
  eval(['n_s' muscles{i} 'full = n_s' muscles{i} ';']);  
  eval(['m_s' muscles{i} 'full = m_s' muscles{i} ';']);
  eval(['sd_s' muscles{i} 'full = sd_s' muscles{i} ';']);  
  eval(['s' muscles{i} 'full = s' muscles{i} ';']);
  eval(['n_a' muscles{i} 'full = n_a' muscles{i} ';']);  
  eval(['m_a' muscles{i} 'full = m_a' muscles{i} ';']);
  eval(['sd_a' muscles{i} 'full = sd_a' muscles{i} ';']);  
  eval(['a' muscles{i} 'full = a' muscles{i} ';']);
end;  
  
eval(['load ''~/Documents/MATLAB/data/' path1 '/processed/EMG/avg_act_thresh/' proc '/' walkspeed '.mat'';']);

for nn = 1:size(muscles,1)
  eval(['figure(' num2str(nn+30) ');']); 
  eval(['niceBars(0:(numpts-1),m_s' muscles{nn} 'full,sd_s' muscles{nn} 'full,''b'',0.5);']);  
  hold on;
  eval(['niceBars(0:(numpts-1),m_s' muscles{nn} ',sd_s' muscles{nn} ',''y'',0.5);']);    
  hold on;
  eval(['plot(0:(numpts-1),m_s' muscles{nn} ',''color'',[1 .5 0],''linewidth'',3);']);
  eval(['ythresh1 = thresh_values.' muscles{nn} ';']); 
  plot([0 100],[ythresh1 ythresh1],'r','linewidth',4);
  grid on;
  eval(['figure(' num2str(nn+50) ');']); 
  eval(['niceBars(0:(numpts-1),m_a' muscles{nn} 'full,sd_a' muscles{nn} 'full,''b'',0.5);']);
  hold on;
  eval(['figure(' num2str(nn+50) ');']); 
  eval(['niceBars(0:(numpts-1),m_a' muscles{nn} ',sd_a' muscles{nn} ',''y'',0.5);']);  
  hold on;
  grid on;
end;

for pp = 1:size(muscles,1)  
  eval(['figure(' num2str(pp+30) ');']); 
  hold on;
  if strcmp(muscles{pp},'BFLH')
    plot(p_sSMEM,'k','linewidth',4);
  elseif strcmp(muscles{pp},'GLU')
    plot(p_sGLU_low,'k','linewidth',4);
    plot(p_sGLU_upp,'y','linewidth',4);        
  elseif strcmp(muscles{pp},'GLUU')
    plot(p_sGLU_upp,'k','linewidth',4);        
  else  
    eval(['plot(p_s' muscles{pp} ',''k'',''linewidth'',4);']);
  end;
  axis tight;
  resp1 = 'n';
  h1 = plot([0 100],[0 0],'r','linewidth',3);
  h2 = plot([0 100],[0 0],'m','linewidth',3);
  while sum(strcmp(resp1,{'Y';'y'})) == 0
    delete(h1);
    delete(h2);
    fprintf('Click vertical level for threshold:\n');
    [xthresh,ythresh] = ginput(1);
    hold on;
    h1 = plot([0 100],[ythresh ythresh],'r','linewidth',4);
    eval(['q1 = m_s' muscles{pp} 'full;']);
    ind = find(q1<=ythresh);
    q1(ind) = 0;
    h2 = plot(0:(numpts-1),q1,'m','linewidth',4);
    resp1 = input('Good enough? (y/n):\n','s');
  end;  
  eval(['thresh_values.' muscles{pp} '=ythresh;']);
end;  

% Perform thresholding:
for ppp = 1:size(muscles,1)
  eval(['num_gc = size(n_a' muscles{ppp} 'full,1);']);
  eval(['full_stim_' muscles{ppp} ' = [];']);
  for qq = 1:num_gc
    eval(['full_stim_' muscles{ppp} ' = cat(1,full_stim_' muscles{ppp} ',s' muscles{ppp} 'full{qq}) ;']);
  end;  
  eval(['ind=find(full_stim_' muscles{ppp} '<=thresh_values.' muscles{ppp} ');']);
  eval(['thresh_stim_' muscles{ppp} ' = full_stim_' muscles{ppp} ';']);
  eval(['thresh_stim_' muscles{ppp} '(ind) = 0;']);  
end;

% Pass remaining stim through Raasch bilinear dynamics:
for rr = 1:size(muscles,1)
  eval(['tt = (0:1/Fs:(length(thresh_stim_' muscles{rr} ')-1)/Fs)'';']);
  simopts = simset('FixedStep',1/Fs,'SrcWorkspace','current');
  if strcmp(muscles{rr},'SOL') == 1
    mfft = 0.20;
  elseif strcmp(muscles{rr},'GAS') == 1
    mfft = 0.50;
  elseif strcmp(muscles{rr},'TA') == 1
    mfft = 0.35;
  elseif strcmp(muscles{rr},'BFSH') == 1    
    mfft = 0.35;
  elseif strcmp(muscles{rr},'BFLH') == 1
    mfft = 0.35;
  elseif strcmp(muscles{rr},'SMEM') == 1
    mfft = 0.35;
  elseif strcmp(muscles{rr},'VAS') == 1
    mfft = 0.50;
  elseif strcmp(muscles{rr},'RF') == 1
    mfft = 0.65;
  elseif strcmp(muscles{rr},'ILL') == 1    
    mfft = 0.50;
  elseif strcmp(muscles{rr},'GLU') == 1
    mfft = 0.45;
  elseif strcmp(muscles{rr},'GMED') == 1
    mfft = 0.50;
  elseif strcmp(muscles{rr},'GLUU') == 1
    mfft = 0.50;  
  elseif strcmp(muscles{rr},'ADDL') == 1
    mfft = 0.35;
  elseif strcmp(muscles{rr},'ADDM') == 1
    mfft = 0.45;
  else
    fprintf('Erroneous muscle label.\n');
    return;
  end;      
  if strcmp(proc(end-1:end),'sh') == 1
    k1 = 1/(.01*8.4/(4.8*(1+1.5*mfft)))-1/(.05*8.4/(4.8*(1+1.5*mfft)));
    k2 = 1/(.05*8.4/(4.8*(1+1.5*mfft)));
  elseif strcmp(proc(end-1:end),'lo') == 1
    k1 = 1/(.08-.047*mfft)-1/(.09-.056*mfft);
    k2 = 1/(.09-.056*mfft);
  else fprintf('Invalid proc.\n');
    return;
  end;
  eval(['x_in = [tt thresh_stim_' muscles{rr} '];']);
  sim('raasch.mdl',tt(end),simopts);
  eval(['a' muscles{rr} '_thresh = {};']);
  eval(['s' muscles{rr} '_thresh = {};']);
  eval(['n_a' muscles{rr} '_thresh = [];']);
  eval(['n_s' muscles{rr} '_thresh = [];']);
  eval(['num_gc = size(n_a' muscles{rr} ',1);']);
  ind1 = 1;
  for qqq = 1:num_gc
    eval(['ind2 = ind1-1+length(s' muscles{rr} '{qqq});']);
    eval(['a' muscles{rr} '_thresh = [a' muscles{rr} '_thresh; {a_out(ind1:ind2)}];']);
    eval(['s' muscles{rr} '_thresh = [s' muscles{rr} '_thresh; {thresh_stim_' muscles{rr} '(ind1:ind2)}];']);
    eval(['x1 = 0:dt:(length(a' muscles{rr} '_thresh{qqq})-1)*dt;']);
    x2 = linspace(0,x1(end),numpts);
    eval(['i_a' muscles{rr} '_thresh=interp1(x1,a' muscles{rr} '_thresh{qqq},x2,''spline'');']);
    %eval(['i_a' muscles{rr} '_thresh=i_a' muscles{rr} '_thresh-min(i_a' muscles{rr} '_thresh);']);
    eval(['n_a' muscles{rr} '_thresh=[n_a' muscles{rr} '_thresh; i_a' muscles{rr} '_thresh];']);
    eval(['i_s' muscles{rr} '_thresh=interp1(x1,s' muscles{rr} '_thresh{qqq},x2,''spline'');']);
    %eval(['i_s' muscles{rr} '_thresh=i_s' muscles{rr} '_thresh-min(i_s' muscles{rr} '_thresh);']);
    eval(['n_s' muscles{rr} '_thresh=[n_s' muscles{rr} '_thresh; i_s' muscles{rr} '_thresh];']);  
    ind1 = ind2 + 1;
  end;
  eval(['m_a' muscles{rr} '_thresh = mean(n_a' muscles{rr} '_thresh,1);']);  
  eval(['m_s' muscles{rr} '_thresh = mean(n_s' muscles{rr} '_thresh,1);']);  
  eval(['sd_a' muscles{rr} '_thresh = std(n_a' muscles{rr} '_thresh,0,1);']);  
  eval(['sd_s' muscles{rr} '_thresh = std(n_s' muscles{rr} '_thresh,0,1);']);  
  % Plot:
  eval(['figure(' num2str(rr+50) ');']);
  hold on; 
  eval(['niceBars(0:(numpts-1),m_a' muscles{rr} '_thresh,sd_a' muscles{rr} '_thresh,''g'',0.5);']);
  axis tight;
  eval(['figure(' num2str(rr+30) ');']);    
  hold on;
  eval(['niceBars(0:(numpts-1),m_s' muscles{rr} '_thresh,sd_s' muscles{rr} '_thresh,''g'',0.5);']);
  axis tight;
end;

eval(['outfile_thr = ''~/Documents/MATLAB/data/' path1 '/processed/EMG/avg_act_thresh/' proc '/' walkspeed '.mat'';']); 

varlist = who('*_thresh');
for jjj = 1:length(varlist)
  eval([varlist{jjj}(1:end-7) '=' varlist{jjj} ';']);
end;
clear *_thresh;
clear *full*
save(outfile_thr,'m_*','n_*','sd_*','HA*in','a*','s*','thresh_values');

if strcmp(proc(1:3),'low') == 1; 
  procr = 'low100';
elseif strcmp(proc(1:4),'sang') == 1; 
  procr = 'sang';
else
  fprintf('Faulty proc input.\n');
  return;
end;
path_out3 = ['~/Documents/MATLAB/data/' path1 '/processed/EMG/thresh_values/' procr ];
outfile3 = [path_out3 '/' walkspeed '.mat'];
save(outfile3,'thresh_values');

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if exist('chngHA') > 0 && chngHA == 1
  fprintf(['Current HAMin = ' num2str(HAMin) '.\n']);
  fprintf(['Current HABin = ' num2str(HABin) '.\n']);  
  HAMin = 2; HABin = 2;
  while ((HAMin==0)+(HAMin==1)) == 0 
    HAMin = input('Enter HAMin value:\n');
  end;
  while ((HABin==0)+(HABin==1)) == 0 
    HABin = input('Enter HABin value:\n');
  end;
  if HAMin == 0
    m_aHAM = m_aSMEM; n_aHAM = n_aSMEM; sd_aHAM = sd_aSMEM; aHAM = aSMEM;
    m_sHAM = m_sSMEM; n_sHAM = n_sSMEM; sd_sHAM = sd_sSMEM; sHAM = sSMEM;   
  elseif HAMin == 1
    m_aHAM = m_aBFLH; n_aHAM = n_aBFLH; sd_aHAM = sd_aBFLH; aHAM = aBFLH;
    m_sHAM = m_sBFLH; n_sHAM = n_sBFLH; sd_sHAM = sd_sBFLH; sHAM = sBFLH;
  else    
  end;
  if HABin == 0
    m_aHAB = m_aGLUU; n_aHAB = n_aGLUU; sd_aHAB = sd_aGLUU; aHAB = aGLUU;
    m_sHAB = m_sGLUU; n_sHAB = n_sGLUU; sd_sHAB = sd_sGLUU; sHAB = sGLUU;   
  elseif HABin == 1
    m_aHAB = m_aGMED; n_aHAB = n_aGMED; sd_aHAB = sd_aGMED; aHAB = aGMED;
    m_sHAB = m_sGMED; n_sHAB = n_sGMED; sd_sHAB = sd_sGMED; sHAB = sGMED;   
  else
  end;
  fprintf(['New HAMin = ' num2str(HAMin) '.\n']);
  fprintf(['New HABin = ' num2str(HABin) '.\n']);    
end;  
