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]:

torch.set_default_dtype(torch.float64)
device = "cpu"

torch.set_default_dtype(torch.float64) device = "cpu"

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)[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)[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(seqs=qp.DigitalNetB2(1,seed=11,randomize="DS")),
    fastgps.FastGPDigitalNetB2(seqs=qp.DigitalNetB2(1,seed=7,randomize="DS"),alpha=4),
    fastgps.FastGPLattice(seqs=qp.Lattice(1,seed=7),alpha=4),
]

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

In [7]:

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

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

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 | loss       | term1      | term2     
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            0.00e+00 | 4.88e+02   | 9.81e+02   | -1.36e+01 
            5.00e+00 | 5.40e+01   | 1.08e+02   | -7.28e+00 
            1.00e+01 | 5.48e+00   | 1.68e+00   | 1.92e+00  
            1.50e+01 | 4.74e+00   | 2.59e+00   | -4.48e-01 
            2.00e+01 | 4.02e+00   | 4.74e+00   | -4.06e+00 
            2.50e+01 | 3.97e+00   | 4.43e+00   | -3.85e+00 
            3.00e+01 | 3.95e+00   | 3.85e+00   | -3.31e+00 
            3.10e+01 | 3.95e+00   | 4.01e+00   | -3.47e+00 
	l2 relative error = 5.4e-01
FastGPDigitalNetB2
     iter of 5.0e+03 | loss       | term1      | term2     
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            0.00e+00 | 5.62e+00   | 1.55e+00   | 2.35e+00  
            5.00e+00 | 5.08e+00   | 4.37e+00   | -1.56e+00 
            1.00e+01 | 4.99e+00   | 3.43e+00   | -7.98e-01 
            1.50e+01 | 4.92e+00   | 3.61e+00   | -1.13e+00 
            2.00e+01 | 4.85e+00   | 3.32e+00   | -9.69e-01 
            2.50e+01 | 4.74e+00   | 3.95e+00   | -1.83e+00 
            3.00e+01 | 4.68e+00   | 3.38e+00   | -1.37e+00 
            3.50e+01 | 4.65e+00   | 3.92e+00   | -1.98e+00 
            4.00e+01 | 4.57e+00   | 3.92e+00   | -2.13e+00 
            4.50e+01 | 4.42e+00   | 3.93e+00   | -2.44e+00 
            5.00e+01 | 4.24e+00   | 3.93e+00   | -2.80e+00 
            5.50e+01 | 4.29e+00   | 2.98e+00   | -1.76e+00 
            6.00e+01 | 4.21e+00   | 3.83e+00   | -2.75e+00 
	l2 relative error = 4.9e-01
FastGPLattice
     iter of 5.0e+03 | loss       | term1      | term2     
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            0.00e+00 | 1.30e+01   | 1.79e+01   | 6.94e-01  
            5.00e+00 | 8.65e+00   | 4.05e+00   | 5.90e+00  
            1.00e+01 | 8.56e+00   | 4.12e+00   | 5.65e+00  
            1.50e+01 | 8.28e+00   | 3.82e+00   | 5.39e+00  
            2.00e+01 | 7.59e+00   | 3.83e+00   | 4.01e+00  
            2.50e+01 | 6.03e+00   | 4.13e+00   | 5.84e-01  
            3.00e+01 | 5.80e+00   | 4.13e+00   | 1.28e-01  
            3.50e+01 | 5.82e+00   | 4.40e+00   | -1.15e-01 
            4.00e+01 | 5.78e+00   | 4.17e+00   | 4.87e-02  
	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 | loss       | term1      | term2     
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            0.00e+00 | 2.04e+05   | 4.09e+05   | -3.73e+01 
            5.00e+00 | 8.32e+03   | 1.67e+04   | -2.68e+01 
            1.00e+01 | 2.01e+01   | 1.31e+01   | 1.25e+01  
            1.50e+01 | 2.07e+01   | 1.59e+01   | 1.07e+01  
            2.00e+01 | 1.95e+01   | 8.70e+00   | 1.55e+01  
            2.50e+01 | 1.90e+01   | 9.77e+00   | 1.35e+01  
            3.00e+01 | 1.89e+01   | 8.93e+00   | 1.42e+01  
            3.50e+01 | 1.89e+01   | 8.21e+00   | 1.49e+01  
            4.00e+01 | 1.89e+01   | 7.89e+00   | 1.52e+01  
            4.50e+01 | 1.89e+01   | 8.31e+00   | 1.47e+01  
            4.80e+01 | 1.89e+01   | 8.00e+00   | 1.50e+01  
	l2 relative error = 2.4e-01
