From 9504abdf05ad44b014ea5e8dd75f1baf44da5ab1 Mon Sep 17 00:00:00 2001 From: noblec04 Date: Wed, 28 Aug 2024 19:43:51 -0700 Subject: [PATCH] added LOO MM --- MatlabGP/+utils/catunique.m | 11 ++ MatlabGP/GP.m | 4 +- .../examples/TestAppliedAviationProblem.asv | 53 ++--- .../examples/TestAppliedAviationProblem.m | 36 ++-- MatlabGP/examples/TestRosenbrockProblem.asv | 183 ++++++++++++++++++ MatlabGP/examples/TestRosenbrockProblem.m | 105 ++++++---- 6 files changed, 318 insertions(+), 74 deletions(-) create mode 100644 MatlabGP/+utils/catunique.m create mode 100644 MatlabGP/examples/TestRosenbrockProblem.asv diff --git a/MatlabGP/+utils/catunique.m b/MatlabGP/+utils/catunique.m new file mode 100644 index 0000000..6a51cc6 --- /dev/null +++ b/MatlabGP/+utils/catunique.m @@ -0,0 +1,11 @@ +function [x,flag] = catunique(x,xn) + +reps = ismembertol(xn,x,1e-6,'ByRows',true); + +if ~reps + x = [x;xn]; +end + +flag = ~reps; + +end \ No newline at end of file diff --git a/MatlabGP/GP.m b/MatlabGP/GP.m index b69362f..79b391a 100644 --- a/MatlabGP/GP.m +++ b/MatlabGP/GP.m @@ -141,7 +141,7 @@ res = obj.Y - obj.mean.eval(obj.X); - kkp = pinv(obj.K,0); + kkp = pinv(obj.K); sigp = sqrt(abs(res'*kkp*res./(size(obj.Y,1)))); @@ -151,7 +151,7 @@ obj.K = obj.K + diag(0*xx(:,1)+obj.kernel.signn); - obj.Kinv = pinv(obj.K,0); + obj.Kinv = pinv(obj.K); obj.alpha = obj.Kinv*(res); diff --git a/MatlabGP/examples/TestAppliedAviationProblem.asv b/MatlabGP/examples/TestAppliedAviationProblem.asv index 7516331..cf29d08 100644 --- a/MatlabGP/examples/TestAppliedAviationProblem.asv +++ b/MatlabGP/examples/TestAppliedAviationProblem.asv @@ -1,5 +1,6 @@ -clear all +clear +close all clc lb = [0.5 0.5 2.5*10^(-3)]; @@ -21,11 +22,11 @@ y{1} = y1; y{2} = y2; %% -% ma = means.linear([1 1 1 1]); -% mb = means.linear([1 1 1]); +ma = means.linear([1 1 1 1]); +mb = means.linear([1 1 1]); -ma = means.zero([1 1 1 1]); -mb = means.zero(); +% ma = means.zero(); +% mb = means.zero(); a = kernels.RQ(2,1,[0.1 0.2 0.1 0.3]);%.periodic(1,10); b = kernels.RQ(2,1,[0.2 0.2 0.2]); @@ -79,13 +80,13 @@ toc %% -% mc = means.linear(ones(1,3));%*means.sine(1,10,0,1); -% c = kernels.RQ(2,1,ones(1,3)); -% c.signn = eps; -% -% LOOZ = GP(mc,c); -% LOOZ = LOOZ.condition(x{1},MF.LOO,lb,ub); -% LOOZ = LOOZ.train(); +mc = means.linear(ones(1,3));%*means.sine(1,10,0,1); +c = kernels.RQ(2,1,ones(1,3)); +c.signn = 0.001; + +LOOZ = GP(mc,c); +LOOZ = LOOZ.condition(x{1},log(abs(MF.LOO)),lb,ub); +LOOZ = LOOZ.train(); %% @@ -115,8 +116,8 @@ C = [50 1]; %% for jj = 1:60 - [xn,Rn] = BO.argmax(@BO.MFSFDelta,MF); - %[xn,Rn] = BO.argmax(@BO.UCB,LOOZ); + %[xn,Rn] = BO.argmax(@BO.MFSFDelta,MF); + [xn,Rn] = BO.argmax(@BO.UCB,LOOZ); sign(1) = Z{1}.eval_var(xn)/C(1); sign(2) = Z{2}.eval_var(xn)/C(2); @@ -124,14 +125,16 @@ for jj = 1:60 [~,in] = max(sign); if in==1 - x{1} = [x{1}; xn]; + [x{1},flag] = utils.catunique(x{1},xn); + if flag + y{1} = [y{1}; testFuncs.StressedPlate(xn,1)]; + end end - x{2} = [x{2}; xn]; - if in==1 - y{1} = [y{1}; testFuncs.StressedPlate(xn,1)]; + [x{2},flag] = utils.catunique(x{2},xn); + if flag + y{2} = [y{2}; testFuncs.StressedPlate(xn,2)]; end - y{2} = [y{2}; testFuncs.StressedPlate(xn,2)]; for ii = 1:2 Z{ii} = Z{ii}.condition(x{ii},y{ii},lb,ub); @@ -140,7 +143,7 @@ for jj = 1:60 MF.GPs = Z; MF = MF.condition(); - %LOOZ = LOOZ.condition(x{1},MF.LOO,lb,ub); + LOOZ = LOOZ.condition(x{1},log(abs(MF.LOO)),lb,ub); pc(jj,1) = size(x{1},1); pc(jj,2) = size(x{2},1); @@ -151,17 +154,19 @@ for jj = 1:60 R2MF(jj) = 1 - mean((yy - MF.eval_mu(xx)).^2)./var(yy); RMAEMF(jj) = max(abs(yy - MF.eval_mu(xx)))./std(yy); + cost(jj) = C(1)*pc(jj,1)+pc(jj,2); + figure(3) clf(3) hold on - plot(pc(1:jj,1),R2z) - plot(pc(1:jj,1),R2MF) + plot(cost,R2z) + plot(cost,R2MF) figure(4) clf(4) hold on - plot(pc(1:jj,1),RMAEz) - plot(pc(1:jj,1),RMAEMF) + plot(cost,RMAEz) + plot(cost,RMAEMF) drawnow diff --git a/MatlabGP/examples/TestAppliedAviationProblem.m b/MatlabGP/examples/TestAppliedAviationProblem.m index 576917a..d839bd6 100644 --- a/MatlabGP/examples/TestAppliedAviationProblem.m +++ b/MatlabGP/examples/TestAppliedAviationProblem.m @@ -1,5 +1,6 @@ -clear all +clear +close all clc lb = [0.5 0.5 2.5*10^(-3)]; @@ -30,7 +31,7 @@ a = kernels.RQ(2,1,[0.1 0.2 0.1 0.3]);%.periodic(1,10); b = kernels.RQ(2,1,[0.2 0.2 0.2]); a.signn = eps; -b.signn = 0.1; +b.signn = 0.01; %% tic @@ -81,10 +82,10 @@ % mc = means.linear(ones(1,3));%*means.sine(1,10,0,1); % c = kernels.RQ(2,1,ones(1,3)); -% c.signn = eps; +% c.signn = 1*10^(-7); % % LOOZ = GP(mc,c); -% LOOZ = LOOZ.condition(x{1},MF.LOO,lb,ub); +% LOOZ = LOOZ.condition(x{1},log(abs(MF.LOO)),lb,ub); % LOOZ = LOOZ.train(); @@ -113,10 +114,11 @@ C = [50 1]; %% -for jj = 1:60 +for jj = 1:100 [xn,Rn] = BO.argmax(@BO.MFSFDelta,MF); %[xn,Rn] = BO.argmax(@BO.UCB,LOOZ); + %[xn,Rn] = BO.argmax(@BO.maxVAR,MF); sign(1) = Z{1}.eval_var(xn)/C(1); sign(2) = Z{2}.eval_var(xn)/C(2); @@ -124,14 +126,16 @@ [~,in] = max(sign); if in==1 - x{1} = [x{1}; xn]; + [x{1},flag] = utils.catunique(x{1},xn); + if flag + y{1} = [y{1}; testFuncs.StressedPlate(xn,1)]; + end end - x{2} = [x{2}; xn]; - if in==1 - y{1} = [y{1}; testFuncs.StressedPlate(xn,1)]; + [x{2},flag] = utils.catunique(x{2},xn); + if flag + y{2} = [y{2}; testFuncs.StressedPlate(xn,2)]; end - y{2} = [y{2}; testFuncs.StressedPlate(xn,2)]; for ii = 1:2 Z{ii} = Z{ii}.condition(x{ii},y{ii},lb,ub); @@ -140,7 +144,7 @@ MF.GPs = Z; MF = MF.condition(); - %LOOZ = LOOZ.condition(x{1},MF.LOO,lb,ub); + %LOOZ = LOOZ.condition(x{1},log(abs(MF.LOO)),lb,ub); pc(jj,1) = size(x{1},1); pc(jj,2) = size(x{2},1); @@ -151,17 +155,19 @@ R2MF(jj) = 1 - mean((yy - MF.eval_mu(xx)).^2)./var(yy); RMAEMF(jj) = max(abs(yy - MF.eval_mu(xx)))./std(yy); + cost(jj) = C(1)*pc(jj,1)+pc(jj,2); + figure(3) clf(3) hold on - plot(C(1)*pc(1:jj,1)+pc(1:jj,2),R2z) - plot(C(1)*pc(1:jj,1)+pc(1:jj,2),R2MF) + plot(cost,R2z) + plot(cost,R2MF) figure(4) clf(4) hold on - plot(C(1)*pc(1:jj,1)+pc(1:jj,2),RMAEz) - plot(C(1)*pc(1:jj,1)+pc(1:jj,2),RMAEMF) + plot(cost,RMAEz) + plot(cost,RMAEMF) drawnow diff --git a/MatlabGP/examples/TestRosenbrockProblem.asv b/MatlabGP/examples/TestRosenbrockProblem.asv new file mode 100644 index 0000000..86529c2 --- /dev/null +++ b/MatlabGP/examples/TestRosenbrockProblem.asv @@ -0,0 +1,183 @@ + +clear all +close all +clc + +D = 3; +nF = 3; + +lb = -2*ones(1,D); +ub = 2*ones(1,D); + +xx = lb + (ub - lb).*lhsdesign(50000,D); +yy = testFuncs.Rosenbrock(xx,1); + +x1 = lb + (ub - lb).*lhsdesign(5,D); +y1 = testFuncs.Rosenbrock(x1,1); + +x2 = [lb + (ub - lb).*lhsdesign(20,D)]; +y2 = testFuncs.Rosenbrock(x2,2); + +x3 = [lb + (ub - lb).*lhsdesign(100,D)]; +y3 = testFuncs.Rosenbrock(x3,3); + +x{1} = x1; +x{2} = x2; +x{3} = x3; + +y{1} = y1; +y{2} = y2; +y{3} = y3; + +%% +ma = means.const(1); +mb = means.linear(ones(1,D)); + +a = kernels.RQ(2,1,ones(1,D+nF-1)); +b = kernels.RQ(2,1,ones(1,D)); +a.signn = 1*10^(-7); +b.signn = 1*10^(-7); + +%% +tic +for i = 1:nF + Z{i} = GP(mb,b); + Z{i} = Z{i}.condition(x{i},y{i},lb,ub); + Z{i} = Z{i}.train(); +end +toc + +%% +tic +MF = NLMFGP(Z,ma,a); +MF = MF.condition(); +MF = MF.train(); +toc + +mc = means.linear(ones(1,D));%*means.sine(1,10,0,1); +c = kernels.RQ(2,1,ones(1,D)); +c.signn = eps; + +LOO = GP(mc,c); + +LOOZ{1} = LOO.condition(x{1},log(abs(Z{1}.LOO)),lb,ub); +LOOZ{1} = LOOZ{1}.train(); +LOOZ{2} = LOO.condition(x{2},log(abs(Z{2}.LOO)),lb,ub); +LOOZ{2} = LOOZ{2}.train(); + +LOOMF = LOO.condition(x{1},log(abs(MF.LOO)),lb,ub); +LOOMF = LOOMF.train(); + + +%% +figure +hold on +utils.plotSurf(Z{1},1,2,'color','r') +utils.plotSurf(Z{2},1,2,'color','b') +%utils.plotSurf(Z{3},1,2,'color','g') + +%% + +figure +hold on +utils.plotSurf(MF,1,2) + +%% + +1 - mean((yy - Z{1}.eval_mu(xx)).^2)./var(yy) +max(abs(yy - Z{1}.eval_mu(xx)))./std(yy) + +1 - mean((yy - MF.eval_mu(xx)).^2)./var(yy) +max(abs(yy - MF.eval_mu(xx)))./std(yy) + +%% + +C = [50 20 1];%20 + +for jj = 1:60 + + %[xn,Rn] = BO.argmax(@BO.UCB,LOOMF); + %[xn,Rn] = BO.argmax(@BO.MFSFDelta,MF); + [xn,Rn] = BO.argmax(@BO.maxVAR,MF); + + siggn(1) = exp((LOOZ{1}.eval(xn)))/(C(1)); + siggn(2) = exp((LOOZ{2}.eval(xn)))/(C(2)); + + % siggn(1) = abs(Z{1}.eval_var(xn))/(C(1)); + % siggn(2) = abs(Z{2}.eval_var(xn))/(C(2)); + %siggn(3) = abs(Z{3}.eval_var(xn))/(C(3)); + + [~,in] = max(siggn); + + if in==1 + [x{1},flag] = utils.catunique(x{1},xn); + if flag + y{1} = [y{1}; testFuncs.Rosenbrock(xn,1)]; + end + end + + if in<=2 + [x{2},flag] = utils.catunique(x{2},xn); + if flag + y{2} = [y{2}; testFuncs.Rosenbrock(xn,2)]; + end + end + + % [x{3},flag] = utils.catunique(x{3},xn); + % if flag + % y{3} = [y{3}; testFuncs.Rosenbrock(xn,3)]; + % end + + for ii = 1:nF + Z{ii} = Z{ii}.condition(x{ii},y{ii},lb,ub); + end + + MF.GPs = Z; + MF = MF.condition(); + + LOOZ{1} = LOOZ{1}.condition(x{1},log(abs(Z{1}.LOO)),lb,ub); + LOOZ{2} = LOOZ{2}.condition(x{2},log(abs(Z{2}.LOO)),lb,ub); + + %LOOMF = LOOMF.condition(x{1},log(abs(MF.LOO)),lb,ub); + + R2z(jj) = 1 - mean((yy - Z{1}.eval_mu(xx)).^2)./var(yy); + RMAEz(jj) = max(abs(yy - Z{1}.eval_mu(xx)))./std(yy); + + R2MF(jj) = 1 - mean((yy - MF.eval_mu(xx)).^2)./var(yy); + RMAEMF(jj) = max(abs(yy - MF.eval_mu(xx)))./std(yy); + + cost(jj) = C(1)*size(x{1},1)+C(2)*size(x{2},1);%+C(3)*size(x{3},1); + + figure(3) + clf(3) + hold on + plot(cost,R2z) + plot(cost,R2MF) + + figure(4) + clf(4) + hold on + plot(cost,RMAEz) + plot(cost,RMAEMF) + + drawnow + + if RMAEMF(jj)<0.05 + break + end +end + + +%% + +% xi = lb + (ub - lb).*lhsdesign(size(x{1},1),D); +% yi = testFuncs.Rosenbrock(xn,1); +% +% Zi = GP(mb,b); +% Zi = Zi.condition(xi,yi,lb,ub); +% Zi = Zi.train(); +% +% %% +% +% 1 - mean((yy - Zi.eval_mu(xx)).^2)./var(yy) +% max(abs(yy - Zi.eval_mu(xx)))./std(yy) \ No newline at end of file diff --git a/MatlabGP/examples/TestRosenbrockProblem.m b/MatlabGP/examples/TestRosenbrockProblem.m index 49fbe4a..db04730 100644 --- a/MatlabGP/examples/TestRosenbrockProblem.m +++ b/MatlabGP/examples/TestRosenbrockProblem.m @@ -1,8 +1,10 @@ clear all +close all clc D = 3; +nF = 3; lb = -2*ones(1,D); ub = 2*ones(1,D); @@ -28,17 +30,17 @@ y{3} = y3; %% -ma = means.const([1]); +ma = means.const(1); mb = means.linear(ones(1,D)); -a = kernels.RQ(2,1,ones(1,D+2)); +a = kernels.RQ(2,1,ones(1,D+nF-1)); b = kernels.RQ(2,1,ones(1,D)); a.signn = eps; -b.signn = eps; +b.signn = 1*10^(-7); %% tic -for i = 1:3 +for i = 1:nF Z{i} = GP(mb,b); Z{i} = Z{i}.condition(x{i},y{i},lb,ub); Z{i} = Z{i}.train(); @@ -54,11 +56,19 @@ mc = means.linear(ones(1,D));%*means.sine(1,10,0,1); c = kernels.RQ(2,1,ones(1,D)); -c.signn = eps; +c.signn = 0.001; -LOOZ = GP(mc,c); -LOOZ = LOOZ.condition(x{1},MF.LOO,lb,ub); -LOOZ = LOOZ.train(); +LOO = GP(mc,c); + +LOOZ{1} = LOO.condition(x{1},log(abs(Z{1}.LOO)),lb,ub); +LOOZ{1} = LOOZ{1}.train(); +LOOZ{2} = LOO.condition(x{2},log(abs(Z{2}.LOO)),lb,ub); +LOOZ{2} = LOOZ{2}.train(); +LOOZ{3} = LOO.condition(x{3},log(abs(Z{3}.LOO)),lb,ub); +LOOZ{3} = LOOZ{3}.train(); + +LOOMF = LOO.condition(x{1},log(abs(MF.LOO)),lb,ub); +LOOMF = LOOMF.train(); %% @@ -83,29 +93,56 @@ max(abs(yy - MF.eval_mu(xx)))./std(yy) %% + +C = [50 20 1];%20 + for jj = 1:60 - [xn,Rn] = BO.argmax(@BO.UCB,LOOZ); + %[xn,Rn] = BO.argmax(@BO.UCB,LOOMF); %[xn,Rn] = BO.argmax(@BO.MFSFDelta,MF); - %[xn,Rn] = BO.argmax(@BO.maxVAR,MF); + [xn,Rn] = BO.argmax(@BO.maxVAR,MF); + siggn(1) = exp((LOOZ{1}.eval(xn)))/(C(1)); + siggn(2) = exp((LOOZ{2}.eval(xn)))/(C(2)); + siggn(3) = exp((LOOZ{3}.eval(xn)))/(C(3)); + + % siggn(1) = abs(Z{1}.eval_var(xn))/(C(1)); + % siggn(2) = abs(Z{2}.eval_var(xn))/(C(2)); + %siggn(3) = abs(Z{3}.eval_var(xn))/(C(3)); + + [~,in] = max(siggn); + + if in==1 + [x{1},flag] = utils.catunique(x{1},xn); + if flag + y{1} = [y{1}; testFuncs.Rosenbrock(xn,1)]; + end + end - x{1} = [x{1}; xn]; - x{2} = [x{2}; xn]; - x{3} = [x{3}; xn]; + if in<=2 + [x{2},flag] = utils.catunique(x{2},xn); + if flag + y{2} = [y{2}; testFuncs.Rosenbrock(xn,2)]; + end + end - y{1} = [y{1}; testFuncs.Rosenbrock(xn,1)]; - y{2} = [y{2}; testFuncs.Rosenbrock(xn,2)]; - y{3} = [y{3}; testFuncs.Rosenbrock(xn,3)]; + [x{3},flag] = utils.catunique(x{3},xn); + if flag + y{3} = [y{3}; testFuncs.Rosenbrock(xn,3)]; + end - for ii = 1:3 + for ii = 1:nF Z{ii} = Z{ii}.condition(x{ii},y{ii},lb,ub); end MF.GPs = Z; MF = MF.condition(); - LOOZ = LOOZ.condition(x{1},MF.LOO,lb,ub); + LOOZ{1} = LOOZ{1}.condition(x{1},log(abs(Z{1}.LOO)),lb,ub); + LOOZ{2} = LOOZ{2}.condition(x{2},log(abs(Z{2}.LOO)),lb,ub); + LOOZ{3} = LOOZ{3}.condition(x{3},log(abs(Z{3}.LOO)),lb,ub); + + %LOOMF = LOOMF.condition(x{1},log(abs(MF.LOO)),lb,ub); R2z(jj) = 1 - mean((yy - Z{1}.eval_mu(xx)).^2)./var(yy); RMAEz(jj) = max(abs(yy - Z{1}.eval_mu(xx)))./std(yy); @@ -113,17 +150,19 @@ R2MF(jj) = 1 - mean((yy - MF.eval_mu(xx)).^2)./var(yy); RMAEMF(jj) = max(abs(yy - MF.eval_mu(xx)))./std(yy); + cost(jj) = C(1)*size(x{1},1)+C(2)*size(x{2},1);%+C(3)*size(x{3},1); + figure(3) clf(3) hold on - plot(1:jj,R2z) - plot(1:jj,R2MF) + plot(cost,R2z) + plot(cost,R2MF) figure(4) clf(4) hold on - plot(1:jj,RMAEz) - plot(1:jj,RMAEMF) + plot(cost,RMAEz) + plot(cost,RMAEMF) drawnow @@ -135,14 +174,14 @@ %% -xi = lb + (ub - lb).*lhsdesign(size(x{1},1),D); -yi = testFuncs.Rosenbrock(xn,1); - -Zi = GP(mb,b); -Zi = Zi.condition(xi,yi,lb,ub); -Zi = Zi.train(); - -%% - -1 - mean((yy - Zi.eval_mu(xx)).^2)./var(yy) -max(abs(yy - Zi.eval_mu(xx)))./std(yy) \ No newline at end of file +% xi = lb + (ub - lb).*lhsdesign(size(x{1},1),D); +% yi = testFuncs.Rosenbrock(xn,1); +% +% Zi = GP(mb,b); +% Zi = Zi.condition(xi,yi,lb,ub); +% Zi = Zi.train(); +% +% %% +% +% 1 - mean((yy - Zi.eval_mu(xx)).^2)./var(yy) +% max(abs(yy - Zi.eval_mu(xx)))./std(yy) \ No newline at end of file