Merge pull request #12 from SpM-lab/spm-v2

SpM v2

CMakeLists.txt

@@ -23,15 +23,27 @@ include_directories(${CMAKE_SOURCE_DIR}/thirdparty/cpplapack/include)
 #  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
 #endif(OPENMP_FOUND)
 
-#Find BLAS
-find_package(BLAS REQUIRED)
-
-#Find LAPACK
-find_package(LAPACK REQUIRED)
-
-#Find FFTW3
-find_package(FFTW REQUIRED)
-include_directories(${FFTW_INCLUDE_DIRS})
+if(USE_MKL)
+  if(NOT MKLROOT)
+    set(MKLROOT $ENV{MKLROOT})
+  endif()
+  if(NOT MKLROOT)
+    message(FATAL "MKLROOT is not set")
+  endif()
+
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -mkl")
+  include_directories(${MKLROOT}/include/fftw)
+else()
+  #Find BLAS
+  find_package(BLAS REQUIRED)
+
+  #Find LAPACK
+  find_package(LAPACK REQUIRED)
+
+  #Find FFTW3
+  find_package(FFTW REQUIRED)
+  include_directories(${FFTW_INCLUDE_DIRS})
+endif()
 
 # Build and enable tests
 # testing setup

README.md

@@ -1,9 +1,6 @@
+# SpM
 
