Compatibility with ApproxFun v0.10

.travis.yml

@@ -5,6 +5,7 @@ os:
   - osx
 julia:
   - 0.7
+  - 1.0
 notifications:
   email: false
 matrix:

README.md

@@ -13,7 +13,7 @@ As an example, take your favourite Markov interval map and give it digital form:
 
 ```julia
 using Poltergeist, ApproxFun
-d = Segment(0,1)
+d = 0..1.
 f1 = 2x+sin(2pi*x)/6; f2(x) = 2-2x
 f = MarkovMap([f1,f2],[0..0.5,0.5..1])
 f(0.25), f'(0.25)
@@ -23,7 +23,7 @@ f(0.25), f'(0.25)
 Similarly, take a circle map, or maps defined by modulo or inverse:
 
 ```julia
-c = CircleMap(x->4x + sin(2pi*x)/2pi,PeriodicInterval(0,1))
+c = CircleMap(x->4x + sin(2pi*x)/2pi,PeriodicSegment(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) # or doubling(0..1)
 ```

REQUIRE

@@ -1,8 +1,8 @@
 julia 0.7
-ApproxFun 0.7
+ApproxFun 0.10
 BandedMatrices 0.2
 Compat 0.31.0
+DomainSets
 DualNumbers 0.3.0
 StaticArrays 0.3
-IntervalSets 0.0.2
 LightGraphs 0.12

images/images.jl

@@ -1,5 +1,5 @@
 using Poltergeist, Plots, ApproxFun; pyplot()
-d = Segment(0,1);
+d = Interval(0,1);
 M1(x) = 2x+sin(2pi*x)/6; M2(x) = 2-2x
 M = MarkovMap([M1,M2],[0..0.5,0.5..1]);
 L = Transfer(M)

src/HGraphs.jl

@@ -59,7 +59,7 @@ eltype(x::HofbauerExtension) = Int # for graph compatibiltiy
 zero(h::HofbauerExtension{I}) where I = HofbauerExtension(0,Vector{HEdge{I}}[],Vector{HEdge{I}}[],h.m,h.hdomains,h.returnto)
 
 HofbauerExtension{I,D}(m::AbstractIntervalMap) where {I,D} = HofbauerExtension(0, Vector{HEdge{I}}[], Vector{HEdge{I}}[],m,D[],Int[])
-HofbauerExtension(m::AbstractIntervalMap) = HofbauerExtension{Int,Segment{Float64}}(m)
+HofbauerExtension(m::AbstractIntervalMap) = HofbauerExtension{Int,Interval{Float64}}(m)
 HofbauerExtension(::Type{T},::Type{D},m::AbstractIntervalMap) where {T,D} = HofbauerExtension{T,D}(m)
 
 edgetype(::HofbauerExtension{T}) where T<: Integer = HEdge{T}

src/Poltergeist.jl

@@ -2,20 +2,18 @@ __precompile__()
 
 module Poltergeist
 
-using Base, ApproxFun, BandedMatrices, Compat, DualNumbers, StaticArrays, IntervalSets, LightGraphs#, PyPlot#, FastTransforms
+using Base, ApproxFun, BandedMatrices, Compat, DualNumbers, DomainSets, StaticArrays, LightGraphs#, PyPlot#, FastTransforms
 using LinearAlgebra
 
 import Base: values,getindex,setindex!,*,.*,+,.+,-,.-,==,<,<=,>,|,
   >=,./,/,.^,^,\,∪,∘,transpose, size, length,  eltype, inv, mod, convert#issymmetric,
 @compat import LinearAlgebra: eigvals, eigvecs
 import ApproxFun: domainspace, rangespace, domain, israggedbelow, RaggedMatrix,
-  resizedata!, colstop, CachedOperator, Infinity, IntervalDomain, UnionDomain,
+  resizedata!, colstop, CachedOperator, Infinity,
   fromcanonicalD, tocanonicalD, fromcanonical, tocanonical, space, eigs, cfstype, prectype
+import DomainSets: AbstractInterval, UnionDomain
 
-export (..), Segment, PeriodicInterval, Fun, eigs
-
-## TEMPORARY PENDING APPROXFUN UPDATE
-temp_in(x,dom) = (x >= dom.a) && (x <= dom.b)
+export (..), Interval, PeriodicSegment, Fun, eigs
 
 include("general.jl")
 include("maps/maps.jl")

src/general.jl

@@ -13,7 +13,7 @@ end
 
 @compat const PureLaurent{D<:Domain} = Hardy{false,D}
 
-@compat const GeneralInterval = Union{Segment,PeriodicInterval}
+@compat const GeneralInterval = Union{Segment,Interval,PeriodicSegment}
 
 # rem()
 
@@ -45,18 +45,29 @@ ApproxFun.prectype(p::InterpolationNode) = prectype(p.sp)
 # Circle domains
 
 interval_mod(x,a,b) = (b-a)*mod((x-a)/(b-a),1)+a
-domain_mod(x,p::PeriodicInterval) = interval_mod(x,p.a,p.b)
+domain_mod(x,p::PeriodicSegment) = interval_mod(x,leftendpoint(p),rightendpoint(p))
 
 # Interval domains
 
-coveringsegment(dsm::AbstractArray{T}) where {T<:Domain} = Segment(minimum(first(convert(Domain,d)) for d in dsm),maximum(last(convert(Domain,d)) for d in dsm))
-coveringsegment(ds::AbstractArray) = coveringsegment([convert(Domain,d) for d in ds])
-mapinterval(f,d::Domain) = Interval(extrema((f(first(d)),f(last(d))))...)
+function approx_in(x,d::GeneralInterval; atol=50eps(arclength(d)), kwargs...)
+    x <= rightendpoint(d)+atol && x >= leftendpoint(d) - atol
+end
+approx_in(x,d) = approx_in(convert(typeof(d),x),d)
+
+function approx_issubset(a,b;kwargs...)
+    acupb = a∪b
+    atol_isapprox(acupb,b;kwargs...)
+end
+
+atol_isapprox(a,b; atol=50eps(max(arclength(a),arclength(b))), kwargs...) = isapprox(a,b;atol=atol,kwargs...)
+coveringinterval(dsm::AbstractArray{T}) where {T<:AbstractInterval} = Interval(minimum(leftendpoint(convert(Domain,d)) for d in dsm),maximum(rightendpoint(convert(Domain,d)) for d in dsm))
+coveringinterval(ds::AbstractArray) = coveringinterval([convert(Domain,d) for d in ds])
+mapinterval(f,d::T) where T<:GeneralInterval = T(extrema((f(leftendpoint(d)),f(rightendpoint(d))))...)
 mapinterval(f,d) = mapinterval(f,convert(Domain,d))
 
 # Newton's method
 
-domain_newton(f,df,y::U,D::Domain,x::U=convert(U,rand(D)),tol=400eps(eltype(U))) where U = basic_newton(f,df,y,x,tol)
+domain_newton(f,df,y::U,D::Domain,x::U=convert(U,ApproxFun.checkpoints(D)[1]),tol=400eps(eltype(U))) where U = basic_newton(f,df,y,x,tol)
 domain_guess(x,dom::Domain,ran::Domain) = rand(ran)
 
 domain_guess(x::SVector{N},dom::ProductDomain,ran::ProductDomain) where N = SVector{N}([interval_guess(x[i],dom.domains[i],ran.domains[i]) for i =1:N])
@@ -131,14 +142,14 @@ function interval_newton(f,df,y::U,da::T,db::T,x::U=(da+db)/2,tol=40eps(max(abs(
   abs(rem) > tol && error("Newton: failure to converge: y = $y, x estimate = $x, rem = $rem")
   x
 end
-function domain_newton(f,df,y::U,D::GeneralInterval,x::U=(D.a+D.b)/2,
-            tol=40eps(max(abs(D.a),abs(D.b)))) where U<:Union{Real,Complex}
+function domain_newton(f,df,y::U,D::GeneralInterval,x::U=(leftendpoint(D)+rightendpoint(D))/2,
+            tol=40eps(max(leftendpoint(D),rightendpoint(D)))) where U<:Union{Real,Complex}
     # x = convert(U,x);
-    interval_newton(f,df,y,D.a,D.b,x,tol)
+    interval_newton(f,df,y,leftendpoint(D),rightendpoint(D),x,tol)
 end
 
 interval_guess(y::Number,dom::Domain,ran::Domain) =
-    (dom.a*ran.b-ran.a*dom.b+y*(dom.b-dom.a))/(ran.b-ran.a)
+    (leftendpoint(dom)*rightendpoint(ran)-leftendpoint(ran)*rightendpoint(dom)+y*(rightendpoint(dom)-leftendpoint(dom)))/(rightendpoint(ran)-leftendpoint(ran))
 domain_guess(y::Number,dom::GeneralInterval,ran::GeneralInterval) =
     interval_guess(y,dom,ran)
 

src/hofbauerextension.jl

@@ -109,9 +109,9 @@ The keyword `maxdepth` says how deep the Hofbauer extension should go, and
 `forcereturn` sets whether a return to given members of `basedomains` should
  forced if possible.
 
- For example, in the case of `f(x) = mod(2x,1)`, if `basedomains` is set to `Segment(0.,0.3)` then
-`Segment(0.,0.5)` might map to `Segment(0.,1.)`` for a given map if `forcereturn=false` but would map to
-Segment(0,0.3)` and `Segment(0.3,1)` if `forcereturn=true`.
+ For example, in the case of `f(x) = mod(2x,1)`, if `basedomains` is set to `Interval(0.,0.3)` then
+`Interval(0.,0.5)` might map to `Interval(0.,1.)`` for a given map if `forcereturn=false` but would map to
+Interval(0,0.3)` and `Interval(0.3,1)` if `forcereturn=true`.
 """
 function hofbauerextension(m::AbstractIntervalMap,basedoms=(maxrangedomain(m),);maxdepth=100,forcereturn=trues(length(basedoms)))
     @assert domain(m) == rangedomain(m)
@@ -133,7 +133,7 @@ function hofbauerextension(m::AbstractIntervalMap,basedoms=(maxrangedomain(m),);
     #     convert(Vector{HofbauerDomain{D}},hdomains(G)))
     G
 end
-hofbauerextension(m::AbstractIntervalMap,basedoms::Union{Domain,IntervalSets.ClosedInterval};maxdepth=100,forcereturn=trues(length(basedoms))) =
+hofbauerextension(m::AbstractIntervalMap,basedoms::Domain;maxdepth=100,forcereturn=trues(length(basedoms))) =
     hofbauerextension(m,(basedoms,);maxdepth=maxdepth,forcereturn=forcereturn)
 
 # walks through domains in HofbauerExtension
@@ -143,30 +143,30 @@ function towerwalk!(G::HofbauerExtension, graphind, hdom)
     for branchind in eachbranchindex(m)
         b = branches(m)[branchind]
         bdom = domain(b)
-        capdom = bdom ∩ hdom.domain
-        if ~isempty(capdom) # hdom is partially contained in the domain of b
+        capdom = bdom ∩ convert(Domain,hdom)
+        if ~isempty(capdom) && arclength(capdom) > 0 # hdom is partially contained in the domain of b
             if isa(b,NeutralBranch) #&& (neutralfixedpoint(b) ∈ capdom) # inducing at a neutral fixed point
                 isnfp = true
-                newdom = rangedomain(b) \ ((bdom ≈ capdom) ? domain(b) : b(capdom))
+                newdom = rangedomain(b) \ (atol_isapprox(bdom,capdom) ? domain(b) : b(capdom))
             else # normal bounce inducing
                 isnfp = false
-                newdom = (bdom ≈ capdom) ? rangedomain(b) : b(capdom)
+                newdom = atol_isapprox(bdom,capdom) ? rangedomain(b) : b(capdom)
             end
             for i in G.returnto
                 forcedom = convert(Domain,hdomains(G)[i])
-                if issubset(forcedom,newdom)
+                if approx_issubset(forcedom,newdom)
                     towerextend!(G, graphind, hdom, forcedom, branchind, isnfp)
                     newdom = newdom \ forcedom
                 end
             end
-            towerextend!(G, graphind, hdom, newdom, branchind, isnfp)
+            arclength(newdom) > 0 && towerextend!(G, graphind, hdom, newdom, branchind, isnfp)
         end
     end
 end
 
 function towerextend!(G::HofbauerExtension, graphind, hdom, newdom, branchind, isnfp)
     arclength(newdom) < 300eps(arclength(convert(Domain,hdom))) && return true
-    newgraphind = findfirst(collect(h.domain≈newdom for h in hdomains(G)))
+    newgraphind = findfirst(collect(atol_isapprox(h.domain,newdom) for h in hdomains(G)))
 
     # add new domain if it doesn't already exist
     if newgraphind isa Nothing

src/maps/CircleMap.jl

@@ -13,7 +13,7 @@
 end
 function FwdCircleMap(f::ff,domd::D,randm::R,dfdx::gg=autodiff(f,domd)) where {D<:PeriodicDomain,R<:PeriodicDomain,ff,gg}
   # @assert isempty(∂(randm)) isempty(∂(domd))
-  fa = f(first(domd)); fb = f(last(domd))
+  fa = f(leftendpoint(domd)); fb = f(rightendpoint(domd))
   cover_est = (fb-fa)/arclength(randm)
   cover_integer = round(Int,cover_est)
   cover_est ≈ cover_integer || error("Circle map lift does not have integer covering number.")
@@ -49,7 +49,7 @@ end
 
 function RevCircleMap(v::ff,domd::D,randm::R,dvdx::gg=autodiff(v,randm), maxcover=10000) where {D<:Domain,R<:Domain,ff,gg}
   # @assert isempty(∂(ran)) isempty(∂(dom))
-  ra = first(randm); va = v(ra)
+  ra = leftendpoint(randm); va = v(ra)
   dr = arclength(randm); dd = arclength(domd)
   cover = 1; vr = v(ra+dr)-va
   while (abs(vr) <= dd && cover < maxcover)
@@ -60,7 +60,7 @@ function RevCircleMap(v::ff,domd::D,randm::R,dvdx::gg=autodiff(v,randm), maxcove
   end
   cover == maxcover && error("Can't get to the end of the inverse lift after $maxcover steps")
 
-  RevCircleMap{D,R,ff,gg,typeof(va)}(v,dvdx,domd,randm,cover,va,v(last(randm)))
+  RevCircleMap{D,R,ff,gg,typeof(va)}(v,dvdx,domd,randm,cover,va,v(rightendpoint(randm)))
 end
 function RevCircleMap(v::ff,dom,ran=dom,dvdx::gg=autodiff(v,ran)) where {ff,gg}
   domd = convert(PeriodicDomain,dom); randm = convert(PeriodicDomain,ran)

src/maps/ComposedMaps.jl

@@ -58,7 +58,7 @@ mapinvD(c::ComposedMarkovMap,b,x) = mapinvP(c,b,x)[2]
 # TODO: map(P,D)(c,b,x)
 
 function getbranchind(m::ComposedMarkovMap,x)
-  temp_in(x,m.domain) || error("DomainError: $x ∉ $(m.domain)")
+  x ∈ m.domain || error("DomainError: $x ∉ $(m.domain)")
   fx = x
   br = getbranchind(m.maps[end],fx)
   for i = complength(c)-1:-1:1

src/maps/ExpandingBranch.jl

@@ -67,15 +67,15 @@ unsafe_mapinvP(b::RevExpandingBranch,x) = (unsafe_mapinv(b,x),unsafe_mapinvD(b,x
 
 for B in (FwdExpandingBranch,RevExpandingBranch)
   for (M,UNS_M) in ((:mapinv,:unsafe_mapinv),(:mapinvD,:unsafe_mapinvD),(:mapinvP,:unsafe_mapinvP))
-    @eval $M(b::$B,x) = begin @assert in(x,rangedomain(b)); $UNS_M(b,x); end
+    @eval $M(b::$B,x) = begin @assert approx_in(x,rangedomain(b)); $UNS_M(b,x); end # TODO: sort out approx_in stuff
   end
 
-  # we are assuming that domains are unoriented
-  @eval (b::$B)(x::Segment) = Segment(sort(b.(∂(x)))...)
+  @eval (b::$B)(x::T) where T<:AbstractInterval = #begin println(leftendpoint(x)); println(rightendpoint(x));
+    T(extrema((b(leftendpoint(x)),b(rightendpoint(x))))...) # end
   #@eval (b::$B)(x::Interval) = Interval(sort(b.(extrema(x)))...)
 end
 
-mapinv(b::ExpandingBranch,x::Segment) = Segment((mapinv.(b,∂(x)))...)
+mapinv(b::ExpandingBranch,x::T) where T<:AbstractInterval = T(extrema((mapinv(b,leftendpoint(x)),mapinv(b,rightendpoint(x))))...)
 # mapinv(b::ExpandingBranch,x::Interval) = Interval(sort(mapinv.(b,extrema(x)))...)
 
 
@@ -101,8 +101,6 @@ function autodiff_dual(f,bi)
   fd
 end
 
-@compat const DomainInput = Union{Domain,IntervalSets.AbstractInterval}
-
 function branch(f,dom,ran,diff=autodiff(f,(dir==Forward ? dom : ran));dir=Forward,
                       ftype=typeof(f),difftype=typeof(diff))
   domd = convert(Domain,dom); randm  = convert(Domain,ran);
@@ -115,9 +113,9 @@ branch(f,dom,ran,diff::Nothing; dir =Forward) = branch(f,dom,ran;dir=dir)
 
 
 for TYP in (:FwdExpandingBranch,:RevExpandingBranch,:NeutralBranch)
-  @eval (b::$TYP)(x) = temp_in(x,b.domain) ? unsafe_call(b,x) : error("DomainError: $x ∉ $(b.domain)")
-  @eval mapD(b::$TYP,x) = temp_in(x,b.domain) ? unsafe_mapD(b,x) : error("DomainError: $x ∉ $(b.domain)")
-  @eval mapP(b::$TYP,x) = temp_in(x,b.domain) ? unsafe_mapP(b,x) : error("DomainError: $x ∉ $(b.domain)")
+  @eval (b::$TYP)(x) = x ∈ b.domain ? unsafe_call(b,x) : error("DomainError: $x ∉ $(b.domain)")
+  @eval mapD(b::$TYP,x) = x ∈ b.domain ? unsafe_mapD(b,x) : error("DomainError: $x ∉ $(b.domain)")
+  @eval mapP(b::$TYP,x) = x ∈ b.domain ? unsafe_mapP(b,x) : error("DomainError: $x ∉ $(b.domain)")
   @eval domain(b::$TYP) = b.domain
   @eval rangedomain(b::$TYP) = b.rangedomain
 end
@@ -126,8 +124,8 @@ end
 
 function transferbranch_int_edges(x,y,b::ExpandingBranch)
   x ∈ rangedomain(b) && y ∈ rangedomain(b) && (return x,y)
-  iv = Segment(x,y)∩rangedomain(b)
-  iv.a, iv.b
+  iv = Interval(x,y)∩rangedomain(b)
+  leftendpoint(iv), rightendpoint(iv)
 end
 
 function transferbranch(x,b::ExpandingBranch,f)

src/maps/IntervalMap.jl

@@ -54,20 +54,20 @@ function MarkovMap(branches::AbstractVector{B},dom,ran) where {B<:ExpandingBranc
 end
 
 """
-    MarkovMap(fs::Vector, ds::Vector, ran = coveringsegment(ds); dir=Forward, diff=autodiff(...)))
+    MarkovMap(fs::Vector, ds::Vector, ran = coveringinterval(ds); dir=Forward, diff=autodiff(...)))
 
 Generate a MarkovMap with branches given by elements of `fs`` defined on subdomains given by `ds`, onto a vector `ran`.
 
 The keyword argument `dir` stipulates whether the elements of `fs` are the branches (`Forward`) or the branches' inverses (`Reverse`).
 
 The keyword argument `diff` provides the derivatives of the `fs`. By default it is the automatic derivatives of `fs`.
 """
-function MarkovMap(fs::AbstractVector,ds::AbstractVector,ran=coveringsegment(ds);dir=Forward,
+function MarkovMap(fs::AbstractVector,ds::AbstractVector,ran=coveringinterval(ds);dir=Forward,
           diff=[autodiff(fs[i],(dir==Forward ? ds[i] : ran)) for i in eachindex(fs)])
   @assert length(fs) == length(ds)
   randm=convert(Domain,ran)
   MarkovMap([branch(fs[i],convert(Domain,ds[i]),randm,diff[i];dir=dir,ftype=eltype(fs),difftype=eltype(diff)) for i in eachindex(fs)],
-        coveringsegment(ds),convert(Domain,ran))
+        coveringinterval(ds),convert(Domain,ran))
 end
 
 
@@ -79,8 +79,8 @@ struct IntervalMap{D<:Domain,R<:Domain,B<:AbstractBranch} <: AbstractIntervalMap
   rangedomain::R
   @compat function IntervalMap{D,R,B}(branches::Vector{B},dom::D,ran::R) where
           {B<:ExpandingBranch, D<:Domain, R<:Domain}
-    @assert all(b->issubset(b.domain,dom),branches)
-    @assert all(b->issubset(b.rangedomain,ran),branches)
+    @assert all(b->approx_issubset(b.domain,dom),branches)
+    @assert all(b->approx_issubset(b.rangedomain,ran),branches)
     new(branches,dom,ran)
   end
 end
@@ -93,13 +93,16 @@ function IntervalMap(branches::AbstractVector{B},dom,ran) where {B<:ExpandingBra
 end
 
 function IntervalMap(fs::AbstractVector,ds::AbstractVector,
-          ran=coveringsegment([mapinterval(fs[i],ds[i]) for i in eachindex(fs)]);dir=Forward,
+          ran=coveringinterval([mapinterval(fs[i],ds[i]) for i in eachindex(fs)]);dir=Forward,
           diff=[autodiff(fs[i],(dir==Forward ? ds[i] : ran)) for i in eachindex(fs)])
   rand = convert(Domain,ran)
   @assert length(fs) == length(ds)
-  @assert all(issubset(mapinterval(fs[i],ds[i]),rand) for i in eachindex(fs))
+  println(rand)
+  println(extrema(mapinterval.(fs[1],ds[1])))
+  println((approx_issubset(mapinterval(fs[i],ds[i]),rand) for i in eachindex(fs))|>collect)
+  @assert all(approx_issubset(mapinterval(fs[i],ds[i]),rand) for i in eachindex(fs))
   IntervalMap([branch(fs[i],convert(Domain,ds[i]),mapinterval(fs[i],ds[i]),diff[i];dir=dir,ftype=eltype(fs),difftype=eltype(diff)) for i in eachindex(fs)],
-        coveringsegment(ds),rand)
+        coveringinterval(ds),rand)
 end
 
 const SimpleBranchedMap{D<:Domain,R<:Domain,B<:AbstractBranch} = Union{MarkovMap{D,R,B},IntervalMap{D,R,B}} #?? ComposedMap
@@ -140,9 +143,9 @@ eachbranchindex(m::SimpleBranchedMap) = 1:nbranches(m)
 
 neutralfixedpoints(m::SimpleBranchedMap) = [nfp(b) for b in branches(m)[isa.(NeutralBranch,branches(m))]]
 nneutral(m::SimpleBranchedMap) = sum([isa(b,NeutralBranch) for b in m.branches])
-getbranchind(m::SimpleBranchedMap,x) = temp_in(x,m.domain) ? findfirst([temp_in(x,domain(b)) for b in m.branches]) : error("DomainError: $x ∉ $(m.domain)")
-getbranchind(m::SimpleBranchedMap,x::IntervalDomain) = issubset(x,m.domain) ? findfirst([issubset(x,domain(b)) for b in m.branches]) : error("DomainError: $x ⊈ $(m.domain)")
-getbranchind(m::SimpleBranchedMap,x::Segment) = getbranchind(m,convert(Domain,x))
+getbranchind(m::SimpleBranchedMap,x) = x ∈ domain(m) ? findfirst([(x ∈ domain(b)) for b in m.branches]) : error("DomainError: $x ∉ $(domain(m))")
+getbranchind(m::SimpleBranchedMap,x::AbstractInterval) = (x ⊆ domain(m)) ? findfirst([(x ⊆ domain(b)) for b in m.branches]) : error("DomainError: $x ⊈ $(domain(m))")
+getbranchind(m::SimpleBranchedMap,x::Interval) = getbranchind(m,convert(Domain,x))
 branchindtype(m::SimpleBranchedMap) = Int
 
 # Transfer function
@@ -176,7 +179,7 @@ end
 
 function forwardmodulomap(f::ff,dom,ran=dom,diff=autodiff(f,dom)) where {ff}
   domd = convert(Domain,dom); randm = convert(Domain,ran)
-  fa = f(first(domd)); fb = f(last(domd))
+  fa = f(leftendpoint(domd)); fb = f(rightendpoint(domd))
   L = arclength(randm)
   nb_est = (fb-fa)/L; nb = round(Int,nb_est) # number of branches
   σ = sign(nb); NB = abs(nb)
@@ -186,12 +189,12 @@ function forwardmodulomap(f::ff,dom,ran=dom,diff=autodiff(f,dom)) where {ff}
   @assert nb != 0
 
   breakpoints = Array{eltype(domd)}(undef,NB+1)
-  breakpoints[1] = first(domd)
+  breakpoints[1] = leftendpoint(domd)
 
   for i = 1:NB-1
     breakpoints[i+1] = domain_newton(f,diff,fa+σ*i*L,domd)
   end
-  breakpoints[end] = last(domd)
+  breakpoints[end] = rightendpoint(domd)
 
   fs = [FwdOffset(f,(i-(σ == 1))*σ*L) for i = 1:NB]
   ds = [Interval(breakpoints[i],breakpoints[i+1]) for i = 1:NB]
@@ -210,9 +213,9 @@ reversemodulomap(args...) = _NI("reversemodulomap")
 # function reversemodulomap{ff}(v::ff,dom,ran=dom,diff=autodiff(v,dom);maxcover=10000)
 #   domd = convert(Domain,dom); randm = convert(Domain,ran)
 #   dr = arclength(randm); dd = arclength(domd)
-#   ra = first(randm); va = v(ra); vb = v()
+#   ra = leftendpoint(randm); va = v(ra); vb = v()
 #   in(va,∂(domd)) || error("v($a) not in boundary of domain")
-#   flip = (va ≈ last(domd))
+#   flip = (va ≈ rightendpoint(domd))
 #   drflip = (-1)^flip * dr
 #
 #   NB = 1; vr = v(ra+drflip)-va