From 3303b0cc22868795b229936d3acbdbc048d3f082 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Fri, 14 Apr 2023 09:14:17 -0500 Subject: [PATCH 1/2] add mac support --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 394f773..f748bfa 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ *.mat *.png *.fig +.DS_Store From 5760f8bcd3f3e183eaf18340e1e3863fd5bfd632 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Fri, 14 Apr 2023 09:14:23 -0500 Subject: [PATCH 2/2] update EnGMF --- src/+datools/+statistical/+ensemble/EnGMF.m | 219 +++++++++++++++++--- 1 file changed, 191 insertions(+), 28 deletions(-) diff --git a/src/+datools/+statistical/+ensemble/EnGMF.m b/src/+datools/+statistical/+ensemble/EnGMF.m index 3d810fe..a43ed28 100644 --- a/src/+datools/+statistical/+ensemble/EnGMF.m +++ b/src/+datools/+statistical/+ensemble/EnGMF.m @@ -1,55 +1,218 @@ classdef EnGMF < datools.statistical.ensemble.EnF + properties + BandwidthScale + BRUFSteps + UseRobustSampling + SamplingType + UseELocalization + UseAdaptiveEnsembleSize = false + RadiusScale + end + + properties (Access = private) + QMCStream + end + + methods + function obj = EnGMF(varargin) + + p = inputParser; + p.KeepUnmatched = true; + + % this is recommended according to the Reich book + + addParameter(p, 'BandwidthScale', 1); + addParameter(p, 'BRUFSteps', 1); + addParameter(p, 'UseRobustSampling', false); + addParameter(p, 'SamplingType', "Gaussian"); + addParameter(p, 'UseELocalization', false); + addParameter(p, 'RadiusScale', 1); + + parse(p, varargin{2:end}); + + s = p.Results; + + kept = p.Unmatched; + + obj@datools.statistical.ensemble.EnF(varargin{1}, kept); + + obj.BandwidthScale = s.BandwidthScale; + obj.BRUFSteps = s.BRUFSteps; + obj.UseRobustSampling = s.UseRobustSampling; + obj.SamplingType = s.SamplingType; + obj.UseELocalization = s.UseELocalization; + obj.RadiusScale = s.RadiusScale; + + end + end + methods function analysis(obj, R, y) tc = obj.Model.TimeSpan(1); - xf = obj.Ensemble; + Xb = obj.Ensemble; ensN = obj.NumEnsemble; - n = size(xf, 1); + n = size(Xb, 1); - beta = ((4/(n + 2))^(2/(n + 4)))*((ensN)^(-2/(n + 4))); + if obj.SamplingType == "QMC" && isempty(obj.QMCStream) + obj.QMCStream = qrandstream('sobol', n + 1, 'Skip', 100); + end - xfm = mean(xf, 2); - Af = xf - repmat(xfm, 1, ensN); - Af = sqrt(beta) * Af/sqrt(ensN - 1); + beta = obj.BandwidthScale*((4/(n + 2))^(2/(n + 4)))*((ensN)^(-2/(n + 4))); + + if obj.UseELocalization + Xbm = mean(Xb, 2); + Ab = Xb - repmat(Xbm, 1, ensN); + Ab = Ab/sqrt(ensN - 1); + + radiusscale = obj.RadiusScale; + epsilon = 1e-4; + % calculate the distances + dist = zeros(ensN, ensN); + wD = zeros(ensN, ensN); + for i = 1:ensN + for j = 1:ensN + dist(i, j) = norm(Xb(:, i) - Xb(:, j)); + end + [~, I] = sort(dist(i, :)); + r = radiusscale*dist(i, I(round(sqrt(ensN)))); + + wD(i, :) = exp(-0.5*(dist(i, :)/r).^2) + epsilon; + end + wD = wD./sum(wD, 2); + + Pb = Ab*Ab.'; + trPb = trace(Pb); + + Btilde = zeros(n, n, ensN); + trPbis = zeros(ensN, 1); + + for i = 1:ensN + wDi = wD(i, :).'; + Xbm = Xb*wDi; + c = 1/(1 - sum(wDi.^2)); + Pbi = c*((Xb - Xbm)*(diag(wDi))*(Xb - Xbm).'); + trPbis(i) = trace(Pbi); + + Btilde(:, :, i) = beta*Pbi; + end + + Btilde = (trPb/mean(trPbis))*Btilde; + + else + Xbm = mean(Xb, 2); + Ab = Xb - repmat(Xbm, 1, ensN); + Ab = sqrt(beta)*Ab/sqrt(ensN - 1); + B = Ab*Ab.'; + Btilde = repmat(B, 1, 1, ensN); + end - Hxf = obj.Observation.observeWithoutError(tc, xf); - Hxfm = mean(Hxf, 2); - HAf = Hxf - repmat(Hxfm, 1, ensN); - HAf = sqrt(beta)*HAf/sqrt(ensN - 1); + Xtilde = Xb; + as = zeros(ensN, 1); + w = obj.Weights.'; + + K = obj.BRUFSteps; + for k = 1:K + wold = w; + HXtilde = obj.Observation.observeWithoutError(tc, Xtilde); + for i = 1:ensN + H = obj.Observation.linearization(tc, Xtilde(:, i)); + + Btildei = Btilde(:, :, i); + + S = (H*Btildei*H.' + K*R); + S = (S + S.')/2; + RS = chol(S); + d = (HXtilde(:, i) - y); + Xtilde(:, i) = Xtilde(:, i) - Btildei*H.'*(S\d); + + % recompute Ptilde + Btildei = (Btildei - Btildei*H.'*(S\(H*Btildei))); + Btilde(:, :, i) = (Btildei + Btildei.')/2; + + a = (RS.'\d); + as(i) = -0.5*(a.'*a) - sum(log(diag(RS))); + end + m = max(as); + w = exp(as-(m + log(sum(exp(as-m))))).'; + w = w/sum(w); + w = w.*wold; + w = w/sum(w); + end - dS = decomposition((HAf * HAf.') + R, 'chol'); + % resample + Xa = Xb; - xtilde = xf - (Af*HAf.')*(dS\(Hxf - y)); - Aa = xtilde - mean(xtilde, 2); - Aa = sqrt(beta)*Aa/sqrt(ensN - 1); + if obj.UseAdaptiveEnsembleSize - t0 = (Hxf - y); + else + ensNa = ensN; + end - as = (-0.5 * sum(t0.*(dS \ t0), 1)).'; - m = max(as); - w = exp(as-(m + log(sum(exp(as-m))))).'; - w = w/sum(w); + switch obj.SamplingType + case "QMC" + for k = 1:ensNa - defensivefactor = 1; - w = defensivefactor*w + (1 - defensivefactor)/ensN; + samples = obj.QMCStream.rand(1, n + 1).'; + ind = find(samples(1) <= cumsum(w), 1, 'first'); - % resample - xa = xf; - for k = 1:ensN - ind = find(rand <= cumsum(w), 1, 'first'); - xa(:, k) = xtilde(:, ind) + Aa*randn(ensN, 1); + sqBa = sqrtm(Btilde(:, :, ind)); + + if obj.UseRobustSampling + Xa(:, k) = Xtilde(:, ind) + sqBa*robustnorminv(samples(2:end)); + else + Xa(:, k) = Xtilde(:, ind) + sqBa*norminv(samples(2:end)); + end + + end + + case "Gaussian" + for k = 1:ensNa + + ind = find(rand <= cumsum(w), 1, 'first'); + + sqBa = sqrtm(Btilde(:, :, ind)); + + if obj.UseRobustSampling + Xa(:, k) = Xtilde(:, ind) + sqBa*robustnorminv(rand(n, 1)); + else + Xa(:, k) = Xtilde(:, ind) + sqBa*randn(n, 1); + end + end + case "None" + obj.Ensemble = Xtilde; + + obj.Weights = w.'; + return end - obj.Ensemble = xa; + obj.Ensemble = Xa; - obj.Weights = ones(ensN, 1) / ensN; + obj.Weights = ones(ensNa, 1) / ensNa; end end end + +function x = robustnorminv(x) + +x(x < 0) = 0; +x(x > 1) = 1; + +I1 = x <= 1/6; +I2 = and(x > 1/6, x <= 5/6); +I3 = x > 5/6; + +x(I1) = (-3)+2.*6.^(1/3).*x(I1).^(1/3); +x(I2) = 2.*3.^(1/2).*sin((1/3).*atan(((-2)+4.*x(I2)).*((-1)+(-16).*((-1)+x(I2)).* ... + x(I2)).^(-1/2))); +x(I3) = 3+(-2).*(6+(-6).*x(I3)).^(1/3); + +x = real(x); + +end