Name That Tune: I Can Code That In One Step

[blackeneth]

A highly fascinating function I’m writing requires, at one point, I create a vector of length nv that contains the result of this calculation:

abs(2*mod(q[i-1]+rand(),1)-1)

That is, each row needs to start with a random uniform value (a new one generated for each row), then you add q[i-1] to it, then take that result mod(,1) [which just takes the decimal result, e.g., mod(1.23,1) = 0.23], multiply it by 2, then subtract 1.

I haven’t been able to figure out how to do that in 1 step with Julia. Below is a snippet of code that illustrates the calculation, but it requires 2 lines of code to compute the x vector.

function tune() 
    q = [1.41421356, 1.73205080, 2.23606798]
    nv = 5

    for i in 2:4
        x=rand(1,nv)
        x=abs.(2*mod.(x .+ q[i-1],1) .- 1)
        # other fascinating code follows
        println(x)
    end 
end

An example of the output is:

 [0.3976719434201268 0.0972997022276254 0.5098824523243595 0.5008470075362657 0.5218571148803717]
 [0.33333003989780785 0.5247038763630942 0.7129593432507679 0.35352136344984864 0.5579598444795035]
 [0.14589717926504786 0.15900475794535218 0.5369653035112769 0.5487666346545677 0.09274380227621393]

Note the two lines that are needed to create the x vector.

How could I code it in 1 step in Julia?

[Sukera]

@blackeneth

x=abs.(2*mod.(x .+ q[i-1],1) .- 1)

x = abs.(2*mod.(rand(1,nv) .+ q[i-1],1) .- 1)

but I’m not sure why this is wrong, I’ve simply moved the first assignment into the second line?

[DNF]

@blackeneth

x=rand(1,nv)
x=abs.(2*mod.(x .+ q[i-1],1) .- 1)

You missed a dot: 2 .* mod. Forgetting this will create an extra allocation. Also, since you are reusing x you can dot the assignment as well. Often it’s best to use the @. macro so you don’t miss any dots:

x = rand(nv)
@. x = abs(2 * mod(x + q[i-1], 1) - 1)

@blackeneth

mod(,1) [which just takes the decimal result, e.g., mod(1.23,1) = 0.23],

For negative numbers this isn’t quite the case:

julia> mod(-0.2, 1)
0.8

On the other hand, rem does this (and is faster, too):

julia> rem(-0.2, 1)
-0.2

Edit: Actually, rem is basically the same speed as mod. Both are slow. If you want fast (and you may not really care), then check out trunc, floor or ceil:

myrem1(x) = x - trunc(x)  # changed name from mymod1 to myrem1

is much faster.

[yha]

@blackeneth

mod(q[i-1]+rand(),1)

Why would you do this calculation? Except for possible floating-point issues, or imperfections of the random number generation, adding q[i-1] has no effect on the (statistical) behavior of this code: the sequences [rand() for _ in q] and [mod(x+rand(),1) for x in q] have the same distribution, regardless of the values in q (both are i.i.d. sequences of samples from the uniform distribution on (0,1)).

(Unless perhaps, you really mean to do rem rather than mod, in which case I’m still curious what the use case is. I’ve never yet seen actual use of rem where mod would be incorrect.)

@DNF

mymod1(x) = x - trunc(x)

That’s computing rem rather than mod

EDIT: actually, the whole thing just computes i.i.d U(0,1) values in a convoluted way. The tune() function as given can be written (up to computational issues such as imperfect RNG)

tune() = foreach(println, rand(1,5) for _=1:3)

[DNF]

@yha

That’s computing rem rather than mod

Yes, I believe I made that clear.

Edit: Ah, the name is bad though. Should call it myrem1.

Edit 2: Actually, there is a function in Base, called modf which does basically the same thing, and perhaps should be called remf?

help?> modf

  modf(x)

  Return a tuple (fpart, ipart) of the fractional and integral parts of a number. Both parts have the same sign as the
  argument.

  Examples
  ≡≡≡≡≡≡≡≡≡≡

  julia> modf(3.5)
  (0.5, 3.0)
  
  julia> modf(-3.5)
  (-0.5, -3.0)

