Skip to content
This repository has been archived by the owner on Nov 11, 2024. It is now read-only.

Commit

Permalink
This is the final release of SMB that is compatible with CADET2.3.2
Browse files Browse the repository at this point in the history 
    figurePlots
        The internal column states are used to plot chromatogram
        Column tensor: row-major and column-major issue

    Sorry for the bugs in figurePlots
  • Loading branch information
KimHe committed May 23, 2017
1 parent 5753aed commit 4e3cd84
Show file tree
Hide file tree
Showing 20 changed files with 2,379 additions and 10,034 deletions.
12 changes: 12 additions & 0 deletions Makefile
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
# MATLAB Makefile

FILE = simulatedMovingBed

all: run clean

run:
matlab -nodesktop -nosplash -r "cd ..; installCADET; cd $(FILE); $(FILE) quit"

clean:
rm *.m~
@echo "all cleaned up"
81 changes: 46 additions & 35 deletions OptAlgorithms.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,13 +6,12 @@
% =============================================================================

properties (Constant)
FUNC = @SMB.simulatedMovingBed; % the objective function
swarm = 20; % the population amount
sample = 200; % the loop count for evolution
dataPoint = 1000; % the amount of data observation
DE_strategy = 5; % the strategy of DE kernel
prior = [];
% prior = load('prior.dat'); % load the prior distribution if possible
prior = []; % default; no prior information
FUNC = @SMB.simulatedMovingBed; % the objective function
end


Expand Down Expand Up @@ -89,7 +88,7 @@
fprintf('****************************************************** \n');
fprintf('Time %5.3g elapsed after %5d Iteration \n', result.optTime, result.Iteration);
fprintf('********* Objective value: %10.3g ********** \n', yValue);
fprintf('%g |', xValue); fprintf('\n');
fprintf(' %g |', xValue); fprintf('\n');

end

Expand Down Expand Up @@ -361,7 +360,7 @@
(ParSwarm, OptSwarm, ToplOptSwarm, k, opt);

% Abstract best information so far from the population and display it
xValue = exp( OptSwarm(opt.swarmSize+1, 1:opt.particleSize) );
xValue = exp(OptSwarm(opt.swarmSize+1, 1:opt.particleSize));
yValue = OptSwarm(opt.swarmSize+1, opt.particleSize+1);

fprintf('Iter = %5d ---------------- Minimum: %10.3g ---------------- \n', k, yValue);
Expand Down Expand Up @@ -395,7 +394,7 @@
fprintf('****************************************************** \n');
fprintf('Time %5.3g elapsed after %5d Iteration \n', result.optTime, result.Iteration);
fprintf('********* Objective value: %10.3g ********** \n', yValue);
fprintf('%g |', xValue); fprintf('\n');
fprintf(' %g |', xValue); fprintf('\n');

end

Expand Down Expand Up @@ -768,7 +767,7 @@

% If the Jacobian matrix cannot be obtained,
% a set of samples is used to generate R
[R, oldpar] = OptAlgorithms.burnInSamples(opt);
[R, oldpar, SS] = OptAlgorithms.burnInSamples(opt);

end

Expand Down Expand Up @@ -799,8 +798,8 @@
newSS = feval( OptAlgorithms.FUNC, exp(newpar) );

% The Metropolis probability
rho12 = exp( -0.5 * (newSS - SS) / sigmaSqu) *...
( OptAlgorithms.priorPDF(newpar) / OptAlgorithms.priorPDF(oldpar) );
rho12 = exp( -0.5 *( (newSS - SS) / sigmaSqu + ...
OptAlgorithms.priorPDF(newpar) - OptAlgorithms.priorPDF(oldpar) ) );

% The new proposal is accepted with Metropolis probability
if rand <= min(1, rho12)
Expand All @@ -827,24 +826,24 @@
newSS2 = feval( OptAlgorithms.FUNC, exp(newpar2) );

% The conventional version of calculation
rho32 = exp( -0.5 * (newSS - newSS2) / sigmaSqu) * ...
( OptAlgorithms.priorPDF(newpar) / OptAlgorithms.priorPDF(newpar2) );
rho32 = exp( -0.5 *( (newSS - newSS2) / sigmaSqu + ...
OptAlgorithms.priorPDF(newpar) - OptAlgorithms.priorPDF(newpar2) ) );

