is there a way to make a range(-max, max, step=n) that is guaranteed to go through zero?, i.e. I want it to be symmetrical/centered around zero
1 Like
From @briochemc on Slack
make a struct
struct SymmetricRange{A,B}
max::A
step::B
end
and define everything you need on it? E.g.,
Base.collect(r::SymmetricRange) = [reverse(collect(range(0, -r.max, step=-r.step)))[1:end-1]; collect(range(0, r.max, step=r.step))]
gives
julia> collect(SymmetricRange(1, 0.32))
7-element Vector{Float64}:
-0.96
-0.64
-0.32
0.0
0.32
0.64
0.96
4 Likes
From @Sukera:
also iterate and you should be set
won’t even need to roll your own collect then and it’ll work in loops just as well
4 Likes
I would not recommend the solutions already given but instead just create a function that returns a carefully assembled range -max:step:max that has such guarantee (i.e., by slightly changing both -max and max). Of course, if you need to deal with floating point numbers this may be a problem.
Why not range(-x, x, length=some_odd_number)? A few quick loops shows it seems to work pretty consistently, but floating point is wild and I don’t trust it fully. How important is the exact zero? How important are the endpoints? How important is the step size? You’re not going to be able to hit all of those requirements exactly no matter what you do.
3 Likes
Yeah, so recompute it to the equivalent length? Again, this gets to your requirements: it’s impossible to answer this question without knowing how exactly you want to hit endpoints, zero, and the step itself.
In MATLAB people do this:
x = [-1*fliplr(dx:dx:max), 0:dx:max];
Which I could do in Julia.e.g.
mysymrange(max, dx) = vcat(range(0,-max,step=-dx) |> reverse, range(0, max, step=dx))
… and get over with, but I would like to have a range 
julia> mysymrange(1.1, 0.2)
12-element Vector{Float64}:
-1.0
-0.8
-0.6
-0.4
-0.2
0.0
0.0
0.2
0.4
0.6
0.8
1.0
I think if you want (i) an exact zero (ii) working with floating points and (iii) not allocating more memory than necessary then the solution with a custom struct implementing iterate is, in fact, the best one. Maybe inherit AbstractRange too.
1 Like
You may be interested in using StepRangeLen directly:
help?> StepRangeLen
search: StepRangeLen StepRange
StepRangeLen{T,R,S}(ref::R, step::S, len, [offset=1]) where {T,R,S}
StepRangeLen( ref::R, step::S, len, [offset=1]) where { R,S}
A range r where r[i] produces values of type T (in the second form, T is deduced automatically), parameterized by a
reference value, a step, and the length. By default ref is the starting value r[1], but alternatively you can supply
it as the value of r[offset] for some other index 1 <= offset <= len. In conjunction with TwicePrecision this can be
used to implement ranges that are free of roundoff error.
julia> StepRangeLen(0, 0.35353, 11, 6)
-1.7676500000000002:0.35353:1.7676500000000002
julia> collect(StepRangeLen(0, 0.35353, 11, 6))
11-element Vector{Float64}:
-1.7676500000000002
-1.41412
-1.06059
-0.70706
-0.35353
0.0
0.35353
0.70706
1.06059
1.41412
1.7676500000000002
Set the reference point to 0 and set the offset to be the middle of your range.
4 Likes
I think it ought to be possible to re-use StepRangeLen for this, with a new constructor. Here’s a naiive one which takes a given step, and (if I understood it right) it measures distances from zero here (not by adding to one end) so should hit that exactly:
julia> function zed(step, length)
isodd(length) || error()
StepRangeLen(zero(step), step, length, 1 + lengthĂ·2)
end;
julia> zed(1/3, 7) |> collect
7-element Vector{Float64}:
-1.0
-0.6666666666666666
-0.3333333333333333
0.0
0.3333333333333333
0.6666666666666666
1.0
A constructor to always hit a given max exactly would probably need similar TwicePrecision tricks to what Base uses.
(Edit: scooped by @mkitti!)
Edit: It looks like range(-m, m; length) often produces exactly this, with ref == 0.0 and offset == (length+1)/2:
julia> m = nextfloat(123e123);
julia> range(-m, m, length=10^6 + 1) |> dump
StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}
ref: Base.TwicePrecision{Float64}
hi: Float64 0.0
lo: Float64 0.0
step: Base.TwicePrecision{Float64}
hi: Float64 2.459999999995291e119
lo: Float64 4.7094854703847126e107
len: Int64 1000001
offset: Int64 500001
But not always, here’s an example:
julia> m = 1.2057884f13; n = 105; r = range(-m, m, length=n)
-1.2057884f13:2.3188238f11:1.2057884f13
julia> r.ref, r.offset
(1.0, 53)
julia> r[[52, 53, 54]]
3-element Vector{Float32}:
-2.3188238f11
1.0
2.3188238f11
2 Likes
Another way, edited to be in line with requested behavior:
function range0(x, dx)
n = Int(floor(x/dx))
(-n:n)*dx
end
x = 1.1
dx = 0.25
julia> range0(x,dx)
-1.0:0.25:1.0
Hmmm… the problem with StepRangeLen is that I have to specify the length which I guess you thought it could me max but it is not an integer
Just divide the max by the step and floor it, double it, and add one. Again, though, you need to carefully consider your requirements that I keep trying to get you to define
4 Likes
I was just asking what is the problem with this:
julia> function srange(dx,xmax)
n = floor(xmax/dx)
smin = -n*dx
smax = n*dx
smin:dx:smax
end
srange (generic function with 1 method)
julia> collect(srange(0.8,5.6))
13-element Vector{Float64}:
-4.800000000000001
-4.000000000000001
-3.2000000000000006
-2.4000000000000004
-1.6000000000000005
-0.8000000000000007
0.0
0.7999999999999998
1.5999999999999996
2.3999999999999995
3.1999999999999993
4.0
4.800000000000001
1 Like
The problem is that we need obey this equation:
max = step * length
and have isinteger(length) == true.
We could do something like this, but as you can see it gets really messy very quickly.
julia> function zed_step_max(step::T, max::T) where T
length = max / step
isinteger(length) || error("step does not evenly divide max")
length = Integer(length)*2 + 1
StepRangeLen{T}(zero(step), step, length, 1 + lengthĂ·2)
end
zed_step_max (generic function with 2 methods)
julia> zed_step_max(0.0001, 0.3535)
ERROR: step does not evenly divide max
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:33
[2] zed_step_max(step::Float64, max::Float64)
@ Main .\REPL[156]:3
[3] top-level scope
@ REPL[157]:1
julia> zed_step_max(0.0001, 0.3536)
-0.3536:0.0001:0.3536
julia> 0.3535 / 0.0001
3534.9999999999995
julia> 0.3536 / 0.0001
3536.0
The problem is IEEE 754 floating point. Clearly though “0.3535” should be divisble by “0.001”, so how do we fix this to make it Just Work™?
To just get the result, we can try to bypass floating point arithmetic by using rationalize:
julia> function zed_step_max_rationalize(step::T, max::T) where T
step, max = rationalize.( (step,max) )
length = max / step
isinteger(length) || error("step does not evenly divide max")
length = Integer(length)*2 + 1
StepRangeLen{T}(zero(step), step, length, 1 + lengthĂ·2)
end
zed_step_max_rationalize (generic function with 1 method)
julia> zed_step_max_rationalize( 0.0001, 0.3535 )
-0.3535:0.0001:0.3535
julia> zed_step_max_rationalize( 0.0001, 0.3536 )
-0.3536:0.0001:0.3536
julia> r = zed_step_max_rationalize(0.0001, 0.3535)
mid = length(r)Ă·2 + 1
@info "" r[ 1] r[ 2] r[mid-1] r[mid ] r[mid+1] r[end-1] r[end]
┌ Info:
│ r[1] = -0.3535
│ r[2] = -0.3534
│ r[mid - 1] = -0.0001
│ r[mid] = 0.0
│ r[mid + 1] = 0.0001
│ r[end - 1] = 0.3534
â”” r[end] = 0.3535
julia> r = zed_step_max_rationalize(0.0001, 0.3536)
mid = length(r)Ă·2 + 1
@info "" r[ 1] r[ 2] r[mid-1] r[mid ] r[mid+1] r[end-1] r[end]
┌ Info:
│ r[1] = -0.3536
│ r[2] = -0.3535
│ r[mid - 1] = -0.0001
│ r[mid] = 0.0
│ r[mid + 1] = 0.0001
│ r[end - 1] = 0.3535
â”” r[end] = 0.3536
Inspired by @mbauman solution I think that this is my solution.
function srange(dx,xmax)
n = ceil(2*xmax/dx) |> Int
n = isodd(n) ? n : n+1
range(-xmax, xmax, length=n)
end
1 Like
Note, though, that won’t always hit zero exactly. Not sure if that’s important.
1 Like