minor changes

Project.toml

@@ -0,0 +1,5 @@
+[deps]
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

README.md

@@ -32,9 +32,10 @@ It uses a two-stage approach based on a two-metric projected truncated Newton me
 
 ## How to Use
 ### Python Wrapper
-The most convenient way of using the program is with the included Python wrapper `pqdts.py` and supports input matrices ($P$,$F$) from numpy (dense, npz) and scipy (sparse, npy). The Python wrapper only supports the OpenMP-parallelized version because the running the MPI-parallelized version usually requires an adaption to the HPC system where is should run, e.g., the MPI-wrapper of the system like `mpirun`, `mpiexec`, `srun` or others.
+The most convenient way of using the program is with the included Python wrapper `pqdts.py` which supports input matrices ($P$,$F$) from numpy (dense, npz) and scipy (sparse, npy). The Python wrapper only supports the OpenMP-parallelized version because running the MPI-parallelized version usually requires an adaption to the HPC system where is should run, e.g., the MPI-wrapper of the system like `mpirun`, `mpiexec`, `srun` or others.
+
+This wrapper has the following options:
 
-This wrapper has the following options
 ```
 usage: pqdts.py [-h] -P PMATRIX [-F FMATRIX] [-D DMAX] [-t THREADS] -p
                 PQDTSPATH [-o OUTPUT] [-e EPSILON] [-g GAMMA] [-m MAXITER]
@@ -64,7 +65,7 @@ optional arguments:
   -b, --benchmark       benchmark mode, don't write output POVMs
   -v, --verbose         be more verbose
 ```
-The dependencies of the wrapper can be installed with `pip install -r requirements.txt`.
+The dependencies of the wrapper can be installed with `pip3 install -r requirements.txt`.
 
 ### Command Line Arguments
 Without the Python wrapper, the command line arguments of `pqdts_omp.x` and `pqdts_omp_mpi.x` are:
@@ -75,7 +76,7 @@ Without the Python wrapper, the command line arguments of `pqdts_omp.x` and `pqd
 5. number of non-zero elements in $P$
 6. computation mode: 2 for two-metric projected truncated Newton method, 1 for projected gradient method 
 7. maxiter: maximal number of iterations in stages
-8. output: 0 for no output of the POVMs, 0 to disable output of POVMs, 1 to enable output of POVMs after every minimization stage
+8. output: 0 to disable output of the POVMs, 1 to enable output of POVMs after every minimization stage
 9. gamma: regularization parameter $\gamma$
 10. epsilon: value of the convergence criterion
 11. index of stage to start with, i.e., 1 or 2
@@ -91,31 +92,31 @@ For $F$ in addition:
 * `data/F_row.bin` row indices of elements in $F$ in Fortran binary format
 * `data/F_col.bin` column indices of elements in $F$ in Fortran binary format
 * `data/F_data.bin` values of elements in $F$ in Fortran binary format
-For details have a look at the routine `read_sparse_from_python` in `pqdts.f90` andthe Python wrapper `pqdts.py`.
+For details have a look at the routine `read_sparse_from_python` in `pqdts.f90` and the Python wrapper `pqdts.py`.
 
 For reading in existing POVMs as a starting point for the minimization:
 * files of the form `"rank_"+rank+"_oiter"+oiter+".dat"` that were written by pqdts.
 
 ## Outputs
-* std-out: progress information
+* stdout: progress information
 * `"rank_"+rank+"_oiter"+oiter+".dat"`: files containing the POVMs in Fortran file format of rank `rank` after stage `oiter`
 * `"timing_"+rank+".out"` file containing timing information of the seperate stages and timed routines
 * `"status_"+rank+".out"` content of `/proc/self/status` at the end of the program execution to retreive memory usage from high-water mark
 
 # Wigner Functions for High Photon Numbers
-The computation of the Wigner function for very high photon numbers is not possible with conventional implementations because it requires the representability of very large floating-point numbers that exceed the range of 64-bit floating-pint numbers. To be able to use arbitrary-precision floating-point numbers we make use of the multiple dispatch feature of the Julia programming language. With this, the Wigner-function implmentation of QuantumOptics.jl (https://qojulia.org/, functions `wigner` and `_clenshaw` from https://github.com/qojulia/QuantumOptics.jl/blob/v1.0.15/src/phasespace.jl) can be generalized to arbitrary-precision floating-point numbers. Julia uses the GNU MPFR library for the BigFloat data type.
+The computation of the Wigner function for very high photon numbers is not possible with conventional implementations because it requires the representability of very large floating-point numbers that exceed the range of 64-bit double-precision. To be able to use arbitrary-precision floating-point numbers, we make use of the multiple dispatch feature of the Julia programming language. With this, the Wigner-function implementation of QuantumOptics.jl (https://qojulia.org/, functions `wigner` and `_clenshaw` from https://github.com/qojulia/QuantumOptics.jl/blob/v1.0.15/src/phasespace.jl) can be generalized to arbitrary-precision floating-point numbers. Julia uses the GNU MPFR library for the BigFloat data type.
 
 We have optimized the implementation for phase-insensitive detectors, i.e., diagonal density matrices and have trivially parallelized the computation because computations with arbitrary-precision floating-point numbers are much more cumbersome than with double-precision floating-point numbers.
 
 ## How to Build
-A Julia installation is required to run the Wigner-programs. See https://julialang.org/downloads/ for details. The Julia packages that the program depends on can be installed with `julia wigner_init.jl`.
+A Julia installation is required to run the Wigner-programs. See https://julialang.org/downloads/ for details. The Julia packages that the program depends on can be installed with `julia --project=. -e 'using Pkg; Pkg.instantiate()'`.
 
 ## How to Use
 ### Variants
 There are three variants:
-* `wigner_fp64.jl`: using conventional 64-bit floating-point numbers and can be used to explore the need the larger arbitrary-precision floating-point numbers
+* `wigner_fp64.jl`: using conventional 64-bit floating-point numbers, useful for exploring the need for larger arbitrary-precision floating-point numbers
 * `wigner_mpfr.jl`: using arbitrary-precision floating-point numbers
-* `wigner_mpfr_mpi.jl`: using arbitrary-precision floating-point numbers and a trivial MPI-parallelization to scale to be able to many CPU-cores/compute nodes. This variant needs to be run as an MPI program, e.g, `mpirun julia wigner_mpfr_mpi.jl ...`. 
+* `wigner_mpfr_mpi.jl`: using arbitrary-precision floating-point numbers and a trivial MPI-parallelization to scale to many CPU-cores/compute nodes. This variant needs to be run as an MPI program, e.g, `mpirun julia --project=. wigner_mpfr_mpi.jl ...`. 
 
 ### Command Line Arguments
 The three variants take the following command line arguments
@@ -126,8 +127,9 @@ The three variants take the following command line arguments
 5. precision: for `wigner_mpfr.jl` and `wigner_mpfr_mpi.jl` the bit size of the mantissa
 
 ### Input HDF5-File
-The input HDF5 file is expected to have a one-dimensional data set named `n` containing the POVM. It can be created with the python script from $\varPi$ 
-```
+The input HDF5 file is expected to have a one-dimensional data set named `n` containing the POVM. It can be created with the following python script from $\varPi$ 
+
+```python
 import h5py
 #povm (Pi) is a MxN numpy matrix 
 mat=np.asarray(povm)
@@ -140,8 +142,9 @@ hf.close()
 ### Output Files 
 The scripts plot the Wigner functions and store the plot as `"n_"+(n)*"_M_"+(M)+"_precision_"+(precision)+".png"`. For plotting with Python the Wigner function is written as an HDF5-file `n_"+(n)+"_M_"+(M)+"_precision_"+(precision)+".h5", "w")` with group name `wigner` and data sets `x` for the points on the x-axis and `w` for the values of the Wigner function on these points.
 
-This HDF5-file can be read in Python with
-```
+This HDF5-file can be read in Python as follows:
+
+```python
 import h5py
 hf = h5py.File(filename, 'r')
 x=hf['wigner']["x"][:]

pqdts.py

