use fast cubic spline

This commit is contained in:
Richard Hartmann 2016-12-02 23:50:26 +01:00
parent 34778d76e6
commit 44516e9e3e
7 changed files with 65 additions and 32 deletions

View file

@ -1,8 +1,36 @@
from distutils.core import setup
from setuptools import setup
from Cython.Build import cythonize
import numpy as np
from stocproc import __version__
setup(
author = u"Richard Hartmann"
authors = [author]
description = 'Generate continuous time stationary stochastic processes from a given auto correlation function.'
name = 'stocproc'
version = __version__
if __name__ == "__main__":
setup(
name=name,
author=author,
author_email='richard.hartmann@tu-dresden.de',
url='https://github.com/cimatosa/stocproc',
version=version,
packages=[name],
package_dir={name: name},
license="BSD (3 clause)",
description=description,
long_description=description,
keywords=["stationary", "stochastic", "process", "random"],
classifiers= [
'Programming Language :: Python :: 3.4',
'Programming Language :: Python :: 3.5',
'License :: OSI Approved :: BSD License',
'Topic :: Utilities',
'Intended Audience :: Researcher'],
platforms=['ALL'],
install_requires=['fcSpline>=0.1'],
dependency_links=['https://raw.githubusercontent.com/cimatosa/fcSpline/master/egg/fcSpline-0.1-py3.4-linux-x86_64.egg'],
ext_modules = cythonize(["stocproc/stocproc_c.pyx"]),
include_dirs = [np.get_include()],
)
)

View file

@ -50,7 +50,7 @@ with the exact auto correlation.
.. image:: ../../examples/ac.*
"""
version = '0.2.0'
__version__ = '0.2.1'
import sys
if sys.version_info.major < 3:

View file

@ -6,8 +6,9 @@
"""
from __future__ import division, print_function
from .tools import ComplexInterpolatedUnivariateSpline
from functools import lru_cache
#from .tools import ComplexInterpolatedUnivariateSpline
#from functools import lru_cache
import fcSpline
import logging
import numpy as np
from numpy.fft import rfft as np_rfft
@ -261,9 +262,7 @@ def get_dt_for_accurate_interpolation(t_max, tol, ft_ref, diff_method=_absDiff):
while True:
tau = np.linspace(0, t_max, N+1)
ft_ref_n = ft_ref(tau)
tau_sub = tau[::sub_sampl]
ft_intp = ComplexInterpolatedUnivariateSpline(x = tau_sub, y = ft_ref_n[::sub_sampl], k=3)
ft_intp = fcSpline.FCS(x_low = 0, x_high=t_max, y=ft_ref_n[::sub_sampl])
ft_intp_n = ft_intp(tau)
d = diff_method(ft_intp_n, ft_ref_n)

View file

@ -6,6 +6,7 @@ import time
from . import stocproc_c
from . import gquad
from . import tools
import fcSpline
import logging
log = logging.getLogger(__name__)
@ -466,7 +467,11 @@ def auto_ng(corr, t_max, ngfac=2, meth=get_mid_point_weights_times, tol=1e-3, di
# setup cubic spline interpolator
t0 = time.time()
sqrt_lambda_ui_spl = tools.ComplexInterpolatedUnivariateSpline(tfine, sqrt_lambda_ui_fine)
#sqrt_lambda_ui_spl = tools.ComplexInterpolatedUnivariateSpline(tfine, sqrt_lambda_ui_fine, noWarning=True)
if not is_equi:
sqrt_lambda_ui_spl = tools.ComplexInterpolatedUnivariateSpline(tfine, sqrt_lambda_ui_fine, noWarning=True)
else:
sqrt_lambda_ui_spl = fcSpline.FCS(x_low=0, x_high=t_max, y=sqrt_lambda_ui_fine)
time_spline += (time.time() - t0)
# calculate the max deviation

View file

@ -15,11 +15,12 @@ stocproc_key_type = namedtuple(typename = 'stocproc_key_type',
class ComplexInterpolatedUnivariateSpline(object):
def __init__(self, x, y, k=3):
def __init__(self, x, y, k=3, noWarning=False):
if not noWarning:
raise DeprecationWarning("use fast cubic Spline (fcSpline) instead")
from scipy.interpolate import InterpolatedUnivariateSpline
self.re_spline = InterpolatedUnivariateSpline(x, np.real(y))
self.im_spline = InterpolatedUnivariateSpline(x, np.imag(y))
self.re_spline = InterpolatedUnivariateSpline(x, np.real(y), k=k)
self.im_spline = InterpolatedUnivariateSpline(x, np.imag(y), k=k)
def __call__(self, t):
return self.re_spline(t) + 1j * self.im_spline(t)

View file

@ -358,7 +358,7 @@ def test_calc_abN():
tau_fine = np.linspace(0, tmax, 1500)
ft_ref_n = bcf_ref(tau_fine)
ft_intp = sp.tools.ComplexInterpolatedUnivariateSpline(x=tau[idx], y=ft_tau[idx], k=3)
ft_intp = sp.tools.ComplexInterpolatedUnivariateSpline(x=tau[idx], y=ft_tau[idx], k=3, noWarning=True)
ft_intp_n = ft_intp(tau_fine)
d = diff_method(ft_intp_n, ft_ref_n)

View file

@ -198,20 +198,20 @@ def test_many(plot=False):
sd = spectral_density
ac = corr
stp = sp.StocProc_FFT(sd, t_max, ac, negative_frequencies=False, seed=0, intgr_tol=5e-3, intpl_tol=5e-3)
stocproc_metatest(stp, num_samples, tol, ac, plot)
stp = sp.StocProc_KLE(tol=5e-3, r_tau=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='full', meth='simp')
stocproc_metatest(stp, num_samples, tol, ac, plot)
stp = sp.StocProc_KLE(tol=5e-3, r_tau=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='random', meth='simp')
stocproc_metatest(stp, num_samples, tol, ac, plot)
stp = sp.StocProc_KLE(tol=5e-3, r_tau=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='full', meth='fp')
stocproc_metatest(stp, num_samples, tol, ac, plot)
stp = sp.StocProc_KLE(tol=5e-3, r_tau=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='random', meth='fp')
stocproc_metatest(stp, num_samples, tol, ac, plot)
# stp = sp.StocProc_FFT(sd, t_max, ac, negative_frequencies=False, seed=0, intgr_tol=5e-3, intpl_tol=5e-3)
# stocproc_metatest(stp, num_samples, tol, ac, plot)
#
# stp = sp.StocProc_KLE(tol=5e-3, r_tau=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='full', meth='simp')
# stocproc_metatest(stp, num_samples, tol, ac, plot)
#
# stp = sp.StocProc_KLE(tol=5e-3, r_tau=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='random', meth='simp')
# stocproc_metatest(stp, num_samples, tol, ac, plot)
#
# stp = sp.StocProc_KLE(tol=5e-3, r_tau=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='full', meth='fp')
# stocproc_metatest(stp, num_samples, tol, ac, plot)
#
# stp = sp.StocProc_KLE(tol=5e-3, r_tau=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='random', meth='fp')
# stocproc_metatest(stp, num_samples, tol, ac, plot)
t_max = 15