Standard GP¶

In [1]:

import fastgps
import torch
import numpy as np

import fastgps import torch import numpy as np

In [2]:

torch.set_default_dtype(torch.float64)

torch.set_default_dtype(torch.float64)

True Function¶

In [3]:

def f_ackley(x, a=20, b=0.2, c=2*np.pi, scaling=32.768):
    # https://www.sfu.ca/~ssurjano/ackley.html
    assert x.ndim==2
    x = 2*scaling*x-scaling
    t1 = a*torch.exp(-b*torch.sqrt(torch.mean(x**2,1)))
    t2 = torch.exp(torch.mean(torch.cos(c*x),1))
    t3 = a+np.exp(1)
    y = -t1-t2+t3
    return y
f_low_fidelity = lambda x: f_ackley(x,c=0)
f_high_fidelity = lambda x: f_ackley(x)
f_cos = lambda x: torch.cos(2*np.pi*x).sum(1)
fs = [f_low_fidelity,f_high_fidelity,f_cos]
d = 1 # dimension
rng = torch.Generator().manual_seed(17)
x = torch.rand((2**7,d),generator=rng) # random testing locations
y = torch.vstack([f(x) for f in fs]) # true values at random testing locations
z = torch.rand((2**8,d),generator=rng) # other random locations at which to evaluate covariance
print("x.shape = %s"%str(tuple(x.shape)))
print("y.shape = %s"%str(tuple(y.shape)))
print("z.shape = %s"%str(tuple(z.shape)))

def f_ackley(x, a=20, b=0.2, c=2*np.pi, scaling=32.768): # https://www.sfu.ca/~ssurjano/ackley.html assert x.ndim==2 x = 2*scaling*x-scaling t1 = a*torch.exp(-b*torch.sqrt(torch.mean(x**2,1))) t2 = torch.exp(torch.mean(torch.cos(c*x),1)) t3 = a+np.exp(1) y = -t1-t2+t3 return y f_low_fidelity = lambda x: f_ackley(x,c=0) f_high_fidelity = lambda x: f_ackley(x) f_cos = lambda x: torch.cos(2*np.pi*x).sum(1) fs = [f_low_fidelity,f_high_fidelity,f_cos] d = 1 # dimension rng = torch.Generator().manual_seed(17) x = torch.rand((2**7,d),generator=rng) # random testing locations y = torch.vstack([f(x) for f in fs]) # true values at random testing locations z = torch.rand((2**8,d),generator=rng) # other random locations at which to evaluate covariance print("x.shape = %s"%str(tuple(x.shape))) print("y.shape = %s"%str(tuple(y.shape))) print("z.shape = %s"%str(tuple(z.shape)))

x.shape = (128, 1)
y.shape = (3, 128)
z.shape = (256, 1)

Construct Fast GP¶

In [4]:

fgp = fastgps.StandardGP(d,seed_for_seq=7,num_tasks=len(fs))
x_next = fgp.get_x_next(n=[2**4,2**2,2**3])
y_next = [fs[i](x_next[i]) for i in range(fgp.num_tasks)]
fgp.add_y_next(y_next)
assert len(x_next)==len(y_next)
for i in range(len(x_next)):
    print("i = %d"%i)
    print("\tx_next[%d].shape = %s"%(i,str(tuple(x_next[i].shape))))
    print("\ty_next[%d].shape = %s"%(i,str(tuple(y_next[i].shape))))

fgp = fastgps.StandardGP(d,seed_for_seq=7,num_tasks=len(fs)) x_next = fgp.get_x_next(n=[2**4,2**2,2**3]) y_next = [fs[i](x_next[i]) for i in range(fgp.num_tasks)] fgp.add_y_next(y_next) assert len(x_next)==len(y_next) for i in range(len(x_next)): print("i = %d"%i) print("\tx_next[%d].shape = %s"%(i,str(tuple(x_next[i].shape)))) print("\ty_next[%d].shape = %s"%(i,str(tuple(y_next[i].shape))))

i = 0
	x_next[0].shape = (16, 1)
	y_next[0].shape = (16,)
i = 1
	x_next[1].shape = (4, 1)
	y_next[1].shape = (4,)
i = 2
	x_next[2].shape = (8, 1)
	y_next[2].shape = (8,)

In [5]:

pmean = fgp.post_mean(x)
print("pmean.shape = %s"%str(tuple(pmean.shape)))
print("l2 relative error =",(torch.linalg.norm(y-pmean,dim=1)/torch.linalg.norm(y,dim=1)))

