# Name That Tune: I Can Code That In One Step

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?

``````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?

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)
``````

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.

1 Like

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.)

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)
``````
1 Like

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)
``````

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.

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.

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.

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

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]`.

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.

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 `Array`s.

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!"`?

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!

Wow, never seen that way of using it before.

Well, I guess you’re off the hook