From 7fbaf66418be67242a62c5e680af7e3c7044b7de Mon Sep 17 00:00:00 2001 From: Usman Rashid Date: Sun, 6 Oct 2019 14:27:15 +1300 Subject: [PATCH] Removed bugs in PSOt integration --- emgEventsDetectTool.m | 12 +- libs/PSOt/pso_Trelea_vectorized.m | 833 +++++++++++++------------ optimisation/autoFindActivations.m | 4 +- optimisation/autoFindActivationsNoTB.m | 26 +- optimisation/psoStatus.m | 2 +- optimisation/tuneDetectionParams.m | 78 ++- optimisation/tuneDetectionParamsNoTB.m | 52 ++ 7 files changed, 550 insertions(+), 457 deletions(-) create mode 100644 optimisation/tuneDetectionParamsNoTB.m diff --git a/emgEventsDetectTool.m b/emgEventsDetectTool.m index ad49dc9..87a05a6 100644 --- a/emgEventsDetectTool.m +++ b/emgEventsDetectTool.m @@ -296,8 +296,18 @@ function autoFind(~, ~) end function autoTune(~, ~) + + dlgOpts.Interpreter='tex'; + x = inputdlg({'Randomise parameters [Y:1, N:0]'},... + 'nOptim', [1 50], {'0'}, dlgOpts); + if(isempty(x)) + return; + end + answer = x; + isRandomParams = str2double(answer{1}); + paramsVector = tuneDetectionParams(vars.channelStream(:, vars.channelNum), vars.fs,... - vars.detectionParams(:, vars.channelNum),... + isRandomParams, vars.detectionParams(:, vars.channelNum),... vars.options.paramLowerBounds, vars.options.paramUpperBounds,... vars.detectionCellarray{vars.channelNum}{end, vars.ONSETS_COLUMN_NUM}, vars.detectionCellarray{vars.channelNum}{end, vars.OFFSETS_COLUMN_NUM},... vars.options.detectionAlgo); diff --git a/libs/PSOt/pso_Trelea_vectorized.m b/libs/PSOt/pso_Trelea_vectorized.m index 1a40825..372a2af 100644 --- a/libs/PSOt/pso_Trelea_vectorized.m +++ b/libs/PSOt/pso_Trelea_vectorized.m @@ -1,6 +1,6 @@ % pso_Trelea_vectorized.m % a generic particle swarm optimizer -% to find the minimum or maximum of any +% to find the minimum or maximum of any % MISO matlab function % % Implements Common, Trelea type 1 and 2, and Clerc's class 1". It will @@ -22,15 +22,15 @@ % Inputs: % functname - string of matlab function to optimize % D - # of inputs to the function (dimension of problem) -% +% % Optional Inputs: % mv - max particle velocity, either a scalar or a vector of length D -% (this allows each component to have it's own max velocity), +% (this allows each component to have it's own max velocity), % default = 4, set if not input or input as NaN % -% VarRange - matrix of ranges for each input variable, +% VarRange - matrix of ranges for each input variable, % default -100 to 100, of form: -% [ min1 max1 +% [ min1 max1 % min2 max2 % ... % minD maxD ] @@ -40,7 +40,7 @@ % = 2, funct is targeted to P(12) (minimizes distance to errgoal) % % PSOparams - PSO parameters -% P(1) - Epochs between updating display, default = 100. if 0, +% P(1) - Epochs between updating display, default = 100. if 0, % no display % P(2) - Maximum number of iterations (epochs) to train, default = 2000. % P(3) - population size, default = 24 @@ -50,10 +50,10 @@ % P(6) - Initial inertia weight, default = 0.9 % P(7) - Final inertia weight, default = 0.4 % P(8) - Epoch when inertial weight at final value, default = 1500 -% P(9)- minimum global error gradient, -% if abs(Gbest(i+1)-Gbest(i)) < gradient over +% P(9)- minimum global error gradient, +% if abs(Gbest(i+1)-Gbest(i)) < gradient over % certain length of epochs, terminate run, default = 1e-25 -% P(10)- epochs before error gradient criterion terminates run, +% P(10)- epochs before error gradient criterion terminates run, % default = 150, if the SSE does not change over 250 epochs % then exit % P(11)- error goal, if NaN then unconstrained min or max, default=NaN @@ -95,77 +95,77 @@ rand('state',sum(100*clock)); if nargin < 2 - error('Not enough arguments.'); + error('Not enough arguments.'); end % PSO PARAMETERS if nargin == 2 % only specified functname and D - VRmin=ones(D,1)*-100; - VRmax=ones(D,1)*100; - VR=[VRmin,VRmax]; - minmax = 0; - P = []; - mv = 4; - plotfcn='goplotpso'; + VRmin=ones(D,1)*-100; + VRmax=ones(D,1)*100; + VR=[VRmin,VRmax]; + minmax = 0; + P = []; + mv = 4; + plotfcn='goplotpso'; elseif nargin == 3 % specified functname, D, and mv - VRmin=ones(D,1)*-100; - VRmax=ones(D,1)*100; - VR=[VRmin,VRmax]; - minmax = 0; - mv=varargin{1}; - if isnan(mv) - mv=4; - end - P = []; - plotfcn='goplotpso'; + VRmin=ones(D,1)*-100; + VRmax=ones(D,1)*100; + VR=[VRmin,VRmax]; + minmax = 0; + mv=varargin{1}; + if isnan(mv) + mv=4; + end + P = []; + plotfcn='goplotpso'; elseif nargin == 4 % specified functname, D, mv, Varrange - mv=varargin{1}; - if isnan(mv) - mv=4; - end - VR=varargin{2}; - minmax = 0; - P = []; - plotfcn='goplotpso'; + mv=varargin{1}; + if isnan(mv) + mv=4; + end + VR=varargin{2}; + minmax = 0; + P = []; + plotfcn='goplotpso'; elseif nargin == 5 % Functname, D, mv, Varrange, and minmax - mv=varargin{1}; - if isnan(mv) - mv=4; - end - VR=varargin{2}; - minmax=varargin{3}; - P = []; - plotfcn='goplotpso'; + mv=varargin{1}; + if isnan(mv) + mv=4; + end + VR=varargin{2}; + minmax=varargin{3}; + P = []; + plotfcn='goplotpso'; elseif nargin == 6 % Functname, D, mv, Varrange, minmax, and psoparams - mv=varargin{1}; - if isnan(mv) - mv=4; - end - VR=varargin{2}; - minmax=varargin{3}; - P = varargin{4}; % psoparams - plotfcn='goplotpso'; + mv=varargin{1}; + if isnan(mv) + mv=4; + end + VR=varargin{2}; + minmax=varargin{3}; + P = varargin{4}; % psoparams + plotfcn='goplotpso'; elseif nargin == 7 % Functname, D, mv, Varrange, minmax, and psoparams, plotfcn - mv=varargin{1}; - if isnan(mv) - mv=4; - end - VR=varargin{2}; - minmax=varargin{3}; - P = varargin{4}; % psoparams - plotfcn = varargin{5}; + mv=varargin{1}; + if isnan(mv) + mv=4; + end + VR=varargin{2}; + minmax=varargin{3}; + P = varargin{4}; % psoparams + plotfcn = varargin{5}; elseif nargin == 8 % Functname, D, mv, Varrange, minmax, and psoparams, plotfcn, PSOseedValue - mv=varargin{1}; - if isnan(mv) - mv=4; - end - VR=varargin{2}; - minmax=varargin{3}; - P = varargin{4}; % psoparams - plotfcn = varargin{5}; - PSOseedValue = varargin{6}; -else - error('Wrong # of input arguments.'); + mv=varargin{1}; + if isnan(mv) + mv=4; + end + VR=varargin{2}; + minmax=varargin{3}; + P = varargin{4}; % psoparams + plotfcn = varargin{5}; + PSOseedValue = varargin{6}; +else + error('Wrong # of input arguments.'); end % sets up default pso params @@ -189,440 +189,445 @@ % used with trainpso, for neural net training if strcmp(functname,'pso_neteval') - net = evalin('caller','net'); + net = evalin('caller','net'); Pd = evalin('caller','Pd'); Tl = evalin('caller','Tl'); Ai = evalin('caller','Ai'); - Q = evalin('caller','Q'); + Q = evalin('caller','Q'); TS = evalin('caller','TS'); end % error checking - if ((minmax==2) & isnan(errgoal)) - error('minmax= 2, errgoal= NaN: choose an error goal or set minmax to 0 or 1'); - end +if ((minmax==2) & isnan(errgoal)) + error('minmax= 2, errgoal= NaN: choose an error goal or set minmax to 0 or 1'); +end - if ( (PSOseed==1) & ~exist('PSOseedValue') ) - error('PSOseed flag set but no PSOseedValue was input'); - end +if ( (PSOseed==1) & ~exist('PSOseedValue') ) + error('PSOseed flag set but no PSOseedValue was input'); +end + +if exist('PSOseedValue') + tmpsz=size(PSOseedValue); + if D < tmpsz(2) + error('PSOseedValue column size must be D or less'); + end + if ps < tmpsz(1) + error('PSOseedValue row length must be # of particles or less'); + end +end - if exist('PSOseedValue') - tmpsz=size(PSOseedValue); - if D < tmpsz(2) - error('PSOseedValue column size must be D or less'); - end - if ps < tmpsz(1) - error('PSOseedValue row length must be # of particles or less'); - end - end - % set plotting flag if (P(1))~=0 - plotflg=1; + plotflg=1; else - plotflg=0; + plotflg=0; end % preallocate variables for speed up - tr = ones(1,me)*NaN; +tr = ones(1,me)*NaN; % take care of setting max velocity and position params here if length(mv)==1 - velmaskmin = -mv*ones(ps,D); % min vel, psXD matrix - velmaskmax = mv*ones(ps,D); % max vel -elseif length(mv)==D - velmaskmin = repmat(forcerow(-mv),ps,1); % min vel - velmaskmax = repmat(forcerow( mv),ps,1); % max vel + velmaskmin = -mv*ones(ps,D); % min vel, psXD matrix + velmaskmax = mv*ones(ps,D); % max vel +elseif length(mv)==D + velmaskmin = repmat(forcerow(-mv),ps,1); % min vel + velmaskmax = repmat(forcerow( mv),ps,1); % max vel else - error('Max vel must be either a scalar or same length as prob dimension D'); + error('Max vel must be either a scalar or same length as prob dimension D'); end posmaskmin = repmat(VR(1:D,1)',ps,1); % min pos, psXD matrix posmaskmax = repmat(VR(1:D,2)',ps,1); % max pos posmaskmeth = 3; % 3=bounce method (see comments below inside epoch loop) % PLOTTING - message = sprintf('PSO: %%g/%g iterations, GBest = %%20.20g.\n',me); - +message = sprintf('PSO: %%g/%g iterations, GBest = %%20.20g.\n',me); + % INITIALIZE INITIALIZE INITIALIZE INITIALIZE INITIALIZE INITIALIZE - + % initialize population of particles and their velocities at time zero, % format of pos= (particle#, dimension) - % construct random population positions bounded by VR - pos(1:ps,1:D) = normmat(rand([ps,D]),VR',1); - - if PSOseed == 1 % initial positions user input, see comments above +% construct random population positions bounded by VR +pos(1:ps,1:D) = normmat(rand([ps,D]),VR',1); + +if PSOseed == 1 % initial positions user input, see comments above tmpsz = size(PSOseedValue); - pos(1:tmpsz(1),1:tmpsz(2)) = PSOseedValue; - end + pos(1:tmpsz(1),1:tmpsz(2)) = PSOseedValue; +end - % construct initial random velocities between -mv,mv - vel(1:ps,1:D) = normmat(rand([ps,D]),... - [forcecol(-mv),forcecol(mv)]',1); +% construct initial random velocities between -mv,mv +vel(1:ps,1:D) = normmat(rand([ps,D]),... + [forcecol(-mv),forcecol(mv)]',1); % initial pbest positions vals - pbest = pos; +pbest = pos; -% VECTORIZE THIS, or at least vectorize cost funct call - out = feval(functname,pos); % returns column of cost values (1 for each particle) +% VECTORIZE THIS, or at least vectorize cost funct call +out = feval(functname,pos); % returns column of cost values (1 for each particle) %--------------------------- - - pbestval=out; % initially, pbest is same as pos + +pbestval=out; % initially, pbest is same as pos % assign initial gbest here also (gbest and gbestval) - if minmax==1 - % this picks gbestval when we want to maximize the function +if minmax==1 + % this picks gbestval when we want to maximize the function [gbestval,idx1] = max(pbestval); - elseif minmax==0 - % this works for straight minimization +elseif minmax==0 + % this works for straight minimization [gbestval,idx1] = min(pbestval); - elseif minmax==2 - % this works when you know target but not direction you need to go - % good for a cost function that returns distance to target that can be either - % negative or positive (direction info) +elseif minmax==2 + % this works when you know target but not direction you need to go + % good for a cost function that returns distance to target that can be either + % negative or positive (direction info) [temp,idx1] = min((pbestval-ones(size(pbestval))*errgoal).^2); gbestval = pbestval(idx1); - end +end + +% preallocate a variable to keep track of gbest for all iters +bestpos = zeros(me,D+1)*NaN; +gbest = pbest(idx1,:); % this is gbest position +% used with trainpso, for neural net training +% assign gbest to net at each iteration, these interim assignments +% are for plotting mostly +if strcmp(functname,'pso_neteval') + net=setx(net,gbest); +end +%tr(1) = gbestval; % save for output +bestpos(1,1:D) = gbest; - % preallocate a variable to keep track of gbest for all iters - bestpos = zeros(me,D+1)*NaN; - gbest = pbest(idx1,:); % this is gbest position - % used with trainpso, for neural net training - % assign gbest to net at each iteration, these interim assignments - % are for plotting mostly - if strcmp(functname,'pso_neteval') - net=setx(net,gbest); - end - %tr(1) = gbestval; % save for output - bestpos(1,1:D) = gbest; - % this part used for implementing Carlisle and Dozier's APSO idea % slightly modified, this tracks the global best as the sentry whereas % their's chooses a different point to act as sentry % see "Tracking Changing Extremea with Adaptive Particle Swarm Optimizer", % part of the WAC 2002 Proceedings, June 9-13, http://wacong.com - sentryval = gbestval; - sentry = gbest; - +sentryval = gbestval; +sentry = gbest; + if (trelea == 3) -% calculate Clerc's constriction coefficient chi to use in his form - kappa = 1; % standard val = 1, change for more or less constriction - if ( (ac1+ac2) <=4 ) - chi = kappa; - else - psi = ac1 + ac2; - chi_den = abs(2-psi-sqrt(psi^2 - 4*psi)); - chi_num = 2*kappa; - chi = chi_num/chi_den; - end + % calculate Clerc's constriction coefficient chi to use in his form + kappa = 1; % standard val = 1, change for more or less constriction + if ( (ac1+ac2) <=4 ) + chi = kappa; + else + psi = ac1 + ac2; + chi_den = abs(2-psi-sqrt(psi^2 - 4*psi)); + chi_num = 2*kappa; + chi = chi_num/chi_den; + end end % INITIALIZE END INITIALIZE END INITIALIZE END INITIALIZE END rstflg = 0; % for dynamic environment checking % start PSO iterative procedures - cnt = 0; % counter used for updating display according to df in the options - cnt2 = 0; % counter used for the stopping subroutine based on error convergence - iwt(1) = iw1; +cnt = 0; % counter used for updating display according to df in the options +cnt2 = 0; % counter used for the stopping subroutine based on error convergence +iwt(1) = iw1; for i=1:me % start epoch loop (iterations) - - out = feval(functname,[pos;gbest]); - outbestval = out(end,:); - out = out(1:end-1,:); - - tr(i+1) = gbestval; % keep track of global best val - te = i; % returns epoch number to calling program when done - bestpos(i,1:D+1) = [gbest,gbestval]; - - %assignin('base','bestpos',bestpos(i,1:D+1)); - %------------------------------------------------------------------------ - % this section does the plots during iterations - if plotflg==1 - if (rem(i,df) == 0 ) | (i==me) | (i==1) - fprintf(message,i,gbestval); - cnt = cnt+1; % count how many times we display (useful for movies) - - eval(plotfcn); % defined at top of script - - end % end update display every df if statement + + out = feval(functname,[pos;gbest]); + outbestval = out(end,:); + out = out(1:end-1,:); + + tr(i+1) = gbestval; % keep track of global best val + te = i; % returns epoch number to calling program when done + bestpos(i,1:D+1) = [gbest,gbestval]; + + %assignin('base','bestpos',bestpos(i,1:D+1)); + %------------------------------------------------------------------------ + % this section does the plots during iterations + if plotflg==1 + if (rem(i,df) == 0 ) | (i==me) | (i==1) + fprintf(message,i,gbestval); + cnt = cnt+1; % count how many times we display (useful for movies) + + eval(plotfcn); % defined at top of script + + end % end update display every df if statement end % end plotflg if statement - + % check for an error space that changes wrt time/iter - % threshold value that determines dynamic environment + % threshold value that determines dynamic environment % sees if the value of gbest changes more than some threshold value % for the same location chkdyn = 1; rstflg = 0; % for dynamic environment checking - + if chkdyn==1 - threshld = 0.05; % percent current best is allowed to change, .05 = 5% etc - letiter = 5; % # of iterations before checking environment, leave at least 3 so PSO has time to converge - outorng = abs( 1- (outbestval/gbestval) ) >= threshld; - samepos = (max( sentry == gbest )); - - if (outorng & samepos) & rem(i,letiter)==0 - rstflg=1; - % disp('New Environment: reset pbest, gbest, and vel'); - %% reset pbest and pbestval if warranted -% outpbestval = feval( functname,[pbest] ); -% Poutorng = abs( 1-(outpbestval./pbestval) ) > threshld; -% pbestval = pbestval.*~Poutorng + outpbestval.*Poutorng; -% pbest = pbest.*repmat(~Poutorng,1,D) + pos.*repmat(Poutorng,1,D); - - pbest = pos; % reset personal bests to current positions - pbestval = out; - vel = vel*10; % agitate particles a little (or a lot) + threshld = 0.05; % percent current best is allowed to change, .05 = 5% etc + letiter = 5; % # of iterations before checking environment, leave at least 3 so PSO has time to converge + outorng = abs( 1- (outbestval/gbestval) ) >= threshld; + samepos = (max( sentry == gbest )); - % recalculate best vals - if minmax == 1 - [gbestval,idx1] = max(pbestval); - elseif minmax==0 - [gbestval,idx1] = min(pbestval); - elseif minmax==2 % this section needs work - [temp,idx1] = min((pbestval-ones(size(pbestval))*errgoal).^2); - gbestval = pbestval(idx1); - end + if (outorng & samepos) & rem(i,letiter)==0 + rstflg=1; + % disp('New Environment: reset pbest, gbest, and vel'); + %% reset pbest and pbestval if warranted + % outpbestval = feval( functname,[pbest] ); + % Poutorng = abs( 1-(outpbestval./pbestval) ) > threshld; + % pbestval = pbestval.*~Poutorng + outpbestval.*Poutorng; + % pbest = pbest.*repmat(~Poutorng,1,D) + pos.*repmat(Poutorng,1,D); + + pbest = pos; % reset personal bests to current positions + pbestval = out; + vel = vel*10; % agitate particles a little (or a lot) + + % recalculate best vals + if minmax == 1 + [gbestval,idx1] = max(pbestval); + elseif minmax==0 + [gbestval,idx1] = min(pbestval); + elseif minmax==2 % this section needs work + [temp,idx1] = min((pbestval-ones(size(pbestval))*errgoal).^2); + gbestval = pbestval(idx1); + end + + gbest = pbest(idx1,:); + + % used with trainpso, for neural net training + % assign gbest to net at each iteration, these interim assignments + % are for plotting mostly + if strcmp(functname,'pso_neteval') + net=setx(net,gbest); + end + end % end if outorng - gbest = pbest(idx1,:); + sentryval = gbestval; + sentry = gbest; - % used with trainpso, for neural net training - % assign gbest to net at each iteration, these interim assignments - % are for plotting mostly - if strcmp(functname,'pso_neteval') - net=setx(net,gbest); - end - end % end if outorng - - sentryval = gbestval; - sentry = gbest; - end % end if chkdyn - % find particles where we have new pbest, depending on minmax choice + % find particles where we have new pbest, depending on minmax choice % then find gbest and gbestval - %[size(out),size(pbestval)] + %[size(out),size(pbestval)] if rstflg == 0 - if minmax == 0 - [tempi] = find(pbestval>=out); % new min pbestvals - pbestval(tempi,1) = out(tempi); % update pbestvals - pbest(tempi,:) = pos(tempi,:); % update pbest positions - - [iterbestval,idx1] = min(pbestval); - - if gbestval >= iterbestval - gbestval = iterbestval; - gbest = pbest(idx1,:); - % used with trainpso, for neural net training - % assign gbest to net at each iteration, these interim assignments - % are for plotting mostly - if strcmp(functname,'pso_neteval') - net=setx(net,gbest); - end - end - elseif minmax == 1 - [tempi,dum] = find(pbestval<=out); % new max pbestvals - pbestval(tempi,1) = out(tempi,1); % update pbestvals - pbest(tempi,:) = pos(tempi,:); % update pbest positions - - [iterbestval,idx1] = max(pbestval); - if gbestval <= iterbestval - gbestval = iterbestval; - gbest = pbest(idx1,:); - % used with trainpso, for neural net training - % assign gbest to net at each iteration, these interim assignments - % are for plotting mostly - if strcmp(functname,'pso_neteval') - net=setx(net,gbest); - end + if minmax == 0 + [tempi] = find(pbestval>=out); % new min pbestvals + pbestval(tempi,1) = out(tempi); % update pbestvals + pbest(tempi,:) = pos(tempi,:); % update pbest positions + + [iterbestval,idx1] = min(pbestval); + + if gbestval >= iterbestval + gbestval = iterbestval; + gbest = pbest(idx1,:); + % used with trainpso, for neural net training + % assign gbest to net at each iteration, these interim assignments + % are for plotting mostly + if strcmp(functname,'pso_neteval') + net=setx(net,gbest); + end + end + elseif minmax == 1 + [tempi,dum] = find(pbestval<=out); % new max pbestvals + pbestval(tempi,1) = out(tempi,1); % update pbestvals + pbest(tempi,:) = pos(tempi,:); % update pbest positions + + [iterbestval,idx1] = max(pbestval); + if gbestval <= iterbestval + gbestval = iterbestval; + gbest = pbest(idx1,:); + % used with trainpso, for neural net training + % assign gbest to net at each iteration, these interim assignments + % are for plotting mostly + if strcmp(functname,'pso_neteval') + net=setx(net,gbest); + end + end + elseif minmax == 2 % this won't work as it is, fix it later + egones = errgoal*ones(ps,1); % vector of errgoals + sqrerr2 = ((pbestval-egones).^2); + sqrerr1 = ((out-egones).^2); + [tempi,dum] = find(sqerr1 <= sqrerr2); % find particles closest to targ + pbestval(tempi,1) = out(tempi,1); % update pbestvals + pbest(tempi,:) = pos(tempi,:); % update pbest positions + + sqrerr = ((pbestval-egones).^2); % need to do this to reflect new pbests + [temp,idx1] = min(sqrerr); + iterbestval = pbestval(idx1); + + if (iterbestval-errgoal)^2 <= (gbestval-errgoal)^2 + gbestval = iterbestval; + gbest = pbest(idx1,:); + % used with trainpso, for neural net training + % assign gbest to net at each iteration, these interim assignments + % are for plotting mostly + if strcmp(functname,'pso_neteval') + net=setx(net,gbest); + end + end end - elseif minmax == 2 % this won't work as it is, fix it later - egones = errgoal*ones(ps,1); % vector of errgoals - sqrerr2 = ((pbestval-egones).^2); - sqrerr1 = ((out-egones).^2); - [tempi,dum] = find(sqerr1 <= sqrerr2); % find particles closest to targ - pbestval(tempi,1) = out(tempi,1); % update pbestvals - pbest(tempi,:) = pos(tempi,:); % update pbest positions - - sqrerr = ((pbestval-egones).^2); % need to do this to reflect new pbests - [temp,idx1] = min(sqrerr); - iterbestval = pbestval(idx1); - - if (iterbestval-errgoal)^2 <= (gbestval-errgoal)^2 - gbestval = iterbestval; - gbest = pbest(idx1,:); - % used with trainpso, for neural net training - % assign gbest to net at each iteration, these interim assignments - % are for plotting mostly - if strcmp(functname,'pso_neteval') - net=setx(net,gbest); - end - end - end end - % % build a simple predictor 10th order, for gbest trajectory - % if i>500 - % for dimcnt=1:D - % pred_coef = polyfit(i-250:i,(bestpos(i-250:i,dimcnt))',20); - % % pred_coef = polyfit(200:i,(bestpos(200:i,dimcnt))',20); - % gbest_pred(i,dimcnt) = polyval(pred_coef,i+1); - % end - % else -% gbest_pred(i,:) = zeros(size(gbest)); -% end - - %gbest_pred(i,:)=gbest; - %assignin('base','gbest_pred',gbest_pred); - - % % convert to non-inertial frame - % gbestoffset = gbest - gbest_pred(i,:); - % gbest = gbest - gbestoffset; - % pos = pos + repmat(gbestoffset,ps,1); - % pbest = pbest + repmat(gbestoffset,ps,1); - - %PSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSO - - % get new velocities, positions (this is the heart of the PSO algorithm) - % each epoch get new set of random numbers - rannum1 = rand([ps,D]); % for Trelea and Clerc types - rannum2 = rand([ps,D]); - if trelea == 2 + % % build a simple predictor 10th order, for gbest trajectory + % if i>500 + % for dimcnt=1:D + % pred_coef = polyfit(i-250:i,(bestpos(i-250:i,dimcnt))',20); + % % pred_coef = polyfit(200:i,(bestpos(200:i,dimcnt))',20); + % gbest_pred(i,dimcnt) = polyval(pred_coef,i+1); + % end + % else + % gbest_pred(i,:) = zeros(size(gbest)); + % end + + %gbest_pred(i,:)=gbest; + %assignin('base','gbest_pred',gbest_pred); + + % % convert to non-inertial frame + % gbestoffset = gbest - gbest_pred(i,:); + % gbest = gbest - gbestoffset; + % pos = pos + repmat(gbestoffset,ps,1); + % pbest = pbest + repmat(gbestoffset,ps,1); + + %PSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSO + + % get new velocities, positions (this is the heart of the PSO algorithm) + % each epoch get new set of random numbers + rannum1 = rand([ps,D]); % for Trelea and Clerc types + rannum2 = rand([ps,D]); + if trelea == 2 % from Trelea's paper, parameter set 2 - vel = 0.729.*vel... % prev vel - +1.494.*rannum1.*(pbest-pos)... % independent - +1.494.*rannum2.*(repmat(gbest,ps,1)-pos); % social - elseif trelea == 1 - % from Trelea's paper, parameter set 1 - vel = 0.600.*vel... % prev vel - +1.700.*rannum1.*(pbest-pos)... % independent - +1.700.*rannum2.*(repmat(gbest,ps,1)-pos); % social - elseif trelea ==3 + vel = 0.729.*vel... % prev vel + +1.494.*rannum1.*(pbest-pos)... % independent + +1.494.*rannum2.*(repmat(gbest,ps,1)-pos); % social + elseif trelea == 1 + % from Trelea's paper, parameter set 1 + vel = 0.600.*vel... % prev vel + +1.700.*rannum1.*(pbest-pos)... % independent + +1.700.*rannum2.*(repmat(gbest,ps,1)-pos); % social + elseif trelea ==3 % Clerc's Type 1" PSO - vel = chi*(vel... % prev vel - +ac1.*rannum1.*(pbest-pos)... % independent - +ac2.*rannum2.*(repmat(gbest,ps,1)-pos)) ; % social - else - % common PSO algo with inertia wt + vel = chi*(vel... % prev vel + +ac1.*rannum1.*(pbest-pos)... % independent + +ac2.*rannum2.*(repmat(gbest,ps,1)-pos)) ; % social + else + % common PSO algo with inertia wt % get inertia weight, just a linear funct w.r.t. epoch parameter iwe - if i<=iwe + if i<=iwe iwt(i) = ((iw2-iw1)/(iwe-1))*(i-1)+iw1; - else + else iwt(i) = iw2; - end + end % random number including acceleration constants - ac11 = rannum1.*ac1; % for common PSO w/inertia - ac22 = rannum2.*ac2; - - vel = iwt(i).*vel... % prev vel - +ac11.*(pbest-pos)... % independent - +ac22.*(repmat(gbest,ps,1)-pos); % social - end - - % limit velocities here using masking - vel = ( (vel <= velmaskmin).*velmaskmin ) + ( (vel > velmaskmin).*vel ); - vel = ( (vel >= velmaskmax).*velmaskmax ) + ( (vel < velmaskmax).*vel ); + ac11 = rannum1.*ac1; % for common PSO w/inertia + ac22 = rannum2.*ac2; - % update new position (PSO algo) - pos = pos + vel; + vel = iwt(i).*vel... % prev vel + +ac11.*(pbest-pos)... % independent + +ac22.*(repmat(gbest,ps,1)-pos); % social + end - % position masking, limits positions to desired search space - % method: 0) no position limiting, 1) saturation at limit, - % 2) wraparound at limit , 3) bounce off limit - minposmask_throwaway = pos <= posmaskmin; % these are psXD matrices - minposmask_keep = pos > posmaskmin; - maxposmask_throwaway = pos >= posmaskmax; - maxposmask_keep = pos < posmaskmax; - - if posmaskmeth == 1 - % this is the saturation method - pos = ( minposmask_throwaway.*posmaskmin ) + ( minposmask_keep.*pos ); - pos = ( maxposmask_throwaway.*posmaskmax ) + ( maxposmask_keep.*pos ); - elseif posmaskmeth == 2 - % this is the wraparound method - pos = ( minposmask_throwaway.*posmaskmax ) + ( minposmask_keep.*pos ); - pos = ( maxposmask_throwaway.*posmaskmin ) + ( maxposmask_keep.*pos ); - elseif posmaskmeth == 3 - % this is the bounce method, particles bounce off the boundaries with -vel - pos = ( minposmask_throwaway.*posmaskmin ) + ( minposmask_keep.*pos ); - pos = ( maxposmask_throwaway.*posmaskmax ) + ( maxposmask_keep.*pos ); - - vel = (vel.*minposmask_keep) + (-vel.*minposmask_throwaway); - vel = (vel.*maxposmask_keep) + (-vel.*maxposmask_throwaway); - else - % no change, this is the original Eberhart, Kennedy method, - % it lets the particles grow beyond bounds if psoparams (P) - % especially Vmax, aren't set correctly, see the literature - end - - %PSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSO -% check for stopping criterion based on speed of convergence to desired - % error + % limit velocities here using masking + vel = ( (vel <= velmaskmin).*velmaskmin ) + ( (vel > velmaskmin).*vel ); + vel = ( (vel >= velmaskmax).*velmaskmax ) + ( (vel < velmaskmax).*vel ); + + % update new position (PSO algo) + pos = pos + vel; + + % position masking, limits positions to desired search space + % method: 0) no position limiting, 1) saturation at limit, + % 2) wraparound at limit , 3) bounce off limit + minposmask_throwaway = pos <= posmaskmin; % these are psXD matrices + minposmask_keep = pos > posmaskmin; + maxposmask_throwaway = pos >= posmaskmax; + maxposmask_keep = pos < posmaskmax; + + if posmaskmeth == 1 + % this is the saturation method + pos = ( minposmask_throwaway.*posmaskmin ) + ( minposmask_keep.*pos ); + pos = ( maxposmask_throwaway.*posmaskmax ) + ( maxposmask_keep.*pos ); + elseif posmaskmeth == 2 + % this is the wraparound method + pos = ( minposmask_throwaway.*posmaskmax ) + ( minposmask_keep.*pos ); + pos = ( maxposmask_throwaway.*posmaskmin ) + ( maxposmask_keep.*pos ); + elseif posmaskmeth == 3 + % this is the bounce method, particles bounce off the boundaries with -vel + pos = ( minposmask_throwaway.*posmaskmin ) + ( minposmask_keep.*pos ); + pos = ( maxposmask_throwaway.*posmaskmax ) + ( maxposmask_keep.*pos ); + + vel = (vel.*minposmask_keep) + (-vel.*minposmask_throwaway); + vel = (vel.*maxposmask_keep) + (-vel.*maxposmask_throwaway); + else + % no change, this is the original Eberhart, Kennedy method, + % it lets the particles grow beyond bounds if psoparams (P) + % especially Vmax, aren't set correctly, see the literature + end + + %PSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSO + % check for stopping criterion based on speed of convergence to desired + % error tmp1 = abs(tr(i) - gbestval); if tmp1 > ergrd - cnt2 = 0; + cnt2 = 0; elseif tmp1 <= ergrd - cnt2 = cnt2+1; - if cnt2 >= ergrdep - if plotflg == 1 - fprintf(message,i,gbestval); - disp(' '); - disp(['--> Solution likely, GBest hasn''t changed by at least ',... - num2str(ergrd),' for ',... - num2str(cnt2),' epochs.']); - eval(plotfcn); - end - break - end + cnt2 = cnt2+1; + if cnt2 >= ergrdep + if plotflg == 1 + fprintf(message,i,gbestval); + disp(' '); + disp(['--> Solution likely, GBest hasn''t changed by at least ',... + num2str(ergrd),' for ',... + num2str(cnt2),' epochs.']); + haltingState = 1; % Added for exiting from outside + eval(plotfcn); + end + break + end end - % this stops if using constrained optimization and goal is reached + % this stops if using constrained optimization and goal is reached if ~isnan(errgoal) - if ((gbestval<=errgoal) & (minmax==0)) | ((gbestval>=errgoal) & (minmax==1)) - - if plotflg == 1 - fprintf(message,i,gbestval); - disp(' '); - disp(['--> Error Goal reached, successful termination!']); - haltingState = 1; % Added for exiting from outside - eval(plotfcn); - end - break - end - - % this is stopping criterion for constrained from both sides - if minmax == 2 - if ((tr(i)=errgoal)) | ((tr(i)>errgoal) ... - & (gbestval <= errgoal)) - if plotflg == 1 - fprintf(message,i,gbestval); - disp(' '); - disp(['--> Error Goal reached, successful termination!']); - haltingState = 1; % Added for exiting from outside - eval(plotfcn); - end - break - end - end % end if minmax==2 + if ((gbestval<=errgoal) & (minmax==0)) | ((gbestval>=errgoal) & (minmax==1)) + + if plotflg == 1 + fprintf(message,i,gbestval); + disp(' '); + disp(['--> Error Goal reached, successful termination!']); + haltingState = 1; % Added for exiting from outside + eval(plotfcn); + end + break + end + + % this is stopping criterion for constrained from both sides + if minmax == 2 + if ((tr(i)=errgoal)) | ((tr(i)>errgoal) ... + & (gbestval <= errgoal)) + if plotflg == 1 + fprintf(message,i,gbestval); + disp(' '); + disp(['--> Error Goal reached, successful termination!']); + haltingState = 1; % Added for exiting from outside + eval(plotfcn); + end + break + end + end % end if minmax==2 end % end ~isnan if - - % % convert back to inertial frame - % pos = pos - repmat(gbestoffset,ps,1); - % pbest = pbest - repmat(gbestoffset,ps,1); - % gbest = gbest + gbestoffset; - - % Added for exiting from outside - if haltingState - break; - end - + + % % convert back to inertial frame + % pos = pos - repmat(gbestoffset,ps,1); + % pbest = pbest - repmat(gbestoffset,ps,1); + % gbest = gbest + gbestoffset; + + % Added for exiting from outside + if haltingState + break; + end + if i==me + haltingState = 1; + eval(plotfcn); + end + end % end epoch loop %% clear temp outputs % evalin('base','clear temp_pso_out temp_te temp_tr;'); % output & return - OUT=[gbest';gbestval]; - varargout{1}=[1:te]; - varargout{2}=[tr(find(~isnan(tr)))]; - - return \ No newline at end of file +OUT=[gbest';gbestval]; +varargout{1}=[1:te]; +varargout{2}=[tr(find(~isnan(tr)))]; + +return \ No newline at end of file diff --git a/optimisation/autoFindActivations.m b/optimisation/autoFindActivations.m index f1dcb48..5d71294 100644 --- a/optimisation/autoFindActivations.m +++ b/optimisation/autoFindActivations.m @@ -65,7 +65,7 @@ 'SwarmSize', PARAM_MULTIPLIER * totalParams,... 'UseParallel', true,... 'OutputFcn', @psOutFunc,... - 'ObjectiveLimit', 1-1e-3,... + 'ObjectiveLimit', 0,... 'FunctionTolerance', 0.1,... 'HybridFcn', {@fmincon, hybridopts}); else @@ -74,7 +74,7 @@ 'UseParallel', true,... 'OutputFcn', @psOutFunc,... 'InitialSwarmMatrix', swarmMatrix,... - 'ObjectiveLimit', 1-1e-3,... + 'ObjectiveLimit', 0,... 'FunctionTolerance', 0.1,... 'HybridFcn', {@fmincon, hybridopts}); end diff --git a/optimisation/autoFindActivationsNoTB.m b/optimisation/autoFindActivationsNoTB.m index b174011..e144bff 100644 --- a/optimisation/autoFindActivationsNoTB.m +++ b/optimisation/autoFindActivationsNoTB.m @@ -23,19 +23,19 @@ % Problem preparation -PSOparams(1) = 1; -PSOparams(2) = 100; -PSOparams(3) = PARAM_MULTIPLIER*totalParams; -PSOparams(4) = 2; -PSOparams(5) = 2; -PSOparams(6) = 0.9; -PSOparams(7) = 0.4; -PSOparams(8) = 1500; -PSOparams(9) = 0.1; -PSOparams(10) = 20; -PSOparams(11) = 1-1e-3; -PSOparams(12) = 1; -PSOparams(13) = ~randomInit; +PSOparams(1) = 1; +PSOparams(2) = 100; +PSOparams(3) = PARAM_MULTIPLIER*totalParams; +PSOparams(4) = 2; +PSOparams(5) = 2; +PSOparams(6) = 0.9; +PSOparams(7) = 0.4; +PSOparams(8) = 1500; +PSOparams(9) = 0.1; +PSOparams(10) = 20; +PSOparams(11) = 0; +PSOparams(12) = 1; +PSOparams(13) = ~randomInit; % Run particle swarm optimiser, first pass [optOUT,~,~]= pso_Trelea_vectorized(f,totalParams,4,[lb ub], 0, PSOparams, 'psoStatus', swarmMatrix); diff --git a/optimisation/psoStatus.m b/optimisation/psoStatus.m index 131e1d7..9fd97e6 100644 --- a/optimisation/psoStatus.m +++ b/optimisation/psoStatus.m @@ -32,7 +32,7 @@ function d = infoDialog(msg) d = dialog('Position', [300 300 250 150],... - 'Name','nOptim',... + 'Name','Optimisation Status',... 'WindowStyle', 'modal',... 'CloseRequestFcn', @haltAndCloseOperation); diff --git a/optimisation/tuneDetectionParams.m b/optimisation/tuneDetectionParams.m index ad96d6e..6942149 100644 --- a/optimisation/tuneDetectionParams.m +++ b/optimisation/tuneDetectionParams.m @@ -1,4 +1,4 @@ -function paramsVector = tuneDetectionParams(singleChannel, fs, initialParams,... +function paramsVector = tuneDetectionParams(singleChannel, fs, randomInit, initialParams,... lowerBounds, upperBounds, uOnsets, uOffsets, detectionAlgo, lambda, funcVer) %tuneDetectionParams maxConcordance optimisation method. % @@ -10,17 +10,20 @@ % Licensed under the MIT License. See LICENSE in the project root for % license information. -if nargin < 9 - lambda = 0; -end - -if nargin < 9 +if nargin < 10 lambda = 0; funcVer = 'v2'; -elseif nargin < 10 +elseif nargin < 11 funcVer = 'v2'; end +if ~(exist("particleswarm", 'file') && exist("optimoptions", 'file') && exist("fmincon", 'file')) + paramsVector = tuneDetectionParamsNoTB(singleChannel, fs, randomInit, initialParams,... + lowerBounds, upperBounds, uOnsets, uOffsets, detectionAlgo, lambda); + warning(sprintf('MATLAB optimisation functions were not found.\nA third party optimisation toolbox was used.\nResults may be better or worse as full testing has not been conducted.\nInstalling Global Optimisation Toolbox and Optimisation Toolbox may solve the problem.')); + return; +end + % Defaults PARAM_MULTIPLIER = 3; LAMBDA = lambda; @@ -35,26 +38,49 @@ swarmMatrix = initialParams'; % M by nvars if strcmp(funcVer, 'v1') - optimOptions = optimoptions('particleswarm',... - 'SwarmSize', PARAM_MULTIPLIER * totalParams,... - 'UseParallel', true,... - 'OutputFcn', @psOutFunc,... - 'InitialSwarmMatrix', swarmMatrix,... - 'ObjectiveLimit', 0,... - 'FunctionTolerance', 1e-6); + if randomInit + optimOptions = optimoptions('particleswarm',... + 'SwarmSize', PARAM_MULTIPLIER * totalParams,... + 'UseParallel', true,... + 'OutputFcn', @psOutFunc,... + 'ObjectiveLimit', 0,... + 'FunctionTolerance', 1e-6); + else + optimOptions = optimoptions('particleswarm',... + 'SwarmSize', PARAM_MULTIPLIER * totalParams,... + 'UseParallel', true,... + 'OutputFcn', @psOutFunc,... + 'InitialSwarmMatrix', swarmMatrix,... + 'ObjectiveLimit', 0,... + 'FunctionTolerance', 1e-6); + end elseif strcmp(funcVer, 'v2') - hybridopts = optimoptions('fmincon',... - 'Algorithm', 'interior-point',... - 'OutputFcn', @ipOutFunc,... - 'FunctionTolerance', 1e-6); - optimOptions = optimoptions('particleswarm',... - 'SwarmSize', PARAM_MULTIPLIER * totalParams,... - 'UseParallel', true,... - 'OutputFcn', @psOutFunc,... - 'InitialSwarmMatrix', swarmMatrix,... - 'ObjectiveLimit', 0,... - 'FunctionTolerance', 1e-6,... - 'HybridFcn', {@fmincon, hybridopts}); + if randomInit + hybridopts = optimoptions('fmincon',... + 'Algorithm', 'interior-point',... + 'OutputFcn', @ipOutFunc,... + 'FunctionTolerance', 1e-6); + optimOptions = optimoptions('particleswarm',... + 'SwarmSize', PARAM_MULTIPLIER * totalParams,... + 'UseParallel', true,... + 'OutputFcn', @psOutFunc,... + 'ObjectiveLimit', 0,... + 'FunctionTolerance', 1e-6,... + 'HybridFcn', {@fmincon, hybridopts}); + else + hybridopts = optimoptions('fmincon',... + 'Algorithm', 'interior-point',... + 'OutputFcn', @ipOutFunc,... + 'FunctionTolerance', 1e-6); + optimOptions = optimoptions('particleswarm',... + 'SwarmSize', PARAM_MULTIPLIER * totalParams,... + 'UseParallel', true,... + 'OutputFcn', @psOutFunc,... + 'InitialSwarmMatrix', swarmMatrix,... + 'ObjectiveLimit', 0,... + 'FunctionTolerance', 1e-6,... + 'HybridFcn', {@fmincon, hybridopts}); + end end % Run particle swarm optimiser diff --git a/optimisation/tuneDetectionParamsNoTB.m b/optimisation/tuneDetectionParamsNoTB.m new file mode 100644 index 0000000..cdc8df2 --- /dev/null +++ b/optimisation/tuneDetectionParamsNoTB.m @@ -0,0 +1,52 @@ +function paramsVector = tuneDetectionParamsNoTB(singleChannel, fs, randomInit, initialParams,... + lowerBounds, upperBounds, uOnsets, uOffsets, detectionAlgo, lambda) +%tuneDetectionParamsNoTB maxConcordance optimisation method when MATLAB's +% optimisation toolboxes are not installed. +% +% Copyright (c) <2019> +% Licensed under the MIT License. See LICENSE in the project root for +% license information. + +% Defaults +PARAM_MULTIPLIER = 3; +LAMBDA = lambda; + +totalParams = length(initialParams); +f = @(P)costFunc(singleChannel, fs, uOnsets, uOffsets, detectionAlgo, P); +lb = lowerBounds; +ub = upperBounds; +swarmMatrix = initialParams'; % Pop. size by numParams +swarmMatrix = repmat(swarmMatrix, PARAM_MULTIPLIER*totalParams, 1); + + +% Problem preparation +PSOparams(1) = 1; +PSOparams(2) = 100; +PSOparams(3) = PARAM_MULTIPLIER*totalParams; +PSOparams(4) = 2; +PSOparams(5) = 2; +PSOparams(6) = 0.9; +PSOparams(7) = 0.4; +PSOparams(8) = 1500; +PSOparams(9) = 1e-6; +PSOparams(10) = 20; +PSOparams(11) = 0; +PSOparams(12) = 1; +PSOparams(13) = ~randomInit; + +% Run particle swarm optimiser, first pass +[optOUT,~,~]= pso_Trelea_vectorized(f,totalParams,4,[lb ub], 0, PSOparams, 'psoStatus', swarmMatrix); +paramsVector = optOUT(1:end-1); + + + function cost = costFunc(singleChannel, fs, uOnsets, uOffsets, detectionAlgo, P) + cost = zeros(size(P, 1), 1); + for i=1:size(P, 1) + [ onSets, offSets ] = detectionAlgo(singleChannel, fs, P(i, :)); + uBurst = createBusts(singleChannel, uOnsets, uOffsets); + algoBurst = createBusts(singleChannel, onSets, offSets); + concordance = sum(uBurst == algoBurst) / length(uBurst); + cost(i) = norm(1-concordance) + LAMBDA*norm(P); + end + end +end \ No newline at end of file