text: # -*- coding: utf-8 -*-
  #
  # Copyright (C) 2008-2014 Jonathan F. Donges
  # Author: Jonathan F. Donges <donges@pik-potsdam.de>
  # URL: <http://www.pik-potsdam.de/members/donges/software>
  
  """
  Computes and plots auto-correlation functions for a COPRA ensemble.
  
  This script provides analyses for this publication:
  
  J.F. Donges, R.V. Donner, N. Marwan, S.F.M. Breitenbach, K. Rehfeld, and J. Kurths,
  Nonlinear regime shifts in Holocene Asian monsoon variability: Potential impacts on cultural change and migratory patterns,
  Climate of the Past 11, 709-741 (2015),
  DOI: 10.5194/cp-11-709-2015
  """
  
  #
  #  Imports
  #
  
  import sys
  
  import numpy as np
  import pylab
  
  from matplotlib import rc
  rc('text', usetex=True)
  
  #  Import progress bar for easy progress bar handling
  import progressbar
  
  from pyunicorn.funcnet import CouplingAnalysis
  
  #
  #  Settings
  #
  
  #  Name of data directory
  DATA_DIR = "../../data/copra_ensemble_data/"
  
  #  Give filename
  FILENAME = "dongge_realisations.txt"
  
  #  Specify symbol used for commenting in data file
  COMMENT_SYMBOL = "%"
  
  #  Maximum delay
  DELAY_TIME_MAX = 1000. # in units of years
  
  #  Settings for detrending
  DETREND = True
  DETRENDING_WINDOW_SIZE = 1000. # measured in [a] / [ka]
  
  #
  #  Functions
  #
  
  def detrend_time_series(data, window_size):
      """
      """
      #  Get length of data array
      n = data.shape[0]
      #  Initialize a local copy of data array
      detrended_data = np.empty(n)
  
      #  Detrend data
      for j in xrange(n):
          #  Get distance of sample from boundaries of time series
          dist = min(j, n - 1 - j)
  
          if window_size / 2 > dist:
              half_size = dist
          else:
              half_size = window_size / 2
  
          detrended_data[j] = data[j] - data[j - half_size:j + half_size + 1].mean()
  
      return detrended_data
  
  #
  #  Main script
  #
  
  #
  #  The main script
  #
  
  print "Analyzing ensemble of time series generated by COPRA"
  print "----------------------------------------------------"
  
  #
  #  Import data
  #
  
  print "Loading data..."
  
  #  Load data
  data = np.loadtxt(DATA_DIR + FILENAME, comments=COMMENT_SYMBOL).transpose()
  
  #  Look at smaller number of 100 realizations for now!
  time = data[0, ::10]
  values = data[1:1001, ::10]
  
  #  Remove NaNs
  is_not_nan = np.all(np.logical_not(np.isnan(values)), axis=0)
  time = time[is_not_nan]
  values = values[:, is_not_nan]
  
  #  Get time series properties
  average_sampling_time = np.diff(time).mean()
  n_realizations = values.shape[0]
  
  DELAY_MAX = int(DELAY_TIME_MAX / average_sampling_time)
  
  if DETREND:
      #  Detrend data!
      print "Detrending data..."
  
      for i in xrange(n_realizations):
          values[i,:] = detrend_time_series(data=values[i,:], window_size=DETRENDING_WINDOW_SIZE / average_sampling_time)
  
  #  Initialize
  acf = np.zeros((n_realizations, DELAY_MAX + 1))
  
  print "Computing auto-correlation functions..."
  
  #  Initialize progress bar
  progress = progressbar.ProgressBar().start()
  
  #  Compute autocorrelation functions
  for i in xrange(n_realizations):
  	#  Reshape
  	series = values[i,:].copy()
  	series.shape = (-1,1)
  
  	#  Initialize instance of CouplingAnalysis class
  	ca = CouplingAnalysis(data=series)
  
  	#  Calculate auto correlation function
  	acf[i,:] = ca.cross_correlation(tau_max=DELAY_MAX, lag_mode="all").flatten()
  
  	#  Update progress bar every step
  	progress.update(int(100 * i / float(n_realizations)))
  
  #  Terminate progress bar
  progress.finish()
  
  #
  #  Plot
  #
  
  #  Set plotting parameters (for Clim. Past paper)
  params = { 'figure.figsize': (5.,4.),
            'axes.labelsize': 12,
            'text.fontsize': 12,
            'xtick.labelsize': 12,
            'ytick.labelsize': 12,
            'legend.fontsize': 10,
            'title.fontsize': 12,
            'text.usetex': False,
            'font': 'Helvetica',
            'mathtext.bf': 'helvetica:bold',
            'xtick.major.pad': 6,
            'ytick.major.pad': 6,
            'xtick.major.size': 5,
            'ytick.major.size': 5,
            'xtick.minor.size': 3,      # minor tick size in points
            'xtick.major.width': 1.,    # major tick width in points
            'xtick.minor.width': 1.,    # minor tick width in points
            'ytick.minor.size': 3,      # minor tick size in points
            'ytick.major.width': 1.,    # major tick width in points
            'ytick.minor.width': 1.,    # minor tick width in points
            'tick.labelsize': 'small'
             }
  
  #pylab.rcParams.update(params)
  
  pylab.figure()
  
  #for i in xrange(n_realizations):
  #	pylab.plot(np.arange(DELAY_MAX + 1) * average_sampling_time, acf[i, :], "k-")
  
  #  Plot min and mac auto-correlation for each lag from ensemble
  pylab.fill_between(np.arange(DELAY_MAX + 1) * average_sampling_time, acf.min(axis=0), acf.max(axis=0), color="0.75")
  
  #  Load Kernel ACF computed using Kira's method
  kernel_acf = pylab.loadtxt("kernel_acf_dongge.txt")
  
  pylab.plot(kernel_acf[:,0], kernel_acf[:,1], "k-")
  
  pylab.xlabel("Time lag (a)")
  pylab.ylabel("Autocorrelation")
  
  pylab.xlim(0, DELAY_TIME_MAX)
  pylab.ylim(-0.4, 1)
  
  #  Enable minor ticks
  pylab.minorticks_on()
  
  #  Set number of ticks and tick labels
  pylab.locator_params(axis="y", tight=True, nbins=5)
  
  pylab.axhline(y=0, color="k", linestyle="dashed")
  
  pylab.savefig("dongge_acf.pdf")
  
  pylab.show()
  