# Optimization of Code for NaN Operations

**URL:** https://discourse.julialang.org/t/optimization-of-code-for-nan-operations/31339
**Category:** Performance
**Created:** [November 21, 2019, 8:20am UTC](https://discourse.julialang.org/t/optimization-of-code-for-nan-operations/31339 "2019-11-21T08:20:03Z")
**Posts on this page:** 4
**Page:** 1

<div class="post-metadata">

### Author: ![natgeo-wong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/natgeo-wong/32/10965_2.png) [@natgeo-wong](https://discourse.julialang.org/u/natgeo-wong)
#### Post date: [November 21, 2019, 8:20am UTC](https://discourse.julialang.org/t/optimization-of-code-for-nan-operations/31339/1 "2019-11-21T08:20:03Z")

</div>

Hi, I’ve been trying to perform operations on Arrays that may have missing values, using the below prototype function:

```julia
function nanop(f::Function,data::Array;dim=1)

    ndim = ndims(data); dsize = size(data);
    if ndim==2
        if dim==1; newdim = [2,1]; data = permutedims(data,newdim); nsize = dsize[2]; end
    elseif ndim>2
        if dim==1; newdim = convert(Array,2:ndim); newdim = vcat(newdim,1);
            data = permutedims(data,newdim);
            nsize = dsize[2:end]
        elseif dim<ndim
            newdim1 = convert(Array,1:dim-1);
            newdim2 = convert(Array,dim+1:ndim);
            newdim = vcat(newdim1,newdim2,dim)
            data = permutedims(data,newdim);
            nsize = tuple(dsize[newdim1]...,dsize[newdim2]...)
        else nsize = dsize[1:end-1]
        end
    else; nsize = dsize[1:end-1]
    end

    data = reshape(data,:,dsize[dim]); l = size(data,1); out = zeros(l);

    for ii = 1 : l; dataii = data[ii,:];
        out[ii] = f(skipmissing(dataii));
    end

    return reshape(out,nsize)

end

```

However, benchmarking using @btime shows that something is seriously slowing down performance for large arrays.

```julia
@btime nanop(mean,rand(480,241,365),dim=3);;
@btime nanop(maximum,rand(480,241,365),dim=3);
@btime nanop(minimum,rand(480,241,365),dim=3);
@btime nanop(std,rand(480,241,365),dim=3);

```

gives me

```julia
  619.931 ms (231370 allocations: 663.69 MiB)
  618.972 ms (115690 allocations: 661.93 MiB)
  619.175 ms (115690 allocations: 661.93 MiB)
  829.410 ms (231370 allocations: 663.69 MiB)

```

Is there any way to reduce memory allocation, because frankly I’m going to be doing operations like

```julia
@btime nanop(mean,rand(480,241,365,4),dim=4);;
@btime nanop(maximum,rand(480,241,365,4),dim=4);
@btime nanop(minimum,rand(480,241,365,4),dim=4);
@btime nanop(std,rand(480,241,365,4),dim=4);

```

which frankly scares the crap out of me since I’ll have to perform the for-loop operation many many more times than in my first few test cases.

---

<div class="post-metadata">

### Author: ![oheil](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oheil/32/220745_2.png) [@oheil](https://discourse.julialang.org/u/oheil)
#### Post date: [November 21, 2019, 10:14am UTC](https://discourse.julialang.org/t/optimization-of-code-for-nan-operations/31339/2 "2019-11-21T10:14:28Z")

</div>

skipmissing brings a performance decrease of about factor 3-4 (edited, was 5 but 5 is too much):

```julia
julia> a=rand(480,241,365,4);

julia> @benchmark [x for x in a[:]]
BenchmarkTools.Trial:
  memory estimate: 2.52 GiB
  allocs estimate: 5
  --------------
  minimum time: 1.651 s (43.07% GC)
  median time: 1.775 s (47.18% GC)
  mean time: 1.783 s (47.36% GC)
  maximum time: 1.922 s (51.20% GC)
  --------------
  samples: 3
  evals/sample: 1

julia> @benchmark [x for x in skipmissing(a[:])]
BenchmarkTools.Trial:
  memory estimate: 2.72 GiB
  allocs estimate: 34
  --------------
  minimum time: 5.952 s (49.85% GC)
  median time: 5.952 s (49.85% GC)
  mean time: 5.952 s (49.85% GC)
  maximum time: 5.952 s (49.85% GC)
  --------------
  samples: 1
  evals/sample: 1

```

You can try to set the missing values to something which is neutral to the operation, e.g. case of mean  
you can do:

```julia
data[:]=[if ismissing(x) zero(x) else x end for x in data[:] ]
for ii = 1 : l;
        out[ii] = f(data[ii,:]);
end

```

(you don’t need to allocate the column into dataii, I didn’t check if it is really allocated, I guess not)

---

<div class="post-metadata">

### Author: ![natgeo-wong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/natgeo-wong/32/10965_2.png) [@natgeo-wong](https://discourse.julialang.org/u/natgeo-wong)
#### Post date: [November 21, 2019, 4:50pm UTC](https://discourse.julialang.org/t/optimization-of-code-for-nan-operations/31339/3 "2019-11-21T16:50:35Z")

</div>

Thanks so much for your help! I will have to look into this later. will put up some benchmarking later once I get the time.

---

<div class="post-metadata">

### Author: ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [November 21, 2019, 8:50pm UTC](https://discourse.julialang.org/t/optimization-of-code-for-nan-operations/31339/4 "2019-11-21T20:50:33Z")

</div>

> [@natgeo-wong](#):
>
> ```julia
> @btime nanop(mean,rand(480,241,365),dim=4);
> 
> ```

First of all, you are benchmarking also the creation of the array `rand(480,241,365)`, which takes time and allocates memory. I assume you already have have the array.

```julia
julia> @btime nanop(mean, rand(480,241,365); dim=3);
  517.498 ms (231370 allocations: 663.69 MiB)

julia> r = rand(480,241,365);

julia> @btime nanop(mean, $r, dim=3);
  214.726 ms (231368 allocations: 341.55 MiB)

```

If you mostly want to reduce memory use, you can do this instead:

```julia
mapslices(mean∘skipmissing, $r; dims=3);`

```

Benchmarking shows poorer timing, but far less memory use:

```julia
julia> @btime mapslices(mean∘skipmissing, $r; dims=3);
  355.505 ms (1388200 allocations: 34.43 MiB)

```