output weight while calculating

[zhyiyu]

I know we could access the weight used in a calculation by calling the random method, however, in the calculation I am doing, some input can produce unphysical values and the function will return zero for those values. For example

def sigma(x):
    pt = 54.5
    w = pt * x[0]
    v = pt * x[1]
    if w < 3.0 or v > 51.0:
        return 0.0
    else:
        return w * v

In this case, not all the weight values produced by the random method are being used in the actual calculation. It is hard to pick out the weight from random that are being used in the calculation. Therefore I would suggest to add the feature such that user could print weight along with calculation, like the following

def sigma(x):
    pt = 54.5
    w = pt * x[0]
    v = pt * x[1]
    if w < 3.0 or v > 51.0:
        return 0.0
    else:
        print weight
        return w * v

[gplepage]

I am not sure what you are asking for here. It would help if I knew why you want to print out the weights. What do you want to do with them?

A common situation is where one wants to normalize the integral:

integral(f(x) with constraint) / integral(1 with constraint)

where the constraint limits the integral to a sub-volume as in your example (w<3 or v>51). If this is what you want to do, a very good way is to modify sigma(x) to do two integrals at the same time:

def sigma(x):
pt = 54.5
w = pt * x[0]
v = pt * x[1]
if w < 3.0 or v > 51.0:
    return [0.0, 0.0]
else:
    return [w * v, 1.0]

Then the integrator will return an array, call it result, and the normalized integral is result[0] / result[1]. Here result[1] is the sum of the weights but only for the region where the function is nonzero.

This is like an example in the documentation: see,

https://vegas.readthedocs.io/en/latest/tutorial.html#multiple-integrands-simultaneously

[zhyiyu]

The reason why I would like to return the weight is, they are later used to evaluate other integrals with equation 7.8.20 of

Numerical Recipes in Fortran 77

And since the integrals I am trying to evaluate with the weight are actually parts in function sigma, I would like to use only the weight that give sigma and physical value.

Please let me know if it is still confusing, I can try to write a sample code for this. Thanks!

[gplepage]

The Numerical Recipes version of vegas is a direct copy of the first version I wrote, back in the 1970s. I set it up to pass the weight to the integrand so multiple integrals could be done, as in Equation 7.8.20. I no longer do this because it is more accurate (sometimes much more accurate) to do multiple integrals the way I describe in

https://vegas.readthedocs.io/en/latest/tutorial.html#multiple-integrands-simultaneously

This would have been hard to implement in Fortran 77 but works easily in Python. Have you tried to set up your multiple integrals this way? My last response shows how to do 2 integrals simultaneously, sharing parts of the sigma(x). It sounds like your problem is of this type.

[zhyiyu]

Thanks for your advice, since I am doing multiple other integrals over and over again with different complex parameters, I may not be able to employ your method.

Here is a sample code of what I am trying to do, functions pdf_generator and amplitude_generator return a float. Pretending we can access weight

import numpy
import vegas

def pdf_generator(scale):
    ...
    return pdf

def amplitude_generator(scale):
    ...
    return amplitude

def sigma(x):
    global pretable
    pt = 54.5
    w = pt * x[0]
    v = pt * x[1]
    if w < 3.0 or v > 51.0:
        return 0.0
    else:
        global_pretable['w'].append(w)
        global_pretable['v'].append(v)
        global_pretable['weight'].append(weight)
        pdf_w = pdf_generator(w)
        pdf_v = pdf_generator(v)
        amplitude_w = amplitude_generator(w)
        amplitude_v = amplitude_generator(v)
        global_pretable['amplitude'].append(amplitude_w * amplitude_v)
        sigma = (pdf_w * amplitude_w) + (pdf_v * amplitude_v)
        return sigma

global pretable
pretable = {'w': [], 'v': [], 'weight': [], 'amplitude': []}

integrator = vegas.Integrator([[0.0, 1.0], [0.0, 1.0]])
integrator(sigma, nitn = 10, neval = 1000)

Now with pretable, I will do the following with complex parameters m and n, which are both one dimensional lists

import numpy

m = ...
n = ...

w_power = numpy.exp(numpy.einsum('i, j', - m, numpy.log(pretable['w'])))
v_power = numpy.exp(numpy.einsum('i, j', - n, numpy.log(pretable['v'])))

table = numpy.einsum('il, jl, l, l', w_power, v_power, pretable['q_qp']['h_12'], weight)

The matrix table will be two dimensional and m and n are the different complex parameters passed to the "multiple other integrals".

To summarize, some input values of x will produce unphysical values for function sigma, therefore not all weight should be used. And then since I would do multiple other integrals with different complex parameters, I may not be able to use your approach.

I would love to know if the approach you described is much more accurate. If so, I would try to rewrite with this approach.

[gplepage]

One easy way to do what you want (adds only 3 extra lines) is the following code. It runs the integrator on sigma(x) to adapt the grid and then does an additional pass using integrator.random() to collect w, v, amplitude, and weight values in pretable, but only when sigma(x) is nonzero. I added nonsense code to pdf_generator and amplitude_generator so I could run the script, to make sure it does what it is supposed to.

import numpy
import vegas

def pdf_generator(scale):
    return numpy.exp(-scale/5)

def amplitude_generator(scale):
    return scale

def sigma(x, weight=None):
    global pretable
    pt = 54.5
    w = pt * x[0]
    v = pt * x[1]
    if w < 3.0 or v > 51.0:
        return 0.0
    else:
        pdf_w = pdf_generator(w)
        pdf_v = pdf_generator(v)
        amplitude_w = amplitude_generator(w)
        amplitude_v = amplitude_generator(v)
        sigma = (pdf_w * amplitude_w) + (pdf_v * amplitude_v)
        if weight is not None:
            pretable['w'].append(w)
            pretable['v'].append(v)
            pretable['weight'].append(weight)
            pretable['amplitude'].append(amplitude_w * amplitude_v)
        return sigma

pretable = {'w': [], 'v': [], 'weight': [], 'amplitude': []}

integrator = vegas.Integrator([[0.0, 1.0], [0.0, 1.0]])
results = integrator(sigma, nitn = 10, neval = 1000)

# Add an extra pass through vegas, now giving the weight to sigma
# so it can decide whether or not to store information in pretable.
#
for x, weight in integrator.random():
    sigma(x, weight)

You might want to experiment with rewriting the code to pull your extra integrals inside sigma. As I say it tends to be more accurate. (This is because vegas takes advantage of the hypercubes, as described later in thevegas documentation on random; see vegas as a Random Number Generator.) It also works out the statistical correlations between the results of the different integrals, which can be quite significant if those integrals are combined arithmetically (see brief discussion in the vegas documentation on Multiple Integrands Simultaneously).

[zhyiyu]

Thanks for your reply, I think I will try to implement your method of doing all integrals all together. I tested a simple code but it seems like complex number is not supported

import numpy
import vegas

f = lambda x: x * (1.0 + 1.0j)

integrator = vegas.Integrator([[0.0, 1.0]])
integral = integrator(f, nitn = 10, neval = 1e5)

print integral

And I get a warning

ComplexWarning: Casting complex values to real discards the imaginary part

Do I have to separate the real and imaginary part manually or the vegas library can support complex number?

[gplepage]

vegas does not support complex numbers directly (it is on the list for some day). For complex integrands I find it easiest to do all the work inside the integrand using complex numbers, separating them into real and imaginary at the very last step. For example, if you want your integrand to return a complex-valued array table you could use code like:

def sigma(x):
    ...
    table = ... complex array ...
    return dict(real=np.real(table), imag=np.imag(table))

The result returned by the integrator would be a dictionary, so you reconstruct complex numbers from it using something like:

result = integrator(sigma, neval=10000)
result = result['real'] + 1j * result['imag']

The dictionary could also contain other integrands, in addition to the real and imaginary parts of table.

[zhyiyu]

Thanks you so much for your help! I will try to implement that.

[zhyiyu]

I tried the method of doing multiple integrals at the same time, the problem is, sometimes it is much more accurate, but sometimes it is not as accurate as the method I used in the beginning. This should because I am trying to do too many other integrals all together, which made the Integrator confusing.

[gplepage]

The Integrator should not get confused, but you need to be careful about what it is adapting to. When there are multiple integrands, the integrator adapts to the first integrand in the result. So in a function like

def sigma(x):
    ... calculate real-valued sigma_value as in your original code ...
    ... calculate complex-valued table ...
    return dict(sigma=sigma_value, real=np.real(table), imag=np.imag(table))

the Integrator will adapt to sigma_value, not to the table values, because sigma_value appears first in the dictionary. (I am assuming you are using Python 3; if not you need to use an ordered dictionary like collections.OrderedDict.) This would be analogous to your original plan; it shouldn't be less accurate than the original method. If there is some other combination of sigma_value and table that would be better for the adaptation, form that combination and make it the first entry in the dictionary.

One other point, if calculating table is expensive in the above example, you might want to make two versions of sigma(x), one with table and one without. The latter version would be used to adapt the grid and then the former version used with the already adapted grid. This can save a bit run time. The set up would be as follows:

import vegas
import numpy as np 

def sigma(x):
    ... calculate sigma_value as before ..
    return sigma_value

def sigma_table(x):
    ... calculate sigma_value ...
    ... calculate complex-valued table ...
    return dict(sigma=sigma_value, real=np.real(table), imag=np.imag(table))

integrator = vegas.Integrator(2 * [(0,1)])
# adapt to sigma
result = integrator(sigma, neval=10000, ...)
# now add table, but with an adapted integrator
result = integrator(sigma_table, neval=10000, ...)
sigma_value = result['sigma_value']
table = result['real'] + 1j * result['imag']

[zhyiyu]

Sorry for late reply. In order to construct the table, there are 72 other values I have to return and that makes adapting the integrator particularly hard. I will try to look into the structure of each individuals and see if this procedure can be done within a reasonable amount of time.

Thank you again for your helpful advice and patient comments!

[gplepage]

I am closing this issue because it has become inactive (and the questions asked were answered). Feel free to contact me if you have further questions.