Feature/package

README.md

@@ -0,0 +1,36 @@
+# lsdensities: Spectral densities from two-point correlation functions
+
+lsdensities is a Python library for operating spectral density reconstruction and
+related studies, by using the HLT (Hansen-Lupo-Tantalo) method or Gaussian processes.
+
+## Authors
+
+Niccolò Forzano, Alessandro Lupo.
+
+## Installation
+
+One can download, build and install the package
+
+```bash
+pip install https://github.com/LupoA/rhos
+```
+
+## Usage
+
+```python
+import lsdensities
+
+# e.g. compute the S matrix
+lsdensities.core.Smatrix_mp(t_max, alpha)
+```
+
+## Contributing
+
+Pull requests are welcome. For major changes, please open an issue first
+to discuss what you would like to change.
+
+Please make sure to update tests as appropriate.
+
+## License
+
+[GPL](https://choosealicense.com/licenses/gpl-3.0/)

exec/test_GPerrorHLT.py → examples/GPerrorHLT.py

@@ -1,7 +1,29 @@
 import sys
 
-sys.path.append("..")
-from importall import *
+
+import rhoUtils as u
+from .rhoUtils import init_precision
+from .rhoUtils import LogMessage
+from .rhoUtils import end
+from .rhoUtils import Obs
+from .rhoUtils import adjust_precision
+from .rhoUtils import Inputs
+from .rhoUtils import *
+from .rhoStat import *
+from .rhoMath import *
+from .core import *
+from .rhoParser import *
+from .transform import *
+from .abw import *
+from .rhoParallelUtils import *
+from .HLT_class import *
+from .GPHLT_class import *
+from .GP_class import *
+from .correlatorUtils import foldPeriodicCorrelator
+from .correlatorUtils import symmetrisePeriodicCorrelator
+from .mpmath import mp, mpf
+from .InverseProblemWrapper import *
+from .plotutils import *
 
 
 def init_variables(args_):
@@ -79,6 +101,11 @@ def main():
     corr.evaluate_covmatrix(plot=False)
     corr.corrmat_from_covmat(plot=False)
 
+    with open(os.path.join(par.outdir,'corrmatrix.txt'), "w") as output:
+        for i in range(par.time_extent):
+            for j in range(par.time_extent):
+                print(i, j, corr.corrmat[i,j], file=output)
+
     #   Make it into a mp sample
     print(LogMessage(), "Converting correlator into mpmath type")
     # mpcorr_sample = mp.matrix(par.num_boot, tmax)
@@ -88,11 +115,10 @@ def main():
     cNorm = mpf(str(corr.central[1] ** 2))
 
     lambdaMax = 1e+4
-
     #   Prepare
     hltParams = AlgorithmParameters(
         alphaA=0,
-        alphaB=-1.99,
+        alphaB=1/2,
         alphaC=+1.99,
         lambdaMax=lambdaMax,
         lambdaStep=lambdaMax/2,
@@ -104,7 +130,7 @@ def main():
     matrix_bundle = MatrixBundle(Bmatrix=corr.mpcov, bnorm=cNorm)
 
     #   Wrapper for the Inverse Problem
-    HLTGP = HLTWGPrapper(
+    HLTGP = HLTGPWrapper(
         par=par, algorithmPar=hltParams, matrix_bundle=matrix_bundle, correlator=corr
     )
     HLTGP.prepareHLT()

exec/testGP.py → examples/onlyGP.py

@@ -1,7 +1,30 @@
 import sys
 
-sys.path.append("..")
-from importall import *
+
+import lsdensities.utils.rhoUtils as u
+from lsdensities.utils.rhoUtils import init_precision
+from lsdensities.utils.rhoUtils import LogMessage
+from lsdensities.utils.rhoUtils import end
+from lsdensities.utils.rhoUtils import Obs
+from lsdensities.utils.rhoUtils import adjust_precision
+from lsdensities.utils.rhoUtils import Inputs
+from lsdensities.utils.rhoUtils import *
+from lsdensities.utils.rhoStat import *
+from lsdensities.utils.rhoMath import *
+from lsdensities.core import *
+from lsdensities.utils.rhoParser import *
+from lsdensities.transform import *
+from lsdensities.abw import *
+from lsdensities.utils.rhoParallelUtils import *
+from lsdensities.HLT_class import *
+from lsdensities.GPHLT_class import *
+from lsdensities.GP_class import *
+from lsdensities.correlator.correlatorUtils import foldPeriodicCorrelator
+from lsdensities.correlator.correlatorUtils import symmetrisePeriodicCorrelator
+from mpmath import mp, mpf
+from lsdensities.InverseProblemWrapper import *
+from lsdensities.plotutils import *
+import lsdensities
 
 read_SIGMA_ = True
 
@@ -11,6 +34,7 @@ def init_variables(args_):
     in_.periodicity = args_.periodicity
     in_.prec = args_.prec
     in_.datapath = args_.datapath
+    in_.kerneltype = args_.kerneltype
     in_.outdir = args_.outdir
     in_.massNorm = args_.mpi
     in_.num_boot = args_.nboot
@@ -90,7 +114,7 @@ def main():
     lambdaMax = 1e+6
 
     #   Prepare
-    hltParams = AlgorithmParameters(
+    hltParams = lsdensities.GP_class.AlgorithmParameters(
         alphaA=0,
         alphaB=1/2,
         alphaC=1.99,
@@ -101,18 +125,18 @@ def main():
         kfactor=0.1,
         lambdaMin=1e-4
     )
-    matrix_bundle = MatrixBundle(Bmatrix=corr.mpcov, bnorm=cNorm)
+    matrix_bundle = lsdensities.GP_class.MatrixBundle(Bmatrix=corr.mpcov, bnorm=cNorm)
 
     #   Wrapper for the Inverse Problem
-    GP = GaussianProcessWrapper(
+    GP = lsdensities.GP_class.GaussianProcessWrapper(
         par=par, algorithmPar=hltParams, matrix_bundle=matrix_bundle, correlator=corr, read_SIGMA=read_SIGMA_
     )
     GP.prepareGP()
 
     #   Run
     GP.run(how_many_alphas=par.Na)
     GP.plotParameterScan(how_many_alphas=par.Na, save_plots=True, plot_live=True)
-    GP.plotRhos(savePlot=True)
+    GP.plothltrho(savePlot=True)
 
     end()
 

exec/testHLT.py → examples/onlyHLT.py

@@ -1,7 +1,30 @@
 import sys
 
-sys.path.append("..")
-from importall import *
+
+import lsdensities.utils.rhoUtils as u
+from lsdensities.utils.rhoUtils import init_precision
+from lsdensities.utils.rhoUtils import LogMessage
+from lsdensities.utils.rhoUtils import end
+from lsdensities.utils.rhoUtils import Obs
+from lsdensities.utils.rhoUtils import adjust_precision
+from lsdensities.utils.rhoUtils import Inputs
+from lsdensities.utils.rhoUtils import *
+from lsdensities.utils.rhoStat import *
+from lsdensities.utils.rhoMath import *
+from lsdensities.core import *
+from lsdensities.utils.rhoParser import *
+from lsdensities.transform import *
+from lsdensities.abw import *
+from lsdensities.utils.rhoParallelUtils import *
+from lsdensities.HLT_class import *
+from lsdensities.GPHLT_class import *
+from lsdensities.GP_class import *
+from lsdensities.correlator.correlatorUtils import foldPeriodicCorrelator
+from lsdensities.correlator.correlatorUtils import symmetrisePeriodicCorrelator
+from mpmath import mp, mpf
+from lsdensities.InverseProblemWrapper import *
+from lsdensities.plotutils import *
+import lsdensities
 
 
 def init_variables(args_):
@@ -10,6 +33,7 @@ def init_variables(args_):
     in_.periodicity = args_.periodicity
     in_.prec = args_.prec
     in_.datapath = args_.datapath
+    in_.kerneltype = args_.kerneltype
     in_.outdir = args_.outdir
     in_.massNorm = args_.mpi
     in_.num_boot = args_.nboot
@@ -95,7 +119,7 @@ def main():
     lambdaMax = 1e+8
 
     #   Prepare
-    hltParams = AlgorithmParameters(
+    hltParams = lsdensities.HLT_class.AlgorithmParameters(
         alphaA=0,
         alphaB=1/2,
         alphaC=+1.99,
@@ -106,18 +130,18 @@ def main():
         kfactor=0.1,
         lambdaMin=1e-6
     )
-    matrix_bundle = MatrixBundle(Bmatrix=corr.mpcov, bnorm=cNorm)
+    matrix_bundle = lsdensities.HLT_class.MatrixBundle(Bmatrix=corr.mpcov, bnorm=cNorm)
 
     #   Wrapper for the Inverse Problem
-    HLT = HLTWrapper(
+    HLT = lsdensities.HLT_class.HLTWrapper(
         par=par, algorithmPar=hltParams, matrix_bundle=matrix_bundle, correlator=corr
     )
     HLT.prepareHLT()
 
     #   Run
     HLT.run(how_many_alphas=par.Na)
     HLT.plotParameterScan(how_many_alphas=par.Na, save_plots=True, plot_live=True)
-    HLT.plotRhos(savePlot=True)
+    HLT.plothltrho(savePlot=True)
 
     end()
 

exec/testHLT_singleAlpha.py → examples/onlyHLT_singleAlpha.py

@@ -1,7 +1,29 @@
 import sys
 
-sys.path.append("..")
-from importall import *
+
+import rhoUtils as u
+from .rhoUtils import init_precision
+from .rhoUtils import LogMessage
+from .rhoUtils import end
+from .rhoUtils import Obs
+from .rhoUtils import adjust_precision
+from .rhoUtils import Inputs
+from .rhoUtils import *
+from .rhoStat import *
+from .rhoMath import *
+from .core import *
+from .rhoParser import *
+from .transform import *
+from .abw import *
+from .rhoParallelUtils import *
+from .HLT_class import *
+from .GPHLT_class import *
+from .GP_class import *
+from .correlatorUtils import foldPeriodicCorrelator
+from .correlatorUtils import symmetrisePeriodicCorrelator
+from .mpmath import mp, mpf
+from .InverseProblemWrapper import *
+from .plotutils import *
 
 eNorm = False
 

exec/runExact.py → examples/runExact.py

@@ -6,10 +6,39 @@
 import matplotlib.pyplot as plt
 import os
 import sys
-sys.path.append("..")
-from importall import *
+
 from mpmath import mpf
 
+
+import numpy as np
+import os
+import matplotlib.pyplot as plt
+
+import lsdensities.utils.rhoUtils as u
+from lsdensities.utils.rhoUtils import init_precision
+from lsdensities.utils.rhoUtils import LogMessage
+from lsdensities.utils.rhoUtils import end
+from lsdensities.utils.rhoUtils import Obs
+from lsdensities.utils.rhoUtils import adjust_precision
+from lsdensities.utils.rhoUtils import Inputs
+from lsdensities.utils.rhoUtils import *
+from lsdensities.utils.rhoStat import *
+from lsdensities.utils.rhoMath import *
+from lsdensities.core import *
+from lsdensities.utils.rhoParser import *
+from lsdensities.transform import *
+from lsdensities.abw import *
+from lsdensities.utils.rhoParallelUtils import *
+from lsdensities.HLT_class import *
+from lsdensities.GPHLT_class import *
+from lsdensities.GP_class import *
+from lsdensities.correlator.correlatorUtils import foldPeriodicCorrelator
+from lsdensities.correlator.correlatorUtils import symmetrisePeriodicCorrelator
+from mpmath import mp, mpf
+from lsdensities.InverseProblemWrapper import *
+from lsdensities.plotutils import *
+import lsdensities
+
 pion_mass = 0.140       # in Gev
 a = 0.4                 # in Gev ^-1 ( 1 fm = 5.068 GeV^-1 )
 a_fm = a / 5.068        # lattice spacing in fm
@@ -26,6 +55,7 @@ def init_variables(args_):
     in_.periodicity = args_.periodicity
     in_.prec = args_.prec
     in_.outdir = args_.outdir
+    in_.kerneltype = args_.kerneltype
     in_.massNorm = args_.mpi
     in_.num_boot = args_.nboot
     in_.sigma = args_.sigma
@@ -70,7 +100,12 @@ def generate(par, seed, espace):
 
     for e_i in range(par.Ne):
         for _n in range(STATES):
-            rhoStrue[e_i] += gauss_fp(peaks_location[_n], espace[e_i], par.sigma, norm="full") * weights[_n]
+            if par.kerneltype == 'FULLNORMGAUSS':
+                rhoStrue[e_i] += gauss_fp(peaks_location[_n], espace[e_i], par.sigma, norm="full") * weights[_n]
+            elif par.kerneltype == 'HALFNORMGAUSS':
+                rhoStrue[e_i] += gauss_fp(peaks_location[_n], espace[e_i], par.sigma, norm="half") * weights[_n]
+            elif par.kerneltype == 'CAUCHY':
+                rhoStrue[e_i] += cauchy(peaks_location[_n], par.sigma, espace[e_i]) * weights[_n]
 
     if par.periodicity == "EXP":
         return exact_correlator, espace, rhoStrue
@@ -119,6 +154,9 @@ def main():
         print(LogMessage(), "Energy [a^-1]", espace[e_i])
         _g_t_estar = h_Et_mp_Eslice(Sinv, par, espace[e_i], alpha_=0)
         rhos[e_i] = y_combine_central_Eslice_mp(_g_t_estar, exact_correlator, par)
+        print('Exact Rho: ', rhoStrue[e_i])
+        print('Reconstructed Rho: ', rhos[e_i])
+
 
     plt.plot(
         espace/a,