From 2c6da6a35a08d287f10579f422b826420ac12106 Mon Sep 17 00:00:00 2001
From: tayral <thomas.ayral@cpht.polytechnique.fr>
Date: Wed, 12 Mar 2014 11:28:25 +0000
Subject: [PATCH] Made autocorrtime calculation from binning faster.

Bin size is doubled at each step to converge faster and reuse previous binned series.
---
 .../triqs/statistics/autocorrelation_time.cpp | 12 ++++++-
 .../statistics/autocorrelation_time.output    |  8 +++--
 triqs/statistics/statistics.hpp               | 33 ++++++++++++++++++-
 3 files changed, 48 insertions(+), 5 deletions(-)

diff --git a/test/triqs/statistics/autocorrelation_time.cpp b/test/triqs/statistics/autocorrelation_time.cpp
index 216740b1..d2264b62 100644
--- a/test/triqs/statistics/autocorrelation_time.cpp
+++ b/test/triqs/statistics/autocorrelation_time.cpp
@@ -2,6 +2,7 @@
 #include <triqs/statistics.hpp>
 #include "./correlated_gaussian.hpp"
 #include <iostream>
+#include <ctime>
 #define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> "<< (X) <<std::endl<<std::endl;
 using namespace triqs::statistics;
 
@@ -31,13 +32,20 @@ void test_1(int argc, char ** argv){
  correlated_gaussian_vector(A, seed, L);
  double intrinsic_variance = 1;
 
+ auto t1 = clock();
  TEST( autocorrelation_time(A));
+ //std::cout << "time = " << double( clock()-t1)/CLOCKS_PER_SEC << std::endl;
+ //t1 = clock();
  TEST( autocorrelation_time_from_binning(A,intrinsic_variance));
  TEST( autocorrelation_time_from_binning(A));
+ //std::cout << "time = " << double( clock()-t1)/CLOCKS_PER_SEC << std::endl;
+ //t1 = clock();
+ //TEST( autocorrelation_time_from_binning2(A));
+ //std::cout << "time = " << double( clock()-t1)/CLOCKS_PER_SEC << std::endl;
 }
 
 void test_2(int argc, char ** argv){
- int N=10000, L=40;
+ int N=100000, L=40;
  if(argc==3){
   N = atoi(argv[1]); //size
   L = atoi(argv[2]); //autocorrelation time
@@ -50,6 +58,8 @@ void test_2(int argc, char ** argv){
  observable<double> V;
  for (auto & x:A) V << x;
  TEST(autocorrelation_time(V));
+ TEST(autocorrelation_time_from_binning(V));
+ //TEST(autocorrelation_time_from_binning2(V));
  TEST(autocorrelation_time(V*V)); 
 }
 
diff --git a/test/triqs/statistics/autocorrelation_time.output b/test/triqs/statistics/autocorrelation_time.output
index 0cd6af3b..9a8d704d 100644
--- a/test/triqs/statistics/autocorrelation_time.output
+++ b/test/triqs/statistics/autocorrelation_time.output
@@ -6,9 +6,11 @@ L = 100
 
 (autocorrelation_time_from_binning(A,intrinsic_variance)) ---> 80.2806
 
-(autocorrelation_time_from_binning(A)) ---> 106.764
+(autocorrelation_time_from_binning(A)) ---> 95.7671
 
-(autocorrelation_time(V)) ---> 41
+(autocorrelation_time(V)) ---> 40
 
-(autocorrelation_time(V*V)) ---> 40
+(autocorrelation_time_from_binning(V)) ---> 38.0031
+
+(autocorrelation_time(V*V)) ---> 37
 
diff --git a/triqs/statistics/statistics.hpp b/triqs/statistics/statistics.hpp
index 73726429..93d3168c 100644
--- a/triqs/statistics/statistics.hpp
+++ b/triqs/statistics/statistics.hpp
@@ -437,7 +437,7 @@ namespace statistics {
 
  // ------  Auto-correlation time from binning --------------------
 
- template <typename TimeSeries> double autocorrelation_time_from_binning(TimeSeries const& A) {
+ template <typename TimeSeries> double autocorrelation_time_from_binning2(TimeSeries const& A) {
 
   auto size = make_immutable_time_series(A).size();
   double var1 = empirical_variance(make_immutable_time_series(A));
@@ -462,5 +462,36 @@ namespace statistics {
   }
   return empirical_average(t);
  }
+
+ template<typename TimeSeries>
+ double t_cor(TimeSeries const & A, int bin_size, double var1){
+   double var = empirical_variance(A);
+   return .5 * var / var1 * bin_size; 
+ }
+
+ template <typename TimeSeries> double autocorrelation_time_from_binning(TimeSeries const& A) {
+
+  auto size = make_immutable_time_series(A).size();
+  double var1 = empirical_variance(make_immutable_time_series(A));
+
+  int B = 2;
+  auto Ab=make_binned_series(A,2);
+  double autocorr_time = t_cor(Ab, B, var1);
+  double slope = 1.;
+  int small_slope_count = 0;
+  std::vector<double> t;
+  while (small_slope_count < 5 && B < size / 10) {
+   B*=2;
+   Ab=make_binned_series(Ab,2);
+   double t_cor_new = t_cor(Ab, B, var1);
+   slope = (std::abs(t_cor_new - autocorr_time) / autocorr_time);
+   if (slope < .5*1e-1) small_slope_count++;
+   if (small_slope_count > 0) t.push_back(t_cor_new);
+   autocorr_time = t_cor_new;
+   //std::cout << B << "\t" << t_cor_new << "\t" << slope << std::endl;
+  }
+  if (t.size()>0)  {return empirical_average(t);}
+  else{ std::cout << "autocorrelation time not converged!!" << std::endl; return autocorr_time;}
+ }
 } // namespace statistics
 } // triqs