Also, it is quite slow:

julia> remf(x) = (a=trunc(x); (x-a, a))
remf (generic function with 1 method)

julia> @btime modf(x) setup=(x=10rand()-5);
  3.193 ns (0 allocations: 0 bytes)

julia> @btime remf(x) setup=(x=10rand()-5);
  1.738 ns (0 allocations: 0 bytes)

[blackeneth]

The purpose isn’t to generate random numbers. The purpose is to add random perturbation to the q[i-1] values.

Just where do those q values come from ?

q = [1.41421356, 1.73205080, 2.23606798]

Do 1.4142 and 1.732 look familiar?

q is constructed from the square root of the prime numbers.

The code is part of constructing Richtmyer sequences for quasi-Monte Carlo integration.

[blackeneth]

Thank you for your answer. I did indeed miss the dot for the “2”; being a numeric constant, it didn’t trigger the “dot impulse” in my brain.

I also didn’t know about the @. macro. I’m sure I’ll be using that a lot. Broadcasting is one of the nicer features of Julia.

[DNF]

@blackeneth

The purpose isn’t to generate random numbers. The purpose is to add random perturbation to the q[i-1] values.

But there’s no way to distinguish between rand() and mod(q + rand(), 1) no matter what q is. They are both random variables with identical distributions, i.e. uniformly distributed between 0 and 1.

So you can just replace one with the other.

[blackeneth]

The original line was before final debugging. The correct line is:

x=transpose(abs.(2 .* mod.(collect(UnitRange(1:nv)) .* q[i-1] .+ rand(), 1) .- 1)) 
# note: the rand() is not broadcast -- it's a single random uniform value added to all elements

[DNF]

OK, previously you made a random vector, rand(nv) (well actually, you made a random matrix, rand(1, nv)), so this is different, if it’s just one single random number to the entire range.

BTW, separately from the previous point: this construct is a bit odd:

collect(UnitRange(1:nv))

So, 1:nv is already a UnitRange. Putting it in a UnitRange constructor is redundant. Then, you convert it to an Array{Int, 1} using collect. Ranges have really nice properties, so I suggest that you just write (1:nv) .* q[i-1] instead of collect(UnitRange(1:nv)) .* q[i-1].

[blackeneth]

Thanks, I was hoping you’d have some tips on this.

collect(UnitRange(1:nv)) came about because it works in the REPL. I look up how to create a sequence from 1…n – find UnitRange(). Also see that collect() explodes the range out.

Consider (I’ll use “3” instead of q[i-1])

julia> x=(1:5)
1:5 
julia> x=(1:5) .* 3 
3:3:15 
julia> x=UnitRange(1:5) .* 3
3:3:15 
julia> x=collect(UnitRange(1:5)) .* 3
5-element Array{Int64,1}:
  3
  6
  9
 12
 15

So that’s how collect(UnitRange(1:nv)) came about.

I’ve now replaced the line with:

x=transpose(abs.(2 .* mod.((1:nv) .* q[i-1] .+ rand(),1) .- 1))

Which works

However, if I pull the guts out of the mod function and plug it into the REPL:

julia> (1:nv) .* q[i-1] .+ rand()
2.681369819129881:1.7320508075688772:577.722237931997

I usually check out the functionality in the REPL before putting it into the program.

[DNF]

Yes, after mod it cannot any longer be expressed as a range, so mod converts ranges to ordinary Array{.,1}s. Note, though, that ranges are arrays, they’re just not Arrays.

@blackeneth

So that’s how collect(UnitRange(1:nv)) came about.

Hmmmm. I still don’t quite get why you put it inside the UnitRange call I mean, you wouldn’t write s = String("Hello!"), right? Just s = "Hello!"?

[blackeneth]

Here’s how I think it happened.

Lookup UnitRange in the documentation.

The example given is:

julia> collect(UnitRange(2.3, 5.2))
3-element Array{Float64,1}:
 2.3
 3.3
 4.3

That’s just what I need! Modify it to use integers, and boom, good to go!

[DNF]

Wow, never seen that way of using it before.

Well, I guess you’re off the hook