How to create a `range` centered around zero?

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
3 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

@mbauman because step is a Real number. Maybe I should have written that step=dx instead of n

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 :slight_smile:

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

@Henrique_Becker I thought so :smiley:

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

When I saw mcabbott starting to type, I knew I better hurry up and push post. :slight_smile:

7 Likes

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 :slight_smile:

4 Likes

:grimacing:

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