# 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`

.

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

```
qmcpy 1.3
```

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

```
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.

```
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.

```
distribution = Lattice(dimension=2, randomize=True, seed=7)
distribution.gen_samples(n_min=0,n_max=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:

where is the unit hypercube in . The integrand is defined everywhere except at 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 and return a real value, but also that it would be able to take a set of sampling points as rows in a Numpy array of size and return an array with values evaluated at each sampling point. The following examples illustrate this point.

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

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

```
x = 0.01
f(x)
```

```
0.6309573444801932
```

What about an array that represents sampling points in a two-dimensional domain, i.e., ?

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

```
1.001650000560437
```

Now, the function should have returned real values that corresponding to each of the sampling points. Let’s debug our Python function.

```
norm(x)
```

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

```
norm(x, axis = 1)
```

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

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

```
array([1. , 0.1 , 0.23784142])
```

It is. Putting everything together, we have:

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

```
array([1. , 0.63095734, 0.50502242])
```

We have got our proper function definition now.

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

```
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)
```

```
MeanVarData (AccumulateData Object)
solution 0.659
error_bound 0.010
n_total 4400
n 3376
levels 1
time_integrate 0.001
CubMCCLT (StoppingCriterion Object)
abs_tol 0.010
rel_tol 0
n_init 2^(10)
n_max 10000000000
inflate 1.200
alpha 0.010
CustomFun (Integrand Object)
Uniform (TrueMeasure Object)
lower_bound 0
upper_bound 1
IIDStdUniform (DiscreteDistribution Object)
d 1
entropy 7
spawn_key ()
```

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

```
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 by simply changing the input parameter value of dimension for QuickConstruct?

```
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)
```

```
MeanVarData (AccumulateData Object)
solution 0.831
error_bound 0.010
n_total 6659
n 5635
levels 1
time_integrate 0.002
CubMCCLT (StoppingCriterion Object)
abs_tol 0.010
rel_tol 0
n_init 2^(10)
n_max 10000000000
inflate 1.200
alpha 0.010
CustomFun (Integrand Object)
Uniform (TrueMeasure Object)
lower_bound 0
upper_bound 1
IIDStdUniform (DiscreteDistribution Object)
d 2^(1)
entropy 7
spawn_key ()
```

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

```
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")
```