diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..cd11344 --- /dev/null +++ b/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" diff --git a/README.md b/README.md index 3df4444..ea95777 100755 --- a/README.md +++ b/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"][:] diff --git a/pqdts.py b/pqdts.py index e133f33..ceeefd9 100644 --- a/pqdts.py +++ b/pqdts.py @@ -1,3 +1,4 @@ +#!/usr/bin/env python3 import numpy as np import scipy as sp import matplotlib.pyplot as plt diff --git a/requirements.txt b/requirements.txt index f23634b..61769bd 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,4 @@ +h5py +matplotlib numpy scipy -matplotlib diff --git a/wigner_fp64.jl b/wigner_fp64.jl index 7c78a93..a52fc0c 100644 --- a/wigner_fp64.jl +++ b/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α) diff --git a/wigner_init.jl b/wigner_init.jl deleted file mode 100644 index 3a7c4cc..0000000 --- a/wigner_init.jl +++ /dev/null @@ -1,5 +0,0 @@ -import Pkg; Pkg.update() -import Pkg; Pkg.add("Plots") -import Pkg; Pkg.add("HDF5") -import Pkg; Pkg.add("MPI") - diff --git a/wigner_mpfr.jl b/wigner_mpfr.jl index f05d268..c0cf761 100644 --- a/wigner_mpfr.jl +++ b/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α) diff --git a/wigner_mpfr_mpi.jl b/wigner_mpfr_mpi.jl index 5ae3080..a933886 100644 --- a/wigner_mpfr_mpi.jl +++ b/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α)