diff --git a/CMakeLists.txt b/CMakeLists.txt index e1cf201..5b533bc 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/README.md b/README.md index c801eb9..58e493b 100755 --- a/README.md +++ b/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). diff --git a/c++/include/Gf.h b/c++/include/Gf.h index 5a200d1..7ff5e2e 100644 --- a/c++/include/Gf.h +++ b/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; diff --git a/c++/include/SVD_matrix.h b/c++/include/SVD_matrix.h index 90f447a..543f3bb 100644 --- a/c++/include/SVD_matrix.h +++ b/c++/include/SVD_matrix.h @@ -40,10 +40,10 @@ class SVD_matrix { void rearrange_col(std::vector &); // 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; diff --git a/c++/include/admm_svd.h b/c++/include/admm_svd.h index f18ecbd..33bf2a9 100644 --- a/c++/include/admm_svd.h +++ b/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 x, xsv, z1, z1sv, z2, z2sv; + std::vector x, xsv, z1, z1sv, z2, z2sv, z3, z3sv; std::vector 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::vectorconst &rho_ref, std::vector const & rho_ref_weight); + void set_omega_coeff(std::vector 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 &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 rho_ref; + std::vector 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 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; }; diff --git a/c++/include/common.h b/c++/include/common.h index 146c5c5..cdf7faf 100644 --- a/c++/include/common.h +++ b/c++/include/common.h @@ -24,6 +24,7 @@ CPPL::dcovector vec2cppl_col(const std::vector &); // vector -> CPPL std::vector cppl2vec(const CPPL::dcovector &); // CPPL -> vector std::vector 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 &ref, + const std::vector &coeff); + // return shrinked matrix CPPL::dgematrix low_rank_matrix(CPPL::dgematrix &mat, int m, int n); diff --git a/c++/include/fft.h b/c++/include/fft.h index 3a95a04..7b1457f 100644 --- a/c++/include/fft.h +++ b/c++/include/fft.h @@ -1,16 +1,16 @@ - - #ifndef _FFT_H #define _FFT_H #include #include +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 &G_tau, std::vector< std::complex > &G_iw, const double beta, const double tail=1.); void fft_fermion_iw2tau(std::vector &G_tau, const std::vector< std::complex > &G_iw, const double beta, const double tail=1.); - void fft_boson_tau2iw(const std::vector &G_tau, std::vector< std::complex > &G_iw, const double beta, const double tail=0.); void fft_boson_iw2tau(std::vector &G_tau, const std::vector< std::complex > &G_iw, const double beta, const double tail=0.); diff --git a/c++/include/set_initial.h b/c++/include/set_initial.h index 39cf76a..220a729 100644 --- a/c++/include/set_initial.h +++ b/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 mapForKeyWordToValue; std::map mapForKeyWordToRead; - template - 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(); }; diff --git a/c++/include/spm_core.h b/c++/include/spm_core.h index 2096c27..757500e 100644 --- a/c++/include/spm_core.h +++ b/c++/include/spm_core.h @@ -39,6 +39,8 @@ class SPM_Core { std::vector omega; std::vector result; std::vector info; + std::string StatisticsType; + double beta; double integrate(std::vector &y, double width); @@ -66,7 +68,8 @@ class SPM_Core { void GetSpectrum(std::vector &_spectrum); - void GetResults(std::vector &_vmse, std::vector &_vmse_full, std::vector &_vl1_norm, + void GetResults(std::vector &_vmse, std::vector &_vmse_full, + std::vector &_vl0_norm, std::vector &_vl1_norm, std::vector &_valid); int SolveEquation( @@ -74,15 +77,19 @@ class SPM_Core { double _Beta, std::vector > &_AIn, std::vector &_Gtau, - std::vector &_lambda, - std::vector &_omega); + std::vector &_Gtau_error, + std::vector &_lambda, + std::vector &_omega); int SolveEquationCore( std::vector > &_AIn, std::vector &_Gtau, std::vector &_omega, - std::vector &_lambda, - const double _sum_G + std::vector &_omega_coeff, + std::vector &_lambda, + const double _sum_G, + std::vector &_ref_rho, + std::vector &_ref_coeff ); }; diff --git a/c++/include/spm_param.h b/c++/include/spm_param.h index cc54e9d..d6936ba 100644 --- a/c++/include/spm_param.h +++ b/c++/include/spm_param.h @@ -16,53 +16,64 @@ #ifndef _SPM_PARAM_HEADER #define _SPM_PARAM_HEADER -class SPM_Param{ - private: - struct Lambda{ - int Nl; - double lbegin; - double lend; - int lvalid; - double dlambda; - Lambda(){ - Nl=1; - lbegin = 1e-1; - lend = 1e+0; - lvalid=0; - dlambda=-1; - } - }; +#include - struct Admm{ - double penalty; - double tolerance; - int max_iter; - bool flag_penalty_auto; - Admm(){ - penalty=10.; - tolerance = 1e-6; - max_iter=1000; - flag_penalty_auto=false; - } - }; +class SPM_Param { +private: + struct Lambda { + int Nl; + double lbegin; + double lend; + int lvalid; + double dlambda; - struct SVD{ - double sv_min; - SVD(){ - sv_min=0; - } - }; - - public: + Lambda() { + Nl = 1; + lbegin = 1e-1; + lend = 1e+0; + lvalid = 0; + dlambda = -1; + } + }; + + struct Admm { + double penalty; + double tolerance; + int max_iter; + bool flag_penalty_auto; + + Admm() { + penalty = 10.; + tolerance = 1e-6; + max_iter = 1000; + flag_penalty_auto = false; + } + }; + + struct SVD { + double sv_min; + SVD() { + sv_min = 0; + } + }; + + struct Pade { + std::string filename; + int nsample; + double eta; + }; + +public: Lambda lambda; Admm admm; SVD svd; - + Pade pade; }; -struct SPM_Flags{ - bool validation; - bool nonnegative; +struct SPM_Flags { + bool validation; + bool nonnegative; + bool refrho; }; - + #endif diff --git a/c++/src/G2spectrum.cpp b/c++/src/G2spectrum.cpp index 587da5a..dd99ecb 100644 --- a/c++/src/G2spectrum.cpp +++ b/c++/src/G2spectrum.cpp @@ -143,26 +143,16 @@ static vector operator*(CPPL::dgematrix &mat, vector &vec) { int main(int argc, char *argv[]) { SetInitial setInitial; int iret = 0; - if (argc > 1) { - if (strcmp(argv[1], "-i") == 0) { - if (argc != 3) iret = -1; - else { - if (setInitial.ReadParam(argv[2]) == false) { + if(argc !=2){ + std::cout << "input file \""< Gtau_samples; // for cross validation vector Gtau_rest; // for cross validation if (!flags.validation) { - Gtau = read_Gtau(setInitial.fileInfo.filein_G.c_str(), setInitial.fileInfo.col); - } else { // cross validation - read_Gtau_postfix(setInitial.fileInfo.filein_G.c_str(), setInitial.fileInfo.col, Gtau_samples); - printf("size=%lu\n", Gtau_samples.size()); - Gtau = average(Gtau_samples); - for (unsigned i = 0; i < Gtau_samples.size(); i++) Gtau_rest.push_back(average_rest(Gtau_samples, i)); + Gtau = read_Gtau(setInitial.fileInfo.filein_G.c_str(), setInitial.fileInfo.col); + } else { + // cross validation + // For Obuchi-Kabashima method + Gtau = read_Gtau(setInitial.fileInfo.filein_G.c_str(), setInitial.fileInfo.col); + + //read_Gtau_postfix(setInitial.fileInfo.filein_G.c_str(), setInitial.fileInfo.col, Gtau_samples); + //printf("size=%lu\n", Gtau_samples.size()); + //Gtau = average(Gtau_samples); + //for (unsigned i = 0; i < Gtau_samples.size(); i++) Gtau_rest.push_back(average_rest(Gtau_samples, i)); + } + + type_gtau Gtau_sigma; + if(setInitial.GetParam().pade.eta != 0.0){ + if(!std::isfinite(setInitial.calcInfo.sigma)){ + if(setInitial.fileInfo.colsigma < 0){ + printf("Neither G_sigma nor column_sigma is set.\n"); + exit(-1); + } + Gtau_sigma = read_Gtau(setInitial.fileInfo.filein_Gsigma.c_str(), + setInitial.fileInfo.colsigma); + }else{ + Gtau_sigma.resize(Gtau.size(), setInitial.calcInfo.sigma); + } + }else{ + Gtau_sigma.resize(Gtau.size(), 0.0); } Kernel kernel; @@ -228,6 +238,7 @@ int main(int argc, char *argv[]) { printf("\nConstructing kernel matrix...\n"); vector > A; vector vmse, vmse_full, l1_norm, valid; + vector l0_norm; int errinfo = 0; //TODO: Make Kernel should be treated as a private function @@ -238,7 +249,7 @@ int main(int argc, char *argv[]) { printf("\nSolving the equation...\n"); int l_valid; - errinfo = spm_core.SolveEquation(statistics, beta, A, Gtau, lambda, omega); + errinfo = spm_core.SolveEquation(statistics, beta, A, Gtau, Gtau_sigma, lambda, omega); if (errinfo != 0) { return errinfo; } @@ -247,14 +258,33 @@ int main(int argc, char *argv[]) { vector spec; spm_core.GetSpectrum(spec); - spm_core.GetResults(vmse, vmse_full, l1_norm, valid); + if (statistics == "boson"){ + for(int i = 0; i < omega.size(); i++){ + double x = exp(fabs(omega[i])*beta); + double factor = (1.0 - x)/ (1.0 + x); + if(fabs(omega[i]) < 1e-15) factor = beta/2.0; + else{ + if (omega[i] > 0) factor *= -1.0; + factor /= omega[i]; + } + spec[i] *= factor; + } + } + spm_core.GetResults(vmse, vmse_full, l0_norm, l1_norm, valid); //--------------------------------------------------------- // Save spectrum { FILE *fp = fopen((outputDir + setInitial.fileInfo.fileout_spec).c_str(), "w"); fprintf(fp, "# lambda=%.3e (l=%d)\n", lambda[l_valid], l_valid); - for (int i = 0; i < setInitial.omegaInfo.NW; i++) { - fprintf(fp, "%.5e %.5e\n", omega[i], spec[i] / d_omega); + if (statistics == "fermion") { + for (int i = 0; i < setInitial.omegaInfo.NW; i++) { + fprintf(fp, "%.5e %.5e\n", omega[i], spec[i] / d_omega); + } + } + else{ + for (int i = 0; i < setInitial.omegaInfo.NW; i++) { + fprintf(fp, "%.5e %.5e %.5e\n", omega[i], spec[i]*omega[i] / d_omega, spec[i] / d_omega); + } } fclose(fp); printf("'%s'\n", setInitial.fileInfo.fileout_spec.c_str()); @@ -266,7 +296,7 @@ int main(int argc, char *argv[]) { string fileout_lambda("lambda_dep.dat"); FILE *fp = fopen((outputDir + fileout_lambda).c_str(), "w"); for (unsigned l = 0; l < lambda.size(); l++) { - fprintf(fp, "%.5e %.5e %.5e %.5e %.5e\n", lambda[l], vmse[l], vmse_full[l], l1_norm[l], valid[l]); + fprintf(fp, "%.5e %.5e %.5e %d %.5e %.5e\n", lambda[l], vmse[l], vmse_full[l], l0_norm[l], l1_norm[l], valid[l]); } fclose(fp); printf("'%s'\n", fileout_lambda.c_str()); diff --git a/c++/src/Gf.cpp b/c++/src/Gf.cpp index 13ae194..3980e2e 100644 --- a/c++/src/Gf.cpp +++ b/c++/src/Gf.cpp @@ -7,11 +7,6 @@ const complex IMAG = complex(0, 1.); void Gf::set_Gtau(const std::vector &G_tau, const statistics stat, const double beta, const double tail) { - if( stat == BOSON ){ - // TODO - cout << "BOSON not yet implemented" << endl; - exit(0); - } gtau = G_tau; // copy N = G_tau.size(); this->stat = stat; @@ -25,17 +20,28 @@ void Gf::set_Gtau(const std::vector &G_tau, const statistics stat, const // initialize by G(iw_n). n=[0:N/2) void Gf::set_Giw(const std::vector< std::complex > &G_iw, const statistics stat, const double beta, const double tail) { - if( stat == BOSON ){ - // TODO - cout << "BOSON not yet implemented" << endl; - exit(0); - } - giw = G_iw; // copy - N = G_iw.size()*2; this->stat = stat; this->beta = beta; this->tail = tail; + if(stat == FERMION){ + int n = G_iw.size(); + N = 2*n; + giw.resize(N); + for(int i=0; i > &G_iw, const statisti void Gf::get_Gtau(std::vector &G_tau) { if( gtau.empty() ) compute_iw2tau(); - G_tau = gtau; // copy + G_tau.assign(gtau.begin(), gtau.end()); } void Gf::get_Giw(std::vector< std::complex > &G_iw) { if( giw.empty() ) compute_tau2iw(); - G_iw = giw; // copy + int n = (stat==FERMION ? N/2 : N/2+1 ); + G_iw.assign(giw.begin(), giw.begin()+n); } // void get_Gomega(std::vector< std::complex > &G_omega, const double omega_min, const double omega_max, const int n) @@ -62,10 +69,16 @@ std::complex Gf::Gomega(const double omega) if( !flag_pade ){ if( giw.empty() ) compute_tau2iw(); - vector< complex > matsubara(N); - // TODO: boson - for(int i=0; i > matsubara(n); + vector< complex > Giw(n); + int offset = stat==BOSON?0:1; + for(int i=0; iget_rank(); // # of SV bases - + svd_num = K; x.resize(K); z1.resize(K); u1.resize(K); z2.resize(N); u2.resize(N); + z3.resize(N); + u3.resize(N); clear_x(); } @@ -62,6 +64,22 @@ void admm_svd::set_sumrule(double sum_x) { flag_sumrule = true; } +void admm_svd::set_rho_ref(std::vector const & rho_ref, std::vector const & rho_ref_weight) { + this->rho_ref.assign(rho_ref.begin(), rho_ref.end()); + this->rho_ref_weight.assign(rho_ref_weight.begin(), rho_ref_weight.end()); + flag_rho_ref = false; + for(std::vector::const_iterator itr=rho_ref_weight.begin(); itr!=rho_ref_weight.end(); ++itr){ + if(*itr > 0.0){ + flag_rho_ref = true; + break; + } + } +} + +void admm_svd::set_omega_coeff(std::vector const & omega_coeff){ + this->omega_coeff = omega_coeff; +} + void admm_svd::set_nonnegative(bool _flag) { flag_nonnegative = _flag; } @@ -78,34 +96,59 @@ void admm_svd::set_fileout_iter(string filename) { void admm_svd::set_print_level(int print_level) { this->print_level = print_level; } -//============================================================================= -void admm_svd::set_coef(double regulariz, double penalty1, double penalty2, bool flag_penalty_auto) { +void admm_svd::set_coef(double regulariz, + double penalty1, double penalty2, double penalty3, + bool flag_penalty_auto) { this->regulariz = regulariz; this->penalty1 = penalty1; - this->penalty2 = flag_nonnegative ? penalty2 : 0; + this->penalty2 = (flag_nonnegative ? penalty2 : 0.0); + // this->penalty2 = (flag_nonnegative || flag_rho_ref) ? penalty2 : 0; + this->penalty3 = (flag_rho_ref ? penalty3 : 0.0); this->flag_penalty_auto = flag_penalty_auto; - pre_update(); } void admm_svd::pre_update() { int K = svd->get_rank(); + int N = A.n; - CPPL::dgbmatrix diag(K, K, 0, 0); // diagonal matrix D^{-1} - for (int i = 0; i < K; i++) diag(i, i) = 1. / (penalty1 + penalty2 + pow(svd->S(i, i), 2) / regulariz); + CPPL::dgematrix A_(K, K); + if(flag_rho_ref){ + CPPL::dgbmatrix eta(N, N, 0, 0); + for(int i=0; iVT * eta * CPPL::t(svd->VT); + }else{ + for(int i=0; i < K; i++) for(int j=0; j < K; j++) A_(i,j) = 0.0; + } + for (int i = 0; i < K; i++) A_(i, i) += penalty1 + penalty2 + pow(svd->S(i, i), 2) / regulariz; + CPPL::dgematrix Ainv = CPPL::i(A_); + // CPPL::dgbmatrix diag(K, K, 0, 0); // diagonal matrix D^{-1} + // for (int i = 0; i < K; i++) diag(i, i) = 1. / (penalty1 + penalty2 + penalty3 + pow(svd->S(i, i), 2) / regulariz); - pre.Y = diag * svd->S / regulariz; - pre.B = diag * penalty1; - pre.C = diag * svd->VT * penalty2; + pre.Y = Ainv * svd->S / regulariz; + pre.B = Ainv * penalty1; + pre.C = Ainv * svd->VT * penalty2; CPPL::dcovector p = dcovector_all1(A.n); - pre.w = diag * svd->VT * p; // D^{-1} * V^t * p + pre.w = Ainv * svd->VT * p; // D^{-1} * V^t * p pre.v_row = CPPL::t(p) * CPPL::t(svd->VT); // p^t * V pre.sum_Vw = pre.v_row * pre.w; // < V*w > = p^t * V * w + // SpM-Pade + CPPL::dgbmatrix W(N,N, 0, 0); + for (int i=0; i < N; i++) W(i,i) = omega_coeff[i]; + pre.Pade_B0 = Ainv * svd->VT * W * penalty3; + for (int i=0; i < N; i++) W(i,i) = 1.0 / (rho_ref_weight[i]/regulariz + penalty3); + pre.Pade_B3 = W * penalty3; + pre.Pade_y = W * vec2cppl_col(rho_ref); + for (int i=0; i < N; i++) pre.Pade_y(i) *= rho_ref_weight[i]/regulariz; + { // check matrix size - int M = A.m, N = A.n; + int M = A.m; + int N = A.n; assert(pre.Y.m == K && pre.Y.n == K); // Y: K*K assert(pre.B.m == K && pre.B.n == K); // B: K*K assert(pre.C.m == K && pre.C.n == N); // C: K*N @@ -127,6 +170,8 @@ void admm_svd::clear_x() { u1.zero(); z2.zero(); u2.zero(); + z3.zero(); + u3.zero(); } //============================================================================= @@ -151,23 +196,37 @@ void admm_svd::update_x() { // update x x = pre.Yy; x += pre.B * (z1 - u1); - if (flag_nonnegative) x += pre.C * (z2 - u2); + if (flag_nonnegative) x += pre.C * (z2 - u2); if (flag_sumrule) { - double nu = (sum_x - pre.v_row * x) / pre.sum_Vw; // pre.v_row*x　=　 + double nu = (sum_x - pre.v_row * x ) / pre.sum_Vw; // pre.v_row*x = x += nu * pre.w; } - Vx = svd->transform_sv2x(x); // V x' + if(flag_rho_ref) x += pre.Pade_B0 * (z3-u3); // update z1 and u1 z1 = soft_threshold(1. / penalty1, x + u1); u1 += x - z1; + Vx = svd->transform_sv2x(x); // V x' + // update z2 and u2 if (flag_nonnegative) { - z2 = positive(Vx + u2); // projection + z2 = positive(Vx + u2); // projection to non-negative u2 += Vx - z2; } + + // update z3 and u3 + if (flag_rho_ref){ + const size_t nw = omega_coeff.size(); + WVx = svd->transform_sv2x(x); // V x' + for(size_t i=0; iS * x); // MSE computed in SV basis + info.l0_norm = norm_l0(z1); info.l1_norm = norm_l1(x); // sum_x_calc = sum_vector(Vx); info.sum_x_calc = pre.v_row * x; // < V*x > @@ -244,6 +310,7 @@ int admm_svd::solve(double tolerance, int max_iter, int min_iter) { if (flag_fileout_iter && iter) { fprintf(fp, "%d %.6e", iter, diff_x); fprintf(fp, " %.6e %.6e %.6e %.6e", info.res1_pri, info.res1_dual, info.res2_pri, info.res2_dual); + fprintf(fp, " %.6e %.6e", info.res3_pri, info.res3_dual); fprintf(fp, " %.6e %.6e %.6e %.6e", info.mse, info.l1_norm, info.sum_x_calc, info.negative_weight); fprintf(fp, "\n"); } @@ -259,6 +326,9 @@ int admm_svd::solve(double tolerance, int max_iter, int min_iter) { if (flag_nonnegative) { flag_update |= update_penalty(info.res2_pri, info.res2_dual, penalty2, u2); } + if (flag_rho_ref) { + flag_update |= update_penalty(info.res3_pri, info.res3_dual, penalty3, u3); + } if (flag_update) { // set_matrix_SVD(); pre_update(); @@ -278,6 +348,7 @@ int admm_svd::solve(double tolerance, int max_iter, int min_iter) { // Ref: Eq.(3.12) in Boyd et al., Foundations and Trends in Machine Learning Vol. 3, No. 1 (2010) 1–122 flag_converge = if_converge(info.res1_pri / sqrt(z1.l), tolerance); flag_converge &= if_converge(info.res2_pri / sqrt(x.l), tolerance); + flag_converge &= if_converge(info.res3_pri / sqrt(x.l), tolerance); //--------------------------------------------------------------------- iter++; @@ -300,6 +371,8 @@ int admm_svd::solve(double tolerance, int max_iter, int min_iter) { printf(" | (dual) = %.3e\n", info.res1_dual); printf(" | x>=0 (pri) = %.3e\n", info.res2_pri); printf(" | (dual) = %.3e\n", info.res2_dual); + printf(" | Pade (pri) = %.3e\n", info.res3_pri); + printf(" | (dual) = %.3e\n", info.res3_dual); printf(" |--- target quantities\n"); printf(" | MSE(y-A.x) = %.3e (computed in SV basis)\n", info.mse); printf(" | MSE(y-A.x) = %.3e (computed in tau-omega basis)\n", info.mse_full); @@ -311,6 +384,7 @@ int admm_svd::solve(double tolerance, int max_iter, int min_iter) { printf(" |--- penalty parameters\n"); printf(" | for L1 regul. = %.3lf\n", penalty1); printf(" | for x>=0 = %.3lf\n", penalty2); + printf(" | for Pade = %.3lf\n", penalty3); printf(" |---\n"); } @@ -321,6 +395,8 @@ int admm_svd::solve(double tolerance, int max_iter, int min_iter) { result.z1sv = cppl2vec(z1); result.z2 = cppl2vec(z2); result.z2sv = cppl2vec(svd->transform_x2sv(z2)); // V^t * z2 + result.z3 = cppl2vec(z3); // z3 + result.z3sv = cppl2vec(svd->transform_x2sv(z3)); result.y = cppl2vec(y); result.ysv = cppl2vec(y_sv); diff --git a/c++/src/common.cpp b/c++/src/common.cpp index 0e70efa..6e0a8cb 100644 --- a/c++/src/common.cpp +++ b/c++/src/common.cpp @@ -14,6 +14,7 @@ /* You should have received a copy of the GNU General Public License */ /* along with this program. If not, see . */ #include +#include #include "common.h" using namespace std; @@ -44,6 +45,17 @@ vector cppl2vec(const CPPL::drovector &u) { return v; } +// L0 norm +int norm_l0(const CPPL::dcovector &v) { + int s = 0; + for (int i = 0; i < v.l; i++) { + if(fabs(v(i)) > pow(10.0, -15)) s++; + } + return s; +} + + + // L1 norm double norm_l1(const CPPL::dcovector &v) { double s = 0; @@ -125,6 +137,21 @@ CPPL::dgematrix positive(const CPPL::dgematrix &A) { return B; } +CPPL::dcovector project_ref(const CPPL::dcovector &v, + const std::vector &ref, + const std::vector &coeff) { + CPPL::dcovector u(v.l); + for (int i = 0; i < v.l; i++){ + /* + u(i) = v(i) + coeff[i]*ref[i]; + u(i) /= 1.0+coeff[i]; + */ + u(i) = (1.0-coeff[i])*v(i) + coeff[i] * ref[i]; + } + return u; +} + + // return shrinked matrix CPPL::dgematrix low_rank_matrix(CPPL::dgematrix &mat, int m, int n) { assert(m <= mat.m); diff --git a/c++/src/fft.cpp b/c++/src/fft.cpp index 1a555f6..bd234bb 100644 --- a/c++/src/fft.cpp +++ b/c++/src/fft.cpp @@ -3,104 +3,180 @@ #include #include #include +#include using namespace std; const complex IMAG = complex(0, 1.); -static void fft_tau2iw(const std::vector &g_tau, std::vector< std::complex > &g_iw, const double beta) +void fft_fermion_tau2iw(const std::vector &G_tau, + std::vector< std::complex > &G_iw, + const double beta, + const double tail) { - int N = g_tau.size(); + // G_tau excludes G(beta) - fftw_complex *in, *out; - fftw_plan p; + const int N = G_tau.size(); + const int N2 = 2*N; + G_iw.resize(N); - in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * 2*N); - out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * 2*N); - p = fftw_plan_dft_1d(2*N, in, out, FFTW_BACKWARD, FFTW_ESTIMATE); + fftw_complex* in = static_cast(fftw_malloc(sizeof(fftw_complex)*N2)); + fftw_complex* out = static_cast(fftw_malloc(sizeof(fftw_complex)*N2)); + fftw_plan p = fftw_plan_dft_1d(N2, in, out, FFTW_BACKWARD, FFTW_ESTIMATE); + // anti-periodic for(int i=0; i(out[n_odd][0], out[n_odd][1]) / double(2*N) * beta; + const int n = 2*i+1; + double iw = n*pib; + G_iw[i] = coeff*std::complex(out[n][0], out[n][1]); + G_iw[i] += tail / (IMAG*iw); + G_iw[N-i-1] = conj(G_iw[i]); + } + if(isodd(N)){ + G_iw[N/2] = coeff*std::complex(out[N][0], out[N][1]); + G_iw[N/2] += tail / (pib*N*IMAG); } - fftw_destroy_plan(p); - fftw_free(in); fftw_free(out); + fftw_free(in); + fftw_destroy_plan(p); } -static void fft_iw2tau(std::vector &g_tau, const std::vector< std::complex > &g_iw, const double beta) +void fft_fermion_iw2tau(std::vector &G_tau, + const std::vector< std::complex > &G_iw, + const double beta, + const double tail) { - int N = g_iw.size()*2; + int N = G_iw.size(); + const int N2 = 2*N; + G_tau.resize(N); - fftw_complex *in, *out; - fftw_plan p; + fftw_complex* in = static_cast(fftw_malloc(sizeof(fftw_complex) * N2)); + fftw_complex* out = static_cast(fftw_malloc(sizeof(fftw_complex) * N2)); + fftw_plan p = fftw_plan_dft_1d(N2, in, out, FFTW_FORWARD, FFTW_ESTIMATE); - in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * 2*N); - out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * 2*N); - p = fftw_plan_dft_1d(2*N, in, out, FFTW_FORWARD, FFTW_ESTIMATE); - for(int i=0; i<2*N; i++) in[i][0] = in[i][1] = 0; + for(int i=0; i &G_tau, std::vector< std::complex > &G_iw, const double beta, const double tail) +static void gen_g0_iw_boson(std::vector< std::complex > &g0_iw, const double beta, const double tail) { - double N = G_tau.size(); - G_iw.resize(N/2); - vector G_tau_dif(N); + int N = g0_iw.size(); + { + // g0_iw[0] = tail * beta / 2.; + g0_iw[0] = 0; + for(int i=1; i &G_tau, + std::vector< std::complex > &G_iw, + const double beta, + const double tail) +{ + // G_tau excludes G(beta) + int N = G_tau.size(); + G_iw.resize(N); + + fftw_complex *in = static_cast(fftw_malloc(sizeof(fftw_complex) * N)); + fftw_complex *out = static_cast(fftw_malloc(sizeof(fftw_complex) * N)); + fftw_plan p = fftw_plan_dft_1d(N, in, out, FFTW_BACKWARD, FFTW_ESTIMATE); for(int i=0; i > g0_iw(N); // tail + gen_g0_iw_boson(g0_iw, beta, tail); + + for(int i=0; i(out[i][0], out[i][1]) * (beta/N) + g0_iw[i]; } + + fftw_destroy_plan(p); + fftw_free(out); + fftw_free(in); } -void fft_fermion_iw2tau(std::vector &G_tau, const std::vector< std::complex > &G_iw, const double beta, const double tail) +void fft_boson_iw2tau(std::vector &G_tau, const std::vector< std::complex > &G_iw, const double beta, const double tail) { - double N = G_iw.size()*2; + int N = G_iw.size(); G_tau.resize(N); - vector< complex > G_iw_dif(N/2); - for(int i=0; i(fftw_malloc(sizeof(fftw_complex) * N)); + fftw_complex *out = static_cast(fftw_malloc(sizeof(fftw_complex) * N)); + fftw_plan p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE); + + std::vector< complex > g0_iw(N); // tail + gen_g0_iw_boson(g0_iw, beta, tail); + + for(int i=0; i > &z, vector > &u) swap(g0, g1); } - delete g0, g1; + delete g1; + delete g0; // truncate the coefficients m_a if m_a[n]=0 for(int i=0; i #include -template -std::string SetInitial::argv_or_defaultvalue(int n, T value) { - std::stringstream ss; - ss << value; - if (argc > n) return vargv[n]; - else return (ss.str()); -} - - -bool SetInitial::InputFromArgs(int _argc, char *_argv[]) { - // init params - this->argc = _argc; - if (argc <= 7) { - printf("too few arguments\n"); - return false; - } - vargv.resize(argc); - for (unsigned int i = 0; i < vargv.size(); i++) { - vargv[i] = _argv[i]; - } - - calcInfo.statistics = vargv[1]; - calcInfo.beta = atof(vargv[2].c_str()); - fileInfo.filein_G = vargv[3]; - fileInfo.col = atoi(vargv[4].c_str()); - fileInfo.fileout_spec = vargv[5]; - - omegaInfo.NW = atoi(vargv[6].c_str()); - omegaInfo.omega_min = atof(vargv[7].c_str()); - omegaInfo.omega_max = atof(vargv[8].c_str()); - - // The inputs above are mandatory - // The followings are optional - std::stringstream ss; - param.svd.sv_min = atof(argv_or_defaultvalue(9, param.svd.sv_min).c_str()); - - // string algorithm( argv_or_default(9, "oneshot") ); - param.lambda.Nl = atoi(argv_or_defaultvalue(10, param.lambda.Nl).c_str()); - param.lambda.lbegin = atof(argv_or_defaultvalue(11, param.lambda.lbegin).c_str()); - param.lambda.lend = atof(argv_or_defaultvalue(12, param.lambda.lend).c_str()); - - param.admm.penalty = atof(argv_or_defaultvalue(13, param.admm.penalty).c_str()); - param.admm.tolerance = atof(argv_or_defaultvalue(14, param.admm.tolerance).c_str()); - param.admm.max_iter = atoi(argv_or_defaultvalue(15, param.admm.max_iter).c_str()); - param.admm.flag_penalty_auto = param.admm.penalty < 0 ? true : false; - fileInfo.print_level = atoi(argv_or_defaultvalue(16, fileInfo.print_level).c_str()); - fileInfo.output_interval=atoi(argv_or_defaultvalue(17, fileInfo.output_interval).c_str()); - flags.validation = argv_or_defaultvalue(17, "y") == "y"; - flags.nonnegative = argv_or_defaultvalue(18, "ON") == "ON"; - return true; -} void SetInitial::PrintInfo() { printf("\nParameters:\n"); @@ -81,7 +30,9 @@ void SetInitial::PrintInfo() { printf(" beta = %lf\n", calcInfo.beta); //printf(" filein_G = '%s' \n", fileInfo.filein_G.c_str(), fileInfo.col); printf(" filein_G = '%s' \n", fileInfo.filein_G.c_str()); + printf(" filein_Gsigma = '%s' \n", fileInfo.filein_Gsigma.c_str()); printf(" fileout_spec = '%s'\n", fileInfo.fileout_spec.c_str()); + printf(" fileout_pade = '%s'\n", param.pade.filename.c_str()); printf(" omega = [ %lf : %lf ] (%d points)\n", omegaInfo.omega_min, omegaInfo.omega_max, omegaInfo.NW); @@ -99,6 +50,12 @@ void SetInitial::PrintInfo() { printf(" OutputInterval = %d\n", fileInfo.output_interval); printf(" FlagNonNegative = %s\n", flags.nonnegative ? "ON": "OFF"); + if (flags.refrho){ + printf("\n SpM-Pade:\n"); + printf(" PadeEta = %lf\n", param.pade.eta); + printf(" NSamplePade = %d\n", param.pade.nsample); + } + if (flags.validation) { printf("\n CROSS VALIDATION\n"); @@ -153,9 +110,13 @@ void SetInitial::RegisterMap(std::string _keyword, std::string _value) { bool SetInitial::SetDefaultValue() { RegisterMap("Statistics", ""); RegisterMap("beta", ""); + RegisterMap("G_sigma", "Inf"); RegisterMap("column", ""); + RegisterMap("column_sigma", "-1"); RegisterMap("filein_G", ""); + RegisterMap("filein_G_sigma", ""); RegisterMap("fileout_spec", "spectrum.out"); + RegisterMap("fileout_pade", "pade.out"); RegisterMap("NOmega", ""); RegisterMap("OmegaMin", ""); RegisterMap("OmegaMax", ""); @@ -172,6 +133,8 @@ bool SetInitial::SetDefaultValue() { RegisterMap("CrossValidation", "n"); RegisterMap("OutputInterval", "1"); RegisterMap("FlagNonNegative", "ON"); + RegisterMap("PadeEta", "0.0"); + RegisterMap("NSamplePade", "30"); /* std::map::iterator it; @@ -251,10 +214,19 @@ void SetInitial::SetInputValue() { calcInfo.statistics = GetValue("Statistics"); transform(calcInfo.statistics.begin(), calcInfo.statistics.end(), calcInfo.statistics.begin(), tolower); calcInfo.beta = std::stod(GetValue("beta")); + calcInfo.sigma = std::stod(GetValue("G_sigma")); fileInfo.col = std::stoi(GetValue("column")); + fileInfo.colsigma = std::stoi(GetValue("column_sigma")); fileInfo.filein_G = GetValue("filein_G"); fileInfo.fileout_spec = GetValue("fileout_spec"); fileInfo.output_interval = std::stoi(GetValue("OutputInterval")); + + { + std::string buf = GetValue("filein_G_sigma"); + if(buf==""){fileInfo.filein_Gsigma = fileInfo.filein_G;} + else{fileInfo.filein_Gsigma = buf;} + } + omegaInfo.NW = std::stoi(GetValue("NOmega")); omegaInfo.omega_min = std::stod(GetValue("OmegaMin")); omegaInfo.omega_max = std::stod(GetValue("OmegaMax")); @@ -290,6 +262,10 @@ void SetInitial::SetInputValue() { } */ + param.pade.filename = GetValue("fileout_pade"); + param.pade.nsample = std::stoi(GetValue("NSamplePade")); + param.pade.eta = std::stod(GetValue("PadeEta")); + flags.refrho = param.pade.eta != 0.0; if (param.lambda.Nl == 0) { if (fabs(param.lambda.dlambda) < pow(10.0, -12)) { diff --git a/c++/src/spm_core.cpp b/c++/src/spm_core.cpp index 57b5aba..5881994 100644 --- a/c++/src/spm_core.cpp +++ b/c++/src/spm_core.cpp @@ -14,9 +14,14 @@ /* You should have received a copy of the GNU General Public License */ /* along with this program. If not, see . */ #include +#include // for debug #include #include "common.h" +#include #include +#include + +#include double SPM_Core::integrate(std::vector &y, double width) { long unsigned int n = y.size(); @@ -75,33 +80,140 @@ int SPM_Core::SolveEquation( double _Beta, std::vector > &_AIn, std::vector &_Gtau, + std::vector &_Gtau_error, std::vector &_lambda, std::vector &_omega) { + + this->StatisticsType = _StatisticsType; + this->beta = _Beta; std::vector y = _Gtau; // input - if (_StatisticsType == "fermion") { + if (this->StatisticsType == "fermion") { for (unsigned int i = 0; i < y.size(); i++) { y[i] *= -1.0; } } + std::vector _omega_coeff(_omega); double sum_G; - if (_StatisticsType == "fermion") { + if (this->StatisticsType == "fermion") { sum_G = -_Gtau[0] - _Gtau[_Gtau.size() - 1]; + std::fill(_omega_coeff.begin(), _omega_coeff.end(), 1.0); } else if (_StatisticsType == "boson") { - sum_G = integrate(_Gtau, _Beta); + sum_G = _Gtau[0] + _Gtau[_Gtau.size() - 1]; + const size_t nw = _omega_coeff.size(); + for(size_t iw = 0; iw < nw; iw++){ + const double bw = _Beta * _omega[iw]; + if(bw > 0.0){ + _omega_coeff[iw] = -std::expm1(-bw) / (1.0 + std::exp(-bw)); + }else{ + _omega_coeff[iw] = std::expm1(bw) / (1.0 + std::exp(bw)); + } + } } else { std::cerr << "Error: StatisticsType must be fermion or boson." << std::endl; return -1; } - return SolveEquationCore(_AIn, y, _omega, _lambda, sum_G); + std::vector rho_ref(_omega.size()); + std::vector rho_ref_weight(_omega.size()); + std::vector rhopade(_omega.size()); + + // Pade from G(tau) + { + std::vector gt2(_Gtau.begin(), _Gtau.end()-1); + Gf gf; + if(_StatisticsType=="fermion"){ + // The tail (4th argument) is not used in SpM + gf.set_Gtau(gt2, Gf::FERMION, _Beta, 1.0); + }else{ + // The tail (4th argument) is not used in SpM + for(int i=0; i randn(0.0,1.0); + const int Nsample = param.pade.nsample; + for(int isample=0; isample gt2(_Gtau.begin(), _Gtau.end()-1); + const int Nt = gt2.size(); + for(int i=0; i > &_AIn, std::vector &_y, std::vector &_omega, + std::vector &_omega_coeff, std::vector &_lambda, - const double _sum_G) { + const double _sum_G, + std::vector &_rho_ref, + std::vector &_rho_ref_weight + ) { SetKernel(_AIn); // Solve the equations lambda = _lambda; @@ -134,7 +246,10 @@ int SPM_Core::SolveEquationCore( //core D.set_sumrule(_sum_G); - D.set_coef(lambda[l], param.admm.penalty, param.admm.penalty, flag_penalty_auto); + + D.set_rho_ref(_rho_ref, _rho_ref_weight); + D.set_omega_coeff(_omega_coeff); + D.set_coef(lambda[l], param.admm.penalty, param.admm.penalty, param.admm.penalty, flag_penalty_auto); D.set_y(y); D.solve(param.admm.tolerance, param.admm.max_iter); @@ -167,8 +282,27 @@ int SPM_Core::SolveEquationCore( // determine lambda param.lambda.lvalid = 0; if (flags.validation) { - param.lambda.lvalid = distance(valid.begin(), - min_element(valid.begin(), valid.end())); // get position of the minimum element + std::vector looc(param.lambda.Nl); + int M = D.get_svd_num(); + char filename[64] = "find_lambda_opt_CV.dat"; + FILE *fp = fopen((outputOrgDir + filename).c_str(), "w"); + fprintf(fp, "# log10(lambda) log10(E_LOCC) log10(lambda-E_LOCC) log10(E_LOCC-lambda) M S(lambda)\n"); + for (int l = 0; l < param.lambda.Nl; l++){ + int Slambda=info[l].l0_norm; + double x = M/static_cast(M-Slambda); + looc[l] = info[l].mse*x*x; + fprintf(fp, " %.15e %.15e %.15e %.15e %d %d\n", + log10(lambda[l]), log10(looc[l]), + log10(lambda[l]-looc[l]), log10(looc[l]-lambda[l]), + M, Slambda); + } + fclose(fp); + printf("'%s'\n", filename); + + // get position of the minimum element + param.lambda.lvalid = std::distance(looc.begin(), std::min_element(looc.begin(), looc.end())); + valid.resize(lambda.size()); + fill(valid.begin(), valid.end(), 0.); } else { std::vector mse, diff, log_f; for (int l = 0; l < param.lambda.Nl; l++) mse.push_back(info[l].mse); @@ -233,8 +367,7 @@ int SPM_Core::find_kink(std::vector &x, std::vector &y, std::vec } void SPM_Core::GetSpectrum(std::vector &_spectrum) { - //_spectrum = result[param.lambda.lvalid].x; - if(flags.nonnegative==true) { + if(flags.nonnegative==true || flags.refrho ) { _spectrum = result[param.lambda.lvalid].z2; } else{ @@ -248,14 +381,17 @@ void SPM_Core::GetLambdaOpt(std::vector &_lambda, int *_opt_l) { } void SPM_Core::GetResults(std::vector &_vmse, std::vector &_vmse_full, - std::vector &_vl1_norm, std::vector &_valid) { + std::vector &_vl0_norm, std::vector &_vl1_norm, + std::vector &_valid) { _vmse.resize(info.size()); _vmse_full.resize(info.size()); + _vl0_norm.resize(info.size()); _vl1_norm.resize(info.size()); _valid.resize(info.size()); for (unsigned i = 0; i < info.size(); i++) { _vmse[i] = info[i].mse; _vmse_full[i] = info[i].mse_full; + _vl0_norm[i] = info[i].l0_norm; _vl1_norm[i] = info[i].l1_norm; _valid[i] = valid[i]; } @@ -269,7 +405,7 @@ void SPM_Core::print_results_admm(admm_result &r, std::vector w, double filename = prefix + "x_sv.dat"; fp = fopen(filename.c_str(), "w"); for (unsigned i = 0; i < r.xsv.size(); i++) { - fprintf(fp, "%d %.5e %.5e %.5e\n", i, r.xsv[i], r.z1sv[i], r.z2sv[i]); + fprintf(fp, "%d %.5e %.5e %.5e %.5e\n", i, r.xsv[i], r.z1sv[i], r.z2sv[i], r.z3sv[i]); // fprintf(fp, "%d %.5e %.5e %.5e\n", i, r.xsv[i]/dw, r.z1sv[i]/dw, r.z2sv[i]/dw); // fprintf(fp, "%d %.5e %.5e %.5e\n", i, r.xsv[i]/sqrt(dw), r.z1sv[i]/sqrt(dw), r.z2sv[i]/sqrt(dw)); } @@ -280,7 +416,7 @@ void SPM_Core::print_results_admm(admm_result &r, std::vector w, double filename = prefix + "x_tw.dat"; fp = fopen(filename.c_str(), "w"); for (unsigned i = 0; i < r.x.size(); i++) { - fprintf(fp, "%.5e %.5e %.5e %.5e\n", w[i], r.x[i] / dw, r.z1[i] / dw, r.z2[i] / dw); + fprintf(fp, "%.5e %.5e %.5e %.5e %.5e\n", w[i], r.x[i] / dw, r.z1[i] / dw, r.z2[i] / dw, r.z3[i] / dw); } fclose(fp); printf("'%s'\n", filename.c_str()); @@ -293,12 +429,11 @@ void SPM_Core::print_results_admm(admm_result &r, std::vector w, double } fclose(fp); printf("'%s'\n", filename.c_str()); - // y in tau-omega basis filename = prefix + "y_tw.dat"; fp = fopen(filename.c_str(), "w"); for (unsigned i = 0; i < r.y.size(); i++) { - fprintf(fp, "%.5e %.5e %.5e %.5e\n", static_cast(i) / r.y.size(), r.y[i], r.y_recovered_x[i], + fprintf(fp, "%.5e %.5e %.5e %.5e\n", static_cast(i) / (r.y.size()-1), r.y[i], r.y_recovered_x[i], r.y_recovered_z1[i]); } fclose(fp); diff --git a/cmake/EnableGtests.cmake b/cmake/EnableGtests.cmake new file mode 100644 index 0000000..ed164e8 --- /dev/null +++ b/cmake/EnableGtests.cmake @@ -0,0 +1,28 @@ +# +# enables testing with google test +# provides add_gtest(test) and add_gtest_release(test) commands +# + +# check xml output +option(TestXMLOutput "Output tests to xml" OFF) + +# custom function to add gtest with xml output +# arg0 - test (assume the source is ${test}.cpp +function(add_gtest test) + if (TestXMLOutput) + set (test_xml_output --gtest_output=xml:${test}.xml) + endif(TestXMLOutput) + + if(${ARGC} EQUAL 2) + set(source "${ARGV1}/${test}") + set(gtest_src "${ARGV1}/gtest_main.cc;${ARGV1}/gtest-all.cc") + else(${ARGC} EQUAL 2) + set(source "${test}") + set(gtest_src "gtest_main.cc;gtest-all.cc") + endif(${ARGC} EQUAL 2) + + add_executable(${test} ${source} ${gtest_src}) + target_link_libraries(${test} ${LINK_ALL}) + add_test(NAME ${test} COMMAND ${test} ${test_xml_output}) +endfunction(add_gtest) + diff --git a/docs/manual/build/html/_images/math/08a277da44635b36e49e97a084c18ddbebeb65c6.png b/docs/manual/build/html/_images/math/08a277da44635b36e49e97a084c18ddbebeb65c6.png deleted file mode 100644 index fdf10a5..0000000 Binary files a/docs/manual/build/html/_images/math/08a277da44635b36e49e97a084c18ddbebeb65c6.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/08ee2ced3e4c74283f24491dbb1e4612e1953b38.png b/docs/manual/build/html/_images/math/08ee2ced3e4c74283f24491dbb1e4612e1953b38.png deleted file mode 100644 index f915d1b..0000000 Binary files a/docs/manual/build/html/_images/math/08ee2ced3e4c74283f24491dbb1e4612e1953b38.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/0949ed70c72f091f7e8ac4a7dd9c7bd9d859e06e.png b/docs/manual/build/html/_images/math/0949ed70c72f091f7e8ac4a7dd9c7bd9d859e06e.png deleted file mode 100644 index bbeeb80..0000000 Binary files a/docs/manual/build/html/_images/math/0949ed70c72f091f7e8ac4a7dd9c7bd9d859e06e.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/26e8299b484849735af4360c90f176656e4c035a.png b/docs/manual/build/html/_images/math/26e8299b484849735af4360c90f176656e4c035a.png deleted file mode 100644 index be35853..0000000 Binary files a/docs/manual/build/html/_images/math/26e8299b484849735af4360c90f176656e4c035a.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/34643d2d999fcc53423d45f1bb2392eb2ca4772b.png b/docs/manual/build/html/_images/math/34643d2d999fcc53423d45f1bb2392eb2ca4772b.png deleted file mode 100644 index 5a1eac4..0000000 Binary files a/docs/manual/build/html/_images/math/34643d2d999fcc53423d45f1bb2392eb2ca4772b.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/359c8be672c286f0451ce49e8c8b7db7f99fad74.png b/docs/manual/build/html/_images/math/359c8be672c286f0451ce49e8c8b7db7f99fad74.png deleted file mode 100644 index 8559b23..0000000 Binary files a/docs/manual/build/html/_images/math/359c8be672c286f0451ce49e8c8b7db7f99fad74.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/379721c563c5940c2c493c638e2bbd7607140312.png b/docs/manual/build/html/_images/math/379721c563c5940c2c493c638e2bbd7607140312.png deleted file mode 100644 index 6be6786..0000000 Binary files a/docs/manual/build/html/_images/math/379721c563c5940c2c493c638e2bbd7607140312.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/3a3d1af28de6ef34516ac9d699c48cbbb23ce998.png b/docs/manual/build/html/_images/math/3a3d1af28de6ef34516ac9d699c48cbbb23ce998.png deleted file mode 100644 index f07dea4..0000000 Binary files a/docs/manual/build/html/_images/math/3a3d1af28de6ef34516ac9d699c48cbbb23ce998.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/3cc42bfbdd8b170705ca942fa5fc4974d1580ec3.png b/docs/manual/build/html/_images/math/3cc42bfbdd8b170705ca942fa5fc4974d1580ec3.png deleted file mode 100644 index 924e69f..0000000 Binary files a/docs/manual/build/html/_images/math/3cc42bfbdd8b170705ca942fa5fc4974d1580ec3.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/3ef352c72eac36439ea87320904d745f331d0af6.png b/docs/manual/build/html/_images/math/3ef352c72eac36439ea87320904d745f331d0af6.png deleted file mode 100644 index 7e7ceaf..0000000 Binary files a/docs/manual/build/html/_images/math/3ef352c72eac36439ea87320904d745f331d0af6.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/3f9b358d1499cd88c763c093b52622a9820128c0.png b/docs/manual/build/html/_images/math/3f9b358d1499cd88c763c093b52622a9820128c0.png deleted file mode 100644 index b7187a8..0000000 Binary files a/docs/manual/build/html/_images/math/3f9b358d1499cd88c763c093b52622a9820128c0.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/410a9d0df9c135dd73b269cba7ef04dcfd932b1f.png b/docs/manual/build/html/_images/math/410a9d0df9c135dd73b269cba7ef04dcfd932b1f.png deleted file mode 100644 index 7416cf4..0000000 Binary files a/docs/manual/build/html/_images/math/410a9d0df9c135dd73b269cba7ef04dcfd932b1f.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/4569268c314a44232c428f1e7bf62b8bd225381a.png b/docs/manual/build/html/_images/math/4569268c314a44232c428f1e7bf62b8bd225381a.png deleted file mode 100644 index 8cd6a7f..0000000 Binary files a/docs/manual/build/html/_images/math/4569268c314a44232c428f1e7bf62b8bd225381a.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/47aa9388b4a90760ec185881759daa4766d51191.png b/docs/manual/build/html/_images/math/47aa9388b4a90760ec185881759daa4766d51191.png deleted file mode 100644 index 9d95db4..0000000 Binary files a/docs/manual/build/html/_images/math/47aa9388b4a90760ec185881759daa4766d51191.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/47cdea4ba963069691bfe84c50e17683025c0979.png b/docs/manual/build/html/_images/math/47cdea4ba963069691bfe84c50e17683025c0979.png deleted file mode 100644 index 6449c86..0000000 Binary files a/docs/manual/build/html/_images/math/47cdea4ba963069691bfe84c50e17683025c0979.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/50c5e67e5e093ef5d995560386cea9511bde7f33.png b/docs/manual/build/html/_images/math/50c5e67e5e093ef5d995560386cea9511bde7f33.png deleted file mode 100644 index 31c0a23..0000000 Binary files a/docs/manual/build/html/_images/math/50c5e67e5e093ef5d995560386cea9511bde7f33.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/518712cb83e744e0d12c45a910b011dcca076a39.png b/docs/manual/build/html/_images/math/518712cb83e744e0d12c45a910b011dcca076a39.png deleted file mode 100644 index 1f9c42c..0000000 Binary files a/docs/manual/build/html/_images/math/518712cb83e744e0d12c45a910b011dcca076a39.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/604b498d61fdb6213317eab6a9f7240af06a9911.png b/docs/manual/build/html/_images/math/604b498d61fdb6213317eab6a9f7240af06a9911.png deleted file mode 100644 index dc6065a..0000000 Binary files a/docs/manual/build/html/_images/math/604b498d61fdb6213317eab6a9f7240af06a9911.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/61e319ecabcb5c23765c66fb100f64ec3aaf1f14.png b/docs/manual/build/html/_images/math/61e319ecabcb5c23765c66fb100f64ec3aaf1f14.png deleted file mode 100644 index b8708fe..0000000 Binary files a/docs/manual/build/html/_images/math/61e319ecabcb5c23765c66fb100f64ec3aaf1f14.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/636d4c832194e64ac08693518831e906009d3161.png b/docs/manual/build/html/_images/math/636d4c832194e64ac08693518831e906009d3161.png deleted file mode 100644 index 333a2da..0000000 Binary files a/docs/manual/build/html/_images/math/636d4c832194e64ac08693518831e906009d3161.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/64666578228826a7dc2a0ecfa07b9e1618527ed5.png b/docs/manual/build/html/_images/math/64666578228826a7dc2a0ecfa07b9e1618527ed5.png deleted file mode 100644 index 583b27f..0000000 Binary files a/docs/manual/build/html/_images/math/64666578228826a7dc2a0ecfa07b9e1618527ed5.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/666d5022792b496e5bc4af1c5fb7361d9e5ed685.png b/docs/manual/build/html/_images/math/666d5022792b496e5bc4af1c5fb7361d9e5ed685.png deleted file mode 100644 index ec8450b..0000000 Binary files a/docs/manual/build/html/_images/math/666d5022792b496e5bc4af1c5fb7361d9e5ed685.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/684381a21cd73ebbf43b63a087d3f7410ee99ce8.png b/docs/manual/build/html/_images/math/684381a21cd73ebbf43b63a087d3f7410ee99ce8.png deleted file mode 100644 index aa57866..0000000 Binary files a/docs/manual/build/html/_images/math/684381a21cd73ebbf43b63a087d3f7410ee99ce8.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/6ed94dac5ae2ab85d25d0c11ed15a62e4b9568d0.png b/docs/manual/build/html/_images/math/6ed94dac5ae2ab85d25d0c11ed15a62e4b9568d0.png deleted file mode 100644 index 0f06994..0000000 Binary files a/docs/manual/build/html/_images/math/6ed94dac5ae2ab85d25d0c11ed15a62e4b9568d0.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/6f30456384f9635601349f6bb0de3ccdd58b0837.png b/docs/manual/build/html/_images/math/6f30456384f9635601349f6bb0de3ccdd58b0837.png deleted file mode 100644 index d9cc617..0000000 Binary files a/docs/manual/build/html/_images/math/6f30456384f9635601349f6bb0de3ccdd58b0837.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/76f1d8ace30435987c01a00ca53a71cba1f40e6c.png b/docs/manual/build/html/_images/math/76f1d8ace30435987c01a00ca53a71cba1f40e6c.png deleted file mode 100644 index e6fcdcc..0000000 Binary files a/docs/manual/build/html/_images/math/76f1d8ace30435987c01a00ca53a71cba1f40e6c.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/800a289af47e2f3052ffb324c8f6987bc496052f.png b/docs/manual/build/html/_images/math/800a289af47e2f3052ffb324c8f6987bc496052f.png deleted file mode 100644 index 9dee3ae..0000000 Binary files a/docs/manual/build/html/_images/math/800a289af47e2f3052ffb324c8f6987bc496052f.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/85c913e5997c47b58a8a2ce9ecac21612aeee15b.png b/docs/manual/build/html/_images/math/85c913e5997c47b58a8a2ce9ecac21612aeee15b.png deleted file mode 100644 index c32dc6b..0000000 Binary files a/docs/manual/build/html/_images/math/85c913e5997c47b58a8a2ce9ecac21612aeee15b.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/8cb351dd3b28ea1c9fb991ca73ab98acd8cfce4a.png b/docs/manual/build/html/_images/math/8cb351dd3b28ea1c9fb991ca73ab98acd8cfce4a.png deleted file mode 100644 index cf8c650..0000000 Binary files a/docs/manual/build/html/_images/math/8cb351dd3b28ea1c9fb991ca73ab98acd8cfce4a.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/8ddb02d4e899000ba2ee76d62ac4ab491fd44a75.png b/docs/manual/build/html/_images/math/8ddb02d4e899000ba2ee76d62ac4ab491fd44a75.png deleted file mode 100644 index 7c68a27..0000000 Binary files a/docs/manual/build/html/_images/math/8ddb02d4e899000ba2ee76d62ac4ab491fd44a75.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/9d49d6f3960ae938ac2efcc52449f7f93bc3b97b.png b/docs/manual/build/html/_images/math/9d49d6f3960ae938ac2efcc52449f7f93bc3b97b.png deleted file mode 100644 index 7d38906..0000000 Binary files a/docs/manual/build/html/_images/math/9d49d6f3960ae938ac2efcc52449f7f93bc3b97b.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/a372fe4adb7d7aa51df8ac102e5b7123b3b8f24c.png b/docs/manual/build/html/_images/math/a372fe4adb7d7aa51df8ac102e5b7123b3b8f24c.png deleted file mode 100644 index 3b7afc5..0000000 Binary files a/docs/manual/build/html/_images/math/a372fe4adb7d7aa51df8ac102e5b7123b3b8f24c.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/a3965ac90c9513644606210339d525386626dcb2.png b/docs/manual/build/html/_images/math/a3965ac90c9513644606210339d525386626dcb2.png deleted file mode 100644 index b048c4e..0000000 Binary files a/docs/manual/build/html/_images/math/a3965ac90c9513644606210339d525386626dcb2.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/a3f47b047edcf70c20b71a148ef35825bb7e2fe1.png b/docs/manual/build/html/_images/math/a3f47b047edcf70c20b71a148ef35825bb7e2fe1.png deleted file mode 100644 index 693dd6b..0000000 Binary files a/docs/manual/build/html/_images/math/a3f47b047edcf70c20b71a148ef35825bb7e2fe1.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/a6d47234c6a3b8d24783aaebc72d9d4359d13c82.png b/docs/manual/build/html/_images/math/a6d47234c6a3b8d24783aaebc72d9d4359d13c82.png deleted file mode 100644 index 1a7959e..0000000 Binary files a/docs/manual/build/html/_images/math/a6d47234c6a3b8d24783aaebc72d9d4359d13c82.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/b147153e026a1d013217d1f1cf8a5e039971806e.png b/docs/manual/build/html/_images/math/b147153e026a1d013217d1f1cf8a5e039971806e.png deleted file mode 100644 index 84c1407..0000000 Binary files a/docs/manual/build/html/_images/math/b147153e026a1d013217d1f1cf8a5e039971806e.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/b172478c2e144fb77ed130c729b5ad404895872f.png b/docs/manual/build/html/_images/math/b172478c2e144fb77ed130c729b5ad404895872f.png deleted file mode 100644 index 65638d9..0000000 Binary files a/docs/manual/build/html/_images/math/b172478c2e144fb77ed130c729b5ad404895872f.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/b33b56c5293081e5b6575bf2f6f1eecc81b16d8c.png b/docs/manual/build/html/_images/math/b33b56c5293081e5b6575bf2f6f1eecc81b16d8c.png deleted file mode 100644 index f721824..0000000 Binary files a/docs/manual/build/html/_images/math/b33b56c5293081e5b6575bf2f6f1eecc81b16d8c.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/b359488b993294ebbc2c7b30ab8f749dcbc6826d.png b/docs/manual/build/html/_images/math/b359488b993294ebbc2c7b30ab8f749dcbc6826d.png deleted file mode 100644 index 70ddbe8..0000000 Binary files a/docs/manual/build/html/_images/math/b359488b993294ebbc2c7b30ab8f749dcbc6826d.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/b699a4b61a18d3471ea57883736739f4d4ce3770.png b/docs/manual/build/html/_images/math/b699a4b61a18d3471ea57883736739f4d4ce3770.png deleted file mode 100644 index 2cc54d6..0000000 Binary files a/docs/manual/build/html/_images/math/b699a4b61a18d3471ea57883736739f4d4ce3770.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/b7d76a97ccf41eda71495f66486eff2f457bf3f3.png b/docs/manual/build/html/_images/math/b7d76a97ccf41eda71495f66486eff2f457bf3f3.png deleted file mode 100644 index bd71af0..0000000 Binary files a/docs/manual/build/html/_images/math/b7d76a97ccf41eda71495f66486eff2f457bf3f3.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/b80430baf1b059d0097bc9081abca3b1e0a74287.png b/docs/manual/build/html/_images/math/b80430baf1b059d0097bc9081abca3b1e0a74287.png deleted file mode 100644 index fa9397b..0000000 Binary files a/docs/manual/build/html/_images/math/b80430baf1b059d0097bc9081abca3b1e0a74287.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/bd30e06c49c5b6e0f1ff7b39d06a6f8292259ddb.png b/docs/manual/build/html/_images/math/bd30e06c49c5b6e0f1ff7b39d06a6f8292259ddb.png deleted file mode 100644 index 374327d..0000000 Binary files a/docs/manual/build/html/_images/math/bd30e06c49c5b6e0f1ff7b39d06a6f8292259ddb.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/c2169fa8df7200e1571d40015f85cc20b800f2f0.png b/docs/manual/build/html/_images/math/c2169fa8df7200e1571d40015f85cc20b800f2f0.png deleted file mode 100644 index a126775..0000000 Binary files a/docs/manual/build/html/_images/math/c2169fa8df7200e1571d40015f85cc20b800f2f0.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/c236cbdfa9836dffda883ffed7a9336e0cb5f87f.png b/docs/manual/build/html/_images/math/c236cbdfa9836dffda883ffed7a9336e0cb5f87f.png deleted file mode 100644 index 2db6200..0000000 Binary files a/docs/manual/build/html/_images/math/c236cbdfa9836dffda883ffed7a9336e0cb5f87f.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/c33c44a9858afa000c1c2dedf203753df4baadc8.png b/docs/manual/build/html/_images/math/c33c44a9858afa000c1c2dedf203753df4baadc8.png deleted file mode 100644 index 5507881..0000000 Binary files a/docs/manual/build/html/_images/math/c33c44a9858afa000c1c2dedf203753df4baadc8.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/d1814eaebc091b796e7d216225889ea1ca054937.png b/docs/manual/build/html/_images/math/d1814eaebc091b796e7d216225889ea1ca054937.png deleted file mode 100644 index f50bb5f..0000000 Binary files a/docs/manual/build/html/_images/math/d1814eaebc091b796e7d216225889ea1ca054937.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/d3abeeec19e04fc372e23a3383efcb60678ecdb8.png b/docs/manual/build/html/_images/math/d3abeeec19e04fc372e23a3383efcb60678ecdb8.png deleted file mode 100644 index 4888d00..0000000 Binary files a/docs/manual/build/html/_images/math/d3abeeec19e04fc372e23a3383efcb60678ecdb8.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/db9038319d8e3a6b6b43071d4d08701d9d70c2ac.png b/docs/manual/build/html/_images/math/db9038319d8e3a6b6b43071d4d08701d9d70c2ac.png deleted file mode 100644 index ba05876..0000000 Binary files a/docs/manual/build/html/_images/math/db9038319d8e3a6b6b43071d4d08701d9d70c2ac.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/dd5601c05723358c40a8cb0fa9b1541eb1def298.png b/docs/manual/build/html/_images/math/dd5601c05723358c40a8cb0fa9b1541eb1def298.png deleted file mode 100644 index a125a95..0000000 Binary files a/docs/manual/build/html/_images/math/dd5601c05723358c40a8cb0fa9b1541eb1def298.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/e2bf38137174bbbc865558597cdacfe137d43d52.png b/docs/manual/build/html/_images/math/e2bf38137174bbbc865558597cdacfe137d43d52.png deleted file mode 100644 index bd4fb22..0000000 Binary files a/docs/manual/build/html/_images/math/e2bf38137174bbbc865558597cdacfe137d43d52.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/ed882223f7d35e6c8ab4844af95f33706aaa0c31.png b/docs/manual/build/html/_images/math/ed882223f7d35e6c8ab4844af95f33706aaa0c31.png deleted file mode 100644 index 43b18d3..0000000 Binary files a/docs/manual/build/html/_images/math/ed882223f7d35e6c8ab4844af95f33706aaa0c31.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/f36101c441bcf511a4c70e9f0fd21507af66d746.png b/docs/manual/build/html/_images/math/f36101c441bcf511a4c70e9f0fd21507af66d746.png deleted file mode 100644 index 0b4c833..0000000 Binary files a/docs/manual/build/html/_images/math/f36101c441bcf511a4c70e9f0fd21507af66d746.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/f94e13b8b141ea338343a4edcb7d14d0283209be.png b/docs/manual/build/html/_images/math/f94e13b8b141ea338343a4edcb7d14d0283209be.png deleted file mode 100644 index f93f0cc..0000000 Binary files a/docs/manual/build/html/_images/math/f94e13b8b141ea338343a4edcb7d14d0283209be.png and /dev/null differ diff --git a/docs/manual/build/html/_images/math/fac78ecd202fb268bd8ad30d6fa2176d9f74c01c.png b/docs/manual/build/html/_images/math/fac78ecd202fb268bd8ad30d6fa2176d9f74c01c.png deleted file mode 100644 index 67b0fa7..0000000 Binary files a/docs/manual/build/html/_images/math/fac78ecd202fb268bd8ad30d6fa2176d9f74c01c.png and /dev/null differ diff --git a/docs/manual/build/html/_sources/docs/inputfile.rst.txt b/docs/manual/build/html/_sources/docs/inputfile.rst.txt index ff56082..d9d90fb 100644 --- a/docs/manual/build/html/_sources/docs/inputfile.rst.txt +++ b/docs/manual/build/html/_sources/docs/inputfile.rst.txt @@ -8,8 +8,24 @@ Input files =============================== 1. Parameter file (param.in) - -The variables for which default value is not given are mandatory. + + - The format of the parameter file is the "name = value" pair format. + + - ``#`` means a comment. In other words, a ``#`` character and the subsequent part in a line will be ignored. + + - A blank line is skipped. + + - The variables for which default value is not given are mandatory. + +.. code-block:: + + # character means comment + + statistics = "fermion" # stats. + beta = 100.0 # inv. temp. + + ### continued ### + * INPUT/OUTPUT @@ -17,13 +33,13 @@ The variables for which default value is not given are mandatory. :header-rows: 1 :widths: 1,1,2,4 - Name, Type, Default value, Description - statistics, String, ---, Choose "fermion" or "boson". - beta, Double, ---, The inverse temperature - filein_G, String, Gtau.in, The name of an input file for Green's function. - column, Integer, 1, The column number where the values of G(tau) are stored in the filein_G file. - fileout_spec, String, spectrum.out, The name of output file of spectrum. - + Name, Type, Default value, Description + statistics, String, ---, Choose "fermion" or "boson". + beta, Double, ---, The inverse temperature + filein_G, String, Gtau.in, The name of an input file for Green's function. + column, Integer, 1, The column number where the values of G(tau) are stored in the ``filein_G file``. + fileout_spec, String, spectrum.out, The name of output file of spectrum. + fileout_pade, String, pade.out, The name of output file of spectrum calculated by the Pade approximant. * OMEGA @@ -66,6 +82,19 @@ The variables for which default value is not given are mandatory. maxiteration, Integer,1000, The maximum number of iterations. printlevel, Integer,2, "0; minimum, 1; moderate, 2; verbose." +* SpM-Pade + + .. csv-table:: + :header-rows: 1 + :widths: 1,1,2,4 + + Name, Type, Default value, Description + PadeEta, Double, 0.0, The weight of the spectrum by the Pade approximant. + filein_Gsigma, String, Gtau.in, The name of input file for errors of Green's function. + column_sigma, Integer, 1, The column number where the errors of G(tau) are stored in the ``filein_Gsigma`` file. + g_sigma, Double, inf, "The errors of G(tau). If this is inf (default value), the errors will be loaded from a file specified by ``filein_Gsigma``." + NSamplePade, Integer, 30, The number of random samples (with noise) will be used to calculate the expected value and the standard deviation of the spectrum by the Pade approximant. + 2. Green's function (Gtau.in) In SPM, the values of Green's function is only used for calculation, diff --git a/docs/manual/build/html/_sources/docs/outputfile.rst.txt b/docs/manual/build/html/_sources/docs/outputfile.rst.txt index d206824..a371715 100644 --- a/docs/manual/build/html/_sources/docs/outputfile.rst.txt +++ b/docs/manual/build/html/_sources/docs/outputfile.rst.txt @@ -17,6 +17,7 @@ List of files output ├── spectrum.dat + ├── pade.dat ├── lambda_dep.dat ├── find_lambda_opt.dat ├── SV.dat @@ -54,7 +55,7 @@ Some comments may be given at the header of files (lines beginning with '#'). * 1st column: :math:`\omega_i` * 2nd column: :math:`\rho(\omega_i)` - ex.) calculated results for fermion sample outputted as ``samles/fermion/spectrum.dat``. + ex.) calculated results for fermion sample outputted as ``samples/fermion/output/spectrum.dat``. :: @@ -64,6 +65,15 @@ Some comments may be given at the header of files (lines beginning with '#'). -3.98400e+00 0.00000e+00 -3.97600e+00 0.00000e+00 +* *spectrum.dat* + + Real frequency spectrum :math:`\rho(\omega)` calculated by the Pade approximant. + + * 1st column: :math:`\omega_i` + * 2nd column: :math:`\rho(\omega_i)` from G(tau) in ``filein_G``. + * 3rd column: The expectation value of :math:`\rho(\omega_i)` over the generated ``NSamplePade`` samples with independent noise. + * 4th column: The sample standard deviation of :math:`\rho(\omega_i)` over the generated ``NsamplePade`` samples with independent noise. + * *lambda_dep.dat* Several quantities as a function of lambda. @@ -71,30 +81,31 @@ Some comments may be given at the header of files (lines beginning with '#'). * 1st column: :math:`\lambda` * 2nd column: Square error in the SV basis, :math:`\| \bm{y}'-S\bm{x}' \|_2^2` (the first term in :math:`F`) * 3rd column: Square error computed in the original basis, :math:`\| \bm{y}-K\bm{x} \|_2^2` - * 4th column: L1 norm of :math:`\bm{x}'`, :math:`\| \bm{x}' \|_1` (the second term in :math:`F`) - * 5th column: not used + * 4th column: L0 norm of :math:`\bm{x}'`, :math:`\| \bm{x}' \|_0` (the number of non-zero components in :math:`\bm{x}'`) + * 5th column: L1 norm of :math:`\bm{x}'`, :math:`\| \bm{x}' \|_1` (the second term in :math:`F`) + * 6th column: not used - ex.) calculated results for fermion sample outputted as ``samles/fermion/lambda_dep.dat``. + ex.) calculated results for fermion sample outputted as ``samples/fermion/output/lambda_dep.dat``. :: - 1.00000e+02 1.09110e-01 1.13080e-01 4.81361e-02 0.00000e+00 - 6.30957e+01 1.08854e-01 1.12824e-01 4.81376e-02 0.00000e+00 - 3.98107e+01 1.08752e-01 1.12722e-01 4.81386e-02 0.00000e+00 - 2.51189e+01 1.08710e-01 1.12679e-01 4.81394e-02 0.00000e+00 - 1.58489e+01 1.08609e-01 1.12578e-01 4.81459e-02 0.00000e+00 + 1.00000e+02 1.09110e-01 1.13080e-01 2 4.81361e-02 0.00000e+00 + 6.30957e+01 1.08854e-01 1.12824e-01 2 4.81376e-02 0.00000e+00 + 3.98107e+01 1.08752e-01 1.12722e-01 2 4.81386e-02 0.00000e+00 + 2.51189e+01 1.08710e-01 1.12679e-01 2 4.81394e-02 0.00000e+00 + 1.58489e+01 1.08609e-01 1.12578e-01 2 4.81459e-02 0.00000e+00 * *find_lambda_opt.dat* - Auxiliary data that are used to determine the optimal value of lambda. See ref.?? for details. + Auxiliary data that are used to determine the optimal value of lambda. See our original paper `Phys. Rev. E 95, 061302(R) (2017) `_ for details. * 1st column: :math:`\log_{10}(\lambda)` * 2nd column: :math:`\log_{10}(f(x)/\chi^2)` * 3rd column: :math:`\log_{10}(\chi^2)` * 4th column: :math:`\log_{10}(f(x))` - ex.) calculated results for fermion sample outputted as ``samles/fermion/find_lambda_opt.dat``. + ex.) calculated results for fermion sample outputted as ``samples/fermion/output/find_lambda_opt.dat``. :: @@ -111,7 +122,7 @@ Some comments may be given at the header of files (lines beginning with '#'). * 1st column: index :math:`l` (starting from 0) * 2nd column: The singular value :math:`S_l` in descending order - ex.) calculated results for fermion sample outputted as ``samles/fermion/find_lambda_opt.dat``. + ex.) calculated results for fermion sample outputted as ``samples/fermion/output/SV.dat``. :: @@ -138,7 +149,7 @@ Some comments may be given at the header of files (lines beginning with '#'). * 9th column: sum of spectrum, :math:`\langle V\bm{x}' \rangle` * 10th column: negative weight in the spectrum :math:`V\bm{x}'` - ex.) calculated results for fermion sample outputted as ``samles/fermion/lambda_opt/iter.dat``. + ex.) calculated results for fermion sample outputted as ``samples/fermion/output/lambda_opt/iter.dat``. :: @@ -157,7 +168,7 @@ Some comments may be given at the header of files (lines beginning with '#'). * 3rd column: :math:`z'_l` (**must be sparse**) * 4th column: :math:`V^{\rm t}\bm{z}` - ex.) calculated results for fermion sample outputted as ``samles/fermion/lambda_opt/x_sv.dat``. + ex.) calculated results for fermion sample outputted as ``samples/fermion/output/lambda_opt/x_sv.dat``. :: @@ -177,7 +188,7 @@ Some comments may be given at the header of files (lines beginning with '#'). * 3rd column: :math:`V\bm{z}'` (final result for a given lambda when non-negativity is not imposed) * 4th column: :math:`\bm{z}` (**final result** for a given lambda when non-negativity is imposed) - ex.) calculated results for fermion sample outputted as ``samles/fermion/lambda_opt/x_tw.dat``. + ex.) calculated results for fermion sample outputted as ``samples/fermion/output/lambda_opt/x_tw.dat``. :: @@ -197,7 +208,7 @@ Some comments may be given at the header of files (lines beginning with '#'). * 3rd column: :math:`S\bm{x}'` * 4th column: :math:`S\bm{z}'` - ex.) calculated results for fermion sample outputted as ``samles/fermion/lambda_opt/y_sv.dat``. + ex.) calculated results for fermion sample outputted as ``samples/fermion/output/lambda_opt/y_sv.dat``. :: @@ -216,7 +227,7 @@ Some comments may be given at the header of files (lines beginning with '#'). * 3rd column: :math:`U S \bm{x}'` * 4th column: :math:`U S \bm{z}'` - ex.) calculated results for fermion sample outputted as ``samles/fermion/lambda_opt/y_tw.dat``. + ex.) calculated results for fermion sample outputted as ``samples/fermion/output/lambda_opt/y_tw.dat``. :: diff --git a/docs/manual/build/html/_sources/docs/tutorials.rst.txt b/docs/manual/build/html/_sources/docs/tutorials.rst.txt index 0bcfb13..ad50101 100644 --- a/docs/manual/build/html/_sources/docs/tutorials.rst.txt +++ b/docs/manual/build/html/_sources/docs/tutorials.rst.txt @@ -10,10 +10,15 @@ Tutorials Sample data are provided both for fermionic and bosonic cases: - - ``samples/fermion`` # a sample for fermionic spectrum (data in the article) + - ``samples/fermion`` # a sample for fermionic spectrum (data in the SpM article) + - ``samples/fermion_twopeak`` # another sample for fermionic spectrum (data in the SpM-Pade article) - ``samples/boson`` # a sample for bosonic spectrum -Here, an explanation is given for the fermionic case. + +Fermionic system +~~~~~~~~~~~~~~~~~ + +Here, an explanation is given for the fermionic case (``samples/fermion``). Script ---------------------------------- @@ -103,7 +108,7 @@ In the directory ``samples/fermion/``, just type the command :: - $ build_directory/src/SpM.out -i param.in + $ build_directory/src/SpM.out param.in Results ------- @@ -167,3 +172,45 @@ Let us look at some graphs below. The light blue points show the input data :math:`G'_l` transformed into the SVD basis, and red circles are data used in computing :math:`\rho(\omega)` above. The blue points show, for comparison, the exact :math:`G'_l` without noise, which is provided in the file ``Gtau.in.sv_basis`` (not output of the ``SpM`` program). + + +Bosonic system +~~~~~~~~~~~~~~~~~~~ + +Sample files are available at ``samples/boson``. + +Users can deal with a bosonic system by only changing ``statistics`` parameter to ``"boson"`` from ``"fermion"`` in the parameter file. + + +SpM-Pade method +~~~~~~~~~~~~~~~~~ + +Sample files are available at ``samples/fermion_twopeak``. + +The SpM-Pade method is another SpM AC method, where unphysical oscillation in the SpM spectrum is reduced by using spectra reconstructed by the Pade AC method. + +The cost function is + +.. math:: + + L_\text{SpM-Pade}(\boldsymbol{\rho}) = L_\text{SpM}(\boldsymbol{\rho}) + \frac{\eta}{2} \sum_i w_i \left(\rho(\omega_i) - \rho^\text{Pade}(\omega_i)\right)^2, + +where + +.. math:: + + w_i = \left(1.0 + \left(\frac{\sigma^\text{Pade}(\omega_i)}{\rho^\text{Pade}(\omega_i)}\right)^2\right)^{-1} + +and :math:`\rho^\text{Pade}(\omega_i)` and :math:`\sigma^\text{Pade}(\omega_i)` are the expectation value and standard deviation of the spectrum reconstructed by the Pade AC method from ``NSamplePade`` Green's functions with independent Gaussian noise. +The strength of the Gaussian noise is specified by ``g_sigma`` or the pair of ``filein_Gsigma`` and ``column_sigma`` parameters (see :ref:`inputfiles` for details). +The coefficient of the Pade weight, :math:`eta`, is specified by the ``PadeEta`` parameter. +If ``PadeEta = 0``, the original SpM method will be used. + +``param_spm.in`` and ``param_spmpade.in`` are input files for the SpM (:math:`\eta=0`) and SpM-Pade (:math:`\eta=1`) method, respectively. +Note that optimal :math:`\lambda` for the SpM and SpM-Pade differ, and hence users may need to change the range of :math:`lambda` (``lambdalogbegin`` and ``lambdalogend``). +The following figure shows the result: + +.. image:: figs/spectrum_spmpade.png + +The blue curve and red curve show the spectrum reconstructed by the SpM method (:math:`\eta=0`) and the SpM-Pade method (:math:`\eta=1`), respectively. +The black dashed curve shows the exact spectrum. diff --git a/docs/manual/build/html/_sources/index.rst.txt b/docs/manual/build/html/_sources/index.rst.txt index 590459d..b62b7f5 100644 --- a/docs/manual/build/html/_sources/index.rst.txt +++ b/docs/manual/build/html/_sources/index.rst.txt @@ -15,13 +15,26 @@ License -------------- This package is distributed under GNU General Public License version 3 (GPL v3). -We kindly ask you to cite the article +When you publish a publication that includes results obtained using this package, +we kindly ask you to cite the following articles - 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) `_. +- SpM method -in publications that includes results obtained using this package. + 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) `_. + +- SpM 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) `_. + +- SpM-Pade method (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). Contents diff --git a/docs/manual/build/html/_static/ajax-loader.gif b/docs/manual/build/html/_static/ajax-loader.gif deleted file mode 100644 index 61faf8c..0000000 Binary files a/docs/manual/build/html/_static/ajax-loader.gif and /dev/null differ diff --git a/docs/manual/build/html/_static/alabaster.css b/docs/manual/build/html/_static/alabaster.css index 1d1dc4d..3cb5df1 100644 --- a/docs/manual/build/html/_static/alabaster.css +++ b/docs/manual/build/html/_static/alabaster.css @@ -1,55 +1,3 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @import url("basic.css"); /* -- page layout ----------------------------------------------------------- */ @@ -159,7 +107,7 @@ div.sphinxsidebarwrapper p.blurb { div.sphinxsidebar h3, div.sphinxsidebar h4 { - font-family: 'Garamond', 'Georgia', serif; + font-family: Georgia; color: #444; font-size: 24px; font-weight: normal; @@ -218,6 +166,19 @@ div.sphinxsidebar hr { width: 50%; } +div.sphinxsidebar .badge { + border-bottom: none; +} + +div.sphinxsidebar .badge:hover { + border-bottom: none; +} + +/* To address an issue with donation coming after search */ +div.sphinxsidebar h3.donation { + margin-top: 10px; +} + /* -- body styles ----------------------------------------------------------- */ a { @@ -236,7 +197,7 @@ div.body h3, div.body h4, div.body h5, div.body h6 { - font-family: 'Garamond', 'Georgia', serif; + font-family: Georgia; font-weight: normal; margin: 30px 0px 10px 0px; padding: 0; @@ -277,7 +238,7 @@ div.admonition tt.xref, div.admonition code.xref, div.admonition a tt { } div.admonition p.admonition-title { - font-family: 'Garamond', 'Georgia', serif; + font-family: Georgia; font-weight: normal; font-size: 24px; margin: 0 0 10px 0; @@ -366,7 +327,7 @@ p.admonition-title:after { } pre, tt, code { - font-family: 'Consolas', 'Menlo', 'Deja Vu Sans Mono', 'Bitstream Vera Sans Mono', monospace; + font-family: 'Consolas', 'Menlo', 'DejaVu Sans Mono', 'Bitstream Vera Sans Mono', monospace; font-size: 0.9em; } @@ -690,4 +651,51 @@ table.docutils.citation, table.docutils.citation td, table.docutils.citation th -moz-box-shadow: none; -webkit-box-shadow: none; box-shadow: none; +} + + +/* relbar */ + +.related { + line-height: 30px; + width: 100%; + font-size: 0.9rem; +} + +.related.top { + border-bottom: 1px solid #EEE; + margin-bottom: 20px; +} + +.related.bottom { + border-top: 1px solid #EEE; +} + +.related ul { + padding: 0; + margin: 0; + list-style: none; +} + +.related li { + display: inline; +} + +nav#rellinks { + float: right; +} + +nav#rellinks li+li:before { + content: "|"; +} + +nav#breadcrumbs li+li:before { + content: "\00BB"; +} + +/* Hide certain items when printing */ +@media print { + div.related { + display: none; + } } \ No newline at end of file diff --git a/docs/manual/build/html/_static/basic.css b/docs/manual/build/html/_static/basic.css index 19ced10..24bc73e 100644 --- a/docs/manual/build/html/_static/basic.css +++ b/docs/manual/build/html/_static/basic.css @@ -4,7 +4,7 @@ * * Sphinx stylesheet -- basic theme. * - * :copyright: Copyright 2007-2018 by the Sphinx team, see AUTHORS. + * :copyright: Copyright 2007-2020 by the Sphinx team, see AUTHORS. * :license: BSD, see LICENSE for details. * */ @@ -15,6 +15,12 @@ div.clearer { clear: both; } +div.section::after { + display: block; + content: ''; + clear: left; +} + /* -- relbar ---------------------------------------------------------------- */ div.related { @@ -81,6 +87,10 @@ div.sphinxsidebar input { font-size: 1em; } +div.sphinxsidebar #searchbox form.search { + overflow: hidden; +} + div.sphinxsidebar #searchbox input[type="text"] { float: left; width: 80%; @@ -227,6 +237,16 @@ a.headerlink { visibility: hidden; } +a.brackets:before, +span.brackets > a:before{ + content: "["; +} + +a.brackets:after, +span.brackets > a:after { + content: "]"; +} + h1:hover > a.headerlink, h2:hover > a.headerlink, h3:hover > a.headerlink, @@ -275,6 +295,12 @@ img.align-center, .figure.align-center, object.align-center { margin-right: auto; } +img.align-default, .figure.align-default { + display: block; + margin-left: auto; + margin-right: auto; +} + .align-left { text-align: left; } @@ -283,6 +309,10 @@ img.align-center, .figure.align-center, object.align-center { text-align: center; } +.align-default { + text-align: center; +} + .align-right { text-align: right; } @@ -292,21 +322,27 @@ img.align-center, .figure.align-center, object.align-center { div.sidebar { margin: 0 0 0.5em 1em; border: 1px solid #ddb; - padding: 7px 7px 0 7px; + padding: 7px; background-color: #ffe; width: 40%; float: right; + clear: right; + overflow-x: auto; } p.sidebar-title { font-weight: bold; } +div.admonition, div.topic, blockquote { + clear: left; +} + /* -- topics ---------------------------------------------------------------- */ div.topic { border: 1px solid #ccc; - padding: 7px 7px 0 7px; + padding: 7px; margin: 10px 0 10px 0; } @@ -328,10 +364,6 @@ div.admonition dt { font-weight: bold; } -div.admonition dl { - margin-bottom: 0; -} - p.admonition-title { margin: 0px 10px 5px 0px; font-weight: bold; @@ -342,9 +374,28 @@ div.body p.centered { margin-top: 25px; } +/* -- content of sidebars/topics/admonitions -------------------------------- */ + +div.sidebar > :last-child, +div.topic > :last-child, +div.admonition > :last-child { + margin-bottom: 0; +} + +div.sidebar::after, +div.topic::after, +div.admonition::after, +blockquote::after { + display: block; + content: ''; + clear: both; +} + /* -- tables ---------------------------------------------------------------- */ table.docutils { + margin-top: 10px; + margin-bottom: 10px; border: 0; border-collapse: collapse; } @@ -354,6 +405,11 @@ table.align-center { margin-right: auto; } +table.align-default { + margin-left: auto; + margin-right: auto; +} + table caption span.caption-number { font-style: italic; } @@ -387,6 +443,16 @@ table.citation td { border-bottom: none; } +th > :first-child, +td > :first-child { + margin-top: 0px; +} + +th > :last-child, +td > :last-child { + margin-bottom: 0px; +} + /* -- figures --------------------------------------------------------------- */ div.figure { @@ -427,6 +493,17 @@ table.field-list td, table.field-list th { hyphens: manual; } +/* -- hlist styles ---------------------------------------------------------- */ + +table.hlist { + margin: 1em 0; +} + +table.hlist td { + vertical-align: top; +} + + /* -- other body styles ----------------------------------------------------- */ ol.arabic { @@ -449,11 +526,78 @@ ol.upperroman { list-style: upper-roman; } +:not(li) > ol > li:first-child > :first-child, +:not(li) > ul > li:first-child > :first-child { + margin-top: 0px; +} + +:not(li) > ol > li:last-child > :last-child, +:not(li) > ul > li:last-child > :last-child { + margin-bottom: 0px; +} + +ol.simple ol p, +ol.simple ul p, +ul.simple ol p, +ul.simple ul p { + margin-top: 0; +} + +ol.simple > li:not(:first-child) > p, +ul.simple > li:not(:first-child) > p { + margin-top: 0; +} + +ol.simple p, +ul.simple p { + margin-bottom: 0; +} + +dl.footnote > dt, +dl.citation > dt { + float: left; + margin-right: 0.5em; +} + +dl.footnote > dd, +dl.citation > dd { + margin-bottom: 0em; +} + +dl.footnote > dd:after, +dl.citation > dd:after { + content: ""; + clear: both; +} + +dl.field-list { + display: grid; + grid-template-columns: fit-content(30%) auto; +} + +dl.field-list > dt { + font-weight: bold; + word-break: break-word; + padding-left: 0.5em; + padding-right: 5px; +} + +dl.field-list > dt:after { + content: ":"; +} + +dl.field-list > dd { + padding-left: 0.5em; + margin-top: 0em; + margin-left: 0em; + margin-bottom: 0em; +} + dl { margin-bottom: 15px; } -dd p { +dd > :first-child { margin-top: 0px; } @@ -467,6 +611,11 @@ dd { margin-left: 30px; } +dl > dd:last-child, +dl > dd:last-child > :last-child { + margin-bottom: 0; +} + dt:target, span.highlighted { background-color: #fbe54e; } @@ -526,6 +675,12 @@ dl.glossary dt { font-style: oblique; } +.classifier:before { + font-style: normal; + margin: 0.5em; + content: ":"; +} + abbr, acronym { border-bottom: dotted 1px; cursor: help; @@ -538,6 +693,10 @@ pre { overflow-y: hidden; /* fixes display issues on Chrome browsers */ } +pre, div[class*="highlight-"] { + clear: both; +} + span.pre { -moz-hyphens: none; -ms-hyphens: none; @@ -545,22 +704,57 @@ span.pre { hyphens: none; } +div[class*="highlight-"] { + margin: 1em 0; +} + td.linenos pre { - padding: 5px 0px; border: 0; background-color: transparent; color: #aaa; } table.highlighttable { - margin-left: 0.5em; + display: block; +} + +table.highlighttable tbody { + display: block; +} + +table.highlighttable tr { + display: flex; } table.highlighttable td { - padding: 0 0.5em 0 0.5em; + margin: 0; + padding: 0; +} + +table.highlighttable td.linenos { + padding-right: 0.5em; +} + +table.highlighttable td.code { + flex: 1; + overflow: hidden; +} + +.highlight .hll { + display: block; +} + +div.highlight pre, +table.highlighttable pre { + margin: 0; +} + +div.code-block-caption + div { + margin-top: 0; } div.code-block-caption { + margin-top: 1em; padding: 2px 5px; font-size: small; } @@ -569,8 +763,9 @@ div.code-block-caption code { background-color: transparent; } -div.code-block-caption + div > div.highlight > pre { - margin-top: 0; +table.highlighttable td.linenos, +div.doctest > div.highlight span.gp { /* gp: Generic.Prompt */ + user-select: none; } div.code-block-caption span.caption-number { @@ -582,11 +777,7 @@ div.code-block-caption span.caption-text { } div.literal-block-wrapper { - padding: 1em 1em 0; -} - -div.literal-block-wrapper div.highlight { - margin: 0; + margin: 1em 0; } code.descname { @@ -637,8 +828,7 @@ span.eqno { } span.eqno a.headerlink { - position: relative; - left: 0px; + position: absolute; z-index: 1; } diff --git a/docs/manual/build/html/_static/comment-bright.png b/docs/manual/build/html/_static/comment-bright.png deleted file mode 100644 index 15e27ed..0000000 Binary files a/docs/manual/build/html/_static/comment-bright.png and /dev/null differ diff --git a/docs/manual/build/html/_static/comment-close.png b/docs/manual/build/html/_static/comment-close.png deleted file mode 100644 index 4d91bcf..0000000 Binary files a/docs/manual/build/html/_static/comment-close.png and /dev/null differ diff --git a/docs/manual/build/html/_static/comment.png b/docs/manual/build/html/_static/comment.png deleted file mode 100644 index dfbc0cb..0000000 Binary files a/docs/manual/build/html/_static/comment.png and /dev/null differ diff --git a/docs/manual/build/html/_static/doctools.js b/docs/manual/build/html/_static/doctools.js index 0c15c00..daccd20 100644 --- a/docs/manual/build/html/_static/doctools.js +++ b/docs/manual/build/html/_static/doctools.js @@ -4,7 +4,7 @@ * * Sphinx JavaScript utilities for all documentation. * - * :copyright: Copyright 2007-2018 by the Sphinx team, see AUTHORS. + * :copyright: Copyright 2007-2020 by the Sphinx team, see AUTHORS. * :license: BSD, see LICENSE for details. * */ @@ -70,7 +70,9 @@ jQuery.fn.highlightText = function(text, className) { if (node.nodeType === 3) { var val = node.nodeValue; var pos = val.toLowerCase().indexOf(text); - if (pos >= 0 && !jQuery(node.parentNode).hasClass(className)) { + if (pos >= 0 && + !jQuery(node.parentNode).hasClass(className) && + !jQuery(node.parentNode).hasClass("nohighlight")) { var span; var isInSVG = jQuery(node).closest("body, svg, foreignObject").is("svg"); if (isInSVG) { @@ -85,14 +87,13 @@ jQuery.fn.highlightText = function(text, className) { node.nextSibling)); node.nodeValue = val.substr(0, pos); if (isInSVG) { - var bbox = span.getBBox(); var rect = document.createElementNS("http://www.w3.org/2000/svg", "rect"); - rect.x.baseVal.value = bbox.x; + var bbox = node.parentElement.getBBox(); + rect.x.baseVal.value = bbox.x; rect.y.baseVal.value = bbox.y; rect.width.baseVal.value = bbox.width; rect.height.baseVal.value = bbox.height; rect.setAttribute('class', className); - var parentOfText = node.parentNode.parentNode; addItems.push({ "parent": node.parentNode, "target": rect}); @@ -148,7 +149,9 @@ var Documentation = { this.fixFirefoxAnchorBug(); this.highlightSearchWords(); this.initIndexTable(); - + if (DOCUMENTATION_OPTIONS.NAVIGATION_WITH_KEYS) { + this.initOnKeyListeners(); + } }, /** @@ -280,10 +283,11 @@ var Documentation = { }, initOnKeyListeners: function() { - $(document).keyup(function(event) { + $(document).keydown(function(event) { var activeElementType = document.activeElement.tagName; // don't navigate when in search box or textarea - if (activeElementType !== 'TEXTAREA' && activeElementType !== 'INPUT' && activeElementType !== 'SELECT') { + if (activeElementType !== 'TEXTAREA' && activeElementType !== 'INPUT' && activeElementType !== 'SELECT' + && !event.altKey && !event.ctrlKey && !event.metaKey && !event.shiftKey) { switch (event.keyCode) { case 37: // left var prevHref = $('link[rel="prev"]').prop('href'); @@ -308,4 +312,4 @@ _ = Documentation.gettext; $(document).ready(function() { Documentation.init(); -}); \ No newline at end of file +}); diff --git a/docs/manual/build/html/_static/down-pressed.png b/docs/manual/build/html/_static/down-pressed.png deleted file mode 100644 index 5756c8c..0000000 Binary files a/docs/manual/build/html/_static/down-pressed.png and /dev/null differ diff --git a/docs/manual/build/html/_static/down.png b/docs/manual/build/html/_static/down.png deleted file mode 100644 index 1b3bdad..0000000 Binary files a/docs/manual/build/html/_static/down.png and /dev/null differ diff --git a/docs/manual/build/html/_static/jquery-3.1.0.js b/docs/manual/build/html/_static/jquery-3.1.0.js deleted file mode 100644 index f2fc274..0000000 --- a/docs/manual/build/html/_static/jquery-3.1.0.js +++ /dev/null @@ -1,10074 +0,0 @@ -/*eslint-disable no-unused-vars*/ -/*! - * jQuery JavaScript Library v3.1.0 - * https://jquery.com/ - * - * Includes Sizzle.js - * https://sizzlejs.com/ - * - * Copyright jQuery Foundation and other contributors - * Released under the MIT license - * https://jquery.org/license - * - * Date: 2016-07-07T21:44Z - */ -( function( global, factory ) { - - "use strict"; - - if ( typeof module === "object" && typeof module.exports === "object" ) { - - // For CommonJS and CommonJS-like environments where a proper `window` - // is present, execute the factory and get jQuery. - // For environments that do not have a `window` with a `document` - // (such as Node.js), expose a factory as module.exports. - // This accentuates the need for the creation of a real `window`. - // e.g. var jQuery = require("jquery")(window); - // See ticket #14549 for more info. - module.exports = global.document ? - factory( global, true ) : - function( w ) { - if ( !w.document ) { - throw new Error( "jQuery requires a window with a document" ); - } - return factory( w ); - }; - } else { - factory( global ); - } - -// Pass this if window is not defined yet -} )( typeof window !== "undefined" ? window : this, function( window, noGlobal ) { - -// Edge <= 12 - 13+, Firefox <=18 - 45+, IE 10 - 11, Safari 5.1 - 9+, iOS 6 - 9.1 -// throw exceptions when non-strict code (e.g., ASP.NET 4.5) accesses strict mode -// arguments.callee.caller (trac-13335). But as of jQuery 3.0 (2016), strict mode should be common -// enough that all such attempts are guarded in a try block. -"use strict"; - -var arr = []; - -var document = window.document; - -var getProto = Object.getPrototypeOf; - -var slice = arr.slice; - -var concat = arr.concat; - -var push = arr.push; - -var indexOf = arr.indexOf; - -var class2type = {}; - -var toString = class2type.toString; - -var hasOwn = class2type.hasOwnProperty; - -var fnToString = hasOwn.toString; - -var ObjectFunctionString = fnToString.call( Object ); - -var support = {}; - - - - function DOMEval( code, doc ) { - doc = doc || document; - - var script = doc.createElement( "script" ); - - script.text = code; - doc.head.appendChild( script ).parentNode.removeChild( script ); - } -/* global Symbol */ -// Defining this global in .eslintrc would create a danger of using the global -// unguarded in another place, it seems safer to define global only for this module - - - -var - version = "3.1.0", - - // Define a local copy of jQuery - jQuery = function( selector, context ) { - - // The jQuery object is actually just the init constructor 'enhanced' - // Need init if jQuery is called (just allow error to be thrown if not included) - return new jQuery.fn.init( selector, context ); - }, - - // Support: Android <=4.0 only - // Make sure we trim BOM and NBSP - rtrim = /^[\s\uFEFF\xA0]+|[\s\uFEFF\xA0]+$/g, - - // Matches dashed string for camelizing - rmsPrefix = /^-ms-/, - rdashAlpha = /-([a-z])/g, - - // Used by jQuery.camelCase as callback to replace() - fcamelCase = function( all, letter ) { - return letter.toUpperCase(); - }; - -jQuery.fn = jQuery.prototype = { - - // The current version of jQuery being used - jquery: version, - - constructor: jQuery, - - // The default length of a jQuery object is 0 - length: 0, - - toArray: function() { - return slice.call( this ); - }, - - // Get the Nth element in the matched element set OR - // Get the whole matched element set as a clean array - get: function( num ) { - return num != null ? - - // Return just the one element from the set - ( num < 0 ? this[ num + this.length ] : this[ num ] ) : - - // Return all the elements in a clean array - slice.call( this ); - }, - - // Take an array of elements and push it onto the stack - // (returning the new matched element set) - pushStack: function( elems ) { - - // Build a new jQuery matched element set - var ret = jQuery.merge( this.constructor(), elems ); - - // Add the old object onto the stack (as a reference) - ret.prevObject = this; - - // Return the newly-formed element set - return ret; - }, - - // Execute a callback for every element in the matched set. - each: function( callback ) { - return jQuery.each( this, callback ); - }, - - map: function( callback ) { - return this.pushStack( jQuery.map( this, function( elem, i ) { - return callback.call( elem, i, elem ); - } ) ); - }, - - slice: function() { - return this.pushStack( slice.apply( this, arguments ) ); - }, - - first: function() { - return this.eq( 0 ); - }, - - last: function() { - return this.eq( -1 ); - }, - - eq: function( i ) { - var len = this.length, - j = +i + ( i < 0 ? len : 0 ); - return this.pushStack( j >= 0 && j < len ? [ this[ j ] ] : [] ); - }, - - end: function() { - return this.prevObject || this.constructor(); - }, - - // For internal use only. - // Behaves like an Array's method, not like a jQuery method. - push: push, - sort: arr.sort, - splice: arr.splice -}; - -jQuery.extend = jQuery.fn.extend = function() { - var options, name, src, copy, copyIsArray, clone, - target = arguments[ 0 ] || {}, - i = 1, - length = arguments.length, - deep = false; - - // Handle a deep copy situation - if ( typeof target === "boolean" ) { - deep = target; - - // Skip the boolean and the target - target = arguments[ i ] || {}; - i++; - } - - // Handle case when target is a string or something (possible in deep copy) - if ( typeof target !== "object" && !jQuery.isFunction( target ) ) { - target = {}; - } - - // Extend jQuery itself if only one argument is passed - if ( i === length ) { - target = this; - i--; - } - - for ( ; i < length; i++ ) { - - // Only deal with non-null/undefined values - if ( ( options = arguments[ i ] ) != null ) { - - // Extend the base object - for ( name in options ) { - src = target[ name ]; - copy = options[ name ]; - - // Prevent never-ending loop - if ( target === copy ) { - continue; - } - - // Recurse if we're merging plain objects or arrays - if ( deep && copy && ( jQuery.isPlainObject( copy ) || - ( copyIsArray = jQuery.isArray( copy ) ) ) ) { - - if ( copyIsArray ) { - copyIsArray = false; - clone = src && jQuery.isArray( src ) ? src : []; - - } else { - clone = src && jQuery.isPlainObject( src ) ? src : {}; - } - - // Never move original objects, clone them - target[ name ] = jQuery.extend( deep, clone, copy ); - - // Don't bring in undefined values - } else if ( copy !== undefined ) { - target[ name ] = copy; - } - } - } - } - - // Return the modified object - return target; -}; - -jQuery.extend( { - - // Unique for each copy of jQuery on the page - expando: "jQuery" + ( version + Math.random() ).replace( /\D/g, "" ), - - // Assume jQuery is ready without the ready module - isReady: true, - - error: function( msg ) { - throw new Error( msg ); - }, - - noop: function() {}, - - isFunction: function( obj ) { - return jQuery.type( obj ) === "function"; - }, - - isArray: Array.isArray, - - isWindow: function( obj ) { - return obj != null && obj === obj.window; - }, - - isNumeric: function( obj ) { - - // As of jQuery 3.0, isNumeric is limited to - // strings and numbers (primitives or objects) - // that can be coerced to finite numbers (gh-2662) - var type = jQuery.type( obj ); - return ( type === "number" || type === "string" ) && - - // parseFloat NaNs numeric-cast false positives ("") - // ...but misinterprets leading-number strings, particularly hex literals ("0x...") - // subtraction forces infinities to NaN - !isNaN( obj - parseFloat( obj ) ); - }, - - isPlainObject: function( obj ) { - var proto, Ctor; - - // Detect obvious negatives - // Use toString instead of jQuery.type to catch host objects - if ( !obj || toString.call( obj ) !== "[object Object]" ) { - return false; - } - - proto = getProto( obj ); - - // Objects with no prototype (e.g., `Object.create( null )`) are plain - if ( !proto ) { - return true; - } - - // Objects with prototype are plain iff they were constructed by a global Object function - Ctor = hasOwn.call( proto, "constructor" ) && proto.constructor; - return typeof Ctor === "function" && fnToString.call( Ctor ) === ObjectFunctionString; - }, - - isEmptyObject: function( obj ) { - - /* eslint-disable no-unused-vars */ - // See https://github.com/eslint/eslint/issues/6125 - var name; - - for ( name in obj ) { - return false; - } - return true; - }, - - type: function( obj ) { - if ( obj == null ) { - return obj + ""; - } - - // Support: Android <=2.3 only (functionish RegExp) - return typeof obj === "object" || typeof obj === "function" ? - class2type[ toString.call( obj ) ] || "object" : - typeof obj; - }, - - // Evaluates a script in a global context - globalEval: function( code ) { - DOMEval( code ); - }, - - // Convert dashed to camelCase; used by the css and data modules - // Support: IE <=9 - 11, Edge 12 - 13 - // Microsoft forgot to hump their vendor prefix (#9572) - camelCase: function( string ) { - return string.replace( rmsPrefix, "ms-" ).replace( rdashAlpha, fcamelCase ); - }, - - nodeName: function( elem, name ) { - return elem.nodeName && elem.nodeName.toLowerCase() === name.toLowerCase(); - }, - - each: function( obj, callback ) { - var length, i = 0; - - if ( isArrayLike( obj ) ) { - length = obj.length; - for ( ; i < length; i++ ) { - if ( callback.call( obj[ i ], i, obj[ i ] ) === false ) { - break; - } - } - } else { - for ( i in obj ) { - if ( callback.call( obj[ i ], i, obj[ i ] ) === false ) { - break; - } - } - } - - return obj; - }, - - // Support: Android <=4.0 only - trim: function( text ) { - return text == null ? - "" : - ( text + "" ).replace( rtrim, "" ); - }, - - // results is for internal usage only - makeArray: function( arr, results ) { - var ret = results || []; - - if ( arr != null ) { - if ( isArrayLike( Object( arr ) ) ) { - jQuery.merge( ret, - typeof arr === "string" ? - [ arr ] : arr - ); - } else { - push.call( ret, arr ); - } - } - - return ret; - }, - - inArray: function( elem, arr, i ) { - return arr == null ? -1 : indexOf.call( arr, elem, i ); - }, - - // Support: Android <=4.0 only, PhantomJS 1 only - // push.apply(_, arraylike) throws on ancient WebKit - merge: function( first, second ) { - var len = +second.length, - j = 0, - i = first.length; - - for ( ; j < len; j++ ) { - first[ i++ ] = second[ j ]; - } - - first.length = i; - - return first; - }, - - grep: function( elems, callback, invert ) { - var callbackInverse, - matches = [], - i = 0, - length = elems.length, - callbackExpect = !invert; - - // Go through the array, only saving the items - // that pass the validator function - for ( ; i < length; i++ ) { - callbackInverse = !callback( elems[ i ], i ); - if ( callbackInverse !== callbackExpect ) { - matches.push( elems[ i ] ); - } - } - - return matches; - }, - - // arg is for internal usage only - map: function( elems, callback, arg ) { - var length, value, - i = 0, - ret = []; - - // Go through the array, translating each of the items to their new values - if ( isArrayLike( elems ) ) { - length = elems.length; - for ( ; i < length; i++ ) { - value = callback( elems[ i ], i, arg ); - - if ( value != null ) { - ret.push( value ); - } - } - - // Go through every key on the object, - } else { - for ( i in elems ) { - value = callback( elems[ i ], i, arg ); - - if ( value != null ) { - ret.push( value ); - } - } - } - - // Flatten any nested arrays - return concat.apply( [], ret ); - }, - - // A global GUID counter for objects - guid: 1, - - // Bind a function to a context, optionally partially applying any - // arguments. - proxy: function( fn, context ) { - var tmp, args, proxy; - - if ( typeof context === "string" ) { - tmp = fn[ context ]; - context = fn; - fn = tmp; - } - - // Quick check to determine if target is callable, in the spec - // this throws a TypeError, but we will just return undefined. - if ( !jQuery.isFunction( fn ) ) { - return undefined; - } - - // Simulated bind - args = slice.call( arguments, 2 ); - proxy = function() { - return fn.apply( context || this, args.concat( slice.call( arguments ) ) ); - }; - - // Set the guid of unique handler to the same of original handler, so it can be removed - proxy.guid = fn.guid = fn.guid || jQuery.guid++; - - return proxy; - }, - - now: Date.now, - - // jQuery.support is not used in Core but other projects attach their - // properties to it so it needs to exist. - support: support -} ); - -if ( typeof Symbol === "function" ) { - jQuery.fn[ Symbol.iterator ] = arr[ Symbol.iterator ]; -} - -// Populate the class2type map -jQuery.each( "Boolean Number String Function Array Date RegExp Object Error Symbol".split( " " ), -function( i, name ) { - class2type[ "[object " + name + "]" ] = name.toLowerCase(); -} ); - -function isArrayLike( obj ) { - - // Support: real iOS 8.2 only (not reproducible in simulator) - // `in` check used to prevent JIT error (gh-2145) - // hasOwn isn't used here due to false negatives - // regarding Nodelist length in IE - var length = !!obj && "length" in obj && obj.length, - type = jQuery.type( obj ); - - if ( type === "function" || jQuery.isWindow( obj ) ) { - return false; - } - - return type === "array" || length === 0 || - typeof length === "number" && length > 0 && ( length - 1 ) in obj; -} -var Sizzle = -/*! - * Sizzle CSS Selector Engine v2.3.0 - * https://sizzlejs.com/ - * - * Copyright jQuery Foundation and other contributors - * Released under the MIT license - * http://jquery.org/license - * - * Date: 2016-01-04 - */ -(function( window ) { - -var i, - support, - Expr, - getText, - isXML, - tokenize, - compile, - select, - outermostContext, - sortInput, - hasDuplicate, - - // Local document vars - setDocument, - document, - docElem, - documentIsHTML, - rbuggyQSA, - rbuggyMatches, - matches, - contains, - - // Instance-specific data - expando = "sizzle" + 1 * new Date(), - preferredDoc = window.document, - dirruns = 0, - done = 0, - classCache = createCache(), - tokenCache = createCache(), - compilerCache = createCache(), - sortOrder = function( a, b ) { - if ( a === b ) { - hasDuplicate = true; - } - return 0; - }, - - // Instance methods - hasOwn = ({}).hasOwnProperty, - arr = [], - pop = arr.pop, - push_native = arr.push, - push = arr.push, - slice = arr.slice, - // Use a stripped-down indexOf as it's faster than native - // https://jsperf.com/thor-indexof-vs-for/5 - indexOf = function( list, elem ) { - var i = 0, - len = list.length; - for ( ; i < len; i++ ) { - if ( list[i] === elem ) { - return i; - } - } - return -1; - }, - - booleans = "checked|selected|async|autofocus|autoplay|controls|defer|disabled|hidden|ismap|loop|multiple|open|readonly|required|scoped", - - // Regular expressions - - // http://www.w3.org/TR/css3-selectors/#whitespace - whitespace = "[\\x20\\t\\r\\n\\f]", - - // http://www.w3.org/TR/CSS21/syndata.html#value-def-identifier - identifier = "(?:\\\\.|[\\w-]|[^\0-\\xa0])+", - - // Attribute selectors: http://www.w3.org/TR/selectors/#attribute-selectors - attributes = "\\[" + whitespace + "*(" + identifier + ")(?:" + whitespace + - // Operator (capture 2) - "*([*^$|!~]?=)" + whitespace + - // "Attribute values must be CSS identifiers [capture 5] or strings [capture 3 or capture 4]" - "*(?:'((?:\\\\.|[^\\\\'])*)'|\"((?:\\\\.|[^\\\\\"])*)\"|(" + identifier + "))|)" + whitespace + - "*\\]", - - pseudos = ":(" + identifier + ")(?:\\((" + - // To reduce the number of selectors needing tokenize in the preFilter, prefer arguments: - // 1. quoted (capture 3; capture 4 or capture 5) - "('((?:\\\\.|[^\\\\'])*)'|\"((?:\\\\.|[^\\\\\"])*)\")|" + - // 2. simple (capture 6) - "((?:\\\\.|[^\\\\()[\\]]|" + attributes + ")*)|" + - // 3. anything else (capture 2) - ".*" + - ")\\)|)", - - // Leading and non-escaped trailing whitespace, capturing some non-whitespace characters preceding the latter - rwhitespace = new RegExp( whitespace + "+", "g" ), - rtrim = new RegExp( "^" + whitespace + "+|((?:^|[^\\\\])(?:\\\\.)*)" + whitespace + "+$", "g" ), - - rcomma = new RegExp( "^" + whitespace + "*," + whitespace + "*" ), - rcombinators = new RegExp( "^" + whitespace + "*([>+~]|" + whitespace + ")" + whitespace + "*" ), - - rattributeQuotes = new RegExp( "=" + whitespace + "*([^\\]'\"]*?)" + whitespace + "*\\]", "g" ), - - rpseudo = new RegExp( pseudos ), - ridentifier = new RegExp( "^" + identifier + "$" ), - - matchExpr = { - "ID": new RegExp( "^#(" + identifier + ")" ), - "CLASS": new RegExp( "^\\.(" + identifier + ")" ), - "TAG": new RegExp( "^(" + identifier + "|[*])" ), - "ATTR": new RegExp( "^" + attributes ), - "PSEUDO": new RegExp( "^" + pseudos ), - "CHILD": new RegExp( "^:(only|first|last|nth|nth-last)-(child|of-type)(?:\\(" + whitespace + - "*(even|odd|(([+-]|)(\\d*)n|)" + whitespace + "*(?:([+-]|)" + whitespace + - "*(\\d+)|))" + whitespace + "*\\)|)", "i" ), - "bool": new RegExp( "^(?:" + booleans + ")$", "i" ), - // For use in libraries implementing .is() - // We use this for POS matching in `select` - "needsContext": new RegExp( "^" + whitespace + "*[>+~]|:(even|odd|eq|gt|lt|nth|first|last)(?:\\(" + - whitespace + "*((?:-\\d)?\\d*)" + whitespace + "*\\)|)(?=[^-]|$)", "i" ) - }, - - rinputs = /^(?:input|select|textarea|button)$/i, - rheader = /^h\d$/i, - - rnative = /^[^{]+\{\s*\[native \w/, - - // Easily-parseable/retrievable ID or TAG or CLASS selectors - rquickExpr = /^(?:#([\w-]+)|(\w+)|\.([\w-]+))$/, - - rsibling = /[+~]/, - - // CSS escapes - // http://www.w3.org/TR/CSS21/syndata.html#escaped-characters - runescape = new RegExp( "\\\\([\\da-f]{1,6}" + whitespace + "?|(" + whitespace + ")|.)", "ig" ), - funescape = function( _, escaped, escapedWhitespace ) { - var high = "0x" + escaped - 0x10000; - // NaN means non-codepoint - // Support: Firefox<24 - // Workaround erroneous numeric interpretation of +"0x" - return high !== high || escapedWhitespace ? - escaped : - high < 0 ? - // BMP codepoint - String.fromCharCode( high + 0x10000 ) : - // Supplemental Plane codepoint (surrogate pair) - String.fromCharCode( high >> 10 | 0xD800, high & 0x3FF | 0xDC00 ); - }, - - // CSS string/identifier serialization - // https://drafts.csswg.org/cssom/#common-serializing-idioms - rcssescape = /([\0-\x1f\x7f]|^-?\d)|^-$|[^\x80-\uFFFF\w-]/g, - fcssescape = function( ch, asCodePoint ) { - if ( asCodePoint ) { - - // U+0000 NULL becomes U+FFFD REPLACEMENT CHARACTER - if ( ch === "\0" ) { - return "\uFFFD"; - } - - // Control characters and (dependent upon position) numbers get escaped as code points - return ch.slice( 0, -1 ) + "\\" + ch.charCodeAt( ch.length - 1 ).toString( 16 ) + " "; - } - - // Other potentially-special ASCII characters get backslash-escaped - return "\\" + ch; - }, - - // Used for iframes - // See setDocument() - // Removing the function wrapper causes a "Permission Denied" - // error in IE - unloadHandler = function() { - setDocument(); - }, - - disabledAncestor = addCombinator( - function( elem ) { - return elem.disabled === true; - }, - { dir: "parentNode", next: "legend" } - ); - -// Optimize for push.apply( _, NodeList ) -try { - push.apply( - (arr = slice.call( preferredDoc.childNodes )), - preferredDoc.childNodes - ); - // Support: Android<4.0 - // Detect silently failing push.apply - arr[ preferredDoc.childNodes.length ].nodeType; -} catch ( e ) { - push = { apply: arr.length ? - - // Leverage slice if possible - function( target, els ) { - push_native.apply( target, slice.call(els) ); - } : - - // Support: IE<9 - // Otherwise append directly - function( target, els ) { - var j = target.length, - i = 0; - // Can't trust NodeList.length - while ( (target[j++] = els[i++]) ) {} - target.length = j - 1; - } - }; -} - -function Sizzle( selector, context, results, seed ) { - var m, i, elem, nid, match, groups, newSelector, - newContext = context && context.ownerDocument, - - // nodeType defaults to 9, since context defaults to document - nodeType = context ? context.nodeType : 9; - - results = results || []; - - // Return early from calls with invalid selector or context - if ( typeof selector !== "string" || !selector || - nodeType !== 1 && nodeType !== 9 && nodeType !== 11 ) { - - return results; - } - - // Try to shortcut find operations (as opposed to filters) in HTML documents - if ( !seed ) { - - if ( ( context ? context.ownerDocument || context : preferredDoc ) !== document ) { - setDocument( context ); - } - context = context || document; - - if ( documentIsHTML ) { - - // If the selector is sufficiently simple, try using a "get*By*" DOM method - // (excepting DocumentFragment context, where the methods don't exist) - if ( nodeType !== 11 && (match = rquickExpr.exec( selector )) ) { - - // ID selector - if ( (m = match[1]) ) { - - // Document context - if ( nodeType === 9 ) { - if ( (elem = context.getElementById( m )) ) { - - // Support: IE, Opera, Webkit - // TODO: identify versions - // getElementById can match elements by name instead of ID - if ( elem.id === m ) { - results.push( elem ); - return results; - } - } else { - return results; - } - - // Element context - } else { - - // Support: IE, Opera, Webkit - // TODO: identify versions - // getElementById can match elements by name instead of ID - if ( newContext && (elem = newContext.getElementById( m )) && - contains( context, elem ) && - elem.id === m ) { - - results.push( elem ); - return results; - } - } - - // Type selector - } else if ( match[2] ) { - push.apply( results, context.getElementsByTagName( selector ) ); - return results; - - // Class selector - } else if ( (m = match[3]) && support.getElementsByClassName && - context.getElementsByClassName ) { - - push.apply( results, context.getElementsByClassName( m ) ); - return results; - } - } - - // Take advantage of querySelectorAll - if ( support.qsa && - !compilerCache[ selector + " " ] && - (!rbuggyQSA || !rbuggyQSA.test( selector )) ) { - - if ( nodeType !== 1 ) { - newContext = context; - newSelector = selector; - - // qSA looks outside Element context, which is not what we want - // Thanks to Andrew Dupont for this workaround technique - // Support: IE <=8 - // Exclude object elements - } else if ( context.nodeName.toLowerCase() !== "object" ) { - - // Capture the context ID, setting it first if necessary - if ( (nid = context.getAttribute( "id" )) ) { - nid = nid.replace( rcssescape, fcssescape ); - } else { - context.setAttribute( "id", (nid = expando) ); - } - - // Prefix every selector in the list - groups = tokenize( selector ); - i = groups.length; - while ( i-- ) { - groups[i] = "#" + nid + " " + toSelector( groups[i] ); - } - newSelector = groups.join( "," ); - - // Expand context for sibling selectors - newContext = rsibling.test( selector ) && testContext( context.parentNode ) || - context; - } - - if ( newSelector ) { - try { - push.apply( results, - newContext.querySelectorAll( newSelector ) - ); - return results; - } catch ( qsaError ) { - } finally { - if ( nid === expando ) { - context.removeAttribute( "id" ); - } - } - } - } - } - } - - // All others - return select( selector.replace( rtrim, "$1" ), context, results, seed ); -} - -/** - * Create key-value caches of limited size - * @returns {function(string, object)} Returns the Object data after storing it on itself with - * property name the (space-suffixed) string and (if the cache is larger than Expr.cacheLength) - * deleting the oldest entry - */ -function createCache() { - var keys = []; - - function cache( key, value ) { - // Use (key + " ") to avoid collision with native prototype properties (see Issue #157) - if ( keys.push( key + " " ) > Expr.cacheLength ) { - // Only keep the most recent entries - delete cache[ keys.shift() ]; - } - return (cache[ key + " " ] = value); - } - return cache; -} - -/** - * Mark a function for special use by Sizzle - * @param {Function} fn The function to mark - */ -function markFunction( fn ) { - fn[ expando ] = true; - return fn; -} - -/** - * Support testing using an element - * @param {Function} fn Passed the created element and returns a boolean result - */ -function assert( fn ) { - var el = document.createElement("fieldset"); - - try { - return !!fn( el ); - } catch (e) { - return false; - } finally { - // Remove from its parent by default - if ( el.parentNode ) { - el.parentNode.removeChild( el ); - } - // release memory in IE - el = null; - } -} - -/** - * Adds the same handler for all of the specified attrs - * @param {String} attrs Pipe-separated list of attributes - * @param {Function} handler The method that will be applied - */ -function addHandle( attrs, handler ) { - var arr = attrs.split("|"), - i = arr.length; - - while ( i-- ) { - Expr.attrHandle[ arr[i] ] = handler; - } -} - -/** - * Checks document order of two siblings - * @param {Element} a - * @param {Element} b - * @returns {Number} Returns less than 0 if a precedes b, greater than 0 if a follows b - */ -function siblingCheck( a, b ) { - var cur = b && a, - diff = cur && a.nodeType === 1 && b.nodeType === 1 && - a.sourceIndex - b.sourceIndex; - - // Use IE sourceIndex if available on both nodes - if ( diff ) { - return diff; - } - - // Check if b follows a - if ( cur ) { - while ( (cur = cur.nextSibling) ) { - if ( cur === b ) { - return -1; - } - } - } - - return a ? 1 : -1; -} - -/** - * Returns a function to use in pseudos for input types - * @param {String} type - */ -function createInputPseudo( type ) { - return function( elem ) { - var name = elem.nodeName.toLowerCase(); - return name === "input" && elem.type === type; - }; -} - -/** - * Returns a function to use in pseudos for buttons - * @param {String} type - */ -function createButtonPseudo( type ) { - return function( elem ) { - var name = elem.nodeName.toLowerCase(); - return (name === "input" || name === "button") && elem.type === type; - }; -} - -/** - * Returns a function to use in pseudos for :enabled/:disabled - * @param {Boolean} disabled true for :disabled; false for :enabled - */ -function createDisabledPseudo( disabled ) { - // Known :disabled false positives: - // IE: *[disabled]:not(button, input, select, textarea, optgroup, option, menuitem, fieldset) - // not IE: fieldset[disabled] > legend:nth-of-type(n+2) :can-disable - return function( elem ) { - - // Check form elements and option elements for explicit disabling - return "label" in elem && elem.disabled === disabled || - "form" in elem && elem.disabled === disabled || - - // Check non-disabled form elements for fieldset[disabled] ancestors - "form" in elem && elem.disabled === false && ( - // Support: IE6-11+ - // Ancestry is covered for us - elem.isDisabled === disabled || - - // Otherwise, assume any non-