pmean = fgp.post_mean(x) print("pmean.shape = %s"%str(tuple(pmean.shape))) print("l2 relative error =",(torch.linalg.norm(y-pmean,dim=1)/torch.linalg.norm(y,dim=1)))

pmean.shape = (3, 128)
l2 relative error = tensor([0.1920, 0.4400, 0.1470])

In [6]:

data = fgp.fit()
list(data.keys())

data = fgp.fit() list(data.keys())

     iter of 5.0e+03 | loss       | term1      | term2     
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            0.00e+00 | 5.46e+05   | 1.09e+06   | -1.66e+02 
            5.00e+00 | 1.17e+05   | 2.35e+05   | -1.23e+02 
            1.00e+01 | 4.24e+03   | 8.41e+03   | 6.12e+00  
            1.50e+01 | 8.75e+01   | 8.74e+00   | 1.15e+02  
            2.00e+01 | 7.38e+01   | 1.94e+01   | 7.67e+01  
            2.50e+01 | 6.72e+01   | 2.81e+01   | 5.47e+01  
            3.00e+01 | 6.31e+01   | 2.63e+01   | 4.85e+01  
            3.50e+01 | 5.81e+01   | 2.83e+01   | 3.64e+01  
            4.00e+01 | 5.06e+01   | 2.94e+01   | 2.03e+01  
            4.50e+01 | 5.11e+01   | 3.02e+01   | 2.04e+01  
            5.00e+01 | 4.90e+01   | 2.69e+01   | 1.96e+01  
            5.50e+01 | 4.73e+01   | 2.69e+01   | 1.64e+01  
            6.00e+01 | 4.61e+01   | 2.64e+01   | 1.43e+01  
            6.50e+01 | 4.59e+01   | 2.49e+01   | 1.55e+01  
            7.00e+01 | 4.59e+01   | 2.48e+01   | 1.56e+01  
            7.30e+01 | 4.59e+01   | 2.46e+01   | 1.59e+01  

Out[6]:

['iterations']

In [7]:

pmean,pvar,q,ci_low,ci_high = fgp.post_ci(x,confidence=0.99)
print("pmean.shape = %s"%str(tuple(pmean.shape)))
print("pvar.shape = %s"%str(tuple(pvar.shape)))
print("q = %.2f"%q)
print("ci_low.shape = %s"%str(tuple(ci_low.shape)))
print("ci_high.shape = %s"%str(tuple(ci_high.shape)))
print("l2 relative error =",(torch.linalg.norm(y-pmean,dim=1)/torch.linalg.norm(y,dim=1)))
pcov = fgp.post_cov(x,x)
print("pcov.shape = %s"%str(tuple(pcov.shape)))
_range0,_rangen1 = torch.arange(pcov.size(0)),torch.arange(pcov.size(-1))
assert torch.allclose(pcov[_range0,_range0][:,_rangen1,_rangen1],pvar) and (pvar>=0).all()
pcov2 = fgp.post_cov(x,z)
print("pcov2.shape = %s"%str(tuple(pcov2.shape)))

pmean,pvar,q,ci_low,ci_high = fgp.post_ci(x,confidence=0.99) print("pmean.shape = %s"%str(tuple(pmean.shape))) print("pvar.shape = %s"%str(tuple(pvar.shape))) print("q = %.2f"%q) print("ci_low.shape = %s"%str(tuple(ci_low.shape))) print("ci_high.shape = %s"%str(tuple(ci_high.shape))) print("l2 relative error =",(torch.linalg.norm(y-pmean,dim=1)/torch.linalg.norm(y,dim=1))) pcov = fgp.post_cov(x,x) print("pcov.shape = %s"%str(tuple(pcov.shape))) _range0,_rangen1 = torch.arange(pcov.size(0)),torch.arange(pcov.size(-1)) assert torch.allclose(pcov[_range0,_range0][:,_rangen1,_rangen1],pvar) and (pvar>=0).all() pcov2 = fgp.post_cov(x,z) print("pcov2.shape = %s"%str(tuple(pcov2.shape)))

pmean.shape = (3, 128)
pvar.shape = (3, 128)
q = 2.58
ci_low.shape = (3, 128)
ci_high.shape = (3, 128)
l2 relative error = tensor([0.0459, 0.0691, 0.0736])
pcov.shape = (3, 3, 128, 128)
pcov2.shape = (3, 3, 128, 256)

