add pade, fft and gf

CMakeLists.txt

@@ -3,7 +3,7 @@ cmake_minimum_required(VERSION 2.8.4 FATAL_ERROR)
 project(sparse CXX C)
 
 # add common modules from ../common/cmake
-list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/common/cmake)
+list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
 
 if(NOT CMAKE_BUILD_TYPE)
   set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Type of build" FORCE)
@@ -16,19 +16,34 @@ set(CMAKE_CXX_FLAGS_DEBUG "-std=c++11 -D_DEBUG -g -Wall -Wextra -pedantic -Wcast
 set(CMAKE_CXX_FLAGS_RELEASE "-std=c++11 -O3")
 
 #include directories
-include_directories(${CMAKE_SOURCE_DIR}/include ${CMAKE_SOURCE_DIR}/thirdparty/cpplapack/include)
+include_directories(${CMAKE_SOURCE_DIR}/thirdparty/cpplapack/include)
 
-find_package(OpenMP)
-if(OPENMP_FOUND)
-  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
-endif(OPENMP_FOUND)
+#find_package(OpenMP)
+#if(OPENMP_FOUND)
+#  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})
+
+# Build and enable tests
+# testing setup
+# enable_testing() must be called in the top-level CMakeLists.txt before any add_subdirectory() is called.
+option(Testing "Enable testing" ON)
+
+if (Testing)
+  enable_testing(test)
+  add_subdirectory(test/c++)
+  add_subdirectory(test/python)
+endif()
+
 #Find SWIG
 #find_package(SWIG REQUIRED)
 #include(${SWIG_USE_FILE})
@@ -52,6 +67,4 @@ find_package(LAPACK REQUIRED)
 # CPPLAPACK
 #include_directories(${CPPLAPACK_INCLUDE_DIR})
 
-add_subdirectory(src)
-
-#add_subdirectory(python)
+add_subdirectory(c++)

c++/CMakeLists.txt

@@ -0,0 +1,3 @@
+#include directories
+include_directories(include)
+add_subdirectory(src)

c++/include/Gf.h

@@ -0,0 +1,56 @@
+
+#ifndef _GF_H
+#define _GF_H
+
+#include <vector>
+#include <complex>
+#include "pade.h"
+
+class Gf{
+public:
+    Gf(): flag_pade(false) {};
+
+    enum statistics{ FERMION, BOSON };
+
+    // initialize by G(tau). tau=[0:beta), N points (N=2^integer is recommended)
+    void set_Gtau(const std::vector<double> &G_tau, const statistics stat, const double beta, const double tail=1.);
+
+    // initialize by G(iw_n). n=[0:N/2)
+    void set_Giw(const std::vector< std::complex<double> > &G_iw, const statistics stat, const double beta, const double tail=1.);
+
+    int get_size() { return N; };
+    double get_tail() { return tail; };
+    double get_beta() { return beta; };
+    statistics get_statistics() { return stat; };
+
+    // get G(tau)
+    void get_Gtau(std::vector<double> &G_tau);
+
+    // get G(iw)
+    void get_Giw(std::vector< std::complex<double> > &G_iw);
+
+    // get G(w). w=[omega_min:omega_max], n-points
+    // void get_Gomega(std::vector< std::complex<double> > &G_omega, const double omega_min, const double omega_max, const int n);
+
+    // return G(omega)
+    std::complex<double> Gomega(const double omega);
+
+    // return rho(omega)
+    double rho(const double omega);
+
+private:
+    int N;
+    double beta;
+    double tail;
+    statistics stat;
+    std::vector<double> gtau;  // size N
+    std::vector< std::complex<double> > giw;  // size N/2
+
+    bool flag_pade;
+    Pade pade;
+
+    void compute_tau2iw();  // G(tau) --> G(iw)
+    void compute_iw2tau();  // G(iw) --> G(tau)
+};
+
+#endif // _GF_H

c++/include/SVD_matrix.h

@@ -0,0 +1,58 @@
+/* SPM  -  Sparse Modeling tool */
+/* Copyright (C) 2017 Junya Otsuki, Kazuyoshi Yoshimi, Hiroshi Shinaoka, Masayuki Ohzeki*/
+
+/* This program is free software: you can redistribute it and/or modify */
+/* it under the terms of the GNU General Public License as published by */
+/* the Free Software Foundation, either version 3 of the License, or */
+/* (at your option) any later version. */
+
+/* This program is distributed in the hope that it will be useful, */
+/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
+/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the */
+/* GNU General Public License for more details. */
+
+/* You should have received a copy of the GNU General Public License */
+/* along with this program.  If not, see <http://www.gnu.org/licenses/>. */
+/*
+ perform SVD, A = U*S*V^t, and store those matrices
+ SVs smaller than 'SV_min' are dropped.
+ All SVs are retained if SV_min=0.
+*/
+
+#ifndef _SVD_MATRIX_H
+#define _SVD_MATRIX_H
+
+#include <vector>
+#include "cpplapack.h"
+
+class SVD_matrix {
+    friend class admm;
+    friend class admm_svd;
+
+public:
+    SVD_matrix(CPPL::dgematrix A, double SV_min = 0);  // copy the matrix A because it is destroyed in SVD
+    int get_rank();
+
+    CPPL::dgematrix At();  // A^t
+    CPPL::dgematrix At_A();  // A^t * A
+    CPPL::dgematrix At_A_inv();  // (A^t * A)^{-1}
+    CPPL::dgematrix At_A_inv_At();  // (A^t * A)^{-1} * A^t
+    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)
+    void OutputSVD(std::string _file_SVD);
+private:
+    CPPL::dcovector S_temp;
+    CPPL::dgematrix U_temp, VT_temp;
+    // K <= M,N
+    CPPL::dgbmatrix S, S_inv;  // singular values, K*K diagonal matrix
+    CPPL::dgematrix U, VT, VT_full;  // M*K, K*N
+    std::string file_SVD;
+    std::string file_U, file_V;
+};
+
+#endif

c++/include/admm_svd.h

@@ -0,0 +1,106 @@
+/* SPM  -  Sparse Modeling tool */
+/* Copyright (C) 2017 Junya Otsuki, Kazuyoshi Yoshimi, Hiroshi Shinaoka, Masayuki Ohzeki*/
+
+/* This program is free software: you can redistribute it and/or modify */
+/* it under the terms of the GNU General Public License as published by */
+/* the Free Software Foundation, either version 3 of the License, or */
+/* (at your option) any later version. */
+
+/* This program is distributed in the hope that it will be useful, */
+/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
+/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the */
+/* GNU General Public License for more details. */
+
+/* You should have received a copy of the GNU General Public License */
+/* along with this program.  If not, see <http://www.gnu.org/licenses/>. */
+/*
+ compute sparse solution x of the equation y = A.x by solving
+   min_x { (1/2/lambda) ||y-A.x||_2^2 + ||V^t.x||_1 }
+     subject to x_i >= 0 & sum_i x_i = 'sum_x'
+ where V^t is defined by SVD of A:  A = U.W.V^t
+
+  x: N-dimensional vector
+  y: M-dimensional vector
+  A: (M * N) matrix
+  B: (K * N) matrix, where K is number of singular values retained
+*/
+
+#ifndef _ADMM_H
+#define _ADMM_H
+
+#include <vector>
+#include "cpplapack.h"
+
+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> 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 mse, mse_full;
+	double l1_norm, sum_x_calc, negative_weight;
+	int iter;
+};
+
+class admm_svd{
+public:
+	admm_svd(CPPL::dgematrix &A, double SV_min=0);  // A
+	~admm_svd();
+
+	// [optional]
+	void set_sumrule(double sum_x);
+	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_y(const std::vector<double> &y);
+
+	// [optional]
+	void clear_x();
+
+	// [required]
+	int solve(double tolerance=1e-6, int max_iter=1000, int min_iter=10);  // return 0 if converged, and 1 if not converged
+
+	struct admm_result result;
+	struct admm_info info;
+
+	double validate(const std::vector<double> &y);  // return MSE
+
+private:
+	const CPPL::dgematrix A;
+	SVD_matrix *svd;
+
+	bool flag_nonnegative;
+	bool flag_sumrule;
+	double sum_x;
+
+	CPPL::dcovector x, Vx, z1, u1, z2, u2;  // 1 for L1 norm, 2 for non-negativity
+	CPPL::dcovector y, y_sv;
+
+	double regulariz, penalty1, penalty2;
+	bool flag_penalty_auto;
+	static const int PENALTY_UPDATE_INTERVAL;
+
+	int print_level;
+	bool flag_fileout_iter;
+	std::string file_iter;
+
+	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();
+
+	// 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::dcovector Yy, w;
+		CPPL::drovector v_row;
+		double sum_Vw;
+	} pre;
+};
+
+#endif

c++/include/common.h

@@ -0,0 +1,48 @@
+/* SPM  -  Sparse Modeling tool */
+/* Copyright (C) 2017 Junya Otsuki, Kazuyoshi Yoshimi, Hiroshi Shinaoka, Masayuki Ohzeki*/
+
+/* This program is free software: you can redistribute it and/or modify */
+/* it under the terms of the GNU General Public License as published by */
+/* the Free Software Foundation, either version 3 of the License, or */
+/* (at your option) any later version. */
+
+/* This program is distributed in the hope that it will be useful, */
+/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
+/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the */
+/* GNU General Public License for more details. */
+
+/* You should have received a copy of the GNU General Public License */
+/* along with this program.  If not, see <http://www.gnu.org/licenses/>. */
+#ifndef _COMMON_H
+#define _COMMON_H
+
+#include <vector>
+#include "cpplapack.h"
+
+// converter
+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
+
+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
+double sum_vector(const CPPL::dcovector &);  // sum of vector elements
+double innerproduct(const CPPL::dcovector &, const CPPL::dcovector &);  //innerproduct
+
+bool if_negative(const CPPL::dcovector &);  // true if at least one component is x_i<0
+bool if_negative(const CPPL::dgematrix &);
+
+CPPL::dcovector dcovector_all1(int n);  // return a covector with elements all being 1
+CPPL::drovector drovector_all1(int n);  // return a rovector with elements all being 1
+
+CPPL::dcovector positive(const CPPL::dcovector &);  // set negative elements at zero
+CPPL::dgematrix positive(const CPPL::dgematrix &);
+
+// return shrinked matrix
+CPPL::dgematrix low_rank_matrix(CPPL::dgematrix &mat, int m, int n);
+
+// return matrix in which the columns are rearranged according to 'index_col' (it can be smaller than the original column-length)
+CPPL::dgematrix arrange_matrix_col(const CPPL::dgematrix &mat_full, const std::vector<int> &index_col);
+
+#endif

c++/include/errorcode.h

@@ -0,0 +1,23 @@
+/* SPM  -  Sparse Modeling tool */
+/* Copyright (C) 2017 Junya Otsuki, Kazuyoshi Yoshimi, Hiroshi Shinaoka, Masayuki Ohzeki*/
+
+/* This program is free software: you can redistribute it and/or modify */
+/* it under the terms of the GNU General Public License as published by */
+/* the Free Software Foundation, either version 3 of the License, or */
+/* (at your option) any later version. */
+
+/* This program is distributed in the hope that it will be useful, */
+/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
+/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the */
+/* GNU General Public License for more details. */
+
+/* You should have received a copy of the GNU General Public License */
+/* along with this program.  If not, see <http://www.gnu.org/licenses/>. */
+#ifndef _ERRORCODE_H
+#define _ERRORCODE_H
+
+#define ERROR_SUCCESS 0
+#define ERROR_INVALID_SIZE 1
+#define ERROR_INVALID_MATRIX_SIZE 2
+
+#endif

c++/include/fft.h

@@ -0,0 +1,18 @@
+
+
+#ifndef _FFT_H
+#define _FFT_H
+
+#include <vector>
+#include <complex>
+
+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.);
+
+#endif // _FFT_H