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
arGradAndHessian.m
% calculate gradients and hessian from sensitivities
function arGradAndHessian
global ar
ar.fit.beta = nan(1, length(ar.p));
ar.fit.alpha = nan(length(ar.p));
ar.fit.beta(ar.fit.qFit==1) = - (ar.fit.res * ar.fit.sres); % (15.5.6)
ar.fit.gradient = -2 * (ar.fit.beta); % (15.5.8)
ar.fit.alpha(ar.fit.qFit==1,ar.fit.qFit==1) = ar.fit.sres' * ar.fit.sres; % (15.5.11)
ar.fit.hessian = 2 * (ar.fit.alpha); % (15.5.8)
% repress warning "Matrix is close to singular or badly scaled."
warn_reset = warning;
warning('off', 'all');
ar.fit.covar = nan(length(ar.p));
ar.fit.covar(ar.fit.qFit==1,ar.fit.qFit==1) = inv(ar.fit.alpha(ar.fit.qFit==1,ar.fit.qFit==1)); % (15.5.15)
% reset warnings
warning(warn_reset);
% eigenvectors and eigenvalues of hessian
ar.fit.hessian_eval = nan(1, length(ar.p));
ar.fit.hessian_evec = nan(length(ar.p));
[evec, eval] = eig(ar.fit.hessian(ar.fit.qFit==1,ar.fit.qFit==1));
ar.fit.hessian_eval(ar.fit.qFit==1) = diag(eval);
ar.fit.hessian_evec(ar.fit.qFit==1,ar.fit.qFit==1) = evec;