@@ -1,3 +1,4 @@
+#!/usr/bin/env python3
 import numpy as np
 import scipy as sp
 import matplotlib.pyplot as plt

requirements.txt

@@ -1,3 +1,4 @@
+h5py
+matplotlib
 numpy
 scipy
-matplotlib

wigner_fp64.jl

@@ -1,16 +1,21 @@
-#Coypright (c) 2024: Paderborn University, Paderborn Center for Parallel Computing, Robert Schade
-#Based wigner _clenshaw from Quantumoptics.jl (https://github.com/qojulia/QuantumOptics.jl/blob/v1.0.15/src/phasespace.jl) Copyright (c) 2016: Sebastian Kraemer.
-#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+# Copyright (c) 2024: Paderborn University, Paderborn Center for Parallel Computing, Robert Schade
+# 
+# Based on the wigner and _clenshaw functions from Quantumoptics.jl
+# (https://github.com/qojulia/QuantumOptics.jl/blob/v1.0.15/src/phasespace.jl)
+# Copyright (c) 2016: Sebastian Kraemer.
+# 
+# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 #
-#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
 #
-#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 
 using LinearAlgebra
 using HDF5
 using Plots
 
-#wigner and _clenshaw are adapted/optimized from https://github.com/qojulia/QuantumOptics.jl/blob/v1.0.15/src/phasespace.jl
+#wigner and _clenshaw0 are adapted/optimized from https://github.com/qojulia/QuantumOptics.jl/blob/v1.0.15/src/phasespace.jl
 function wigner(rho,N, x)
     _2α = x*sqrt(2)
     abs2_2α = abs2(_2α)

wigner_mpfr.jl

@@ -1,16 +1,21 @@
-#Coypright (c) 2024: Paderborn University, Paderborn Center for Parallel Computing, Robert Schade
-#Based wigner _clenshaw from Quantumoptics.jl (https://github.com/qojulia/QuantumOptics.jl/blob/v1.0.15/src/phasespace.jl) Copyright (c) 2016: Sebastian Kraemer.
-#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+# Copyright (c) 2024: Paderborn University, Paderborn Center for Parallel Computing, Robert Schade
+# 
+# Based on the wigner and _clenshaw functions from Quantumoptics.jl
+# (https://github.com/qojulia/QuantumOptics.jl/blob/v1.0.15/src/phasespace.jl)
+# Copyright (c) 2016: Sebastian Kraemer.
+# 
+# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 #
-#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
 #
-#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 
 using LinearAlgebra
 using HDF5
 using Plots
 
-#wigner and _clenshaw are adapted/optimized from https://github.com/qojulia/QuantumOptics.jl/blob/v1.0.15/src/phasespace.jl
+#wigner and _clenshaw0 are adapted/optimized from https://github.com/qojulia/QuantumOptics.jl/blob/v1.0.15/src/phasespace.jl
 function wigner(rho,N, x)
     _2α = x*sqrt(2)
     abs2_2α = abs2(_2α)

wigner_mpfr_mpi.jl

@@ -1,16 +1,21 @@
-#Coypright (c) 2024: Paderborn University, Paderborn Center for Parallel Computing, Robert Schade
-#Based wigner _clenshaw from Quantumoptics.jl (https://github.com/qojulia/QuantumOptics.jl/blob/v1.0.15/src/phasespace.jl) Copyright (c) 2016: Sebastian Kraemer.
-#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+# Copyright (c) 2024: Paderborn University, Paderborn Center for Parallel Computing, Robert Schade
+# 
+# Based on the wigner and _clenshaw functions from Quantumoptics.jl
+# (https://github.com/qojulia/QuantumOptics.jl/blob/v1.0.15/src/phasespace.jl)
+# Copyright (c) 2016: Sebastian Kraemer.
+# 
+# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 #
-#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
 #
-#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 
 using LinearAlgebra
 using HDF5
 using Plots
 using MPI
 
+#wigner and _clenshaw0 are adapted/optimized from https://github.com/qojulia/QuantumOptics.jl/blob/v1.0.15/src/phasespace.jl
 function wigner(rho,N, x)
     _2α = x*sqrt(2)
     abs2_2α = abs2(_2α)