FastGPDigitalNetB2
     iter of 5.0e+03 | loss       | term1      | term2     
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            0.00e+00 | 2.74e+02   | 5.28e+02   | 4.80e+00  
            5.00e+00 | 1.22e+02   | 2.18e+02   | 1.09e+01  
            1.00e+01 | 3.70e+01   | 3.38e+01   | 2.55e+01  
            1.50e+01 | 2.68e+01   | 1.34e+00   | 3.76e+01  
            2.00e+01 | 2.97e+01   | 8.95e+00   | 3.58e+01  

            2.50e+01 | 2.97e+01   | 6.99e+00   | 3.77e+01  
            2.60e+01 | 2.97e+01   | 6.61e+00   | 3.81e+01  
	l2 relative error = 6.0e+01
FastGPLattice
     iter of 5.0e+03 | loss       | term1      | term2     
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            0.00e+00 | 2.43e+02   | 4.73e+02   | -2.82e+00 
            5.00e+00 | 6.49e+01   | 1.07e+02   | 8.34e+00  
            1.00e+01 | 2.56e+01   | 4.41e+00   | 3.20e+01  
            1.50e+01 | 2.51e+01   | 6.06e+00   | 2.94e+01  
            2.00e+01 | 2.49e+01   | 7.97e+00   | 2.72e+01  
            2.50e+01 | 2.47e+01   | 7.97e+00   | 2.68e+01  
            3.00e+01 | 2.42e+01   | 7.97e+00   | 2.58e+01  
            3.50e+01 | 2.29e+01   | 7.97e+00   | 2.32e+01  
            4.00e+01 | 2.21e+01   | 3.48e+00   | 2.60e+01  
            4.50e+01 | 2.13e+01   | 1.05e+01   | 1.75e+01  
            5.00e+01 | 2.11e+01   | 7.80e+00   | 1.97e+01  
            5.30e+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],yticks,color="k",linewidth=LINEWIDTH)
    ax[i,0].scatter(gp.x[:,0],gp.y,color="k",s=MARKERSIZE)
    ax[i,0].plot(xticks[:,0],pmeans[i],color=colors[i],linewidth=LINEWIDTH)
    ax[i,0].fill_between(xticks[:,0],pci_lows[i],pci_highs[i],color=colors[i],alpha=_alpha)
for i,gp in enumerate(gps_grad):
    ax[i,1].plot(xticks[:,0],yticks,color="k",linewidth=LINEWIDTH)
    ax[i,2].plot(xticks[:,0],ypticks,color="k",linewidth=LINEWIDTH)
    ax[i,1].scatter(gp.x[0][:,0],gp.y[0],color="k",s=MARKERSIZE)
    ax[i,2].scatter(gp.x[1][:,0],gp.y[1],color="k",s=MARKERSIZE)
    ax[i,1].plot(xticks[:,0],pmeans_grad[i][0],color=colors[i],linewidth=LINEWIDTH)
    ax[i,2].plot(xticks[:,0],pmeans_grad[i][1],color=colors[i],linewidth=LINEWIDTH)
    ax[i,1].fill_between(xticks[:,0],pci_lows_grad[i][0],pci_highs_grad[i][0],color=colors[i],alpha=_alpha)
    ax[i,2].fill_between(xticks[:,0],pci_lows_grad[i][1],pci_highs_grad[i][1],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],yticks,color="k",linewidth=LINEWIDTH) ax[i,0].scatter(gp.x[:,0],gp.y,color="k",s=MARKERSIZE) ax[i,0].plot(xticks[:,0],pmeans[i],color=colors[i],linewidth=LINEWIDTH) ax[i,0].fill_between(xticks[:,0],pci_lows[i],pci_highs[i],color=colors[i],alpha=_alpha) for i,gp in enumerate(gps_grad): ax[i,1].plot(xticks[:,0],yticks,color="k",linewidth=LINEWIDTH) ax[i,2].plot(xticks[:,0],ypticks,color="k",linewidth=LINEWIDTH) ax[i,1].scatter(gp.x[0][:,0],gp.y[0],color="k",s=MARKERSIZE) ax[i,2].scatter(gp.x[1][:,0],gp.y[1],color="k",s=MARKERSIZE) ax[i,1].plot(xticks[:,0],pmeans_grad[i][0],color=colors[i],linewidth=LINEWIDTH) ax[i,2].plot(xticks[:,0],pmeans_grad[i][1],color=colors[i],linewidth=LINEWIDTH) ax[i,1].fill_between(xticks[:,0],pci_lows_grad[i][0],pci_highs_grad[i][0],color=colors[i],alpha=_alpha) ax[i,2].fill_between(xticks[:,0],pci_lows_grad[i][1],pci_highs_grad[i][1],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 [ ]: