Compare GPs + Plot¶

In [1]:

import fastgps
import qmcpy as qp
import numpy as np
import torch
import pandas as pd
from matplotlib import pyplot

import fastgps import qmcpy as qp import numpy as np import torch import pandas as pd from matplotlib import pyplot

In [2]:

device = "cpu"
if device!="mps":
    torch.set_default_dtype(torch.float64)

device = "cpu" if device!="mps": torch.set_default_dtype(torch.float64)

In [3]:

colors = ["xkcd:"+color[:-1] for color in pd.read_csv("../../../xkcd_colors.txt",comment="#").iloc[:,0].tolist()][::-1]
_alpha = 0.25
WIDTH = 2*(500/72)
LINEWIDTH = 3
MARKERSIZE = 100

colors = ["xkcd:"+color[:-1] for color in pd.read_csv("../../../xkcd_colors.txt",comment="#").iloc[:,0].tolist()][::-1] _alpha = 0.25 WIDTH = 2*(500/72) LINEWIDTH = 3 MARKERSIZE = 100

True Function¶

In [4]:

d = 1 
f_smooth = lambda x: 15*(x[:,0]-1/2)**2*torch.sin(2*torch.pi*x[:,0])
def f_ackley(x, a=20, b=0.2, c=2*np.pi, scaling=32.768): # https://www.sfu.ca/~ssurjano/ackley.html
    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 = f_smooth
def fp(x):
    xg = x.clone().requires_grad_(True)
    yg = f(xg)
    yp = torch.autograd.grad(yg,xg,grad_outputs=torch.ones_like(yg))[0]
    return yp[:,0].detach()

d = 1 f_smooth = lambda x: 15*(x[:,0]-1/2)**2*torch.sin(2*torch.pi*x[:,0]) def f_ackley(x, a=20, b=0.2, c=2*np.pi, scaling=32.768): # https://www.sfu.ca/~ssurjano/ackley.html 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 = f_smooth def fp(x): xg = x.clone().requires_grad_(True) yg = f(xg) yp = torch.autograd.grad(yg,xg,grad_outputs=torch.ones_like(yg))[0] return yp[:,0].detach()

In [5]:

xticks = torch.linspace(0,1,252,device=device)[1:-1,None]
yticks = f(xticks)
ypticks = fp(xticks)
print("xticks.shape = %s"%str(tuple(xticks.shape)))
print("yticks.shape = %s"%str(tuple(yticks.shape)))
print("ypticks.shape = %s"%str(tuple(ypticks.shape)))

xticks = torch.linspace(0,1,252,device=device)[1:-1,None] yticks = f(xticks) ypticks = fp(xticks) print("xticks.shape = %s"%str(tuple(xticks.shape))) print("yticks.shape = %s"%str(tuple(yticks.shape))) print("ypticks.shape = %s"%str(tuple(ypticks.shape)))

xticks.shape = (250, 1)
yticks.shape = (250,)
ypticks.shape = (250,)

Fast GP Construction¶

In [6]:

gps = [
    fastgps.StandardGP(qp.KernelGaussian(1,torchify=True,device=device),seqs=qp.DigitalNetB2(1,seed=11,randomize="DS")),
    fastgps.FastGPDigitalNetB2(qp.KernelDigShiftInvar(1,torchify=True,device=device),seqs=qp.DigitalNetB2(1,seed=7,randomize="DS"),alpha=4),
    fastgps.FastGPLattice(qp.KernelShiftInvar(1,torchify=True,device=device),seqs=qp.Lattice(1,seed=7),alpha=4),
]

gps = [ fastgps.StandardGP(qp.KernelGaussian(1,torchify=True,device=device),seqs=qp.DigitalNetB2(1,seed=11,randomize="DS")), fastgps.FastGPDigitalNetB2(qp.KernelDigShiftInvar(1,torchify=True,device=device),seqs=qp.DigitalNetB2(1,seed=7,randomize="DS"),alpha=4), fastgps.FastGPLattice(qp.KernelShiftInvar(1,torchify=True,device=device),seqs=qp.Lattice(1,seed=7),alpha=4), ]

In [7]:

