diff --git a/.travis.yml b/.travis.yml
index d93d96e..ee9667c 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -5,6 +5,7 @@ os:
   - osx
 julia:
   - 0.7
+  - 1.0
 notifications:
   email: false
 matrix:
diff --git a/README.md b/README.md
index 90c1ccb..ac03c24 100644
--- a/README.md
+++ b/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)
 ```
diff --git a/REQUIRE b/REQUIRE
index aa101d1..20ccc62 100644
--- a/REQUIRE
+++ b/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
diff --git a/images/images.jl b/images/images.jl
index 283e8cf..0b3e21e 100644
--- a/images/images.jl
+++ b/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)
diff --git a/src/HGraphs.jl b/src/HGraphs.jl
index d99b4e0..62da17b 100644
--- a/src/HGraphs.jl
+++ b/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}
diff --git a/src/Poltergeist.jl b/src/Poltergeist.jl
index 87ca1af..0f7c18d 100644
--- a/src/Poltergeist.jl
+++ b/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")
diff --git a/src/general.jl b/src/general.jl
index 1776984..3d8a6d6 100644
--- a/src/general.jl
+++ b/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)
 
diff --git a/src/hofbauerextension.jl b/src/hofbauerextension.jl
index c1ebdf5..3c88de4 100644
--- a/src/hofbauerextension.jl
+++ b/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
diff --git a/src/maps/CircleMap.jl b/src/maps/CircleMap.jl
index d536c33..da18a61 100644
--- a/src/maps/CircleMap.jl
+++ b/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)
diff --git a/src/maps/ComposedMaps.jl b/src/maps/ComposedMaps.jl
index 3dadcaa..160dd86 100644
--- a/src/maps/ComposedMaps.jl
+++ b/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
diff --git a/src/maps/ExpandingBranch.jl b/src/maps/ExpandingBranch.jl
index 6a58aef..b5eb812 100644
--- a/src/maps/ExpandingBranch.jl
+++ b/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)
diff --git a/src/maps/IntervalMap.jl b/src/maps/IntervalMap.jl
index 84a3c82..4c23790 100644
--- a/src/maps/IntervalMap.jl
+++ b/src/maps/IntervalMap.jl
@@ -54,7 +54,7 @@ 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`.
 
@@ -62,12 +62,12 @@ The keyword argument `dir` stipulates whether the elements of `fs` are the branc
 
 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
diff --git a/src/maps/conjugacy.jl b/src/maps/conjugacy.jl
index 84e2a05..f940aa6 100644
--- a/src/maps/conjugacy.jl
+++ b/src/maps/conjugacy.jl
@@ -1,6 +1,6 @@
 # MarkovMaps
 function inv(m::MarkovMap)
-  ~isa(domain(m),IntervalDomain) && error("Non-interval domains not supported")
+  ~isa(domain(m),AbstractInterval) && error("Non-interval domains not supported")
   nbranches(m) > 1 && error("Cannot invert map with multiple branches")
   domain(m.branches[1]) != domain(m) && error("Cannot invert map that isn't defined everywhere")
   unsafe_inv(m)
