https://github.com/Data2Dynamics/d2d
Tip revision: d72d4c223fb8adc0e3e44a1d3d1eafce18b995ae authored by Andreas Raue on 06 March 2015, 00:17:07 UTC
Increase C version code, after Joep's submits today to prevent compatibility problems
Increase C version code, after Joep's submits today to prevent compatibility problems
Tip revision: d72d4c2
arPLErolling.m
function arPLErolling(jk, n)
global ar
if(~exist('jk','var') || isempty(jk))
jk = find(ar.qFit==1);
end
if(length(jk)>1)
for j=1:length(jk)
arPLErolling(jk(j));
end
return;
end
if(~exist('n','var'))
n = 100;
end
if(~isfield(ar, 'ple'))
ar.ple.chi2s = {};
ar.ple.errors = {};
ar.ple.lambdas = {};
ar.ple.ps = {};
ar.ple.run = zeros(size(ar.p));
ar.ple.alpha = 0.05;
ar.ple.ndof = 1;
end
arChi2(true);
chi2Reset = ar.chi2fit;
pReset = ar.p;
ar.ple.pStart = ar.p;
fprintf('PLE #%i for %s...\n', jk, ar.pLabel{jk});
arWaitbar(0);
tic;
[chi2sup, psup, errorsup, lambdasup] = ple_task(jk, n, 1, chi2Reset, pReset);
[chi2sdown, psdown, errorsdown, lambdasdown] = ple_task(jk, n, -1, chi2Reset, pReset);
ar.ple.chi2s{jk} = [fliplr(chi2sdown) chi2Reset chi2sup];
ar.ple.ps{jk} = [flipud(psdown); pReset; psup];
ar.ple.errors{jk} = [fliplr(errorsdown) nan errorsup];
ar.ple.lambdas{jk} = [fliplr(lambdasdown); nan(size(ar.p)); lambdasup];
ar.ple.run(jk) = 1;
arWaitbar(-1);
fprintf('PLE #%i %s elapse time\n', jk, secToHMS(toc));
ar.p = pReset;
arChi2(false);
function [chi2s, ps, errors, lambdas] = ple_task(jk, n, direction, chi2Reset, pReset)
global ar
chi2s = nan(1,n);
ps = nan(n,length(pReset));
lambdas = nan(n,length(pReset));
errors = nan(1,n);
dchi2 = chi2inv(1-ar.ple.alpha, ar.ple.ndof);
dp = 1;
addres = zeros(1,length(pReset));
addres(jk) = -dp*direction;
ar.p = pReset;
pLast = pReset;
quadprog_optims = optimset('Display', 'off');
warning('off','optim:lsqlin:LinConstraints');
warning('off','optim:quadprog:SwitchToMedScale');
for j=1:n
if(direction>0)
arWaitbar(j,2*n, sprintf('PLE up for %s', strrep(ar.pLabel{jk},'_','\_')));
else
arWaitbar(j+n,2*n, sprintf('PLE down for %s', strrep(ar.pLabel{jk},'_','\_')));
end
% deltap = 1;
% while(norm(deltap)>1e-3)
arChi2(true);
qFit = ar.qFit==1;
res = ar.res;
sres = ar.sres;
beta = -ar.res*ar.sres;
beta(jk)
res = [res addres(qFit)]; %#ok<AGROW>
sres = [sres(:,qFit);eye(sum(qFit))]; %#ok<AGROW>
[deltap,~,~,~,~,lambda] = lsqlin(-sres, res, [], [], [], [], ...
ar.lb(qFit) - pLast(qFit), ar.ub(qFit) - pLast(qFit), [], quadprog_optims);
% disp([sum(lambda.lower~=0) sum(lambda.upper~=0)]);
pLast(qFit) = pLast(qFit) + deltap';
ar.p = pLast;
% end
ps(j,:) = ar.p;
chi2s(j) = ar.chi2fit;
lambdas(j,:) = lambda.upper~=0 | lambda.lower~=0;
if(ar.chi2fit - chi2Reset > 2*dchi2)
fprintf('PLE #%i reached confidence limit\n', jk);
break;
end
if(lambda.upper(jk)~=0)
fprintf('PLE #%i reached upper parameter bound\n', jk);
break;
end
if(lambda.lower(jk)~=0)
fprintf('PLE #%i reached lower parameter bound\n', jk);
break;
end
% if(norm(deltap) <= 1e-3)
% dp = dp * 1.1;
% addres(jk) = -dp*direction;
% end
end
ar.p = pReset;
% function [chi2s, ps, errors, lambdas] = ple_task(jk, n, direction, chi2Reset, pReset)
%
% global ar
%
% chi2s = nan(1,n);
% ps = nan(n,length(pReset));
% lambdas = nan(n,length(pReset));
% errors = nan(1,n);
%
% dchi2 = chi2inv(1-ar.ple.alpha, ar.ple.ndof);
%
% dp = 1;
% addres = zeros(1,length(pReset));
% addres(jk) = -dp*direction;
%
% ar.p = pReset;
% pLast = pReset;
%
% quadprog_optims = optimset('Display', 'off');
% warning('off','optim:lsqlin:LinConstraints');
% warning('off','optim:quadprog:SwitchToMedScale');
%
% for j=1:n
% if(direction>0)
% arWaitbar(j,2*n, sprintf('PLE up for %s', strrep(ar.pLabel{jk},'_','\_')));
% else
% arWaitbar(j+n,2*n, sprintf('PLE down for %s', strrep(ar.pLabel{jk},'_','\_')));
% end
%
% % deltap = 1;
% % while(norm(deltap)>1e-3)
% arChi2(true);
%
% qFit = ar.qFit==1;
%
% res = ar.res;
% sres = ar.sres;
%
% beta = -ar.res*ar.sres;
% beta(jk)
%
% res = [res addres(qFit)]; %#ok<AGROW>
% sres = [sres(:,qFit);eye(sum(qFit))]; %#ok<AGROW>
%
% [deltap,~,~,~,~,lambda] = lsqlin(-sres, res, [], [], [], [], ...
% ar.lb(qFit) - pLast(qFit), ar.ub(qFit) - pLast(qFit), [], quadprog_optims);
%
% % disp([sum(lambda.lower~=0) sum(lambda.upper~=0)]);
%
% pLast(qFit) = pLast(qFit) + deltap';
% ar.p = pLast;
% % end
%
% ps(j,:) = ar.p;
% chi2s(j) = ar.chi2fit;
% lambdas(j,:) = lambda.upper~=0 | lambda.lower~=0;
%
% if(ar.chi2fit - chi2Reset > 2*dchi2)
% fprintf('PLE #%i reached confidence limit\n', jk);
% break;
% end
%
% if(lambda.upper(jk)~=0)
% fprintf('PLE #%i reached upper parameter bound\n', jk);
% break;
% end
% if(lambda.lower(jk)~=0)
% fprintf('PLE #%i reached lower parameter bound\n', jk);
% break;
% end
%
% % if(norm(deltap) <= 1e-3)
% % dp = dp * 1.1;
% % addres(jk) = -dp*direction;
% % end
% end
%
% ar.p = pReset;
%