gps_grad = [
    fastgps.StandardGP(
        qp.KernelMultiTaskDerivs(qp.KernelGaussian(1,torchify=True,device=device),num_tasks=2),
        seqs = [qp.DigitalNetB2(1,seed=7,randomize="DS"),qp.DigitalNetB2(1,seed=11,randomize="DS")],
        derivatives = [torch.tensor([0],device=device),torch.tensor([1],device=device)],
    ),
    fastgps.FastGPDigitalNetB2(
        qp.KernelMultiTaskDerivs(qp.KernelDigShiftInvar(1,torchify=True,alpha=4,device=device),num_tasks=2),
        seqs = [qp.DigitalNetB2(1,seed=7,randomize="DS"),qp.DigitalNetB2(1,seed=11,randomize="DS")],
        derivatives = [torch.tensor([0],device=device),torch.tensor([1],device=device)],
    ),
    fastgps.FastGPLattice(
        qp.KernelMultiTaskDerivs(qp.KernelShiftInvar(1,torchify=True,alpha=4,device=device),num_tasks=2),
        seqs = [qp.Lattice(1,seed=7),qp.Lattice(1,seed=11)],
        derivatives = [torch.tensor([0],device=device),torch.tensor([1],device=device)],
    ),
]

gps_grad = [ fastgps.StandardGP( qp.KernelMultiTaskDerivs(qp.KernelGaussian(1,torchify=True,device=device),num_tasks=2), seqs = [qp.DigitalNetB2(1,seed=7,randomize="DS"),qp.DigitalNetB2(1,seed=11,randomize="DS")], derivatives = [torch.tensor([0],device=device),torch.tensor([1],device=device)], ), fastgps.FastGPDigitalNetB2( qp.KernelMultiTaskDerivs(qp.KernelDigShiftInvar(1,torchify=True,alpha=4,device=device),num_tasks=2), seqs = [qp.DigitalNetB2(1,seed=7,randomize="DS"),qp.DigitalNetB2(1,seed=11,randomize="DS")], derivatives = [torch.tensor([0],device=device),torch.tensor([1],device=device)], ), fastgps.FastGPLattice( qp.KernelMultiTaskDerivs(qp.KernelShiftInvar(1,torchify=True,alpha=4,device=device),num_tasks=2), seqs = [qp.Lattice(1,seed=7),qp.Lattice(1,seed=11)], derivatives = [torch.tensor([0],device=device),torch.tensor([1],device=device)], ), ]

GP Fitting¶

In [8]:

pmeans = [None]*len(gps)
pci_lows = [None]*len(gps)
pci_highs = [None]*len(gps)
for i,gp in enumerate(gps):
    print(type(gp).__name__)
    x_next = gp.get_x_next(n=2**2)
    gp.add_y_next(f(x_next))
    gp.fit()
    pmeans[i],_,_,pci_lows[i],pci_highs[i] = gp.post_ci(xticks,confidence=0.95)
    print("\tl2 relative error = %.1e"%(torch.linalg.norm(yticks-pmeans[i])/torch.linalg.norm(yticks)))

pmeans = [None]*len(gps) pci_lows = [None]*len(gps) pci_highs = [None]*len(gps) for i,gp in enumerate(gps): print(type(gp).__name__) x_next = gp.get_x_next(n=2**2) gp.add_y_next(f(x_next)) gp.fit() pmeans[i],_,_,pci_lows[i],pci_highs[i] = gp.post_ci(xticks,confidence=0.95) print("\tl2 relative error = %.1e"%(torch.linalg.norm(yticks-pmeans[i])/torch.linalg.norm(yticks)))

StandardGP
     iter of 5.0e+03 | best loss  | loss       | term1      | term2     
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            0.00e+00 | 4.88e+02   | 4.88e+02   | 9.81e+02   | -1.36e+01 
            5.00e+00 | 8.89e+00   | 8.89e+00   | 1.40e+01   | -3.56e+00 
            1.00e+01 | 4.62e+00   | 4.62e+00   | 2.93e+00   | -1.04e+00 
            1.50e+01 | 3.95e+00   | 3.96e+00   | 3.58e+00   | -3.01e+00 
            2.00e+01 | 3.94e+00   | 3.95e+00   | 3.60e+00   | -3.04e+00 
            2.40e+01 | 3.94e+00   | 3.94e+00   | 4.02e+00   | -3.48e+00 
	l2 relative error = 5.5e-01
