function [pca_eigenval]=PCAgenerator_1el
disp('Multi-parameter variation in XPP from Matlab, with calculating principal components')
disp(' ')
disp('The parameters varied are wights w1, stimulation strength Kn, gsyn and iapp')
disp(' ')
% user specifications
odeFileName = 'tenStim1electr.ode';
x1range = 0.3:0.1:0.3;
x2range = 0.0:2.0:75.0;
x3range = 0.5:0.1:1.4;
x4range = 4.0:1.0:8.0;
verboseTog = false; % two levels of notification verbosity to user
% initializations
x1_ix = 0;
lenx1 = length(x1range);
lenx2 = length(x2range);
lenx3 = length(x3range);
lenx4 = length(x4range);
X2 = x2range;
X3 = x3range;
X4 = x4range;
numpoints = lenx1 * lenx2 * lenx3 * lenx4;
fprintf('Calculating %d points\n',numpoints)
pointnum = 0;
count = 0; % line print count, for backspacing purposes
newpars(1).name='w1';
newpars(1).type='PAR';
newpars(2).name='Kn';
newpars(2).type='PAR';
newpars(3).name='gsyn';
newpars(3).type='PAR';
newpars(4).name='iapp';
newpars(4).type='PAR';
for i=x1range
x1_ix = x1_ix+1;
pca_eigenval = cell(lenx2,lenx3,lenx4);
x2_ix = 0;
for j = x2range
x2_ix = x2_ix+1;
x3_ix = 0;
for k = x3range
x3_ix = x3_ix+1;
x4_ix = 0;
for l = x4range
x4_ix = x4_ix+1;
pointnum = pointnum + 1;
if ~verboseTog % then print current loop counter
bspstr = '';
if count > 0
for bix=1:count
bspstr = [bspstr '\b'];
end
end
fprintf(bspstr);
count = fprintf('%d',pointnum);
end
newpars(1).val = i;
newpars(2).val = j;
newpars(3).val = k;
newpars(4).val = l;
if verboseTog
fprintf(' **** Point %d: Variable pt = %.4f, c = %.4f\n',pointnum,newpars(1).val,newpars(2).val)
end
success = ChangeXPPodeFile(odeFileName,newpars); % type `help ChangeXPPodeFile` in Matlab
if success==0
disp('Change XPP ode file failed')
return
end
success = RunXPP(odeFileName,'',true,'C:\xppall\xppaut'); % type `help RunXPP` in Matlab
if success==0
disp('Run XPP failed')
return
end
format = '%f32%*f32%*f32%*f32%*f32%*f32%f32%f32%f32%*f32%*f32%*f32%*f32%f32%f32%f32%f32%f32%f32%f32%f32%f32%f32%*[^\n]';
fid = fopen('output.dat', 'r'); % change this if your ODE file explicitly names a data file
data = textscan(fid, format);
fclose(fid);
k_v4 = 0;
k_v5 = 0;
k_v6 = 0;
k_v = 0;
for m = 120001:139981 %count # of spikes during the last minute of time series
if (data{2}(m-1) < -40) && (data{2}(m) >= -40)
k_v4 = k_v4+1;
end
if (data{3}(m-1) < -40) && (data{3}(m) >= -40)
k_v5 = k_v5+1;
end
if (data{4}(m-1) < -40) && (data{4}(m) >= -40)
k_v6 = k_v6+1;
end
if (data{3}(m) < -200) || (data{3}(m) > 200)
k_v=k_v+1;
end
end
if (min([k_v4 k_v5 k_v6]) > 10) && (max([k_v4 k_v5 k_v6]) < 100) && (k_v==0)
pca_eigenval{x2_ix,x3_ix,x4_ix}=pca([data{5}(40001:139981) data{6}(40001:139981) data{7}(40001:139981) data{8}(40001:139981) data{9}(40001:139981) data{10}(40001:139981) data{11}(40001:139981) data{12}(40001:139981) data{13}(40001:139981) data{14}(40001:139981)]);
%t1 = 6000:0.5:6990; left for plotting purposes
% v5 = data{3}(12001:13981);
% fOut1 = sprintf('Volt_tr_%s%1.1f_%s%2.1f_%s%1.1f_%s%1.1f.mat','w',i,'Kn',j,'gsyn',k,'iapp',l);
% save(fOut1, 'v5')
% clear v5
end
clear data
end
end
end
fOut2 = sprintf('PCA_eigenval_%s%1.1f.mat','w',i);
save(fOut2,'X4','X3','X2','pca_eigenval')
clear pca_eigenval
end
fprintf('\n')
beep
return
end