% q2 = exp( -0.5 * (newSS2 - SS) / sigmaSqu ) * ...
% ( OptAlgorithms.priorPDF(newpar2) / OptAlgorithms.priorPDF(oldpar) );
% q1 = exp( -0.5 * (norm((newpar2 - newpar) * inv(R))^2 - norm((oldpar - newpar) * inv(R))^2) );
% q2 = exp( -0.5 *( (newSS2 - SS) / sigmaSqu + ...
% OptAlgorithms.priorPDF(newpar2) - OptAlgorithms.priorPDF(oldpar) ) );
% q1 = exp( -0.5 * (norm((newpar2 - newpar) * inv(R))^2 - norm((oldpar - newpar) * inv(R))^2) );

% if rho32 == Inf
% rho13 = 0;
% else
% rho13 = q1 * q2 * (1-rho32) / (1-rho12);
% end
% if rho32 == Inf
% rho13 = 0;
% else
% rho13 = q1 * q2 * (1-rho32) / (1-rho12);
% end

% The speed-up version of above calculation
q1q2 = exp( -0.5 * ( (newSS2 - SS) / sigmaSqu + ...
( (newpar2 - newpar) * (R \ (R' \ (newpar2' - newpar')))...
- (oldpar - newpar) * (R \ (R' \ (oldpar' - newpar'))) )) ) * ...
( OptAlgorithms.priorPDF(newpar2) / OptAlgorithms.priorPDF(oldpar) );
q1q2 = exp( -0.5 *( (newSS2 - SS) / sigmaSqu + ...
( (newpar2 - newpar) * (R \ (R' \ (newpar2' - newpar'))) - ...
(oldpar - newpar) * (R \ (R' \ (oldpar' - newpar'))) ) + ...
OptAlgorithms.priorPDF(newpar2) - OptAlgorithms.priorPDF(oldpar) ) );

rho13 = q1q2 * (1 - rho32) / (1 - rho12);

Expand Down Expand Up @@ -920,7 +919,7 @@
OptAlgorithms.FigurePlot(Population, opt);

[yValue, row] = min(Population(:, opt.Nparams+1));
xValue = exp( Population(row, 1:opt.Nparams) );
xValue = exp(Population(row, 1:opt.Nparams));

% Gather some useful information and store them
result.optTime = etime(clock, startTime) / 3600;
Expand Down Expand Up @@ -975,7 +974,7 @@
opt.criterionTol = 0.0001;
opt.burn_in = 0;
opt.convergInt = 100;
opt.rejectValue = 10000;
opt.rejectValue = 1e5;
opt.nDataPoint = OptAlgorithms.dataPoint;

opt.Jacobian = false; % set it false when you do not need Jacobian matrix
Expand Down Expand Up @@ -1020,7 +1019,7 @@

end

function [R, oldpar] = burnInSamples(opt)
function [R, oldpar, SS] = burnInSamples(opt)
%------------------------------------------------------------------------------
% The routine that is used for generating samples in cases that Jocabian matrix
% is not available
Expand All @@ -1045,10 +1044,15 @@
ParSwarm(nullRow, :) = [];
[~, minRow] = min(ParSwarm(:, opt.Nparams+1));
oldpar = ParSwarm(minRow, 1:opt.Nparams);
SS = ParSwarm(minRow, opt.Nparams+1);

% Calculate the covariance matrix
[chaincov, ~, ~] = OptAlgorithms.covUpdate(ParSwarm(:, 1:opt.Nparams), 1, [], [], []);

if isempty(chaincov)
error('OptAlgorithms.burnInSamples:\n %g randomly generated samples are all rejected with opt.rejectValue = %g \n Please change your rejection value in OptAlgorithms.getOptions_MCMC function', Swarm, opt.rejectValue);
end

% Cholesky decomposition
R = chol( chaincov + eye(opt.Nparams) * 1e-7 );

Expand Down Expand Up @@ -1319,7 +1323,7 @@
% Abstract best information so far from the population and display it
[minValue, row] = min(states(1:opt.Nchain, opt.Nparams+1));
fprintf('---------------- Minimum: %3g ---------------- \n', minValue);
fprintf('%10.3g | ', exp( states(row, 1:opt.Nparams) )); fprintf('\n');
fprintf('%10.3g | ', exp(states(row, 1:opt.Nparams)) ); fprintf('\n');
end

% PRML: Evolution of the chains
Expand Down Expand Up @@ -1384,7 +1388,7 @@
OptAlgorithms.FigurePlot(Population, opt)

[yValue, row] = min(Population(:,opt.Nparams+1));
xValue = exp( Population(row,1:opt.Nparams) );
xValue = exp(Population(row,1:opt.Nparams));

% Gather some useful information and store them
result.optTime = etime(clock, startTime)/3600;
Expand Down Expand Up @@ -1558,8 +1562,9 @@
SS = states(j,opt.Nparams+1);

% The Metropolis probability
rho = ( exp( -0.5*(newSS - SS) / sigmaSqu(j)) )^Beta(j) * ...
( OptAlgorithms.priorPDF(proposal) / OptAlgorithms.priorPDF(states(j,1:opt.Nparams)) );
rho = ( exp( -0.5 *( (newSS - SS) / sigmaSqu(j) + ...
OptAlgorithms.priorPDF(proposal) - ...
OptAlgorithms.priorPDF(states(j, 1:opt.Nparams)) ) ) )^Beta(j);

% If the proposal is accepted
if rand <= min(1, rho)
Expand Down Expand Up @@ -1833,7 +1838,7 @@
OptAlgorithms.FigurePlot(Population, opt);

[yValue, row] = min(Population(:,opt.Nparams+1));
xValue = exp( Population(row,1:opt.Nparams) );
xValue = exp(Population(row,1:opt.Nparams));

% Gather some useful information and store them
result.optTime = etime(clock,startTime) / 3600;
Expand Down Expand Up @@ -1958,7 +1963,7 @@

% Abstract best information so far from the population and display it
fprintf('---------------- Minimum: %g ---------------- \n', OptPopul(opt.Nchain+1, opt.Nparams+1));
fprintf('%10.3g | ', exp( OptPopul(1, 1:opt.Nparams) )); fprintf('\n');
fprintf('%10.3g | ', exp(OptPopul(1, 1:opt.Nparams))); fprintf('\n');

% In each chain, the proposal point is accepted in terms of the Metropolis probability
for j = 1: opt.Nchain
Expand All @@ -1973,8 +1978,8 @@
newSS = feval( OptAlgorithms.FUNC, exp(proposal) );
end

rho = exp( -0.5*(newSS - SS) / sigmaSqu(j)) * ...
( OptAlgorithms.priorPDF(proposal) / OptAlgorithms.priorPDF(states(j, 1:opt.Nparams)) );
rho = exp( -0.5 *( (newSS - SS) / sigmaSqu(j) + ...
OptAlgorithms.priorPDF(proposal) - OptAlgorithms.priorPDF(states(j, 1:opt.Nparams)) ) );

if rand <= min(1, rho)
states(j, 1:opt.Nparams) = proposal;
Expand Down Expand Up @@ -2237,12 +2242,17 @@
idx = floor(0.5 * opt.nsamples);
end

if idx < length(chainData_1)
idx = 0;
end

for k = 1:opt.Nchain
eval(sprintf('chainData_%d(1:idx, :) = [];', k));
Population = [Population; eval(sprintf('chainData_%d', k))];
end

save('population.dat', 'Population', '-ascii', '-append');

end

end % MADE
Expand Down Expand Up @@ -2323,7 +2333,7 @@ function FigurePlot(Population, opt)
return;
end

binNum = 20;
binNum = 15;
x = zeros(binNum, 3);
[r, c] = size(OptAlgorithms.prior);

Expand Down Expand Up @@ -2365,6 +2375,7 @@ function checkOptDimension(opt, LEN)
if length(opt.paramBound) ~= LEN
error('The setup of the initial boundary in SMBOptimiztion is incorrent \n');
end

end

function tickLabelFormat(hAxes, axName, format)
Expand Down
Loading

0 comments on commit 4e3cd84

Please sign in to comment.