From 42db28565c1a44e9eb73eb9b1b61be64d3eb9ac9 Mon Sep 17 00:00:00 2001 From: Stefan Karpinski Date: Tue, 10 Mar 2015 17:15:59 -0400 Subject: [PATCH] LinSpace: make linspace return a range-like type. --- base/array.jl | 45 -------------------------------- base/deprecated.jl | 2 ++ base/exports.jl | 1 + base/range.jl | 65 ++++++++++++++++++++++++++++++++++++++++------ test/ranges.jl | 40 ++++++++++++++-------------- 5 files changed, 80 insertions(+), 73 deletions(-) diff --git a/base/array.jl b/base/array.jl index 56926ba21d1c8..b61035ac137dc 100644 --- a/base/array.jl +++ b/base/array.jl @@ -246,51 +246,6 @@ function one{T}(x::AbstractMatrix{T}) eye(T, m) end -linspace(start::Integer, stop::Integer, n::Integer) = - linspace(float(start), float(stop), n) -function linspace(start::Real, stop::Real, n::Integer) - (start, stop) = promote(start, stop) - T = typeof(start) - a = Array(T, Int(n)) - if n == 1 - a[1] = start - return a - end - n -= 1 - S = promote_type(T, Float64) - for i = 0:n - a[i+1] = start*(convert(S, (n-i))/n) + stop*(convert(S, i)/n) - end - a -end -linspace(start::Real, stop::Real) = linspace(start, stop, 100) - -function linspace{T<:FloatingPoint}(start::T, stop::T, n::Int) - n == 1 && return [start] - n -= 1 - a0, b = rat(start) - a = convert(T,a0) - if a/convert(T,b) == start - c0, d = rat(stop) - c = convert(T,c0) - if c/convert(T,d) == stop - e = lcm(b,d) - a *= div(e,b) - c *= div(e,d) - ne = convert(T,n*e) - if a*n/ne == start && c*n/ne == stop - return [ (a*(n-k) + c*k)/ne for k=0:n ] - end - end - end - return [ start*((n-k)/n) + stop*(k/n) for k=0:n ] -end -linspace(start::FloatingPoint, stop::FloatingPoint, n::Integer) = - linspace(promote(start, stop)..., Int(n)) - -logspace(start::Real, stop::Real, n::Integer) = 10.^linspace(start, stop, n) -logspace(start::Real, stop::Real) = logspace(start, stop, 50) - ## Conversions ## convert{T,n}(::Type{Array{T}}, x::Array{T,n}) = x diff --git a/base/deprecated.jl b/base/deprecated.jl index 4b01ea608c8ee..d64330fe820e5 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -522,3 +522,5 @@ export float32_isvalid, float64_isvalid @deprecate parseint(s,base) parse(Int, s, base) @deprecate parseint(T::Type, s) parse(T, s) @deprecate parseint(T::Type, s, base) parse(T, s, base) + +@deprecate linrange linspace diff --git a/base/exports.jl b/base/exports.jl index de9a9548cadb9..ae6f6511f15bb 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -64,6 +64,7 @@ export IO, IOBuffer, IOStream, + LinSpace, LocalProcess, LowerTriangular, MathConst, diff --git a/base/range.jl b/base/range.jl index 0f92c8ed42ed3..904b9c82405b2 100644 --- a/base/range.jl +++ b/base/range.jl @@ -163,10 +163,44 @@ range(a::FloatingPoint, st::FloatingPoint, len::Integer) = FloatRange(a,st,len,o range(a::Real, st::FloatingPoint, len::Integer) = FloatRange(float(a), st, len, one(st)) range(a::FloatingPoint, st::Real, len::Integer) = FloatRange(a, float(st), len, one(a)) -linrange(a::Real, b::Real, len::Integer) = - len >= 2 ? range(a, (b-a)/(len-1), len) : - len == 1 && a == b ? range(a, zero((b-a)/(len-1)), 1) : - throw(ArgumentError("invalid range length")) +## linspace and logspace + +immutable LinSpace{T<:FloatingPoint} <: Range{T} + start::T + stop::T + len::Int + divisor::T +end +LinSpace(start::FloatingPoint, stop::FloatingPoint, len::Real, divisor::Real) = + LinSpace{promote_type(typeof(start),typeof(stop),typeof(divisor))}(start,stop,len,divisor) + +function linspace{T<:FloatingPoint}(start::T, stop::T, len::Int) + 0 <= len || error("invalid length: linspace($start, $stop, $len)") + 1 != len || start == stop || error("start != stop with len == 1: linspace($start, $stop, $len)") + n = len - 1 + a0, b = rat(start) + a = convert(T,a0) + if 0 < n && a/convert(T,b) == start + c0, d = rat(stop) + c = convert(T,c0) + if c/convert(T,d) == stop + e = lcm(b,d) + a *= div(e,b) + c *= div(e,d) + ne = convert(T,n*e) + if a*n/ne == start && c*n/ne == stop + return LinSpace(a, c, len, ne) + end + end + end + return LinSpace(start, stop, len, max(1,n)) +end +linspace(start::Real, stop::Real, len::Real) = + linspace(promote(FloatingPoint(start), FloatingPoint(stop))..., Int(len)) + +show(io::IO, r::LinSpace) = print(io, "linspace($(first(r)),$(last(r)),$(length(r)))") + +logspace(start::Real, stop::Real, n::Integer=50) = 10.^linspace(start, stop, n) ## interface implementations @@ -178,11 +212,13 @@ size(r::Range) = (length(r),) isempty(r::StepRange) = (r.start != r.stop) & ((r.step > zero(r.step)) != (r.stop > r.start)) isempty(r::UnitRange) = r.start > r.stop -isempty(r::FloatRange) = length(r)==0 +isempty(r::FloatRange) = length(r) == 0 +isempty(r::LinSpace) = length(r) == 0 step(r::StepRange) = r.step step(r::UnitRange) = 1 step(r::FloatRange) = r.step/r.divisor +step(r::LinSpace) = (r.stop - r.start)/r.divisor function length(r::StepRange) n = Integer(div(r.stop+r.step - r.start, r.step)) @@ -190,6 +226,7 @@ function length(r::StepRange) end length(r::UnitRange) = Integer(r.stop - r.start + 1) length(r::FloatRange) = Integer(r.len) +length(r::LinSpace) = r.len function length{T<:Union(Int,UInt,Int64,UInt64)}(r::StepRange{T}) isempty(r) && return zero(T) @@ -220,11 +257,13 @@ let smallint = (Int === Int64 ? end first{T}(r::OrdinalRange{T}) = convert(T, r.start) -first(r::FloatRange) = r.start/r.divisor +first{T}(r::FloatRange{T}) = convert(T, r.start/r.divisor) +first{T}(r::LinSpace{T}) = convert(T, (r.len-1)*r.start/r.divisor) last{T}(r::StepRange{T}) = r.stop last(r::UnitRange) = r.stop last{T}(r::FloatRange{T}) = convert(T, (r.start + (r.len-1)*r.step)/r.divisor) +last{T}(r::LinSpace{T}) = convert(T, (r.len-1)*r.stop/r.divisor) minimum(r::UnitRange) = isempty(r) ? throw(ArgumentError("range must be non-empty")) : first(r) maximum(r::UnitRange) = isempty(r) ? throw(ArgumentError("range must be non-empty")) : last(r) @@ -241,8 +280,14 @@ copy(r::Range) = r ## iteration start(r::FloatRange) = 0 -next{T}(r::FloatRange{T}, i::Int) = (convert(T, (r.start + i*r.step)/r.divisor), i+1) -done(r::FloatRange, i::Int) = (length(r) <= i) +done(r::FloatRange, i::Int) = length(r) <= i +next{T}(r::FloatRange{T}, i::Int) = + (convert(T, (r.start + i*r.step)/r.divisor), i+1) + +start(r::LinSpace) = 1 +done(r::LinSpace, i::Int) = length(r) < i +next{T}(r::LinSpace{T}, i::Int) = + (convert(T, ((r.len-i)*r.start + (i-1)*r.stop)/r.divisor), i+1) # NOTE: For ordinal ranges, we assume start+step might be from a # lifted domain (e.g. Int8+Int8 => Int); use that for iterating. @@ -268,6 +313,10 @@ function getindex{T}(r::FloatRange{T}, i::Integer) 1 <= i <= length(r) || throw(BoundsError()) convert(T, (r.start + (i-1)*r.step)/r.divisor) end +function getindex{T}(r::LinSpace{T}, i::Integer) + 1 <= i <= length(r) || throw(BoundsError()) + convert(T, ((r.len-i)*r.start + (i-1)*r.stop)/r.divisor) +end function check_indexingrange(s, r) sl = length(s) diff --git a/test/ranges.jl b/test/ranges.jl index 1982de1fbc594..1e64ff00363c4 100644 --- a/test/ranges.jl +++ b/test/ranges.jl @@ -232,22 +232,22 @@ end # tricky floating-point ranges -@test [0.1:0.1:0.3;] == linspace(0.1,0.3,3) == [1:3;]./10 -@test [0.0:0.1:0.3;] == linspace(0.0,0.3,4) == [0:3;]./10 -@test [0.3:-0.1:-0.1;] == linspace(0.3,-0.1,5) == [3:-1:-1;]./10 -@test [0.1:-0.1:-0.3;] == linspace(0.1,-0.3,5) == [1:-1:-3;]./10 -@test [0.0:0.1:1.0;] == linspace(0.0,1.0,11) == [0:10;]./10 -@test [0.0:-0.1:1.0;] == linspace(0.0,1.0,0) == [] -@test [0.0:0.1:-1.0;] == linspace(0.0,-1.0,0) == [] -@test [0.0:-0.1:-1.0;] == linspace(0.0,-1.0,11) == [0:-1:-10;]./10 -@test [1.0:1/49:27.0;] == linspace(1.0,27.0,1275) == [49:1323;]./49 -@test [0.0:0.7:2.1;] == linspace(0.0,2.1,4) == [0:7:21;]./10 -@test [0.0:1.1:3.3;] == linspace(0.0,3.3,4) == [0:11:33;]./10 -@test [0.1:1.1:3.4;] == linspace(0.1,3.4,4) == [1:11:34;]./10 -@test [0.0:1.3:3.9;] == linspace(0.0,3.9,4) == [0:13:39;]./10 -@test [0.1:1.3:4.0;] == linspace(0.1,4.0,4) == [1:13:40;]./10 -@test [1.1:1.1:3.3;] == linspace(1.1,3.3,3) == [11:11:33;]./10 -@test [0.3:0.1:1.1;] == linspace(0.3,1.1,9) == [3:1:11;]./10 +@test [0.1:0.1:0.3;] == [linspace(0.1,0.3,3);] == [1:3;]./10 +@test [0.0:0.1:0.3;] == [linspace(0.0,0.3,4);] == [0:3;]./10 +@test [0.3:-0.1:-0.1;] == [linspace(0.3,-0.1,5);] == [3:-1:-1;]./10 +@test [0.1:-0.1:-0.3;] == [linspace(0.1,-0.3,5);] == [1:-1:-3;]./10 +@test [0.0:0.1:1.0;] == [linspace(0.0,1.0,11);] == [0:10;]./10 +@test [0.0:-0.1:1.0;] == [linspace(0.0,1.0,0);] == [] +@test [0.0:0.1:-1.0;] == [linspace(0.0,-1.0,0);] == [] +@test [0.0:-0.1:-1.0;] == [linspace(0.0,-1.0,11);] == [0:-1:-10;]./10 +@test [1.0:1/49:27.0;] == [linspace(1.0,27.0,1275);] == [49:1323;]./49 +@test [0.0:0.7:2.1;] == [linspace(0.0,2.1,4);] == [0:7:21;]./10 +@test [0.0:1.1:3.3;] == [linspace(0.0,3.3,4);] == [0:11:33;]./10 +@test [0.1:1.1:3.4;] == [linspace(0.1,3.4,4);] == [1:11:34;]./10 +@test [0.0:1.3:3.9;] == [linspace(0.0,3.9,4);] == [0:13:39;]./10 +@test [0.1:1.3:4.0;] == [linspace(0.1,4.0,4);] == [1:13:40;]./10 +@test [1.1:1.1:3.3;] == [linspace(1.1,3.3,3);] == [11:11:33;]./10 +@test [0.3:0.1:1.1;] == [linspace(0.3,1.1,9);] == [3:1:11;]./10 @test [0.0:1.0:5.5;] == [0:10:55;]./10 @test [0.0:-1.0:0.5;] == [] @@ -275,7 +275,7 @@ for T = (Float32, Float64,),# BigFloat), vals = T[a:s:a+(n-1)*s;]./den r = start:step:stop @test [r;] == vals - @test linspace(start, stop, n) == vals + n == 1 || @test [linspace(start, stop, length(r));] == vals # issue #7420 n = length(r) @test [r[1:n];] == [r;] @@ -353,13 +353,13 @@ r = -0.004532318104333742:1.2597349521122731e-5:0.008065031416788989 @test_throws BoundsError r[0:10] @test_throws BoundsError r[1:10000] -r = linrange(1/3,5/7,6) +r = linspace(1/3,5/7,6) @test length(r) == 6 @test r[1] == 1/3 @test abs(r[end] - 5/7) <= eps(5/7) -r = linrange(0.25,0.25,1) +r = linspace(0.25,0.25,1) @test length(r) == 1 -@test_throws Exception linrange(0.25,0.5,1) +@test_throws Exception linspace(0.25,0.5,1) # issue #7426 @test [typemax(Int):1:typemax(Int);] == [typemax(Int)]