-SPM
-====
-[![Build Status](https://travis-ci.org/SpM-lab/SpM.svg?branch=master)](https://travis-ci.org/SpM-lab/SpM)
-
-Sparse Modeling tool for analytical continuation.
+Sparse Modeling tool for analytical continuation from the imaginary-time Green's function to real-requency spectral function.
 
 The algorithm is presented in the article
 
@@ -15,68 +12,65 @@ The algorithm is presented in the article
 
 This package is distributed under GNU General Public License version 3 ([GPL v3](http://www.gnu.org/licenses/gpl-3.0.en.html)).
 
-We kindly ask you to cite the article above
+We kindly ask you to cite the below articles
 in publications that include results obtained using this package.
 
-## Author
-Junya Otsuki, Kazuyoshi Yoshimi, Hiroshi Shinaoka, Masayuki Ohzeki
+* The article about original algorithm
+    * J. Otsuki, M. Ohzeki, H. Shinaoka, K. Yoshimi,  
+    "Sparse modeling approach to analytical continuation of imaginary-time quantum Monte Carlo data"  
+    [Phys. Rev. E 95, 061302(R) (2017).](https://doi.org/10.1103/PhysRevE.95.061302)
+* The article about this package
+    * K. Yoshimi, J. Otsuki, Y. Motoyama, M. Ohzeki, and H. Shinaoka, "SpM: Sparse modeling tool for analytic continuation of imaginary-time Green's function" [Comput. Phys. Commun. 244, 319-323 (2019)](https://www.sciencedirect.com/science/article/pii/S0010465519302103).
+* The article about the SpM-Pade method (please cite if you use)
+    * Y. Motoyama, K. Yoshimi, and J. Otsuki, "Robust analytic continuation combining the advantages of the sparse modeling approach and Pade approximation" [arXiv:2109.08370](https://arxiv.org/abs/2109.08370)
+
+## Authors
+
+Junya Otsuki, Kazuyoshi Yoshimi, Yuichi Motoyama, Hiroshi Shinaoka, Masayuki Ohzeki
 
 ## Requirement
 
 * LAPACK, BLAS
+* FFTW3
 * cpplapack (included in this package)
 
-
 ## How to build
 
 ### Getting the source codes
 
 Download the latest source codes by
 
-	$ git clone https://github.com/j-otsuki/SpM.git spm.src
+    $ git clone https://github.com/j-otsuki/SpM.git spm.src
 
-Then, the source codes are downloaded in the directory ``spm.src``.
+Then, the source codes are downloaded in the directory `spm.src`.
 
 ### Using Cmake
 
-Build with cmake command is done in a separate directory, e.g. ``spm.build``.
+Build with cmake command is done in a separate directory, e.g. `spm.build`.
 Type the following commands:
 
-	$ mkdir spm.build && cd spm.build
-	$ cmake ../spm.src
-	$ make
+    $ mkdir spm.build && cd spm.build
+    $ cmake ../spm.src
+    $ make
 
-Then, the executable file ``SpM.out`` is created in directory ``spm.build/src``.
+Then, the executable file `SpM.out` is created in directory `spm.build/src`.
 
 ## Sample scripts
 
-Some sample data are provided in ``samples`` directory:
+Some sample data are provided in `samples` directory:
 
-* ``samples/fermion``  # sample for fermionic spectrum (data in the article)
-* ``samples/boson``  # sample for bosonic spectrum
+* `samples/fermion`  # sample for fermionic spectrum (data in the article)
+* `samples/boson`  # sample for bosonic spectrum
 
 A script file, `run.sh`, is provided to run through the program.
 Enter into the directory `samples/fermion`, and execute the script by
 
     $ ./run.sh
 
-You may need to change the path to ``SpM.out`` in the script.
-If succeeded, results including graphs in eps format are created in ``output`` directory.
+You may need to change the path to `SpM.out` in the script.
+If succeeded, results including graphs in eps format are created in `output` directory.
 For details, see the document linked below.
 
-
-## Directories
-The configuration of the directories is shown below.
-
-```
-.
-├── c++         # source files of the SpM program
-├── cmake       # cmake files
-├── docs        # source files of the document page
-├── samples     # sample scripts
-├── test        # files for `make test`
-└── thirdparty  # CPPLAPACK headers
-```
-
 ## Official page
+
 The official page of the SpM is [here](https://spm-lab.github.io/SpM/manual/build/html/index.html).  

c++/include/Gf.h

@@ -1,4 +1,3 @@
-
 #ifndef _GF_H
 #define _GF_H
 
@@ -39,7 +38,7 @@ class Gf{
     double rho(const double omega);
 
 private:
-    int N;
+    int N; // == beta/dt
     double beta;
     double tail;
     statistics stat;

c++/include/SVD_matrix.h

@@ -40,10 +40,10 @@ class SVD_matrix {
     void rearrange_col(std::vector<int> &);  // rearrange column of A (actaully VT)
     void print_basis(std::string _file_U, std::string _file_V, int);
 
-    CPPL::dcovector transform_x2sv(const CPPL::dcovector &);  // V^t * x  (to original basis)
-    CPPL::dcovector transform_y2sv(const CPPL::dcovector &);  // U^t * y  (to original basis)
-    CPPL::dcovector transform_sv2x(const CPPL::dcovector &);  // V * x'  (to SV basis)
-    CPPL::dcovector transform_sv2y(const CPPL::dcovector &);  // U * y'  (to SV basis)
+    CPPL::dcovector transform_x2sv(const CPPL::dcovector &);  // V^t * x  (to SV basis)
+    CPPL::dcovector transform_y2sv(const CPPL::dcovector &);  // U^t * y  (to SV basis)
+    CPPL::dcovector transform_sv2x(const CPPL::dcovector &);  // V * x'  (to original basis)
+    CPPL::dcovector transform_sv2y(const CPPL::dcovector &);  // U * y'  (to original basis)
     void OutputSVD(std::string _file_SVD);
 private:
     CPPL::dcovector S_temp;

c++/include/admm_svd.h

@@ -35,12 +35,14 @@ class SVD_matrix;
 struct admm_result{
 	// w/ sv  : results in SV basis
 	// w/o sv : results in omega-tau basis
-	std::vector<double> x, xsv, z1, z1sv, z2, z2sv;
+	std::vector<double> x, xsv, z1, z1sv, z2, z2sv, z3, z3sv;
 	std::vector<double> y, ysv, y_recovered_x, y_recovered_z1, ysv_recovered_x, ysv_recovered_z1;
 };
 struct admm_info{
 	double res1_pri, res1_dual, res2_pri, res2_dual;  // residual errors
+	double res3_pri, res3_dual; // spmpade
 	double mse, mse_full;
+	int l0_norm;
 	double l1_norm, sum_x_calc, negative_weight;
 	int iter;
 };
@@ -52,13 +54,16 @@ class admm_svd{
 
 	// [optional]
 	void set_sumrule(double sum_x);
+	void set_rho_ref(std::vector<double>const &rho_ref, std::vector<double> const & rho_ref_weight);
+	void set_omega_coeff(std::vector<double> const& omega_coeff);
 	void set_nonnegative(bool _flag);
 	void set_fileout_iter(const std::string filename);  // output convergence in a file.  unset if filename=""
 	void set_print_level(int);  // 0: none,  1: results,  2: verbose
 
 	// [required]
-	void set_coef(double lambda, double penalty1=1.0, double penalty2=1.0, bool flag_penalty_auto=false);
+	void set_coef(double lambda, double penalty1=1.0, double penalty2=1.0, double penalty3=1.0, bool flag_penalty_auto=false);
 	void set_y(const std::vector<double> &y);
+	int get_svd_num(){return svd_num;}
 
 	// [optional]
 	void clear_x();
@@ -77,28 +82,42 @@ class admm_svd{
 
 	bool flag_nonnegative;
 	bool flag_sumrule;
+	bool flag_rho_ref;
 	double sum_x;
+	std::vector<double> rho_ref;
+	std::vector<double> rho_ref_weight;
 
-	CPPL::dcovector x, Vx, z1, u1, z2, u2;  // 1 for L1 norm, 2 for non-negativity
+	CPPL::dcovector x, Vx, WVx, z1, u1, z2, u2;  // 1 for L1 norm, 2 for non-negativity
+	CPPL::dcovector z3, u3; // SpM-Pade
 	CPPL::dcovector y, y_sv;
 
 	double regulariz, penalty1, penalty2;
+	double penalty3; // SpMPade
 	bool flag_penalty_auto;
 	static const int PENALTY_UPDATE_INTERVAL;
 
+	std::vector<double> omega_coeff; // omega
+
 	int print_level;
 	bool flag_fileout_iter;
 	std::string file_iter;
 
+	int svd_num;
 	void pre_update();  // This function must be called before calling update_x, and must be recalled when one of values of lambda, penalty1, penalty2 are changed
 	void update_x();
 
+	CPPL::dcovector transform_x2sv(CPPL::dcovector const& v);
+	CPPL::dcovector transform_sv2x(CPPL::dcovector const& v);
+
 	// quantities used in functions update_x and set_y (set in function pre_update)
 	struct quantities_for_update{
-		CPPL::dgbmatrix Y, B;  // diagonal matrix
-		CPPL::dgematrix C;
+		CPPL::dgematrix Y, B;  // diagonal matrix
+		CPPL::dgbmatrix Pade_B3;
+		CPPL::dgematrix Pade_B0, C;
 		CPPL::dcovector Yy, w;
+		CPPL::dcovector Pade_y;
 		CPPL::drovector v_row;
+        CPPL::dcovector rho_w;
 		double sum_Vw;
 	} pre;
 };

c++/include/common.h

@@ -24,6 +24,7 @@ CPPL::dcovector vec2cppl_col(const std::vector<double> &);  // vector -> CPPL
 std::vector<double> cppl2vec(const CPPL::dcovector &);  // CPPL -> vector
 std::vector<double> cppl2vec(const CPPL::drovector &);  // CPPL -> vector
 
+int norm_l0(const CPPL::dcovector &);  // L0 norm
 double norm_l1(const CPPL::dcovector &);  // L1 norm
 double norm_l2(const CPPL::dcovector &);  // L2 norm
 double norm_l2_sq(const CPPL::dcovector &);  // square of L2 norm
@@ -39,6 +40,11 @@ CPPL::drovector drovector_all1(int n);  // return a rovector with elements all b
 CPPL::dcovector positive(const CPPL::dcovector &);  // set negative elements at zero
 CPPL::dgematrix positive(const CPPL::dgematrix &);
 
+// return (v + coeff*ref) / (1+coeff)
+CPPL::dcovector project_ref(const CPPL::dcovector &v,
+                            const std::vector<double> &ref,
+                            const std::vector<double> &coeff); 
+
 // return shrinked matrix
 CPPL::dgematrix low_rank_matrix(CPPL::dgematrix &mat, int m, int n);
 

c++/include/fft.h

@@ -1,16 +1,16 @@
-
-
 #ifndef _FFT_H
 #define _FFT_H
 
 #include <vector>
 #include <complex>
 
+inline bool isodd(int n) { return (n/2)*2 != n; }
+inline bool iseven(int n){ return !isodd(n); }
+
 void fft_fermion_tau2iw(const std::vector<double> &G_tau, std::vector< std::complex<double> > &G_iw, const double beta, const double tail=1.);
 
 void fft_fermion_iw2tau(std::vector<double> &G_tau, const std::vector< std::complex<double> > &G_iw, const double beta, const double tail=1.);
 
-
 void fft_boson_tau2iw(const std::vector<double> &G_tau, std::vector< std::complex<double> > &G_iw, const double beta, const double tail=0.);
 
 void fft_boson_iw2tau(std::vector<double> &G_tau, const std::vector< std::complex<double> > &G_iw, const double beta, const double tail=0.);

c++/include/set_initial.h

@@ -30,7 +30,9 @@ class SetInitial {
     };
     struct FileInfo {
         std::string filein_G;
+        std::string filein_Gsigma;
         int col;
+        int colsigma;
         int output_interval;
         std::string fileout_spec;
         int print_level;
@@ -43,6 +45,7 @@ class SetInitial {
     struct CalcInfo {
         std::string statistics;
         double beta;
+        double sigma;
     };
 
     int argc;
@@ -53,8 +56,6 @@ class SetInitial {
     std::map<std::string, std::string> mapForKeyWordToValue;
     std::map<std::string, bool> mapForKeyWordToRead;
 
-    template<class T>
-    std::string argv_or_defaultvalue(int n, T value);
     bool SetDefaultValue();
     void SetInputValue();
     void RegisterMap(std::string _keyword, std::string _value);
@@ -79,7 +80,6 @@ class SetInitial {
     CalcInfo calcInfo;
     bool AddDefaulValueMap(char* _filename);
     bool ReadParam(char* _filename);
-    bool InputFromArgs(int argc, char *argv[]);
     void PrintInfo();
 };
 

c++/include/spm_core.h

@@ -39,6 +39,8 @@ class SPM_Core {
 		std::vector<double> omega;
 		std::vector<admm_result> result;
 		std::vector<admm_info> info;
+        std::string StatisticsType;
+        double beta;
 
 		double integrate(std::vector<double> &y, double width);
 
@@ -66,23 +68,28 @@ class SPM_Core {
 
 		void GetSpectrum(std::vector<double> &_spectrum);
 
-		void GetResults(std::vector<double> &_vmse, std::vector<double> &_vmse_full, std::vector<double> &_vl1_norm,
+		void GetResults(std::vector<double> &_vmse, std::vector<double> &_vmse_full,
+										std::vector<int> &_vl0_norm, std::vector<double> &_vl1_norm,
 										std::vector<double> &_valid);
 
 		int SolveEquation(
 						std::string _StatisticsType,
 						double _Beta,
 						std::vector<std::vector<double> > &_AIn,
 						std::vector<double> &_Gtau,
-            std::vector<double> &_lambda,
-            std::vector<double> &_omega);
+						std::vector<double> &_Gtau_error,
+						std::vector<double> &_lambda,
+						std::vector<double> &_omega);
 
 		int SolveEquationCore(
 						std::vector<std::vector<double> > &_AIn,
 						std::vector<double> &_Gtau,
 						std::vector<double> &_omega,
-            std::vector<double> &_lambda,
-						const double _sum_G
+						std::vector<double> &_omega_coeff,
+						std::vector<double> &_lambda,
+						const double _sum_G,
+						std::vector<double> &_ref_rho,
+						std::vector<double> &_ref_coeff
 		);
 };