Genz Function in Dakota and QMCPy

A QMCPy implementation and comparison of Dakota’s Genz function

from numpy import *
from qmcpy import *
import pandas as pd
from matplotlib import pyplot
import tempfile
import os
import subprocess
import numpy as np
pyplot.style.use('../qmcpy/qmcpy.mplstyle')
%matplotlib inline

kinds_func = ['oscilatory','corner-peak']
kinds_coeff = [1,2,3]
ds = 2**arange(8)
ns = 2**arange(7,19)
ds

array([  1,   2,   4,   8,  16,  32,  64, 128])

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

logging: oscilatory.1, oscilatory.2, oscilatory.3, corner-peak.1, corner-peak.2, corner-peak.3,

oscilatory.1 oscilatory.2 oscilatory.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.259072e-02 0.025678 0.038375
64 -0.617228 -0.590764 -0.556348 5.692751e-04 0.000879 0.001606
128 -0.622677 -0.599348 -0.591170 3.614889e-07 0.000001 0.000003
# 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_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:
#         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)

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 (not random, not general)': Halton(d_max,randomize=False,generalize=False).gen_samples(n_max,warn=False),
'Halton (not random, general)': Halton(d_max,randomize=False,generalize=True).gen_samples(n_max,warn=False),
'Halton (random, not general)': Halton(d_max,randomize=True,generalize=False).gen_samples(n_max),
'Halton (random, general)': Halton(d_max,randomize=True,generalize=True).gen_samples(n_max),
'Halton (Dakota)': x_full_dakota[:n_max,:d_max]
}

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 (not random, not general) d: 1, 2, 4, 8, 16, 32, 64, 128,
Halton (not random, general) d: 1, 2, 4, 8, 16, 32, 64, 128,
Halton (random, not general) d: 1, 2, 4, 8, 16, 32, 64, 128,
Halton (random, general) d: 1, 2, 4, 8, 16, 32, 64, 128,
Halton (Dakota) d: 1, 2, 4, 8, 16, 32, 64, 128,