diff --git a/src/maps/examples.jl b/src/maps/examples.jl
index 5621b9b..b7c8284 100644
--- a/src/maps/examples.jl
+++ b/src/maps/examples.jl
@@ -18,21 +18,20 @@ lanford(T=Float64) = MarkovMap([lanford_v0,lanford_v1],[0..lanford_brk(T),lanfor
 
 # tupling map
 """
-    tupling(k::Int, d = Segment(0,1.))
+    tupling(k::Int, d = Interval(0,1.))
 
 Returns the full-branch interval map on domain d with k equally-sized branches.
 
 See also: [`doubling`](@ref), [`lanford`](@ref)
 """
-tupling(k::Int,d) = tupling(k,convert(Domain,d))
-tupling(k::Int,d::Segment=Segment(0,1.)) = modulomap(x->k*(x-d.a)+d.a,d)
-tupling(k::Int,d::PeriodicInterval) = RevCircleMap(x->(x-d.a)/k+d.a,d,d,x->one(x)/k)
+tupling(k::Int,d=Interval(0,1.)) = modulomap(x->k*(x-leftendpoint(d))+leftendpoint(d),d)
+tupling(k::Int,d::PeriodicSegment) = RevCircleMap(x->(x-leftendpoint(d))/k+leftendpoint(d),d,d,x->one(x)/k)
 
 """
-    doubling(d = Segment(0,1.))
+    doubling(d = Interval(0,1.))
 
 Returns the full-branch interval map on domain d with 2 equally-sized branches.
 
 See also: [`tupling`](@ref), [`lanford`](@ref)
 """
-doubling(d=Segment(0,1.)) = tupling(2,d)
+doubling(d=Interval(0,1.)) = tupling(2,d)
diff --git a/src/poetry.jl b/src/poetry.jl
index 86b6bd6..4bc8e32 100644
--- a/src/poetry.jl
+++ b/src/poetry.jl
@@ -35,7 +35,7 @@ Base.adjoint(b::AbstractMarkovMap) = MarkovMapDerivative(b)
 Construct a self-map on domain d: x ↦ x + ϵ X(x)
 """
 perturb(d,X,ϵ) = perturb(convert(Domain,d),X,ϵ)
-perturb(d::IntervalDomain,X,ϵ) = MarkovMap([x->x+ϵ*X(x)],[d],d)
+perturb(d::AbstractInterval,X,ϵ) = MarkovMap([x->x+ϵ*X(x)],[d],d)
 perturb(d::PeriodicDomain,X,ϵ) = FwdCircleMap([x->x+ϵ*X(x)],d)
 
 """
diff --git a/src/spectral/Transfer.jl b/src/spectral/Transfer.jl
index edfffde..bbc97c6 100644
--- a/src/spectral/Transfer.jl
+++ b/src/spectral/Transfer.jl
@@ -85,8 +85,8 @@ transferfunction_nodes(L::ConcreteTransfer{TT,D,R,M},n::Integer,kk) where {TT,D,
 
 
 function transfer_getindex(L::ConcreteTransfer{T},
-    jdat::Tuple{Integer,Integer,Union{Integer,Infinity{Bool}}},
-    k::AbstractRange,padding::Bool=false) where T
+    jdat::Tuple{Integer,Integer,Union{Integer,Infinity}},
+    k::OrdinalRange,padding::Bool=false) where T
   Tr = real(T)
   @compat dat = Array{T}(undef,0)
   @compat cols = Array{eltype(k)}(undef,Base.length(k)+1)
@@ -159,7 +159,7 @@ function transfer_getindex(L::ConcreteTransfer{T},
       abs(coeffs[i])<tol*maxabsc*log2(lcfc)/10 && (coeffs[i] = 0)
     end
 
-    cutcfc = coeffs[(jdat[1]:jdat[2]:min(lcfc,jdat[3]))::AbstractRange]
+    cutcfc = coeffs[(jdat[1]:jdat[2]:min(lcfc,jdat[3]))::OrdinalRange]
 
     K = max(K,length(cutcfc))
     padding && ApproxFun.pad!(cutcfc,K)
@@ -172,15 +172,14 @@ function transfer_getindex(L::ConcreteTransfer{T},
   RaggedMatrix(dat,cols,K)
 end
 
-Base.getindex(L::ConcreteTransfer,j::AbstractRange,k::AbstractRange) = transfer_getindex(L,(start(j),step(j),last(j)),k)
-Base.getindex(L::ConcreteTransfer,j::Colon,k::AbstractRange) = transfer_getindex(L,(1,1,ApproxFun.∞),k)
-Base.getindex(L::ConcreteTransfer,j::ApproxFun.AbstractCount,k::AbstractRange) = transfer_getindex(L,(start(j),step(j),ApproxFun.∞),k)
+Base.getindex(L::ConcreteTransfer,j::OrdinalRange,k::OrdinalRange) = transfer_getindex(L,(start(j),step(j),last(j)),k)
+Base.getindex(L::ConcreteTransfer,j::Colon,k::OrdinalRange) = transfer_getindex(L,(1,1,ApproxFun.∞),k)
 Base.getindex(L::ConcreteTransfer,j::Integer,k::Integer) = Base.getindex(L,j:j,k:k)[1,1]
-Base.getindex(L::ConcreteTransfer,j::Integer,k::AbstractRange) = Base.getindex(L,j:j,k).data
-Base.getindex(L::ConcreteTransfer,j::AbstractRange,k::Integer) = Base.getindex(L,j,k:k).data
+Base.getindex(L::ConcreteTransfer,j::Integer,k::OrdinalRange) = Base.getindex(L,j:j,k).data
+Base.getindex(L::ConcreteTransfer,j::OrdinalRange,k::Integer) = Base.getindex(L,j,k:k).data
 
 Base.convert(::Type{RaggedMatrix},S::ApproxFun.SubOperator{T,LL,Tuple{R1,R2}}) where
-        {T,LL<:AbstractTransfer,R1<:Union{AbstractRange,ApproxFun.AbstractCount},R2<:AbstractRange} =
+        {T,LL<:AbstractTransfer,R1<:OrdinalRange,R2<:OrdinalRange} =
   Base.getindex(parent(S),parentindexes(S)[1],parentindexes(S)[2])
 
 
@@ -216,7 +215,7 @@ function ApproxFun.colstop(co::ApproxFun.CachedOperator{T,DM,CT},n::Integer) whe
 end
 
 # # fast colstop for cached operators
-# function ApproxFun.colstop{T,M<:ConcreteTransfer,DS<:Space,RS<:Space}(B::ApproxFun.CachedOperator{T,RaggedMatrix{T},M,DS,RS,Tuple{Infinity{Bool}}},k::Integer)
+# function ApproxFun.colstop{T,M<:ConcreteTransfer,DS<:Space,RS<:Space}(B::ApproxFun.CachedOperator{T,RaggedMatrix{T},M,DS,RS,Tuple{Infinity},k::Integer)
 #   ApproxFun.resizedata!(B,:,k)
 #   ApproxFun.colstop(co.op,k)
 # end
diff --git a/src/timeseries.jl b/src/timeseries.jl
index d7ef1c2..e456ff4 100644
--- a/src/timeseries.jl
+++ b/src/timeseries.jl
@@ -19,14 +19,15 @@ function timeseries_initialvalue(m::AbstractMarkovMap,n::Integer,xinit::T) where
   x_ts
 end
 timeseries(m::AbstractMarkovMap,n::Integer;x0=map_n(m,rand(domain(m)),10^4)) = timeseries_initialvalue(m,n,x0)
-samplable(ρ::Fun) = isa(domain(ρ),IntervalDomain) ? ρ : Fun(ρ,Segment(domain(ρ)))
+samplable(ρ::Fun) = isa(domain(ρ),AbstractInterval) ? ρ : Fun(ρ,Space(Interval(domain(ρ))))
+        # probably Interval(..) should be a convert routine but that's how ApproxFun has it...
 timeseries(m::AbstractMarkovMap,n::Integer,ρ::Fun) = timeseries_initialvalue(m,n,sample(samplable(ρ)))
 
 
 function timehist_initialvalue(m::AbstractMarkovMap,n::Integer,nbins::Integer,xinit::T) where T
   @assert domain(m) == rangedomain(m)
-  bin_start = domain(m).a
-  bin_step = (domain(m).b-domain(m).a)/nbins
+  bin_start = leftendpoint(domain(m))
+  bin_step = (rightendpoint(domain(m))-bin_start)/nbins
   x_hist = zeros(Int,nbins)
   x = copy(xinit)
   for i = 1:n
@@ -34,7 +35,7 @@ function timehist_initialvalue(m::AbstractMarkovMap,n::Integer,nbins::Integer,xi
     x_hist[min(nbins,max(1,convert(Int,ceil((x-bin_start)/bin_step))))] += 1
   end
 
-  bin_start:bin_step:domain(m).b,x_hist
+  range(bin_start,stop=rightendpoint(domain(m)),length=nbins+1),x_hist
 end
 timehist(m::AbstractMarkovMap,n::Integer,nbins::Integer;x0=map_n(m,rand(domain(m)),10^4)) = timehist_initialvalue(m,n,nbins,x0)
 timehist(m::AbstractMarkovMap,n::Integer,nbins::Integer,ρ::Fun) = timehist_initialvalue(m,n,nbins,sample(samplable(ρ)))
diff --git a/test/runtests.jl b/test/runtests.jl
index 048da7b..d3ef623 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,9 +10,9 @@ f1d(x)=2+cos(2pi*x)/4; f2d = f1d
 fv1(x) = x/2+sin(2pi*x)/8pi; fv2(x) = x/2+1/2+sin(2pi*x)/8pi
 fv1d(x) = 1/2+cos(2pi*x)/4; fv2d = fv1d
 
-#Periodic domain
+# Periodic domain
 println("Fourier tests 🌚🌞")
-d1 = PeriodicInterval(0,1.)
+d1 = PeriodicSegment(0,1.)
 M1b = CircleMap(fv1,d1,dir="rev",diff=fv1d)
 L1b = Transfer(M1b)
 @time L1b = Transfer(M1b)
@@ -32,8 +32,8 @@ println("Should all be ≤0.3s")
 
 # Non-periodic domain
 println("Chebyshev tests 🌳")
-d2 = Segment(0..1.)
-@test Poltergeist.coveringsegment([0..0.5,0.5..1]) == d2
+d2 = (0..1.)::Interval
+@test Poltergeist.coveringinterval([0..0.5,0.5..1]) == d2
 M2b = MarkovMap([fv1,fv2],[0..0.5,0.5..1],dir=Reverse,diff=[fv1d,fv2d]);
 acim(M2b)
 @time ρ2b = acim(M2b)
@@ -60,14 +60,14 @@ println("Lanford map test")
 lan = lanford()
 K = SolutionInv(lan);
 rho = acim(K);
-l_exp = sum(Fun(x->log(abs(lan'(x))),0..1) * rho)
-sigmasq_A = birkhoffvar(K,Fun(x->x^2,0..1))
+l_exp = sum(Fun(x->log(abs(lan'(x))),0..1.) * rho)
+sigmasq_A = birkhoffvar(K,Fun(x->x^2,0..1.))
 L_lan = Transfer(lan)
 K = SolutionInv(L_lan);
 @time rho = acim(K);
 @time l_exp = lyapunov(K)
-@time l_exp2 = sum(Fun(x->log(abs(lan'(x))),0..1) * rho)
-@time sigmasq_A = birkhoffvar(K,Fun(x->x^2,0..1))
+@time l_exp2 = sum(Fun(x->log(abs(lan'(x))),0..1.) * rho)
+@time sigmasq_A = birkhoffvar(K,Fun(x->x^2,0..1.))
 
 @test l_exp ≈ 0.657661780006597677541582
 @test l_exp2 ≈ 0.657661780006597677541582
@@ -76,17 +76,17 @@ K = SolutionInv(L_lan);
 # modulomap test
 println("Modulomap and examples test")
 lan_lift(x) = 5x/2 - x^2/2
-lan = modulomap(lan_lift,0..1);
+lan = modulomap(lan_lift,0..1.);
 @test Transfer(lan)[1:100,1:100] ≈ L_lan[1:100,1:100]
-@test diag(Transfer(doubling(PeriodicInterval(6.,7.)))[1:10,1:10]) ≈ [1.;zeros(9)]
+@test diag(Transfer(doubling(PeriodicSegment(6.,7.)))[1:10,1:10]) ≈ [1.;zeros(9)]
 @test diag(Transfer(tupling(-4,0..4.))[1:10,1:10]) ≈ (-1/4).^(0:9)
 
 # @test diag(Transfer(modulomap(x->1-x/5,0..1,dir=Reverse))[1:10,1:10]) .≈ (-0.2).^(0:9)
 
 # Composing test
 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)
+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.)
 @time doublelan = lan ∘ lanshift ∘ shiftmap
 println("Should be ≤0.07s")
 doubleK = SolutionInv(doublelan)
@@ -96,7 +96,7 @@ println("Should be ≤0.4s")
 @test doublerho ≈ rho
 @test lyapunov(doubleK) ≈ 2l_exp
 
-lanpet = 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")
@@ -120,7 +120,7 @@ cs1f = correlationsum(M1f,A1)
 @test maximum(abs.(cs1f.(pts)-correlationsum(M2f,A2).(pts))) .< 2000eps(1.)
 
 println("Correlation function test")
-A = Fun(x->x^2,0..1); B = Fun(sin,0..1)
+A = Fun(x->x^2,0..1.); B = Fun(sin,0..1.)
 cA,cB = covariancefunction(lan,A,B)
 @test sum(cA)+sum(cB[2:end]) ≈ birkhoffcov(lan,A,B)
 @time covariancefunction(lan,A,B)
@@ -134,7 +134,7 @@ covariancefunction(lan,A,100)
 
 # Calling
 println("Newton's method test ☏")
-test_f = range(d2.a,stop=d2.b,length=20)[1:end-1] # map boundaries are dodgy because multivalued
+test_f = range(leftendpoint(d2),stop=rightendpoint(d2),length=20)[1:end-1] # map boundaries are dodgy because multivalued
 test_x = [Poltergeist.mapinv(M2b,1,tf) for tf in test_f]
  @test M2b.(test_x) ≈ test_f
  @test M1b.(test_x) ≈ test_f
@@ -142,9 +142,11 @@ test_x = [Poltergeist.mapinv(M2b,1,tf) for tf in test_f]
 
 # #Inducing
 println("Inducing tests 🐴")
-f = IntervalMap([x->φ*x,x->φ*x-1],[0..φ-1,φ-1..1],0..1)
+f = IntervalMap([x->φ*x,x->φ*x-1],[0..φ-1,φ-1..1],0..1.)
 # r = e-2
-he = hofbauerextension(f,Segment(φ-1..1),forcereturn=true)
+he = hofbauerextension(f,Interval(φ-1..1),forcereturn=true)
+println(he)
+println(he.hdomains)
 @test nv(he) == 2
 @test ne(he) == 3
 fi = InducedMap(he)
@@ -152,7 +154,7 @@ Lfi = Transfer(fi)
 @time Lfi[:,10]
 println("Should be ≤ ? s")
 # println(diag(Lfi[1:10,1:10]), 1 ./ (φ.^(1:10) - 1) ./ φ.^(1:10))
-@test all(diag(Lfi[1:10,1:10]) .≈ 1 ./ (φ.^(1:10) - 1) ./ φ.^(1:10))
+@test all(diag(Lfi[1:10,1:10]) .≈ 1 ./ (φ.^(1:10) .- 1) ./ φ.^(1:10))
 # M2bd = MarkovMap([fv1,fv2],[0..0.5,0.5..1],d2,dir=Reverse,diff=[fv1d,fv2d]);
 # M2bi = induce(M2bd,1)
 # # acim(M2bi)
@@ -182,8 +184,8 @@ println("Should be ≤2s")
 #   println("α = $α")
 #   @time b = NeutralBranch(x->1+2^α*x,x->2^α,α,0.6/2^α,Interval(0,0.5),Interval(0,1))
 #   # @time b = NeutralBranch(x->1+2^α*x,x->2^α,α,0.6/2^α,Interval(0,0.5),Interval(0,1))
-#   b2 = branch(x->(x+1)/2,x->0.5,Segment(0.5,1.),Segment(0.,1.),dir="rev")
-#   Mint = MarkovMap(Segment(0.,1),Segment(0.,1),[b,b2])
+#   b2 = branch(x->(x+1)/2,x->0.5,Interval(0.5,1.),Interval(0.,1.),dir="rev")
+#   Mint = MarkovMap(Interval(0.,1),Interval(0.,1),[b,b2])
 #
 #   Mint_I = induce(Mint,1)
 #   ρint = acim(Transfer(Mint_I))
@@ -201,7 +203,7 @@ println("Should be ≤2s")
 # using StaticArrays
 # standardmap_inv_lift(x::SVector) = SVector(x[1] - 0.1*sin(x[2] - x[1]),x[2]-x[1]);
 # standardmap_inv_diff(x::SVector) = SMatrix{2,2}(1,0,0,1); # As only determinant is important...
-# dom = PeriodicInterval()^2
+# dom = PeriodicSegment()^2
 #  # binv = branch(standardmap_inv_lift,standardmap_inv_diff,dom,dom,dir="rev"); # deprecated
 # binv= branch(standardmap_inv_lift,dom,dom,standardmap_inv_diff,dir=Reverse)
 # standardmap = MarkovMap([binv],dom,dom)