FastGPDigitalNetB2
     iter of 5.0e+03 | best loss  | loss       | term1      | term2     
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            0.00e+00 | 5.76e+00   | 5.76e+00   | 1.39e+00   | 2.78e+00  
            5.00e+00 | 5.17e+00   | 5.17e+00   | 4.06e+00   | -1.07e+00 
            1.00e+01 | 5.09e+00   | 5.09e+00   | 3.21e+00   | -3.91e-01 
            1.50e+01 | 4.91e+00   | 4.91e+00   | 4.05e+00   | -1.57e+00 
            2.00e+01 | 4.75e+00   | 4.75e+00   | 3.80e+00   | -1.66e+00 
            2.50e+01 | 4.44e+00   | 4.44e+00   | 3.80e+00   | -2.28e+00 
            3.00e+01 | 4.27e+00   | 4.39e+00   | 5.70e+00   | -4.26e+00 
            3.50e+01 | 4.24e+00   | 4.27e+00   | 3.40e+00   | -2.22e+00 
            4.00e+01 | 4.24e+00   | 4.24e+00   | 3.90e+00   | -2.76e+00 
            4.30e+01 | 4.24e+00   | 4.24e+00   | 4.00e+00   | -2.87e+00 
	l2 relative error = 4.8e-01
FastGPLattice
     iter of 5.0e+03 | best loss  | loss       | term1      | term2     
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            0.00e+00 | 6.23e+00   | 6.23e+00   | 1.59e+00   | 3.52e+00  
            5.00e+00 | 5.75e+00   | 5.75e+00   | 4.48e+00   | -3.37e-01 
            1.00e+01 | 5.61e+00   | 5.61e+00   | 3.11e+00   | 7.55e-01  
            1.50e+01 | 5.33e+00   | 5.33e+00   | 3.11e+00   | 2.06e-01  
            2.00e+01 | 4.68e+00   | 4.68e+00   | 3.14e+00   | -1.13e+00 
            2.50e+01 | 4.14e+00   | 4.20e+00   | 2.47e+00   | -1.43e+00 
            3.00e+01 | 3.99e+00   | 3.99e+00   | 3.87e+00   | -3.24e+00 
            3.50e+01 | 3.77e+00   | 3.77e+00   | 4.02e+00   | -3.82e+00 
            4.00e+01 | 3.73e+00   | 3.73e+00   | 3.96e+00   | -3.84e+00 
            4.50e+01 | 3.70e+00   | 3.70e+00   | 3.80e+00   | -3.75e+00 
            5.00e+01 | 3.64e+00   | 3.64e+00   | 3.89e+00   | -3.96e+00 
            5.50e+01 | 3.59e+00   | 3.72e+00   | 2.71e+00   | -2.63e+00 
            6.00e+01 | 3.58e+00   | 3.58e+00   | 3.90e+00   | -4.08e+00 
            6.20e+01 | 3.58e+00   | 3.59e+00   | 4.30e+00   | -4.47e+00 
	l2 relative error = 4.9e-01

In [9]:

pmeans_grad = [None]*len(gps_grad)
pci_lows_grad = [None]*len(gps_grad)
pci_highs_grad = [None]*len(gps_grad)
for i,gp in enumerate(gps_grad):
    print(type(gp).__name__)
    x_next = gp.get_x_next(n=[2**2,2**2])
    gp.add_y_next([f(x_next[0]),fp(x_next[1])])
    gp.fit()
    pmeans_grad[i],_,_,pci_lows_grad[i],pci_highs_grad[i] = gp.post_ci(xticks,confidence=0.95)
    print("\tl2 relative error = %.1e"%(torch.linalg.norm(yticks-pmeans_grad[i][0])/torch.linalg.norm(yticks)))

pmeans_grad = [None]*len(gps_grad) pci_lows_grad = [None]*len(gps_grad) pci_highs_grad = [None]*len(gps_grad) for i,gp in enumerate(gps_grad): print(type(gp).__name__) x_next = gp.get_x_next(n=[2**2,2**2]) gp.add_y_next([f(x_next[0]),fp(x_next[1])]) gp.fit() pmeans_grad[i],_,_,pci_lows_grad[i],pci_highs_grad[i] = gp.post_ci(xticks,confidence=0.95) print("\tl2 relative error = %.1e"%(torch.linalg.norm(yticks-pmeans_grad[i][0])/torch.linalg.norm(yticks)))

