# Reducing verbosity in array calculations

**URL:** <https://discourse.julialang.org/t/reducing-verbosity-in-array-calculations/8778>\
**Category:** General Usage\
**Tags:** question, statistics, array, syntax\
**Created:** [February 2, 2018, 2:10pm UTC](https://discourse.julialang.org/t/reducing-verbosity-in-array-calculations/8778 "2018-02-02T14:10:29Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![mvhulten](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mvhulten/32/2444_2.png) [@mvhulten](https://discourse.julialang.org/u/mvhulten)\
**Post date:** [February 2, 2018, 2:10pm UTC](https://discourse.julialang.org/t/reducing-verbosity-in-array-calculations/8778/1 "2018-02-02T14:10:29Z")

</div>

I wrote simple code, but it looks messy because I am removing NaN and I usually want to apply element-by-element. To illustrate, I show my calculation of the Pearson correlation coefficient:

```julia
function r(obs, mod)
    om = (obs - mean(obs[!isnan.(obs)])) .* (mod - mean(mod[!isnan.(mod)]))
    oo = (obs - mean(obs[!isnan.(obs)])).^2
    mm = (mod - mean(mod[!isnan.(mod)])).^2
    sum(om[!isnan.(om)]) ./ sqrt.(sum(oo[!isnan.(oo)]) .* sum(mm[!isnan.(mm)]))
end

```

where `obs` and `mod` are arrays of floats of arbitrary dimension. It looks messy because of the additions of:

- `[!isnan.()]` to the argument of _every_ use of `sum()` and `mean()`; and
- `.` for the functions `isnan()` and `sqrt()`, but of course _not_ for `sum()` and `mean()`.

There must be a less verbose, or more elegant, way to write this?

If you have comments on the functional part of the code, like bugs, those are welcome as well! My intention was to codify the equation for r, given by [Stow et al., 2009](https://doi.org/10.1016/j.jmarsys.2008.03.011) ([alternative source](https://dacemirror.sci-hub.hk/journal-article/a829a09a3626e0f60b6d77129cbd00e4/stow2009.pdf)).

---

<div class="post-metadata">

**Author:** ![y4lu](https://avatars.discourse-cdn.com/v4/letter/y/47e85d/32.png) [@y4lu](https://discourse.julialang.org/u/y4lu)\
**Post date:** [February 2, 2018, 2:23pm UTC](https://discourse.julialang.org/t/reducing-verbosity-in-array-calculations/8778/2 "2018-02-02T14:23:23Z")

</div>

Should give the same result i think

```julia
function r(obs, mod)
    obs2 = obs[!isnan.(obs)];
    mod2 = mod[!isnan.(mod)]; 
    om = (obs2 - mean(obs2)) .* (mod2 - mean(mod2))
    oo = (obs2 - mean(obs2)).^2
    mm = (mod2 - mean(mod2)).^2
    sum(om) ./ sqrt.(sum(oo) .* sum(mm))
end

```

---

<div class="post-metadata">

**Author:** ![tshort](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tshort/32/43_2.png) [@tshort](https://discourse.julialang.org/u/tshort)\
**Post date:** [February 2, 2018, 2:35pm UTC](https://discourse.julialang.org/t/reducing-verbosity-in-array-calculations/8778/3 "2018-02-02T14:35:30Z")

</div>

For NaN’s, the following might also help:

> **[GitHub - JuliaMath/NaNMath.jl: Julia math built-ins which return NaN and...](https://github.com/JuliaMath/NaNMath.jl)**
>
> Julia math built-ins which return NaN and accumulator functions which ignore NaN - GitHub - JuliaMath/NaNMath.jl: Julia math built-ins which return NaN and accumulator functions which ignore NaN

---

<div class="post-metadata">

**Author:** ![mvhulten](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mvhulten/32/2444_2.png) [@mvhulten](https://discourse.julialang.org/u/mvhulten)\
**Post date:** [February 2, 2018, 2:35pm UTC](https://discourse.julialang.org/t/reducing-verbosity-in-array-calculations/8778/4 "2018-02-02T14:35:32Z")

</div>

It looks cleaner, but then I tested it and found out that it does not do the same. With the original function I got a number \in [0, 1], but with your code I get `NaN`.

The `NaN` entries are, generally, at different indices in the two arrays `obs` and `mod`. Subsequently, `obs2` and `mod2` generally have different dimensions, and the indices of `NaN` don’t match.

Minimal example:

```julia
julia> a = [1 2 3]
1×3 Array{Int64,2}:
 1 2 3
julia> b = [1 NaN 3]
1×3 Array{Float64,2}:
 1.0 NaN 3.0
julia> a.*b
1×3 Array{Float64,2}:
 1.0 NaN 9.0
julia> aa = a[!isnan(a)]
3-element Array{Int64,1}:
 1
 2
 3
julia> bb = b[!isnan(b)]
2-element Array{Float64,1}:
 1.0
 3.0
julia> aa.*bb
ERROR: DimensionMismatch("arrays could not be broadcast to a common size")
Stacktrace:
...

```

This new code is equivalent to the original:

```julia
function r2(obs, mod)
    obs2 = obs[!isnan.(obs)];
    mod2 = mod[!isnan.(mod)];
    om = (obs - mean(obs2)) .* (mod - mean(mod2))
    oo = (obs - mean(obs2)).^2
    mm = (mod - mean(mod2)).^2
    sum(om[!isnan.(om)]) ./ sqrt.(sum(oo[!isnan.(oo)]) .* sum(mm[!isnan.(mm)]))
end

```

Sadly, it does not as nice as @y4lu’s code.

---

<div class="post-metadata">

**Author:** ![mvhulten](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mvhulten/32/2444_2.png) [@mvhulten](https://discourse.julialang.org/u/mvhulten)\
**Post date:** [February 2, 2018, 2:43pm UTC](https://discourse.julialang.org/t/reducing-verbosity-in-array-calculations/8778/5 "2018-02-02T14:43:08Z")

</div>

Yes, I’ve seen this package. Can be very useful. Its main use appears to be _creating_ `NaN`s, but I already have those in my data! 🙂

Of course, I can then type `NaNMath.mean(x)` instead of `mean(x,!isnan(x))`, and one can shorten by first defining `nm=NaNMath`. It looks a bit cleaner, but one depends on an extra package.

It’s a good option, but I am not fully satisfied, well, maybe this is it and I should be.

**edit:** I’m reconsidering this; this looks quite good:

```julia
function r3(obs, mod)
    #using NaNMath
    nm = NaNMath
    om = (obs - nm.mean(obs)) .* (mod - nm.mean(mod))
    oo = (obs - nm.mean(obs)).^2
    mm = (mod - nm.mean(mod)).^2
    nm.sum(om) ./ sqrt.(nm.sum(oo) .* nm.sum(mm))
end

```

`using` in a function is not accepted, probably for good reasons; but that’s fine.

---

<div class="post-metadata">

**Author:** ![piever](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/piever/32/1815_2.png) [@piever](https://discourse.julialang.org/u/piever)\
**Post date:** [February 2, 2018, 3:59pm UTC](https://discourse.julialang.org/t/reducing-verbosity-in-array-calculations/8778/6 "2018-02-02T15:59:24Z")

</div>

A technical point, but it is generally safer to use `const` when renaming a module. For example:

```julia
import NaNMath
const nm = NaNMath

function r3(obs, mod)
...
end

```

---

<div class="post-metadata">

**Author:** ![y4lu](https://avatars.discourse-cdn.com/v4/letter/y/47e85d/32.png) [@y4lu](https://discourse.julialang.org/u/y4lu)\
**Post date:** [February 3, 2018, 12:42am UTC](https://discourse.julialang.org/t/reducing-verbosity-in-array-calculations/8778/7 "2018-02-03T00:42:11Z")

</div>

This doesn’t look too bad still (Switched a few operators to / from elementwise for clarity)

```julia
function r(obs, mod)
  v1 = !isnan.(obs)
  v2 = !isnan.(mod)
  v3 = v1 & v2;
  obs3 = obs[v3] .- mean(obs[v1])
  mod3 = mod[v3] .- mean(mod[v2])
  sum(obs3 .* mod3) / sqrt.( sum(obs3 .^ 2) * sum(mod3 .^ 2))
end
```

---

<div class="post-metadata">

**Author:** ![mvhulten](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mvhulten/32/2444_2.png) [@mvhulten](https://discourse.julialang.org/u/mvhulten)\
**Post date:** [February 3, 2018, 5:52pm UTC](https://discourse.julialang.org/t/reducing-verbosity-in-array-calculations/8778/8 "2018-02-03T17:52:47Z")

</div>

That does not do the same as my code: there must be a `!isnan.()` criterion for the arguments in each `sum()`. But it gave me some ideas, and decided that this function is probably equivalent to my original code:

```julia
using NaNMath
function r32(obs, mod)
    const nm = NaNMath
    dobs = obs - nm.mean(obs)
    dmod = mod - nm.mean(mod)
    nm.sum(dobs .* dmod) ./ sqrt.(nm.sum(dobs .^ 2) .* nm.sum(dmod .^ 2))
end

```

I said _probably_ because it is based on superficial code inspection and one simple input test, namely `a=Float64.([1,2,3]); b=Float64.([1,2,NaN])` (no integers because of [issue](https://github.com/mlubin/NaNMath.jl/issues/26) in `NaNMath`).

But I just found out that all above code is wrong. They _appear_ to be implementing r from Stow et al. (2009), but they don’t. I made a thinking error in selecting my indices. Every time I use an `obs`, I need to mask the `NaN` from both `obs` _and_ `mod`. Same for `mod`. If I adjust the latest code here, the whole `NaNMath` is not useful anymore and the code becomes verbose.

To make a long story short, I believe that this is the correct code, and very readable as well:

```julia
function r23(obs, mod)
  bis = !isnan.(obs) .& !isnan.(mod); # == !isnan.(obs.*mod)
  dobs = obs[bis] - mean(obs[bis])
  dmod = mod[bis] - mean(mod[bis])  
  sum((dobs .* dmod)) / sqrt( sum((dobs .^ 2)) * sum((dmod .^ 2)) )
end

```

where `isobs` means _there is an observation defined at that array index_, `ismod` similarly, and `bis` stands for _both are_ or _twice_ (quite clever, I thought 😄).

---

<div class="post-metadata">

**Author:** ![ScottPJones](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/scottpjones/32/146_2.png) [@ScottPJones](https://discourse.julialang.org/u/ScottPJones)\
**Post date:** [February 5, 2018, 7:06pm UTC](https://discourse.julialang.org/t/reducing-verbosity-in-array-calculations/8778/9 "2018-02-05T19:06:48Z")

</div>

There are a few problems in the code for `r23`, it gets a deprecation warning on v0.6, and doesn’t work at all on v0.7 (master).

Also, if you need decent performance, while the above is correct and readable, it is about twice as slow as a version that does no allocations.

I’ve placed them in the following Gist: [Code and Benchmark functions](https://gist.github.com/ScottPJones/1dbcd665bc442fb76729b62a0f3ab457)

Running them on my MacBookPro, my version was about twice as fast on v0.6.2, and 2.5x as fast on master [v0.7.0-DEV.3716 (2018-02-05 03:16 UTC)].

Note: there is also a significant performance regression between v0.6.2 and v0.7 here, on v0.7 my version was about 10% slower, and yours was about 36% slower than on v0.6.2. That’s something that needs to be addressed, IMO.