Genz Function in Dakota and QMCPy¶
A QMCPy implementation and comparison of Dakota's Genz function
In [1]:
Copied!
from numpy import *
from qmcpy import DigitalNetB2, Genz, Halton, IIDStdUniform, Lattice
import pandas as pd
from matplotlib import pyplot
import tempfile
import os
import subprocess
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({
"text.usetex": True,
"font.family": "serif",
"text.latex.preamble": r"\usepackage{amsmath}\usepackage{amssymb}\newcommand{\bx}{\boldsymbol{x}}"
})
%matplotlib inline
from numpy import *
from qmcpy import DigitalNetB2, Genz, Halton, IIDStdUniform, Lattice
import pandas as pd
from matplotlib import pyplot
import tempfile
import os
import subprocess
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({
"text.usetex": True,
"font.family": "serif",
"text.latex.preamble": r"\usepackage{amsmath}\usepackage{amssymb}\newcommand{\bx}{\boldsymbol{x}}"
})
%matplotlib inline
In [2]:
Copied!
kinds_func = ['oscillatory','corner-peak']
kinds_coeff = [1,2,3]
ds = 2**arange(8)
ns = 2**arange(7,19)
ds
kinds_func = ['oscillatory','corner-peak']
kinds_coeff = [1,2,3]
ds = 2**arange(8)
ns = 2**arange(7,19)
ds
Out[2]:
array([ 1, 2, 4, 8, 16, 32, 64, 128])
In [ ]:
Copied!
# takes about 5.5 min to run for me
ref_sols = {}
print('logging: ',end='',flush=True)
x_full = DigitalNetB2(ds.max(),seed=7).gen_samples(2**22)
for kind_func in kinds_func:
for kind_coeff in kinds_coeff:
tag = '%s.%d'%(kind_func,kind_coeff)
print('%s, '%tag,end='',flush=True)
mu_hats = zeros(len(ds),dtype=float)
for j,d in enumerate(ds):
genz = Genz(IIDStdUniform(d),kind_func=kind_func,kind_coeff=kind_coeff)
y = genz.f(x_full[:,:d])
mu_hats[j] = y.mean()
ref_sols[tag] = mu_hats
print()
ref_sols = pd.DataFrame(ref_sols)
ref_sols['d'] = ds
ref_sols.set_index('d',inplace=True)
ref_sols
# takes about 5.5 min to run for me
ref_sols = {}
print('logging: ',end='',flush=True)
x_full = DigitalNetB2(ds.max(),seed=7).gen_samples(2**22)
for kind_func in kinds_func:
for kind_coeff in kinds_coeff:
tag = '%s.%d'%(kind_func,kind_coeff)
print('%s, '%tag,end='',flush=True)
mu_hats = zeros(len(ds),dtype=float)
for j,d in enumerate(ds):
genz = Genz(IIDStdUniform(d),kind_func=kind_func,kind_coeff=kind_coeff)
y = genz.f(x_full[:,:d])
mu_hats[j] = y.mean()
ref_sols[tag] = mu_hats
print()
ref_sols = pd.DataFrame(ref_sols)
ref_sols['d'] = ds
ref_sols.set_index('d',inplace=True)
ref_sols
Out[ ]:
| oscillatory.1 | oscillatory.2 | oscillatory.3 | corner-peak.1 | corner-peak.2 | corner-peak.3 | |
|---|---|---|---|---|---|---|
| d | ||||||
| 1 | -0.217229 | -0.217229 | -0.217229 | 8.000000e-01 | 0.800000 | 0.800000 |
| 2 | -0.350528 | -0.379658 | -0.217290 | 7.133127e-01 | 0.712088 | 0.719996 |
| 4 | -0.472868 | -0.472147 | -0.223335 | 5.663464e-01 | 0.566334 | 0.589674 |
| 8 | -0.545580 | -0.526053 | -0.276837 | 3.573127e-01 | 0.360014 | 0.402641 |
| 16 | -0.585497 | -0.558159 | -0.390467 | 1.423258e-01 | 0.147353 | 0.185776 |
| 32 | -0.606470 | -0.577987 | -0.492894 | 2.259076e-02 | 0.025678 | 0.038375 |
| 64 | -0.617227 | -0.590764 | -0.556348 | 5.692663e-04 | 0.000879 | 0.001606 |
| 128 | -0.622677 | -0.599348 | -0.591170 | 3.615787e-07 | 0.000001 | 0.000003 |
In [ ]:
Copied!
if os.path.isfile("x_full_dakota.txt"):
x_full_dakota = np.loadtxt("x_full_dakota.txt")
else:
print("please download Dakota's Halton points from https://drive.google.com/uc?id=1ljmpq3w5L4OjjdinAMSLhXBeWGW6EJ3U")
# with tempfile.TemporaryDirectory() as tmp:
# with open(os.path.join(tmp, "dakota.in"), "w") as io:
# io.write(f"environment\
# \ttabular_data\n\
# method\
# \tfsu_quasi_mc halton\
# \t\tsamples = {ns.max()}\
# \toutput silent\n\
# variables\
# \tcontinuous_design = {ds.max()}\
# \tlower_bounds = {' '.join(['0.0' for _ in range(ds.max())])}\
# \tupper_bounds = {' '.join(['1.0' for _ in range(ds.max())])}\n\
# interface\
# \tfork\
# \t\tanalysis_driver = 'dummy'\
# \tbatch\
# \twork_directory named 'work'\n\
# responses\
# \tobjective_functions = 1\
# \tno_gradients\
# \tno_hessians"
# )
# subprocess.run(["dakota", "dakota.in"], cwd=tmp, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
# file = os.listdir(os.path.join(tmp, "work"))[0]
# with open(os.path.join(tmp, "work", file), "r") as io:
# lines = io.readlines()
# x_full_dakota = []
# for n, line in enumerate(lines):
# if f"{ds.max()} variables" in line:
# x_full_dakota.append([float(lines[n + 1 + j].split()[0]) for j in range(ds.max())])
# x_full_dakota = np.vstack(x_full_dakota)
# np.savetxt("data/x_full_dakota.txt",x_full_dakota)
print(x_full_dakota.shape)
if os.path.isfile("x_full_dakota.txt"):
x_full_dakota = np.loadtxt("x_full_dakota.txt")
else:
print("please download Dakota's Halton points from https://drive.google.com/uc?id=1ljmpq3w5L4OjjdinAMSLhXBeWGW6EJ3U")
# with tempfile.TemporaryDirectory() as tmp:
# with open(os.path.join(tmp, "dakota.in"), "w") as io:
# io.write(f"environment\
# \ttabular_data\n\
# method\
# \tfsu_quasi_mc halton\
# \t\tsamples = {ns.max()}\
# \toutput silent\n\
# variables\
# \tcontinuous_design = {ds.max()}\
# \tlower_bounds = {' '.join(['0.0' for _ in range(ds.max())])}\
# \tupper_bounds = {' '.join(['1.0' for _ in range(ds.max())])}\n\
# interface\
# \tfork\
# \t\tanalysis_driver = 'dummy'\
# \tbatch\
# \twork_directory named 'work'\n\
# responses\
# \tobjective_functions = 1\
# \tno_gradients\
# \tno_hessians"
# )
# subprocess.run(["dakota", "dakota.in"], cwd=tmp, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
# file = os.listdir(os.path.join(tmp, "work"))[0]
# with open(os.path.join(tmp, "work", file), "r") as io:
# lines = io.readlines()
# x_full_dakota = []
# for n, line in enumerate(lines):
# if f"{ds.max()} variables" in line:
# x_full_dakota.append([float(lines[n + 1 + j].split()[0]) for j in range(ds.max())])
# x_full_dakota = np.vstack(x_full_dakota)
# np.savetxt("data/x_full_dakota.txt",x_full_dakota)
print(x_full_dakota.shape)
(262144, 128)
In [19]:
Copied!
n_max,d_max = ns.max(),ds.max()
pts = {
'IID Standard Uniform': IIDStdUniform(d_max).gen_samples(n_max),
'Lattice (random shift)': Lattice(d_max).gen_samples(n_max),
'Digital Net (random scramble + shift)': DigitalNetB2(d_max).gen_samples(n_max),
'Halton (QMCPy)': Halton(d_max).gen_samples(n_max,warn=False),
'Halton (Dakota)': x_full_dakota[:n_max,:d_max]
}
n_max,d_max = ns.max(),ds.max()
pts = {
'IID Standard Uniform': IIDStdUniform(d_max).gen_samples(n_max),
'Lattice (random shift)': Lattice(d_max).gen_samples(n_max),
'Digital Net (random scramble + shift)': DigitalNetB2(d_max).gen_samples(n_max),
'Halton (QMCPy)': Halton(d_max).gen_samples(n_max,warn=False),
'Halton (Dakota)': x_full_dakota[:n_max,:d_max]
}
In [20]:
Copied!
nrows = len(ds)
ncols = len(kinds_func)*len(kinds_coeff)
print('logging')
fig,ax = pyplot.subplots(nrows=nrows,ncols=ncols,figsize=(ncols*5,nrows*5),sharey=True,sharex=True)
ax = ax.reshape(nrows,ncols)
colors = pyplot.rcParams['axes.prop_cycle'].by_key()['color'] + ["indigo"]
for v,(name,x_full) in enumerate(pts.items()):
print('%20s d: '%name,end='',flush=True)
for j,d in enumerate(ds):
print('%d, '%d,end='',flush=True)
for i1,kind_func in enumerate(kinds_func):
for i2,kind_coeff in enumerate(kinds_coeff):
i = len(kinds_coeff)*i1+i2
tag = '%s.%d'%(kind_func,kind_coeff)
genz = Genz(IIDStdUniform(d),kind_func=kind_func,kind_coeff=kind_coeff)
y_full = genz.f(x_full[:,:d])
mu_hats = array([y_full[:n].mean() for n in ns],dtype=float)
error = abs(mu_hats-ref_sols.loc[d,tag])
ax[j,i].plot(ns,error,label=name, color=colors[v])
if v==(len(pts)-1): ax[j,i].legend(loc='lower left')
if v>0: continue
ax[j,i].set_xscale('log',base=2)
ax[j,i].set_yscale('log',base=10)
if i==0: ax[j,i].set_ylabel(r'$d=%d$\\$\varepsilon = \lvert \mu - \hat{\mu} \rvert$'%d)
if j==0: ax[j,i].set_title(tag)
if j==(len(ds)-1):
ax[j,i].set_xlabel(r'$n$')
ax[j,i].set_xticks(ns)
ax[j,i].set_xlim([ns.min(),ns.max()])
print()
nrows = len(ds)
ncols = len(kinds_func)*len(kinds_coeff)
print('logging')
fig,ax = pyplot.subplots(nrows=nrows,ncols=ncols,figsize=(ncols*5,nrows*5),sharey=True,sharex=True)
ax = ax.reshape(nrows,ncols)
colors = pyplot.rcParams['axes.prop_cycle'].by_key()['color'] + ["indigo"]
for v,(name,x_full) in enumerate(pts.items()):
print('%20s d: '%name,end='',flush=True)
for j,d in enumerate(ds):
print('%d, '%d,end='',flush=True)
for i1,kind_func in enumerate(kinds_func):
for i2,kind_coeff in enumerate(kinds_coeff):
i = len(kinds_coeff)*i1+i2
tag = '%s.%d'%(kind_func,kind_coeff)
genz = Genz(IIDStdUniform(d),kind_func=kind_func,kind_coeff=kind_coeff)
y_full = genz.f(x_full[:,:d])
mu_hats = array([y_full[:n].mean() for n in ns],dtype=float)
error = abs(mu_hats-ref_sols.loc[d,tag])
ax[j,i].plot(ns,error,label=name, color=colors[v])
if v==(len(pts)-1): ax[j,i].legend(loc='lower left')
if v>0: continue
ax[j,i].set_xscale('log',base=2)
ax[j,i].set_yscale('log',base=10)
if i==0: ax[j,i].set_ylabel(r'$d=%d$\\$\varepsilon = \lvert \mu - \hat{\mu} \rvert$'%d)
if j==0: ax[j,i].set_title(tag)
if j==(len(ds)-1):
ax[j,i].set_xlabel(r'$n$')
ax[j,i].set_xticks(ns)
ax[j,i].set_xlim([ns.min(),ns.max()])
print()
logging
IID Standard Uniform d: 1, 2, 4, 8, 16, 32, 64, 128,
Lattice (random shift) d: 1, 2, 4, 8, 16, 32, 64, 128,
Digital Net (random scramble + shift) d: 1, 2, 4, 8, 16, 32, 64, 128,
Halton (QMCPy) d: 1, 2, 4, 8, 16, 32, 64, 128,
Halton (Dakota) d: 1, 2, 4, 8, 16, 32, 64, 128,