StandardGP
     iter of 5.0e+03 | best loss  | loss       | term1      | term2     
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            0.00e+00 | 2.04e+05   | 2.04e+05   | 4.09e+05   | -3.73e+01 
            5.00e+00 | 1.03e+03   | 1.03e+03   | 2.07e+03   | -1.92e+01 
            1.00e+01 | 1.99e+01   | 2.17e+01   | 1.85e+01   | 1.03e+01  
            1.50e+01 | 1.99e+01   | 1.99e+01   | 8.60e+00   | 1.65e+01  
            1.80e+01 | 1.99e+01   | 1.99e+01   | 7.50e+00   | 1.75e+01  
	l2 relative error = 2.2e-01
FastGPDigitalNetB2
     iter of 5.0e+03 | best loss  | loss       | term1      | term2     
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            0.00e+00 | 4.65e+02   | 4.65e+02   | 9.12e+02   | 2.97e+00  
            5.00e+00 | 1.99e+02   | 1.99e+02   | 3.74e+02   | 8.95e+00  
            1.00e+01 | 4.80e+01   | 4.80e+01   | 5.79e+01   | 2.35e+01  
            1.50e+01 | 3.08e+01   | 3.10e+01   | 1.12e+01   | 3.61e+01  
            2.00e+01 | 3.07e+01   | 3.07e+01   | 7.57e+00   | 3.91e+01  
            2.50e+01 | 3.07e+01   | 3.07e+01   | 8.08e+00   | 3.85e+01  
            2.70e+01 | 3.07e+01   | 3.07e+01   | 7.93e+00   | 3.87e+01  
	l2 relative error = 1.9e+01
FastGPLattice
     iter of 5.0e+03 | best loss  | loss       | term1      | term2     
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            0.00e+00 | 2.43e+02   | 2.43e+02   | 4.73e+02   | -2.82e+00 
            5.00e+00 | 6.49e+01   | 6.49e+01   | 1.07e+02   | 8.34e+00  
            1.00e+01 | 2.52e+01   | 2.56e+01   | 4.41e+00   | 3.20e+01  
            1.50e+01 | 2.51e+01   | 2.51e+01   | 6.06e+00   | 2.94e+01  
            2.00e+01 | 2.49e+01   | 2.49e+01   | 7.97e+00   | 2.72e+01  
            2.50e+01 | 2.47e+01   | 2.47e+01   | 7.97e+00   | 2.68e+01  
            3.00e+01 | 2.42e+01   | 2.42e+01   | 7.97e+00   | 2.58e+01  
            3.50e+01 | 2.29e+01   | 2.29e+01   | 7.97e+00   | 2.32e+01  
            4.00e+01 | 2.15e+01   | 2.21e+01   | 3.48e+00   | 2.60e+01  
            4.50e+01 | 2.11e+01   | 2.13e+01   | 1.05e+01   | 1.75e+01  
            5.00e+01 | 2.11e+01   | 2.11e+01   | 7.80e+00   | 1.97e+01  
            5.30e+01 | 2.11e+01   | 2.11e+01   | 7.60e+00   | 1.99e+01  
	l2 relative error = 1.2e-01

Plot¶

In [10]:

fig,ax = pyplot.subplots(nrows=len(gps),ncols=3,sharex=True,sharey=False,figsize=(WIDTH*1.5,WIDTH/len(gps)*3))
ax = ax.reshape((len(gps),3))
for i,gp in enumerate(gps):
    ax[i,0].set_ylabel(type(gp).__name__,fontsize="xx-large")
    ax[i,0].plot(xticks[:,0].cpu(),yticks.cpu(),color="k",linewidth=LINEWIDTH)
    ax[i,0].scatter(gp.x[:,0].cpu(),gp.y.cpu(),color="k",s=MARKERSIZE)
    ax[i,0].plot(xticks[:,0].cpu(),pmeans[i].cpu(),color=colors[i],linewidth=LINEWIDTH)
    ax[i,0].fill_between(xticks[:,0].cpu(),pci_lows[i].cpu(),pci_highs[i].cpu(),color=colors[i],alpha=_alpha)
for i,gp in enumerate(gps_grad):
    ax[i,1].plot(xticks[:,0].cpu(),yticks.cpu(),color="k",linewidth=LINEWIDTH)
    ax[i,2].plot(xticks[:,0].cpu(),ypticks.cpu(),color="k",linewidth=LINEWIDTH)
    ax[i,1].scatter(gp.x[0][:,0].cpu(),gp.y[0].cpu(),color="k",s=MARKERSIZE)
    ax[i,2].scatter(gp.x[1][:,0].cpu(),gp.y[1].cpu(),color="k",s=MARKERSIZE)
    ax[i,1].plot(xticks[:,0].cpu(),pmeans_grad[i][0].cpu(),color=colors[i],linewidth=LINEWIDTH)
    ax[i,2].plot(xticks[:,0].cpu(),pmeans_grad[i][1].cpu(),color=colors[i],linewidth=LINEWIDTH)
    ax[i,1].fill_between(xticks[:,0].cpu(),pci_lows_grad[i][0].cpu(),pci_highs_grad[i][0].cpu(),color=colors[i],alpha=_alpha)
    ax[i,2].fill_between(xticks[:,0].cpu(),pci_lows_grad[i][1].cpu(),pci_highs_grad[i][1].cpu(),color=colors[i],alpha=_alpha)
ax[0,0].set_title(r"$f$ no grad",fontsize="xx-large")
ax[0,1].set_title(r"$f$ with grad",fontsize="xx-large")
ax[0,2].set_title(r"$\mathrm{d} f / \mathrm{d} x$ with grad",fontsize="xx-large")
for j in range(3):
    ax[-1,j].set_xlabel(r"$x$",fontsize="xx-large")
fig.savefig("./gps_deriv.pdf",bbox_inches="tight")

fig,ax = pyplot.subplots(nrows=len(gps),ncols=3,sharex=True,sharey=False,figsize=(WIDTH*1.5,WIDTH/len(gps)*3)) ax = ax.reshape((len(gps),3)) for i,gp in enumerate(gps): ax[i,0].set_ylabel(type(gp).__name__,fontsize="xx-large") ax[i,0].plot(xticks[:,0].cpu(),yticks.cpu(),color="k",linewidth=LINEWIDTH) ax[i,0].scatter(gp.x[:,0].cpu(),gp.y.cpu(),color="k",s=MARKERSIZE) ax[i,0].plot(xticks[:,0].cpu(),pmeans[i].cpu(),color=colors[i],linewidth=LINEWIDTH) ax[i,0].fill_between(xticks[:,0].cpu(),pci_lows[i].cpu(),pci_highs[i].cpu(),color=colors[i],alpha=_alpha) for i,gp in enumerate(gps_grad): ax[i,1].plot(xticks[:,0].cpu(),yticks.cpu(),color="k",linewidth=LINEWIDTH) ax[i,2].plot(xticks[:,0].cpu(),ypticks.cpu(),color="k",linewidth=LINEWIDTH) ax[i,1].scatter(gp.x[0][:,0].cpu(),gp.y[0].cpu(),color="k",s=MARKERSIZE) ax[i,2].scatter(gp.x[1][:,0].cpu(),gp.y[1].cpu(),color="k",s=MARKERSIZE) ax[i,1].plot(xticks[:,0].cpu(),pmeans_grad[i][0].cpu(),color=colors[i],linewidth=LINEWIDTH) ax[i,2].plot(xticks[:,0].cpu(),pmeans_grad[i][1].cpu(),color=colors[i],linewidth=LINEWIDTH) ax[i,1].fill_between(xticks[:,0].cpu(),pci_lows_grad[i][0].cpu(),pci_highs_grad[i][0].cpu(),color=colors[i],alpha=_alpha) ax[i,2].fill_between(xticks[:,0].cpu(),pci_lows_grad[i][1].cpu(),pci_highs_grad[i][1].cpu(),color=colors[i],alpha=_alpha) ax[0,0].set_title(r"$f$ no grad",fontsize="xx-large") ax[0,1].set_title(r"$f$ with grad",fontsize="xx-large") ax[0,2].set_title(r"$\mathrm{d} f / \mathrm{d} x$ with grad",fontsize="xx-large") for j in range(3): ax[-1,j].set_xlabel(r"$x$",fontsize="xx-large") fig.savefig("./gps_deriv.pdf",bbox_inches="tight")

In [ ]: