# Copyright (c) 2016-2018, The University of Texas at Austin
# & University of California--Merced.
# Copyright (c) 2019-2020, The University of Texas at Austin
# University of California--Merced, Washington University in St. Louis.
#
# All Rights reserved.
# See file COPYRIGHT for details.
#
# This file is part of the hIPPYlib library. For more information and source code
# availability see https://hippylib.github.io.
#
# hIPPYlib is free software; you can redistribute it and/or modify it under the
# terms of the GNU General Public License (as published by the Free
# Software Foundation) version 2.0 dated June 1991.
import numpy as np
[docs]def _acorr(mean_free_samples, lag, norm = 1):
#http://stackoverflow.com/questions/14297012/estimate-autocorrelation-using-python
return (mean_free_samples[:mean_free_samples.size-lag]*mean_free_samples[lag:]).ravel().mean() / norm
[docs]def _acorr_vs_lag(samples, max_lag):
mean = samples.mean()
mean_free_samples = samples - mean
norm = _acorr(mean_free_samples, 0)
lags = np.arange(0,max_lag+1)
acorrs = np.ones(max_lag+1)
for lag in lags[1:]:
acorrs[lag] = _acorr(mean_free_samples, lag, norm)
return lags, acorrs
[docs]def integratedAutocorrelationTime(samples, max_lag = None):
assert len(samples.shape) == 1
if max_lag is None:
max_lag = samples.shape[0] // 10
lags, acorrs = _acorr_vs_lag(samples, max_lag)
iact = 1. + 2.* np.max( acorrs.cumsum() )
return iact, lags, acorrs