In [8]:

pcmean,pcvar,q,cci_low,cci_high = fgp.post_cubature_ci(confidence=0.99)
print("pcmean =",pcmean)
print("pcvar =",pcvar)
print("cci_low =",cci_low)
print("cci_high",cci_high)

pcmean,pcvar,q,cci_low,cci_high = fgp.post_cubature_ci(confidence=0.99) print("pcmean =",pcmean) print("pcvar =",pcvar) print("cci_low =",cci_low) print("cci_high",cci_high)

pcmean = tensor([1.6813e+01, 1.8687e+01, 1.4150e-02])
pcvar = tensor([0.0033, 0.0041, 0.0009])
cci_low = tensor([16.6657, 18.5230, -0.0652])
cci_high tensor([16.9608, 18.8510,  0.0935])

Project and Increase Sample Size¶

In [9]:

n_new = fgp.n*torch.tensor([4,2,8])
pcov_future = fgp.post_cov(x,z,n=n_new)
pvar_future = fgp.post_var(x,n=n_new)
pcvar_future = fgp.post_cubature_var(n=n_new)

n_new = fgp.n*torch.tensor([4,2,8]) pcov_future = fgp.post_cov(x,z,n=n_new) pvar_future = fgp.post_var(x,n=n_new) pcvar_future = fgp.post_cubature_var(n=n_new)

In [10]:

x_next = fgp.get_x_next(n_new)
y_next = [fs[i](x_next[i]) for i in range(fgp.num_tasks)]
for _y in y_next:
    print(_y.shape)
fgp.add_y_next(y_next)
print("l2 relative error =",(torch.linalg.norm(y-fgp.post_mean(x),dim=1)/torch.linalg.norm(y,dim=1)))
assert torch.allclose(fgp.post_cov(x,z),pcov_future)
assert torch.allclose(fgp.post_var(x),pvar_future)
assert torch.allclose(fgp.post_cubature_var(),pcvar_future)

x_next = fgp.get_x_next(n_new) y_next = [fs[i](x_next[i]) for i in range(fgp.num_tasks)] for _y in y_next: print(_y.shape) fgp.add_y_next(y_next) print("l2 relative error =",(torch.linalg.norm(y-fgp.post_mean(x),dim=1)/torch.linalg.norm(y,dim=1))) assert torch.allclose(fgp.post_cov(x,z),pcov_future) assert torch.allclose(fgp.post_var(x),pvar_future) assert torch.allclose(fgp.post_cubature_var(),pcvar_future)

torch.Size([48])
torch.Size([4])
torch.Size([56])
l2 relative error = tensor([0.0309, 0.0606, 0.0021])

In [11]:

data = fgp.fit(verbose=False)
print("l2 relative error =",(torch.linalg.norm(y-fgp.post_mean(x),dim=1)/torch.linalg.norm(y,dim=1)))

data = fgp.fit(verbose=False) print("l2 relative error =",(torch.linalg.norm(y-fgp.post_mean(x),dim=1)/torch.linalg.norm(y,dim=1)))

l2 relative error = tensor([5.6070e-02, 7.2271e-02, 8.8313e-09])

In [12]:

n_new = fgp.n*torch.tensor([4,8,2])
pcov_new = fgp.post_cov(x,z,n=n_new)
pvar_new = fgp.post_var(x,n=n_new)
pcvar_new = fgp.post_cubature_var(n=n_new)
x_next = fgp.get_x_next(n_new)
y_next = [fs[i](x_next[i]) for i in range(fgp.num_tasks)]
fgp.add_y_next(y_next)
assert torch.allclose(fgp.post_cov(x,z),pcov_new)
assert torch.allclose(fgp.post_var(x),pvar_new)
assert torch.allclose(fgp.post_cubature_var(),pcvar_new)

n_new = fgp.n*torch.tensor([4,8,2]) pcov_new = fgp.post_cov(x,z,n=n_new) pvar_new = fgp.post_var(x,n=n_new) pcvar_new = fgp.post_cubature_var(n=n_new) x_next = fgp.get_x_next(n_new) y_next = [fs[i](x_next[i]) for i in range(fgp.num_tasks)] fgp.add_y_next(y_next) assert torch.allclose(fgp.post_cov(x,z),pcov_new) assert torch.allclose(fgp.post_var(x),pvar_new) assert torch.allclose(fgp.post_cubature_var(),pcvar_new)