Skip to content
Snippets Groups Projects
Commit 73291b0b authored by Marco Ratto's avatar Marco Ratto
Browse files

before issuing F singularity, check with rescaled F matrix: this spares lots...

before issuing F singularity, check with rescaled F matrix: this spares lots of computing time when singularity only happens in the first KF step.
parent 09be021d
Branches
Tags
No related merge requests found
Pipeline #2971 passed
......@@ -158,6 +158,7 @@ else
LIKK={likk,dlikk};
end
rescale_prediction_error_covariance0=rescale_prediction_error_covariance;
while notsteady && t<=last
s = t-start+1;
if Zflag
......@@ -175,7 +176,12 @@ while notsteady && t<=last
end
else
if rcond(F)<kalman_tol
badly_conditioned_F = true;
sig=sqrt(diag(F));
if any(diag(F)<kalman_tol) || rcond(F./(sig*sig'))<kalman_tol
badly_conditioned_F = true;
else
rescale_prediction_error_covariance=1;
end
end
end
if badly_conditioned_F
......@@ -191,6 +197,7 @@ while notsteady && t<=last
if rescale_prediction_error_covariance
log_dF = log(det(F./(sig*sig')))+2*sum(log(sig));
iF = inv(F./(sig*sig'))./(sig*sig');
rescale_prediction_error_covariance=rescale_prediction_error_covariance0;
else
log_dF = log(det(F));
iF = inv(F);
......
......@@ -84,6 +84,7 @@ oldK = Inf;
notsteady = 1;
F_singular = true;
s = 0;
rescale_prediction_error_covariance0=rescale_prediction_error_covariance;
while notsteady && t<=last
s = t-start+1;
......@@ -110,7 +111,13 @@ while notsteady && t<=last
end
else
if rcond(F)<kalman_tol
badly_conditioned_F = true;
sig=sqrt(diag(F));
if any(diag(F)<kalman_tol) || rcond(F./(sig*sig'))<kalman_tol
badly_conditioned_F = true;
else
rescale_prediction_error_covariance=1;
end
% badly_conditioned_F = true;
end
end
if badly_conditioned_F
......@@ -126,6 +133,7 @@ while notsteady && t<=last
if rescale_prediction_error_covariance
log_dF = log(det(F./(sig*sig')))+2*sum(log(sig));
iF = inv(F./(sig*sig'))./(sig*sig');
rescale_prediction_error_covariance=rescale_prediction_error_covariance0;
else
log_dF = log(det(F));
iF = inv(F);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment