Setting up lsdensities¶
In your directory of choice, clone the repository,
git clone https://github.com/LupoA/lsdensities
You can then install by running:
cd lsdensities
pip install -e .
Checking that it works¶
Moving to the example directory
cd lsdensities/examples/
Run the minimal example,
python minimal_example.py
Alternatively, in the test directory, you can run the tests, which should take seconds.
cd lsdensities/tests/
pytest
Usage examples¶
There are two ways of using this library, depending on how the user wants to interface with the optimisation of the parameter lambda. One way is to use the library to compute the smeared spectral density at a fixed value of lambda. One can loop over different values of lambda and choose the value that is deemed correct. Alternatively, the optimisation of lambda is automatised in the InverseProblemWrapper class. This is the recommended way.
If the class is used, only a handful of inputs and the data must be passed to the library, and the rest is handled automatically. The user still has a large degree of freedom by tuning the input parameters of the class. Moreover, the InverseProblemWrapper class will output the smeared spectral density both from a Bayesian and a frequentist algorithm.
To immediately look the workflow without the InverseProblemWrapper, look at examples/minimal_example.py.
To immediately look at the workflow using the InverseProblemWrapper, look at examples/runInverseProblem.py (requires a datafile).
We can start by looking at some common, basic features. First, an Inputs object must be created. This object specifies some basic inputs, such as the type of the smearing kernel, its smearing radius, or the energies at which we want to obtain the solution. The Input object is ubiquitous in lsdensities and it should be well understood. While its attributed have default values, they should be set manually or via input files. For instance:
parameters = Inputs()
parameters.periodicity = "EXP" # "EXP" or "COSH" depending
# on the lattice having open
# or periodic boundary conditions in time
parameters.sigma = 0.25 # smearing radius in given energy units
parameters.prec = 128 # working precision, in this 128 decimal digits
Some inputs are dimensionful quantities, and the user should take care of using them consistently (e.g. everything is expressed in GeV). Once the desired attributed for the Inputs object are specify, a function needs to called in order to fill some internal attributed based on the values provided.
parameters.assign_values()
An important variable set by this function is parameters.tmax which specifies the number of datapoints that will be actually used.
If unspecified, it uses the maximum values, which is inferred by parameters.time_extent and parameter.periodicity.
The assign_values() function will also fix the working precision for the mpmath library based on parameters.prec. If needed, this can be changed:
num_digits = 64
init_precision(num_digits)
Another way to set up certain parameters is by command line
par = parse_inputs()
par.assign_values()
Warning
running parameters.assign_values() is mandatory, and applications may not work if this function is not called. The user should be very careful modifying attributes of Inputs after assign_values() is called.
A first look: solving against synthetic correlators (no datafile required) at a single value of \(\lambda\)¶
In the following we create a single sample for a correlator and we extract the smeared spectral density at a given energy. This basic application of the library is intended to familiarise the user with its basic feature, but it does not contain all the information that is necessary to run against a real dataset.
First, initialise the precision and the Input object. Since the data will be synthetic, we have to additionally specify the time extent of the lattice, which is normally red from the datafile.
parameters = Inputs()
parameters.time_extent = 32
parameters.kerneltype = "FULLNORMGAUSS" # Kernel smearing spectral density
parameters.periodicity = "EXP" # EXP / COSH for open / periodic boundary conditions
parameters.sigma = 0.25 # smearing radius in given energy units
parameters.prec = 128
parameters.assign_values() # assigns internal variables based on given inputs
We then create the synthetic correlator to serve as input. This must be an mp.matrix type.
The number of rows is the number of data points at which the correlator is computed (parameters.tmax).
The number of columns is the number of samples for the correlator, in this case one. The associated covariance matrix should have matching size.
lattice_correlator = mp.matrix(
parameters.tmax, 1
)
lattice_covariance = mp.matrix(
parameters.tmax
)
We fill the correlator with a simple exponential function decaying according to a values MASS. We also populate the covariance matrix artificially.
MASS = 1 # in the given energy units
for t in range(parameters.tmax): # mock data
lattice_correlator[t] = mp.exp(-mpf(t + 1) * mpf(str(MASS)))
lattice_covariance[t, t] = lattice_correlator[t] * 0.02 # for this quick example
# we set the covariance to diagonal,
# with a value of e.g. 2% of the correlator
Notice that the argument of the exponential defining lattice_correlator is shifted by one, because the correlator at \(t=0\) cannot be used.
The library then provides function to compute the smeared spectral density as a linear combination of the correlators,
The coefficients \(g_t(E,\sigma)\) are compute through a linear system. To this end, we shall define the appropriate matrix:
ill_conditioned_matrix = hlt_matrix(parameters.tmax, alpha=0)
This has to be in general regularised. In this example it would not be necessary since we created our correlator was defined with mpmath to be exact up to 64 digits, but we do it for pedagogical reasons
regularising_parameter = mpf(str(1e-6))
regularised_matrix = ill_conditioned_matrix + (
regularising_parameter * lattice_covariance
)
matrix_inverse = regularised_matrix ** (-1)
The value of regularising_parameter regularises the solution but introduces a bias.
Its effect accounted for with a high degree of automation by the InverseProblemWrapper class.
In this first example, we simply set it to a small value.
Having computed and inverted the appropriate matrix, we can obtain the coefficients by using the coefficients_ssd function
energy = 0.5 # the energy at which we compute the smeared spectral density
coeff = coefficients_ssd(
matrix_inverse,
parameters,
energy,
alpha=0,
)
The result is then computed with the get_ssd_scalar function
result = get_ssd_scalar(
coeff, # linear combination of data
lattice_correlator,
parameters,
)
In this example, the derived solution can be compared with the true value
true_value = gauss_fp(peak, energy, parameters.sigma, norm="full")
The derived solution should approach the true value as regularising_parameter is reduced towards zero and par.time_extent is increased.
Warning
The time argument of the correlator must be shifted by one unit. This is because the correlator in zero cannot enter the reconstruction process. When using a datafile, you do NOT need to remove the correlator at \(t=0\) from it, because this is done automatically by the library.
Scan over \(\lambda\) and automatised workflow¶
under construction. See examples/runInverseProblem.py for an example.
Example applications in the examples directory¶
The file minimal_example.py contains a slightly less verbose of the code reported above.
The file runExact.py contains a similar example where the smeared spectral density is computed at a large number of energies.
The file runInverseProblem.py contains an example of automated workflow using the InverseProblemWrapper class.