Welcome to QMCPy¶

Importing QMCPy¶

Here we show three different ways to import QMCPy in a Python environment. First, we can import the package qmcpy under the alias qp.

In [1]:

import qmcpy as qp
print(qp.name, qp.__version__)

import qmcpy as qp print(qp.name, qp.__version__)

qmcpy 1.6.3c

Alternatively, we can import individual objects from 'qmcpy' as shown below.

In [2]:

from qmcpy.integrand import *
from qmcpy.true_measure import *
from qmcpy.discrete_distribution import *
from qmcpy.stopping_criterion import *

from qmcpy.integrand import * from qmcpy.true_measure import * from qmcpy.discrete_distribution import * from qmcpy.stopping_criterion import *

Lastly, we can import all objects from the package using an asterisk.

In [3]:

from qmcpy import *

from qmcpy import *

Important Notes¶

IID vs LDS¶

Low discrepancy (LD) sequences such as lattice and Sobol' are not independent like IID (independent identically distributed) points.

The code below generates 4 Sobol samples of 2 dimensions.

In [4]:

distribution = Lattice(dimension=2, randomize=True, seed=7)
distribution.gen_samples(n_min=0,n_max=4)

distribution = Lattice(dimension=2, randomize=True, seed=7) distribution.gen_samples(n_min=0,n_max=4)

Out[4]:

array([[0.04386058, 0.58727432],
       [0.54386058, 0.08727432],
       [0.29386058, 0.33727432],
       [0.79386058, 0.83727432]])

Multi-Dimensional Inputs¶

Suppose we want to create an integrand in QMCPy for evaluating the following integral:

$$\int_{[0,1]^d} \|x\|_2^{\|x\|_2^{1/2}} dx,$$

where $[0,1]^d$ is the unit hypercube in $\mathbb{R}^d$.

The integrand is defined everywhere except at $x=0$ and hence the definite integral is also defined.

The key in defining a Python function of an integrand in the QMCPy framework is that not only it should be able to take one point $x \in \mathbb{R}^d$ and return a real value, but also that it would be able to take a set of $n$ sampling points as rows in a Numpy array of size $n \times d$ and return an array with $n$ values evaluated at each sampling point. The following examples illustrate this point.

In [5]:

from numpy.linalg import norm as norm
from numpy import sqrt, array

from numpy.linalg import norm as norm from numpy import sqrt, array

Our first attempt maybe to create the integrand as a Python function as follows:

In [6]:

def f(x): return norm(x) ** sqrt(norm(x))

def f(x): return norm(x) ** sqrt(norm(x))

It looks reasonable except that maybe the Numpy function norm is executed twice. It's okay for now. Let us quickly test if the function behaves as expected at a point value:

In [7]:

x = 0.01
f(x)

x = 0.01 f(x)

Out[7]:

0.6309573444801932

What about an array that represents $n=3$ sampling points in a two-dimensional domain, i.e., $d=2$?

In [8]:

x = array([[1., 0.], 
           [0., 0.01],
           [0.04, 0.04]])
f(x)

x = array([[1., 0.], [0., 0.01], [0.04, 0.04]]) f(x)

Out[8]:

1.001650000560437

Now, the function should have returned $n=3$ real values that corresponding to each of the sampling points. Let's debug our Python function.

In [9]:

norm(x)

norm(x)

Out[9]:

1.0016486409914407

Numpy's norm(x) is obviously a matrix norm, but we want it to be vector 2-norm that acts on each row of x. To that end, let's add an axis argument to the function:

In [10]:

norm(x, axis = 1)

norm(x, axis = 1)

Out[10]:

array([1.        , 0.01      , 0.05656854])

Now it's working! Let's make sure that the sqrt function is acting on each element of the vector norm results:

In [11]:

sqrt(norm(x, axis = 1))

sqrt(norm(x, axis = 1))

Out[11]:

array([1.        , 0.1       , 0.23784142])

It is. Putting everything together, we have:

In [12]:

norm(x, axis = 1) ** sqrt(norm(x, axis = 1))

norm(x, axis = 1) ** sqrt(norm(x, axis = 1))

Out[12]:

array([1.        , 0.63095734, 0.50502242])

We have got our proper function definition now.

In [13]:

def myfunc(x):
    x_norms = norm(x, axis = 1)
    return x_norms ** sqrt(x_norms)

def myfunc(x): x_norms = norm(x, axis = 1) return x_norms ** sqrt(x_norms)

We can now create an integrand instance with our QuickConstruct class in QMCPy and then invoke QMCPy's integrate function:

In [14]:

dim = 1
abs_tol = .01
integrand = CustomFun(
    true_measure = Uniform(IIDStdUniform(dimension=dim, seed=7)),
    g=myfunc)
solution,data = CubMCCLT(integrand,abs_tol=abs_tol,rel_tol=0).integrate()
print(data)

dim = 1 abs_tol = .01 integrand = CustomFun( true_measure = Uniform(IIDStdUniform(dimension=dim, seed=7)), g=myfunc) solution,data = CubMCCLT(integrand,abs_tol=abs_tol,rel_tol=0).integrate() print(data)

Data (Data)
    solution        0.657
    bound_low       0.648
    bound_high      0.667
    bound_diff      0.019
    n_total         3378
    time_integrate  0.002
CubMCCLT (AbstractStoppingCriterion)
    abs_tol         0.010
    rel_tol         0
    n_init          2^(10)
    n_limit         2^(30)
    inflate         1.200
    alpha           0.010
CustomFun (AbstractIntegrand)
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     1
IIDStdUniform (AbstractIIDDiscreteDistribution)
    d               1
    replications    1
    entropy         7

For our integral, we know the true value. Let's check if QMCPy's solution is accurate enough:

In [15]:

true_sol = 0.658582  # In WolframAlpha: Integral[x**Sqrt[x], {x,0,1}]
abs_tol = data.stopping_crit.abs_tol
qmcpy_error = abs(true_sol - solution)
if qmcpy_error > abs_tol: raise Exception("Error not within bounds")

true_sol = 0.658582 # In WolframAlpha: Integral[x**Sqrt[x], {x,0,1}] abs_tol = data.stopping_crit.abs_tol qmcpy_error = abs(true_sol - solution) if qmcpy_error > abs_tol: raise Exception("Error not within bounds")

It's good. Shall we test the function with $d=2$ by simply changing the input parameter value of dimension for QuickConstruct?

In [16]:

dim = 2
integrand = CustomFun(
    true_measure = Uniform(IIDStdUniform(dimension=dim, seed=7)),
    g = myfunc)
solution2,data2 = CubMCCLT(integrand,abs_tol=abs_tol,rel_tol=0).integrate()
print(data2)

dim = 2 integrand = CustomFun( true_measure = Uniform(IIDStdUniform(dimension=dim, seed=7)), g = myfunc) solution2,data2 = CubMCCLT(integrand,abs_tol=abs_tol,rel_tol=0).integrate() print(data2)

Data (Data)
    solution        0.830
    bound_low       0.820
    bound_high      0.840
    bound_diff      0.020
    n_total         5640
    time_integrate  0.002
CubMCCLT (AbstractStoppingCriterion)
    abs_tol         0.010
    rel_tol         0
    n_init          2^(10)
    n_limit         2^(30)
    inflate         1.200
    alpha           0.010
CustomFun (AbstractIntegrand)
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     1
IIDStdUniform (AbstractIIDDiscreteDistribution)
    d               2^(1)
    replications    1
    entropy         7

Once again, we could test for accuracy of QMCPy with respect to the true value:

In [17]:

true_sol2 = 0.827606  # In WolframAlpha: Integral[Sqrt[x**2+y**2])**Sqrt[Sqrt[x**2+y**2]], {x,0,1}, {y,0,1}]
abs_tol2 = data2.stopping_crit.abs_tol
qmcpy_error2 = abs(true_sol2 - solution2)
if qmcpy_error2 > abs_tol2: raise Exception("Error not within bounds")

true_sol2 = 0.827606 # In WolframAlpha: Integral[Sqrt[x**2+y**2])**Sqrt[Sqrt[x**2+y**2]], {x,0,1}, {y,0,1}] abs_tol2 = data2.stopping_crit.abs_tol qmcpy_error2 = abs(true_sol2 - solution2) if qmcpy_error2 > abs_tol2: raise Exception("Error not within bounds")

In [ ]: