From 4353418dc3a82c637288fd64b8f6b0a713f01081 Mon Sep 17 00:00:00 2001 From: David O'Gara Date: Tue, 7 Oct 2025 11:49:18 -0500 Subject: [PATCH 1/7] change assert == to np.allclose --- tests/unit_tests/test_hetGP.py | 32 ++++++++++++++++---------------- tests/unit_tests/test_homGP.py | 32 ++++++++++++++++---------------- 2 files changed, 32 insertions(+), 32 deletions(-) diff --git a/tests/unit_tests/test_hetGP.py b/tests/unit_tests/test_hetGP.py index 427e714..c538738 100644 --- a/tests/unit_tests/test_hetGP.py +++ b/tests/unit_tests/test_hetGP.py @@ -30,10 +30,10 @@ def test_fit(): test_model = emulator(x=X, theta=np.array([0]), f=Y, method="hetGP") test_model.fit() - assert test_model._info["ll"] == reference_model["ll"] - assert test_model._info["theta"] == reference_model["theta"] - assert test_model._info["g"] == reference_model["g"] - assert test_model._info["beta0"] == reference_model["beta0"] + assert np.allclose(test_model._info["ll"], reference_model["ll"]) + assert np.allclose(test_model._info["theta"], reference_model["theta"]) + assert np.allclose(test_model._info["g"], reference_model["g"]) + assert np.allclose(test_model._info["beta0"], reference_model["beta0"]) assert np.allclose(test_model._info["Delta"], reference_model["Delta"]) assert np.allclose(test_model._info["Lambda"], reference_model["Lambda"]) @@ -56,10 +56,10 @@ def test_matern(): ) test_model.fit() - assert test_model._info["ll"] == reference_model["ll"] - assert test_model._info["theta"] == reference_model["theta"] - assert test_model._info["g"] == reference_model["g"] - assert test_model._info["beta0"] == reference_model["beta0"] + assert np.allclose(test_model._info["ll"], reference_model["ll"]) + assert np.allclose(test_model._info["theta"], reference_model["theta"]) + assert np.allclose(test_model._info["g"], reference_model["g"]) + assert np.allclose(test_model._info["beta0"], reference_model["beta0"]) assert np.allclose(test_model._info["Delta"], reference_model["Delta"]) assert np.allclose(test_model._info["Lambda"], reference_model["Lambda"]) @@ -114,10 +114,10 @@ def test_hetGP_update_kriging_believer(): test_model.fit() test_model.update(Xnew) - assert test_model._info["ll"] == reference_model["ll"] - assert test_model._info["theta"] == reference_model["theta"] - assert test_model._info["g"] == reference_model["g"] - assert test_model._info["beta0"] == reference_model["beta0"] + assert np.allclose(test_model._info["ll"], reference_model["ll"]) + assert np.allclose(test_model._info["theta"], reference_model["theta"]) + assert np.allclose(test_model._info["g"], reference_model["g"]) + assert np.allclose(test_model._info["beta0"], reference_model["beta0"]) def test_hetGP_update(): @@ -132,10 +132,10 @@ def test_hetGP_update(): test_model.fit() test_model.update(x=Xnew, Y=Ypred) - assert test_model._info["ll"] == reference_model["ll"] - assert test_model._info["theta"] == reference_model["theta"] - assert test_model._info["g"] == reference_model["g"] - assert test_model._info["beta0"] == reference_model["beta0"] + assert np.allclose(test_model._info["ll"], reference_model["ll"]) + assert np.allclose(test_model._info["theta"], reference_model["theta"]) + assert np.allclose(test_model._info["g"], reference_model["g"]) + assert np.allclose(test_model._info["beta0"], reference_model["beta0"]) def test_conversion_to_homGP(): diff --git a/tests/unit_tests/test_homGP.py b/tests/unit_tests/test_homGP.py index f5ce2d1..ca72d09 100644 --- a/tests/unit_tests/test_homGP.py +++ b/tests/unit_tests/test_homGP.py @@ -31,10 +31,10 @@ def test_fit(): test_model = emulator(x=X, theta=np.array([0]), f=Y, method="homGP") test_model.fit() - assert test_model._info["ll"] == reference_model["ll"] - assert test_model._info["theta"] == reference_model["theta"] - assert test_model._info["g"] == reference_model["g"] - assert test_model._info["beta0"] == reference_model["beta0"] + assert np.allclose(test_model._info["ll"], reference_model["ll"]) + assert np.allclose(test_model._info["theta"], reference_model["theta"]) + assert np.allclose(test_model._info["g"], reference_model["g"]) + assert np.allclose(test_model._info["beta0"], reference_model["beta0"]) def test_predict(): @@ -93,10 +93,10 @@ def test_matern(): ) test_model.fit() - assert test_model._info["ll"] == reference_model["ll"] - assert test_model._info["theta"] == reference_model["theta"] - assert test_model._info["g"] == reference_model["g"] - assert test_model._info["beta0"] == reference_model["beta0"] + assert np.allclose(test_model._info["ll"], reference_model["ll"]) + assert np.allclose(test_model._info["theta"], reference_model["theta"]) + assert np.allclose(test_model._info["g"], reference_model["g"]) + assert np.allclose(test_model._info["beta0"], reference_model["beta0"]) def test_homGP_update_kriging_believer(): @@ -111,10 +111,10 @@ def test_homGP_update_kriging_believer(): test_model.fit() test_model.update(Xnew) - assert test_model._info["ll"] == reference_model["ll"] - assert test_model._info["theta"] == reference_model["theta"] - assert test_model._info["g"] == reference_model["g"] - assert test_model._info["beta0"] == reference_model["beta0"] + assert np.allclose(test_model._info["ll"], reference_model["ll"]) + assert np.allclose(test_model._info["theta"], reference_model["theta"]) + assert np.allclose(test_model._info["g"], reference_model["g"]) + assert np.allclose(test_model._info["beta0"], reference_model["beta0"]) def test_homGP_update(): @@ -129,10 +129,10 @@ def test_homGP_update(): test_model.fit() test_model.update(x=Xnew, Y=Ypred) - assert test_model._info["ll"] == reference_model["ll"] - assert test_model._info["theta"] == reference_model["theta"] - assert test_model._info["g"] == reference_model["g"] - assert test_model._info["beta0"] == reference_model["beta0"] + assert np.allclose(test_model._info["ll"], reference_model["ll"]) + assert np.allclose(test_model._info["theta"], reference_model["theta"]) + assert np.allclose(test_model._info["g"], reference_model["g"]) + assert np.allclose(test_model._info["beta0"], reference_model["beta0"]) if __name__ == "__main__": From 8713e129747560ca21847222a966f52a694b626a Mon Sep 17 00:00:00 2001 From: David O'Gara Date: Tue, 7 Oct 2025 11:50:07 -0500 Subject: [PATCH 2/7] multihetGP tests failing for some reason --- tests/unit_tests/test_multihetGP.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/unit_tests/test_multihetGP.py b/tests/unit_tests/test_multihetGP.py index 15ee070..54b5954 100644 --- a/tests/unit_tests/test_multihetGP.py +++ b/tests/unit_tests/test_multihetGP.py @@ -12,8 +12,8 @@ lhs = qmc.LatinHypercube(d=2, rng=rng) X = lhs.random(n=100) -reps = rng.choice(len(X), size=100, replace=True) -X = X[reps,] +mult = rng.choice(np.arange(1,6),size = X.shape[0], replace = True) +X = np.vstack(np.repeat(X,mult,axis=0)) # prediction grid Xp = lhs.random(n=100) Y = 0 * X @@ -21,7 +21,7 @@ Y[:, 1] = np.cos(X[:, 1]) # varying noise field noise = rng.normal(size=Y.shape) * np.exp(-(X**2)) -noise *= 2 +noise *= 0.2 Y += noise @@ -68,7 +68,7 @@ def test_fit(): "noiseControl": NOISECONTROL, }, ) - # multiGP.fit() + multiGP.fit() emulator_keys = ["ll", "theta", "beta0", "Delta"] for i in range(len(GPlist)): @@ -182,4 +182,4 @@ def test_update(): if __name__ == "__main__": - test_update() + test_fit() From 953affcdd422223a890eafe82a2cd5229c913192 Mon Sep 17 00:00:00 2001 From: David O'Gara Date: Tue, 7 Oct 2025 12:43:11 -0500 Subject: [PATCH 3/7] revert back to previous version of multihetGP --- tests/unit_tests/test_multihetGP.py | 179 ++++++++++------------------ 1 file changed, 65 insertions(+), 114 deletions(-) diff --git a/tests/unit_tests/test_multihetGP.py b/tests/unit_tests/test_multihetGP.py index 54b5954..9b194d4 100644 --- a/tests/unit_tests/test_multihetGP.py +++ b/tests/unit_tests/test_multihetGP.py @@ -1,185 +1,136 @@ import sys - -sys.path.append("./") +sys.path.append('./') import numpy as np from scipy.stats import qmc from PUQ.surrogate import emulator from hetgpy import hetGP from copy import deepcopy - # create some noise 2D data rng = np.random.default_rng(1) -lhs = qmc.LatinHypercube(d=2, rng=rng) +lhs = qmc.LatinHypercube(d=2,rng=rng) X = lhs.random(n=100) -mult = rng.choice(np.arange(1,6),size = X.shape[0], replace = True) -X = np.vstack(np.repeat(X,mult,axis=0)) +reps = rng.choice(len(X),size=100,replace=True) +X = X[reps,] # prediction grid Xp = lhs.random(n=100) -Y = 0 * X -Y[:, 0] = np.sin(X[:, 0]) -Y[:, 1] = np.cos(X[:, 1]) +Y = 0*X +Y[:,0] = np.sin(X[:,0]) +Y[:,1] = np.cos(X[:,1]) # varying noise field -noise = rng.normal(size=Y.shape) * np.exp(-(X**2)) -noise *= 0.2 +noise = rng.normal(size=Y.shape) * np.exp(-X**2) +noise *=2 -Y += noise +Y+= noise COVTYPE = "Matern5_2" # module wide settings -SETTINGS = { - "linkThetas": "joint", - "logN": True, - "initStrategy": "residuals", - "checkHom": True, - "penalty": True, - "trace": 0, - "return.matrices": True, - "return.hom": False, - "factr": 1e9, -} -NOISECONTROL = {"k_theta_g_bounds": (1, 100), "g_max": 1e2, "g_bounds": (1e-6, 1)} - +SETTINGS = {"linkThetas": 'joint', "logN": True, "initStrategy": 'residuals', + "checkHom": True, "penalty": True, "trace": 0, "return.matrices": True, + "return.hom": False, "factr": 1e9} +NOISECONTROL = {'k_theta_g_bounds': (1, 100), 'g_max': 1e2, 'g_bounds': (1e-6, 1)} + MAXIT = 100 - - def test_fit(): GPlist = [ - hetGP().mle( - X=X, - Z=Y[:, i], - covtype=COVTYPE, - maxit=MAXIT, - settings=SETTINGS, - noiseControl=NOISECONTROL, - ) + hetGP().mle(X=X,Z=Y[:,i],covtype=COVTYPE,maxit=MAXIT, + settings=SETTINGS,noiseControl=NOISECONTROL) for i in range(Y.shape[1]) ] - multiGP = emulator( - x=X, - theta=np.array([0, 1]).reshape(2, 1), - f=Y, - method="multihetGP", - args={ - "covtype": COVTYPE, - "maxit": MAXIT, - "settings": SETTINGS, - "noiseControl": NOISECONTROL, - }, - ) - multiGP.fit() - - emulator_keys = ["ll", "theta", "beta0", "Delta"] + multiGP = emulator(x=X,theta=np.array([0,1]).reshape(2,1),f=Y, + method='multihetGP',args={'covtype':COVTYPE,'maxit':MAXIT,'settings':SETTINGS,'noiseControl':NOISECONTROL}) + #multiGP.fit() + + emulator_keys = ['ll','theta','beta0','Delta'] for i in range(len(GPlist)): GP = GPlist[i] for key in emulator_keys: - assert np.allclose(multiGP._info["emulist"][i]._info[key], GP[key]) - + assert np.allclose(multiGP._info['emulist'][i]._info[key],GP[key]) + return - def test_predict(): - - multiGP = emulator( - x=X, - theta=np.array([0, 1]).reshape(2, 1), - f=Y, - method="multihetGP", - args={"covtype": COVTYPE, "maxit": MAXIT, "settings": SETTINGS}, - ) + + multiGP = emulator(x=X,theta=np.array([0,1]).reshape(2,1),f=Y, + method='multihetGP',args={'covtype':COVTYPE,'maxit':MAXIT,'settings':SETTINGS}) multiGP.fit() - preds = multiGP.predict(x=Xp, thetaprime=Xp) + preds = multiGP.predict(x=Xp,thetaprime=Xp) GPlist = [ - hetGP() - .mle(X=X, Z=Y[:, i], covtype=COVTYPE, maxit=MAXIT, settings=SETTINGS) - .predict(Xp, xprime=Xp) + hetGP().mle(X=X,Z=Y[:,i],covtype=COVTYPE, + maxit=MAXIT,settings=SETTINGS + ).predict(Xp,xprime=Xp) for i in range(Y.shape[1]) ] - em2GPkey = {"mean": "mean", "var": "sd2", "nugs": "nugs"} + em2GPkey = {'mean':'mean','var':'sd2','nugs':'nugs'} for i in range(len(GPlist)): GP = GPlist[i] for em_key, GP_key in em2GPkey.items(): - assert np.allclose(preds._info[em_key][i, :], GP[GP_key]) - + assert np.allclose(preds._info[em_key][i,:],GP[GP_key]) + # test cov and covmat for i in range(len(GPlist)): - assert np.allclose(preds._info["covmat"][i, :, :], GPlist[i]["cov"]) + assert np.allclose(preds._info['covmat'][i,:,:],GPlist[i]['cov']) return - def test_update_kriging_believer(): - - multiGP = emulator( - x=X, - theta=np.array([0, 1]).reshape(2, 1), - f=Y, - method="multihetGP", - args={"covtype": COVTYPE, "maxit": MAXIT, "settings": SETTINGS}, - ) + + multiGP = emulator(x=X,theta=np.array([0,1]).reshape(2,1),f=Y, + method='multihetGP',args={'covtype':COVTYPE,'maxit':MAXIT,'settings':SETTINGS}) multiGP.fit() - initial_lls = [ - multiGP._info["emulist"][i]._info["ll"] for i in range(multiGP._info["numGPs"]) - ] + initial_lls = [multiGP._info['emulist'][i]._info['ll'] for i in range(multiGP._info['numGPs'])] initial_lls = deepcopy(initial_lls) GPlist = [ - hetGP().mle(X=X, Z=Y[:, i], covtype=COVTYPE, maxit=MAXIT, settings=SETTINGS) + hetGP().mle(X=X,Z=Y[:,i],covtype=COVTYPE,maxit=MAXIT, + settings=SETTINGS) for i in range(Y.shape[1]) ] GPlist_initial = deepcopy(GPlist) - Xnew = X.mean(axis=0).reshape(-1, X.shape[1]) + Xnew = X.mean(axis=0).reshape(-1,X.shape[1]) multiGP.update(x=Xnew) - em_keys = ["ll", "theta", "beta0", "Delta"] + em_keys = ['ll','theta','beta0','Delta'] for GP in GPlist: - Ynew = GP.predict(Xnew)["mean"] - GP.update(Xnew, Ynew, maxit=0) - for i in range(multiGP._info["numGPs"]): + Ynew = GP.predict(Xnew)['mean'] + GP.update(Xnew,Ynew,maxit=0) + for i in range(multiGP._info['numGPs']): GP = GPlist[i] for key in em_keys: - assert np.allclose(multiGP._info["emulist"][i]._info[key], GP[key]) - + assert np.allclose(multiGP._info['emulist'][i]._info[key],GP[key]) def test_update(): - """Test update and that log likelihood changes""" - - multiGP = emulator( - x=X, - theta=np.array([0, 1]).reshape(2, 1), - f=Y, - method="multihetGP", - args={"covtype": COVTYPE, "maxit": MAXIT, "settings": SETTINGS}, - ) + '''Test update and that log likelihood changes''' + + multiGP = emulator(x=X,theta=np.array([0,1]).reshape(2,1),f=Y, + method='multihetGP',args={'covtype':COVTYPE,'maxit':MAXIT,'settings':SETTINGS}) multiGP.fit() - initial_lls = [ - multiGP._info["emulist"][i]._info["ll"] for i in range(multiGP._info["numGPs"]) - ] + initial_lls = [multiGP._info['emulist'][i]._info['ll'] for i in range(multiGP._info['numGPs'])] initial_lls = deepcopy(initial_lls) - + GPlist = [ - hetGP().mle(X=X, Z=Y[:, i], covtype=COVTYPE, maxit=MAXIT, settings=SETTINGS) + hetGP().mle(X=X,Z=Y[:,i],covtype=COVTYPE,maxit=MAXIT, + settings=SETTINGS) for i in range(Y.shape[1]) ] GPlist_initial = deepcopy(GPlist) - Xnew = X.mean(axis=0).reshape(-1, X.shape[1]) + Xnew = X.mean(axis=0).reshape(-1,X.shape[1]) mpreds = multiGP.predict(x=Xnew) - Ynew = mpreds._info["mean"].T - multiGP.update(x=Xnew, Y=Ynew) - em_keys = ["ll", "theta", "beta0", "Delta"] + Ynew = mpreds._info['mean'].T + multiGP.update(x=Xnew,Y=Ynew) + em_keys = ['ll','theta','beta0','Delta'] for GP in GPlist: - Ynew = GP.predict(Xnew)["mean"] - GP.update(Xnew, Ynew, maxit=MAXIT) - for i in range(multiGP._info["numGPs"]): + Ynew = GP.predict(Xnew)['mean'] + GP.update(Xnew,Ynew,maxit=MAXIT) + for i in range(multiGP._info['numGPs']): GP = GPlist[i] for key in em_keys: - assert np.allclose(multiGP._info["emulist"][i]._info[key], GP[key]) + assert np.allclose(multiGP._info['emulist'][i]._info[key],GP[key]) # also check that likelihoods changed after update for i in range(len(initial_lls)): old_ll = initial_lls[i] - assert not np.allclose(old_ll, multiGP._info["emulist"][i]._info["ll"]) - + assert not np.allclose(old_ll,multiGP._info['emulist'][i]._info['ll']) if __name__ == "__main__": - test_fit() + test_update() \ No newline at end of file From d4c8924fd395fb38f74094cc1781c5c12229d7d8 Mon Sep 17 00:00:00 2001 From: David O'Gara Date: Tue, 7 Oct 2025 13:02:55 -0500 Subject: [PATCH 4/7] add more numerically stable homGP conversion example --- tests/unit_tests/test_hetGP.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/tests/unit_tests/test_hetGP.py b/tests/unit_tests/test_hetGP.py index c538738..c221cc2 100644 --- a/tests/unit_tests/test_hetGP.py +++ b/tests/unit_tests/test_hetGP.py @@ -141,7 +141,13 @@ def test_hetGP_update(): def test_conversion_to_homGP(): # should return homoskedastic GP X = np.linspace(0, 1, 20).reshape(-1, 1) - Y = X + rand = np.random.default_rng(2) + reps = rand.choice(len(X),size=len(X)) + X = X[reps,:] + Y = np.sin(X).squeeze() + noise = 0.2 * rand.normal(size=Y.shape[0]) + Y += noise + Y = Y.reshape(-1,1) model = emulator(x=X, theta=np.array([0]), f=Y, method="hetGP") Xp = np.linspace(X.min(), X.max(), 100).reshape(-1, 1) preds = model.predict(Xp) From 30b78e3dd43b0e8dd81b71238231beba21d2709f Mon Sep 17 00:00:00 2001 From: David O'Gara Date: Tue, 7 Oct 2025 13:05:26 -0500 Subject: [PATCH 5/7] restore black formatting --- tests/unit_tests/test_multihetGP.py | 178 ++++++++++++++++++---------- 1 file changed, 114 insertions(+), 64 deletions(-) diff --git a/tests/unit_tests/test_multihetGP.py b/tests/unit_tests/test_multihetGP.py index 9b194d4..07e1b64 100644 --- a/tests/unit_tests/test_multihetGP.py +++ b/tests/unit_tests/test_multihetGP.py @@ -1,136 +1,186 @@ + import sys -sys.path.append('./') + +sys.path.append("./") import numpy as np from scipy.stats import qmc from PUQ.surrogate import emulator from hetgpy import hetGP from copy import deepcopy + # create some noise 2D data rng = np.random.default_rng(1) -lhs = qmc.LatinHypercube(d=2,rng=rng) +lhs = qmc.LatinHypercube(d=2, rng=rng) X = lhs.random(n=100) -reps = rng.choice(len(X),size=100,replace=True) +reps = rng.choice(len(X), size=100, replace=True) X = X[reps,] # prediction grid Xp = lhs.random(n=100) -Y = 0*X -Y[:,0] = np.sin(X[:,0]) -Y[:,1] = np.cos(X[:,1]) +Y = 0 * X +Y[:, 0] = np.sin(X[:, 0]) +Y[:, 1] = np.cos(X[:, 1]) # varying noise field -noise = rng.normal(size=Y.shape) * np.exp(-X**2) -noise *=2 +noise = rng.normal(size=Y.shape) * np.exp(-(X**2)) +noise *= 2 -Y+= noise +Y += noise COVTYPE = "Matern5_2" # module wide settings -SETTINGS = {"linkThetas": 'joint', "logN": True, "initStrategy": 'residuals', - "checkHom": True, "penalty": True, "trace": 0, "return.matrices": True, - "return.hom": False, "factr": 1e9} -NOISECONTROL = {'k_theta_g_bounds': (1, 100), 'g_max': 1e2, 'g_bounds': (1e-6, 1)} - +SETTINGS = { + "linkThetas": "joint", + "logN": True, + "initStrategy": "residuals", + "checkHom": True, + "penalty": True, + "trace": 0, + "return.matrices": True, + "return.hom": False, + "factr": 1e9, +} +NOISECONTROL = {"k_theta_g_bounds": (1, 100), "g_max": 1e2, "g_bounds": (1e-6, 1)} + MAXIT = 100 + + def test_fit(): GPlist = [ - hetGP().mle(X=X,Z=Y[:,i],covtype=COVTYPE,maxit=MAXIT, - settings=SETTINGS,noiseControl=NOISECONTROL) + hetGP().mle( + X=X, + Z=Y[:, i], + covtype=COVTYPE, + maxit=MAXIT, + settings=SETTINGS, + noiseControl=NOISECONTROL, + ) for i in range(Y.shape[1]) ] - multiGP = emulator(x=X,theta=np.array([0,1]).reshape(2,1),f=Y, - method='multihetGP',args={'covtype':COVTYPE,'maxit':MAXIT,'settings':SETTINGS,'noiseControl':NOISECONTROL}) - #multiGP.fit() - - emulator_keys = ['ll','theta','beta0','Delta'] + multiGP = emulator( + x=X, + theta=np.array([0, 1]).reshape(2, 1), + f=Y, + method="multihetGP", + args={ + "covtype": COVTYPE, + "maxit": MAXIT, + "settings": SETTINGS, + "noiseControl": NOISECONTROL, + }, + ) + # multiGP.fit() + + emulator_keys = ["ll", "theta", "beta0", "Delta"] for i in range(len(GPlist)): GP = GPlist[i] for key in emulator_keys: - assert np.allclose(multiGP._info['emulist'][i]._info[key],GP[key]) - + assert np.allclose(multiGP._info["emulist"][i]._info[key], GP[key]) + return + def test_predict(): - - multiGP = emulator(x=X,theta=np.array([0,1]).reshape(2,1),f=Y, - method='multihetGP',args={'covtype':COVTYPE,'maxit':MAXIT,'settings':SETTINGS}) + + multiGP = emulator( + x=X, + theta=np.array([0, 1]).reshape(2, 1), + f=Y, + method="multihetGP", + args={"covtype": COVTYPE, "maxit": MAXIT, "settings": SETTINGS}, + ) multiGP.fit() - preds = multiGP.predict(x=Xp,thetaprime=Xp) + preds = multiGP.predict(x=Xp, thetaprime=Xp) GPlist = [ - hetGP().mle(X=X,Z=Y[:,i],covtype=COVTYPE, - maxit=MAXIT,settings=SETTINGS - ).predict(Xp,xprime=Xp) + hetGP() + .mle(X=X, Z=Y[:, i], covtype=COVTYPE, maxit=MAXIT, settings=SETTINGS) + .predict(Xp, xprime=Xp) for i in range(Y.shape[1]) ] - em2GPkey = {'mean':'mean','var':'sd2','nugs':'nugs'} + em2GPkey = {"mean": "mean", "var": "sd2", "nugs": "nugs"} for i in range(len(GPlist)): GP = GPlist[i] for em_key, GP_key in em2GPkey.items(): - assert np.allclose(preds._info[em_key][i,:],GP[GP_key]) - + assert np.allclose(preds._info[em_key][i, :], GP[GP_key]) + # test cov and covmat for i in range(len(GPlist)): - assert np.allclose(preds._info['covmat'][i,:,:],GPlist[i]['cov']) + assert np.allclose(preds._info["covmat"][i, :, :], GPlist[i]["cov"]) return + def test_update_kriging_believer(): - - multiGP = emulator(x=X,theta=np.array([0,1]).reshape(2,1),f=Y, - method='multihetGP',args={'covtype':COVTYPE,'maxit':MAXIT,'settings':SETTINGS}) + + multiGP = emulator( + x=X, + theta=np.array([0, 1]).reshape(2, 1), + f=Y, + method="multihetGP", + args={"covtype": COVTYPE, "maxit": MAXIT, "settings": SETTINGS}, + ) multiGP.fit() - initial_lls = [multiGP._info['emulist'][i]._info['ll'] for i in range(multiGP._info['numGPs'])] + initial_lls = [ + multiGP._info["emulist"][i]._info["ll"] for i in range(multiGP._info["numGPs"]) + ] initial_lls = deepcopy(initial_lls) GPlist = [ - hetGP().mle(X=X,Z=Y[:,i],covtype=COVTYPE,maxit=MAXIT, - settings=SETTINGS) + hetGP().mle(X=X, Z=Y[:, i], covtype=COVTYPE, maxit=MAXIT, settings=SETTINGS) for i in range(Y.shape[1]) ] GPlist_initial = deepcopy(GPlist) - Xnew = X.mean(axis=0).reshape(-1,X.shape[1]) + Xnew = X.mean(axis=0).reshape(-1, X.shape[1]) multiGP.update(x=Xnew) - em_keys = ['ll','theta','beta0','Delta'] + em_keys = ["ll", "theta", "beta0", "Delta"] for GP in GPlist: - Ynew = GP.predict(Xnew)['mean'] - GP.update(Xnew,Ynew,maxit=0) - for i in range(multiGP._info['numGPs']): + Ynew = GP.predict(Xnew)["mean"] + GP.update(Xnew, Ynew, maxit=0) + for i in range(multiGP._info["numGPs"]): GP = GPlist[i] for key in em_keys: - assert np.allclose(multiGP._info['emulist'][i]._info[key],GP[key]) + assert np.allclose(multiGP._info["emulist"][i]._info[key], GP[key]) + def test_update(): - '''Test update and that log likelihood changes''' - - multiGP = emulator(x=X,theta=np.array([0,1]).reshape(2,1),f=Y, - method='multihetGP',args={'covtype':COVTYPE,'maxit':MAXIT,'settings':SETTINGS}) + """Test update and that log likelihood changes""" + + multiGP = emulator( + x=X, + theta=np.array([0, 1]).reshape(2, 1), + f=Y, + method="multihetGP", + args={"covtype": COVTYPE, "maxit": MAXIT, "settings": SETTINGS}, + ) multiGP.fit() - initial_lls = [multiGP._info['emulist'][i]._info['ll'] for i in range(multiGP._info['numGPs'])] + initial_lls = [ + multiGP._info["emulist"][i]._info["ll"] for i in range(multiGP._info["numGPs"]) + ] initial_lls = deepcopy(initial_lls) - + GPlist = [ - hetGP().mle(X=X,Z=Y[:,i],covtype=COVTYPE,maxit=MAXIT, - settings=SETTINGS) + hetGP().mle(X=X, Z=Y[:, i], covtype=COVTYPE, maxit=MAXIT, settings=SETTINGS) for i in range(Y.shape[1]) ] GPlist_initial = deepcopy(GPlist) - Xnew = X.mean(axis=0).reshape(-1,X.shape[1]) + Xnew = X.mean(axis=0).reshape(-1, X.shape[1]) mpreds = multiGP.predict(x=Xnew) - Ynew = mpreds._info['mean'].T - multiGP.update(x=Xnew,Y=Ynew) - em_keys = ['ll','theta','beta0','Delta'] + Ynew = mpreds._info["mean"].T + multiGP.update(x=Xnew, Y=Ynew) + em_keys = ["ll", "theta", "beta0", "Delta"] for GP in GPlist: - Ynew = GP.predict(Xnew)['mean'] - GP.update(Xnew,Ynew,maxit=MAXIT) - for i in range(multiGP._info['numGPs']): + Ynew = GP.predict(Xnew)["mean"] + GP.update(Xnew, Ynew, maxit=MAXIT) + for i in range(multiGP._info["numGPs"]): GP = GPlist[i] for key in em_keys: - assert np.allclose(multiGP._info['emulist'][i]._info[key],GP[key]) + assert np.allclose(multiGP._info["emulist"][i]._info[key], GP[key]) # also check that likelihoods changed after update for i in range(len(initial_lls)): old_ll = initial_lls[i] - assert not np.allclose(old_ll,multiGP._info['emulist'][i]._info['ll']) + assert not np.allclose(old_ll, multiGP._info["emulist"][i]._info["ll"]) + if __name__ == "__main__": - test_update() \ No newline at end of file + test_update() From 2dfa27cdbef2891671e5e98d93284618f93f41e1 Mon Sep 17 00:00:00 2001 From: David O'Gara Date: Tue, 7 Oct 2025 14:16:01 -0500 Subject: [PATCH 6/7] add more reps to stabilize example --- tests/unit_tests/test_multihetGP.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/unit_tests/test_multihetGP.py b/tests/unit_tests/test_multihetGP.py index 07e1b64..1666823 100644 --- a/tests/unit_tests/test_multihetGP.py +++ b/tests/unit_tests/test_multihetGP.py @@ -9,11 +9,12 @@ from copy import deepcopy # create some noise 2D data +np.random.seed(1) rng = np.random.default_rng(1) lhs = qmc.LatinHypercube(d=2, rng=rng) X = lhs.random(n=100) -reps = rng.choice(len(X), size=100, replace=True) +reps = rng.choice(len(X), size=5000, replace=True) X = X[reps,] # prediction grid Xp = lhs.random(n=100) @@ -21,7 +22,7 @@ Y[:, 0] = np.sin(X[:, 0]) Y[:, 1] = np.cos(X[:, 1]) # varying noise field -noise = rng.normal(size=Y.shape) * np.exp(-(X**2)) +noise = rng.normal(size=Y.shape) noise *= 2 Y += noise @@ -42,7 +43,7 @@ } NOISECONTROL = {"k_theta_g_bounds": (1, 100), "g_max": 1e2, "g_bounds": (1e-6, 1)} -MAXIT = 100 +MAXIT = 200 def test_fit(): From 4aa8d8f25802a29d1ffbb11cb65b070812568614 Mon Sep 17 00:00:00 2001 From: David O'Gara Date: Tue, 7 Oct 2025 15:17:27 -0500 Subject: [PATCH 7/7] black --- tests/unit_tests/test_hetGP.py | 6 +++--- tests/unit_tests/test_multihetGP.py | 1 - 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/tests/unit_tests/test_hetGP.py b/tests/unit_tests/test_hetGP.py index c221cc2..c06560d 100644 --- a/tests/unit_tests/test_hetGP.py +++ b/tests/unit_tests/test_hetGP.py @@ -142,12 +142,12 @@ def test_conversion_to_homGP(): # should return homoskedastic GP X = np.linspace(0, 1, 20).reshape(-1, 1) rand = np.random.default_rng(2) - reps = rand.choice(len(X),size=len(X)) - X = X[reps,:] + reps = rand.choice(len(X), size=len(X)) + X = X[reps, :] Y = np.sin(X).squeeze() noise = 0.2 * rand.normal(size=Y.shape[0]) Y += noise - Y = Y.reshape(-1,1) + Y = Y.reshape(-1, 1) model = emulator(x=X, theta=np.array([0]), f=Y, method="hetGP") Xp = np.linspace(X.min(), X.max(), 100).reshape(-1, 1) preds = model.predict(Xp) diff --git a/tests/unit_tests/test_multihetGP.py b/tests/unit_tests/test_multihetGP.py index 1666823..61c9dbd 100644 --- a/tests/unit_tests/test_multihetGP.py +++ b/tests/unit_tests/test_multihetGP.py @@ -1,4 +1,3 @@ - import sys sys.path.append("./")