Fixed ComposedMaps performance issues

README.md

@@ -1,9 +1,9 @@
-# Poltergeist.jl 
+# Poltergeist.jl
 
 [![Build Status](https://travis-ci.org/wormell/Poltergeist.jl.svg?branch=master)](https://travis-ci.org/wormell/Poltergeist.jl)
 [![Coverage Status](https://coveralls.io/repos/github/wormell/Poltergeist.jl/badge.svg?branch=master)](https://coveralls.io/github/wormell/Poltergeist.jl?branch=master)
 
-Poltergeist is a package for quick, accurate and abstract approximation of statistical properties of one-dimensional chaotic dynamical systems. 
+Poltergeist is a package for quick, accurate and abstract approximation of statistical properties of one-dimensional chaotic dynamical systems.
 
 It treats chaotic systems through the framework of spectral methods (i.e. transfer operator-based approaches to dynamical systems) and is numerically implemented via spectral methods (e.g. Fourier and Chebyshev). For the latter, Poltergeist relies on and closely interfaces with the adaptive function approximation package [ApproxFun](https://github.com/ApproxFun/ApproxFun.jl).  
 
@@ -24,7 +24,7 @@ Similarly, take a circle map, or maps defined by modulo or inverse:
 
 ```julia
 c = CircleMap(x->4x + sin(2pi*x)/2pi,PeriodicInterval(0,1))
-lanford = modulomap(x->2x+x*(1-x)/2,0..1)
+lanford = modulomap(x->2x+x*(1-x)/2,0..1) # Or call lanford()
 doubling = MarkovMap([x->x/2,x->(x+1)/2],[0..0.5,0.5..1],dir=Reverse)
 ```
 
@@ -40,9 +40,11 @@ L = Transfer(M)
 f0 = Fun(x->sin(3pi*x),d) #ApproxFun function
 f1 = L*f0
 g = ((2I-L)\f0)'
-eigvals(L,30)
 det(I-4L) # Fredholm determinant
-``` 
+
+using Plots
+scatter(eigvals(L,80))
+```
 
 In particular, you can solve for many statistical properties, many of which Poltergeist has built-in commands for. Most of these commands allow you to use the ```MarkovMap``` directly (bad, zero caching between uses), transfer operator (caches transfer operator entries, usually the slowest step), or the ```SolutionInv``` operator (caches QR factorisation as well).
 
@@ -54,12 +56,10 @@ birkhoffvar(K,Fun(x->x^2,d))
 birkhoffcov(K,Fun(x->x^2,d),Fun(identity,d))
 dρ = linearresponse(K,Fun(sinpi,d))
 
-using Plots
 plot(ρ)
 ε = 0.05
 plot!(ρ + ε*dρ,title="Linear response")
-pertε(func) = x-> func(x) + ε*sinpi(func(x))
-plot!(acim(MarkovMap([pertε(f1),pertε(f2)],[0..0.5,0.5..1])))
+plot!(acim(perturb(M,sinpi,ϵ)))
 ```
 <!--- TODO: plot!(linearresponse(L,Fun(x->x*(1-x),d))) --->
 <img src=https://github.com/johnwormell/Poltergeist.jl/raw/master/images/acim.png width=500>
@@ -74,10 +74,10 @@ This package is based on academic work. If you find this package useful in your
 
 
 <!---
-* J. P. Wormell (2017 in preparation), Fast numerical methods for intermittent systems 
+* J. P. Wormell (2017 in preparation), Fast numerical methods for intermittent systems
 
 
 _____________
 * A map f on [-1,1] is C-expanding if acos◦f◦cos is uniformly expanding. Most uniformly expanding maps are C-expanding, or if not there is always a conjugate or iterate that is.
 
---->
+--->

images/images.jl

@@ -2,12 +2,20 @@ using Poltergeist, Plots, ApproxFun; pyplot()
 d = Segment(0,1);
 M1(x) = 2x+sin(2pi*x)/6; M2(x) = 2-2x
 M = MarkovMap([M1,M2],[0..0.5,0.5..1]);
-K = SolutionInv(M);
+L = Transfer(M)
+
+scatter(eigvals(L,80),legend=true,grid=false,legend=false,
+title = "Eigenvalues of transfer operator")
+xlims!(-1.1,1.1)
+ylims!(-0.5,0.5)
+savefig("eigvals.pdf")
+run(`sips -s format png eigvals.pdf --out eigvals.png`)
+
+K = SolutionInv(L);
 ρ = acim(K);
 dρ = linearresponse(K,sinpi)
 ϵ = 0.05
-ϵpert(x) = x + ϵ*sinpi(x)
-Mϵ = MarkovMap([ϵpert∘M1,ϵpert∘M2],[0..0.5,0.5..1])
+Mϵ = perturb(M,sinpi,ϵ)
 plot(ρ,legend=true,grid=false,label="\$\\rho\_\{0\}\$",
 title = "Linear response");
 plot!(ρ+ϵ*dρ,label="\$\\rho\_\{\\epsilon\}\$ estimate")
@@ -16,4 +24,4 @@ xlabel!("\$x\$"); ylabel!("\$\\rho(x)\$")
 # xlabel!("x"); ylabel!("ρ(x)")
 ylims!(0.,1.8)
 savefig("acim.pdf")
-println("First image done")
+run(`sips -s format png acim.pdf --out acim.png`)

src/ComposedMaps.jl

@@ -1,55 +1,52 @@
-struct ComposedMarkovMap{M<:AbstractMarkovMap,D<:Domain,R<:Domain,N} <: AbstractMarkovMap{D,R}
-  maps::NTuple{N,M}
+struct ComposedMarkovMap{T<:Tuple,D<:Domain,R<:Domain} <: AbstractMarkovMap{D,R}
+  maps::T
   domain::D
   rangedomain::R
 end
 (∘)(f::AbstractMarkovMap,g::AbstractMarkovMap) = ComposedMarkovMap(f,g)
 
 function ComposedMarkovMap(maps...)
-  N = length(maps)
   ran = rangedomain(maps[1])
   mps = isa(maps[1],ComposedMarkovMap) ? maps[1].maps : (maps[1],)
-  for i = 1:N-1
+  for i = 1:length(maps)-1
     @assert rangedomain(maps[i+1]) == domain(maps[i])
     mps = isa(maps[i+1],ComposedMarkovMap) ? (mps...,maps[i+1].maps...) : (mps...,maps[i+1])
   end
   dom = domain(maps[end])
 
-  Nf = length(mps)
-  T = typejoin(map(typeof,mps)...)
-  ComposedMarkovMap{T,typeof(dom),typeof(ran),Nf}(convert(NTuple{Nf,T},mps),dom,ran)
+  ComposedMarkovMap(mps,dom,ran)
 end
 
-complength{M,D,R,N}(::ComposedMarkovMap{M,D,R,N}) = N
+complength(c::ComposedMarkovMap) = length(c.maps)
 
-function (c::ComposedMarkovMap{M,D,R,N})(x) where {M,D,R,N}
+function (c::ComposedMarkovMap)(x)
   f = c.maps[end](x)
-  for i = N-1:-1:1
+  for i = complength(c)-1:-1:1
     f = c.maps[i](f)
   end
   f
 end
 
-function mapP(c::ComposedMarkovMap{M,D,R,N},x) where {M,D,R,N}
+function mapP(c::ComposedMarkovMap,x)
   f,dfdx = mapP(c.maps[end],x)
-  for i = N-1:-1:1
+  for i = complength(c)-1:-1:1
     f,dfdxs = mapP(c.maps[i],f)
     dfdx *= dfdxs
   end
   f,dfdx
 end
 mapD(c::ComposedMarkovMap,x) = mapP(c,x)[2]
 
-function mapinv(c::ComposedMarkovMap{M,D,R,N},b,x) where {M,D,R,N}
+function mapinv(c::ComposedMarkovMap,b,x)
   v = mapinv(c.maps[1],b[1],x)
-  for i = 2:N
+  for i = 2:complength(c)
     v = mapinv(c.maps[i],b[i],v)
   end
   v
 end
-function mapinvP(c::ComposedMarkovMap{M,D,R,N},b,x) where {M,D,R,N}
+function mapinvP(c::ComposedMarkovMap,b,x)
   v,dvdx = mapinvP(c.maps[1],b[1],x)
-  for i = 2:N
+  for i = 2:complength(c)
     v,dvdxs = mapinvP(c.maps[i],b[i],v)
     dvdx *= dvdxs
   end
@@ -59,11 +56,11 @@ mapinvD(c::ComposedMarkovMap,b,x) = mapinvP(c,b,x)[2]
 
 # TODO: map(P,D)(c,b,x)
 
-function getbranch(m::ComposedMarkovMap{M,D,R,N},x) where {M,D,R,N}
+function getbranch(m::ComposedMarkovMap,x)
   temp_in(x,m.domain) || error("DomainError: $x ∉ $(m.domain)")
   fx = x
-  br = getbranch(m.maps[N],fx)
-  for i = N-1:-1:1
+  br = getbranch(m.maps[end],fx)
+  for i = complength(c)-1:-1:1
     fx = m.maps[i+1](fx)
     br = (getbranch(m.maps[i],fx),br...)
   end
@@ -74,14 +71,24 @@ nbranches(C::ComposedMarkovMap) = prod(nbranches(mm for mm in C.maps))
 eachbranchindex(C::ComposedMarkovMap) = product(eachbranchindex(mm) for mm in C.maps)
 
 #TODO: must be faster??
-function transferfunction{M,D,R,N}(x,m::ComposedMarkovMap{M,D,R,N},f,T)
-  m2 = ComposedMarkovMap(m.maps[2:end]...)
-  transferfunction(x,m.maps[1],x->transferfunction(x,m2,f,T),T)
-end
-function transferfunction{M,D,R}(x,m::ComposedMarkovMap{M,D,R,1},f,T)
+function transferfunction{M<:AbstractMarkovMap}(x,m::ComposedMarkovMap{Tuple{M}},f,T)
   transferfunction(x,m.maps[1],f,T)
 end
 
+struct TransferCall{M<:AbstractMarkovMap,ff,T}
+  m::M
+  f::ff
+  t::Type{T}
+end
+# TransferCall(m,f,T) = TransferCall{typeof(m),typeof(f),T}(m,f)
+(t::TransferCall)(x) = transferfunction(x,t.m,t.f,t.t)
+
+function transferfunction(x,m::ComposedMarkovMap,f,T)
+  m2 = ComposedMarkovMap(m.maps[2:end],domain(m.maps[end]),rangedomain(m.maps[2]))
+  transferfunction(x,m.maps[1],#x->transferfunction(x,m2,f,T),T)
+    TransferCall(m2,f,T),T)
+end
+
 # TODO: transferfunction_int
 
 

src/Poltergeist.jl

@@ -6,7 +6,7 @@ import Base: values,getindex,setindex!,*,.*,+,.+,-,.-,==,<,<=,>,|,
 import ApproxFun: domainspace, rangespace, domain, israggedbelow, RaggedMatrix, resizedata!, colstop, CachedOperator, Infinity,
 fromcanonicalD, tocanonicalD, fromcanonical, tocanonical, space
 
-export (..)
+export (..), Segment, PeriodicInterval, Fun
 
 ## TEMPORARY PENDING APPROXFUN UPDATE
 temp_in(x,dom) = (x >= dom.a) && (x <= dom.b)

src/Transfer.jl

@@ -91,14 +91,14 @@ function transfer_getindex{T}(L::ConcreteTransfer{T},jdat::Tuple{Integer,Integer
     kind > 1 && (mc = max(mc,maximum(L.colstops[k[kind-1]+1:kk])))
 
     #    f(x) = transferfunction(x,L,kk,T)
-    tol =T==Any?20eps():20eps(T)
+    tol =T==Any?200eps():200eps(T)
 
     if L.colstops[kk] >= 1
       coeffs = ApproxFun.transform(rs,transferfunction_nodes(L,max(16,nextpow2(L.colstops[kk])),kk,T))[1:L.colstops[kk]]
       maxabsc = max(maximum(abs.(coeffs)),one(T))
-      chop!(coeffs,tol*maxabsc*log2(length(coeffs))/10)
+      chop!(coeffs,tol*maxabsc*log2(length(coeffs)))
     elseif L.colstops[kk]  == 0
-      coeffs = [0.]
+      coeffs = zeros(T,1)
     else
 
       #       if mc ≤ 2^4
@@ -118,7 +118,7 @@ function transfer_getindex{T}(L::ConcreteTransfer{T},jdat::Tuple{Integer,Integer
 
         maxabsc = max(one(T),maximum(abs.(coeffs)))
         if maxabsc == 0 && maxabsfr == 0
-          coeffs = [0.]
+          coeffs = zeros(T,1)
           break
         else
           maxabsfr = max(maxabsfr,one(T))

src/poetry.jl

@@ -34,6 +34,9 @@ perturb(d::IntervalDomain,X,ϵ) = MarkovMap([x->x+ϵ*X(x)],[d],d)
 perturb(d::PeriodicDomain,X,ϵ) = FwdCircleMap([x->x+ϵ*X(x)],d)
 perturb(m::AbstractMarkovMap,X,ϵ) = perturb(rangedomain(m),X,ϵ)∘m
 
+# Eigvals overloads
+eigvals(m::AbstractMarkovMap,n) = eigvals(Transfer(m),n)
+
 # # plotting
 # function plot(m::MarkovMap)
 #   pts = eltype(m)[]

src/statistics.jl

@@ -139,8 +139,6 @@ function covariancefunction(L::Operator,A::Fun,B::Fun;r=acim(L),tol=eps(norm(coe
 end
 
 # Lyapunov exponent
-# IMPORTANT TODO: check that this works in general cause the function to be transferred
-# isn't necessarily continuous
 function lyapunov(f::AbstractMarkovMap,r=acim(f),sp=Space(rangedomain(f)))
   sum(Fun(x->transferfunction(x,f,x->log(abs(f'(x)))*r(x),eltype(f)),sp))
 end
@@ -150,7 +148,22 @@ function lyapunov(K::SolutionInvWrapper)
 end
 lyapunov(m::Operator) = lyapunov(markovmap(m),acim(m),rangespace(m))
 
-#TODO: lyapunov exponent of ComposedMarkovMap
+struct LyapContainer{M<:AbstractMarkovMap,rr}
+  m::M
+  tr::rr
+end
+(lc::LyapContainer)(x) = log(abs(lc.m'(x)))*lc.tr(x)
+
+function lyapunov(f::ComposedMarkovMap,r=acim(f))#,sp=Space(rangedomain(f)))
+  T = eltype(f)
+  tr = copy(r)
+  lyap = sum(Fun(x->transferfunction(x,f.maps[end],LyapContainer(f.maps[end],tr),T),rangedomain(f.maps[end])))
+  for k = complength(f)-1:-1:1
+    tr = Fun(x->transferfunction(x,f.maps[k+1],tr,T),rangedomain(f.maps[k+1]))
+    lyap += sum(Fun(x->transferfunction(x,f.maps[k],LyapContainer(f.maps[k],tr),T),rangedomain(f.maps[k])))
+  end
+  lyap
+end
 
 # Gaussian process parametrisation
 function logMA_process(L,A::Fun)

test/runtests.jl

@@ -47,9 +47,9 @@ acim(M2f)
 println("Should be ≤0.12s")
 
 pts = [points(space(ρ1b),100);points(space(ρ2b),100)]
-@test maximum(abs.(ρ1f.(pts) - ρ2f.(pts))) < 1000eps(1.)
-@test maximum(abs.(ρ1b.(pts) - ρ2b.(pts))) < 1000eps(1.)
-@test maximum(abs.(ρ2b.(pts) - ρ2ba.(pts))) < 1000eps(1.)
+@test maximum(abs.(ρ1f.(pts) - ρ2f.(pts))) < 2000eps(1.)
+@test maximum(abs.(ρ1b.(pts) - ρ2b.(pts))) < 2000eps(1.)
+@test maximum(abs.(ρ2b.(pts) - ρ2ba.(pts))) < 2000eps(1.)
 
 # # Transfer
 # @test transfer(M1f,x->Fun(Fourier(d1),[0.,1.])(x),0.28531) == Poltergeist.transferfunction(0.28531,M1f,Poltergeist.BasisFun(Fourier(d1),2),Float64)
@@ -76,16 +76,21 @@ K = SolutionInv(lan);
 println("Composition test 🎼")
 shiftmap = modulomap(x->5x+30,0..1,30..35)
 lanshift = modulomap(x->5(x/5-6)/2-(x/5-6)^2/2,30..35,0..1)
-doublelan = lan ∘ lanshift ∘ shiftmap
+@time doublelan = lan ∘ lanshift ∘ shiftmap
+println("Should be ≤0.07s")
 doubleK = SolutionInv(doublelan)
 doublerho = acim(doubleK)
-#TODO: why is this so slow:
-@time doublerho = acim(lan)
+@time doublerho = acim(doublelan)
+println("Should be ≤0.4s")
 @test doublerho ≈ rho
 @test lyapunov(doubleK) ≈ 2l_exp
 
-@test lyapunov(perturb(lan,sinpi,-0.1)∘inv(perturb(0..1,sinpi,-0.1))) ≈ l_exp
-@time acim(perturb(lan,sinpi,-0.1)∘inv(perturb(0..1,sinpi,-0.1)))
+lanpet = perturb(lan,sinpi,-0.1)∘inv(perturb(0..1,sinpi,-0.1))
+@test lyapunov(lanpet) ≈ l_exp
+@time lyapunov(lanpet)
+println("Should be ≤0.01s")
+# @time lanpet = acim(lanpet)
+# @time lyapunov(lanpet)
 
 # Correlation sums
 println("Correlation sum test")