Comprehension (generator) with statement IF and the number of true

[Martin Florek]

Hi all,

I'm new in Julia and I'm doing refactoring. I have the following function:

function mapeBase_v1(A::Vector{Float64}, F::Vector{Float64})
  s = 0.0
  count = 0
  for i in 1:length(A)
    if(A[i] != 0.0)
      s += abs( (A[i] - F[i]) / A[i])
      count += 1
    end
  end

  s, count 

end   

I'm looking for a simpler variant which is as follows:

function mapeBase_v2(A::Vector{Float64}, F::Vector{Float64})
# A - actual target values
# F - forecasts (model estimations)

  s = sumabs((x - y) / x for (x, y) in zip(A, F) if x != 0) # Generator

  count = length(A) # ???
  s, count
end

However with this variant can not determine the number of non-zero elements. I found option with length(A[A .!= 0.0]), but it has a large allocation. Please, someone knows a solution with generator, or variant v1 is very good choice?

Thanks in advance,
Martin

[Jeffrey Sarnoff]

This may do what you want.

function mapeBase_v3(actuals::Vector{Float64}, forecasts::Vector{Float64})
# actuals - actual target values
# forecasts - forecasts (model estimations)

  sum_reldiffs = sumabs((x - y) / x for (x, y) in zip(actuals, forecasts) if x != 0.0)  # Generator

  count_zeros = sum( map(x->(x==0.0), actuals) )
  count_nonzeros = length(actuals) - count_zeros
  sum_reldiffs, count_nonzeros
end

[Martin Florek]

Thnaks, It is true, but when I apply @benchmark v3 is 6 times slower as v1, also has a large allocation and I do not want it. For me speed is important and v3 is not without visual noise, too. Any more thoughts?

ben1 = @benchmark mapeBase_v1(a,f)
BenchmarkTools.Trial: 
  samples:          848
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%
  memory estimate:  32.00 bytes
  allocs estimate:  1
  minimum time:     4.35 ms (0.00% GC)
  median time:      5.87 ms (0.00% GC)
  mean time:        5.89 ms (0.00% GC)
  maximum time:     7.57 ms (0.00% GC)

ben2 = @benchmark mapeBase_v3(a,f)
BenchmarkTools.Trial: 
  samples:          145
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%
  memory estimate:  977.03 kb
  allocs estimate:  14
  minimum time:     32.69 ms (0.00% GC)
  median time:      33.91 ms (0.00% GC)
  mean time:        34.55 ms (0.10% GC)
  maximum time:     49.03 ms (3.25% GC)

[Jeffrey Sarnoff]

(you may like that visual noise if the code finds more than 1 year's use)

Which matters more to you saving time or saving space?

[Jeffrey Sarnoff]

As given above, the time is in sum_reldiffs, the space is in the count_zeros.

[Jeffrey Sarnoff]

this change saves some space: `count_zeros = length( find(x->x==0.0, actuals) )`

[Gregory Salvan]

maybe with reduce ?

```
function mapeBase_v4(A::Vector{Float64}, F::Vector{Float64})
  function abscount(prev, current)
    current[2] == 0.0 && return prev
    index, item = current
    (prev[1] + abs( 1 - F[index] / item), prev[2] + 1)
  end
  reduce(abscount, (0.0, 0), enumerate(A))
end

```

[DNF]

Could you explain why you are unhappy with v1? It seems like a near-optimal implementation to me. It is idiomatic Julia, clean, clear, and fast. Are you certain that you're not looking for a one-liner out of force of habit?

[Cedric St-Jean]

count_zeros = sum( map(x->(x==0.0), actuals) 

)

# can be rewritten

count_zeros = sum(x==0.0 for x in actuals)

This shouldn't allocate, because it's using a generator.

I agree with DNF that _v1 is already straight-forward, as far as performance-sensitive code goes.

[Martin Florek]

Thank you everyone. v1 is very nice, as it turned out. I was looking for the magic of language Julia especially for the generator. 

[Gregory Salvan]

Sorry reduce was a bad idea, even if syntax is nice, it's really slow and greedy.

V1 can take advantage of @inbounds and @simd optimizations: http://docs.julialang.org/en/release-0.5/manual/performance-tips/#performance-annotations

I hope reduce will be optimized in future because here it expresses rightly the problem we solve and has a nice syntax:

mapeBase_v4(A::Vector{Float64}, F::Vector{Float64}) = reduce((0.0, 0), enumerate(A)) do result::Tuple{Float64, Int64}, current::Tuple{Int64, Float64}
    current[2] == 0.0 && return result # guard clause
    result[1] + abs( 1 - F[current[1]] / current[2]), result[2] + 1
end

unfortunately it's unusable ;)

[DNF]

I don't have the impression that reduce is slow. The reduce function that you're using is complicated and may have features that preclude optimizations, such as vectorization.

But perhaps more importantly, the reduce version, while probably very clever, is almost completely impossible to understand. I know what it's supposed to do, and I still cannot decipher it, while the straight loop is crystal clear and easy to understand. And almost as concise!

[Steven G. Johnson]

If you want to spend your effort on making this code cleaner, the more Julian thing would be to focus on making it more type-generic, so that it can handle arguments of more types.  For example (requires Julia 0.5):

function mapeBase{T<:Number,S<:Number}(A::AbstractArray{T}, F::AbstractArray{S})
  indices(A) != indices(F) && error("arguments must be arrays of the same shape")

  s = zero(float(real(promote_type(T,S))))
  count = 0

  for i in eachindex(A)
    @inbounds if A[i] != 0
      s += abs((A[i] - F[i]) / A[i])

      count += 1
    end
  end

  return s, count 

end   

There is no performance penalty to declaring more generic function argument types.  When you pass in Vector{Float64} arrays, a version of mapeBase is specialized to that type and compiled, just as if you had declared those types.

Note, by the way, that this way of comparing two arrays is very susceptible to floating-point roundoff errors — think about what happens if A[i] is supposed to be 0.0, but is actually 1e-15 due to roundoff.    I would normally recommend something like vecnorm(A - F, 1) / vecnorm(A, 1) instead -- i.e. take the sum of the |A[i] - F[i]| and divide by the sum of |A[i]|.

I agree with other posters that this is a case in which a loop is much cleaner.  It is possible to solve this problem efficiently with reduce in 0.5 (where higher-order functions are now fast), but it would require much more careful coding and the resulting code would be much more obscure and unmaintainable, nor would a reduce-based implementation be much shorter as you have seen.

Higher-order functions are great, and it is great to use them where they simplify your code.  But sometimes a loop is cleaner, and (unlike in other dynamic languages) there is no reason not to use loops in Julia.

[Gregory Salvan]

2016-10-26 13:09 GMT+02:00 DNF <oyv...@gmail.com>:

I don't have the impression that reduce is slow. The reduce function that you're using is complicated and may have features that preclude optimizations, such as vectorization.

I don't know exactly why, the difference is bigger than what I was expected. for example with A and F = rand(100_000_000) I have:

mapeBase_v1
BenchmarkTools.Trial:
  samples:          26

  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%
  memory estimate:  32.00 bytes
  allocs estimate:  1

  minimum time:     194.51 ms (0.00% GC)
  median time:      197.77 ms (0.00% GC)
  mean time:        198.77 ms (0.00% GC)
  maximum time:     220.32 ms (0.00% GC)
mapeBase_v4
BenchmarkTools.Trial:
  samples:          1

  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%

  memory estimate:  37.25 gb
  allocs estimate:  1799998401
  minimum time:     58.77 s (3.04% GC)
  median time:      58.77 s (3.04% GC)
  mean time:        58.77 s (3.04% GC)
  maximum time:     58.77 s (3.04% GC)

If you have tips to help me understand this huge difference and eventually optimize (I've tried to look at llvm ir but it's quite hard for my level, I've just noted a lot of "store" but don't understand why it's done this way)

 

But perhaps more importantly, the reduce version, while probably very clever, is almost completely impossible to understand. I know what it's supposed to do, and I still cannot decipher it, while the straight loop is crystal clear and easy to understand. And almost as concise!

"crystal clear"... depends on your background and habits.
V4 is more natural for me than V1, probably because when I need a single result (sum, abs...) from a list of values my first option is always reduce, and I assimilate for loop with kind of "repeat something for each of these values"
Then I also avoid nested blocks and use a lot guard clauses at the point I read them as a part of function signature (by default I would naturally extract v1 for loop content to a function with a guard clause), so my eyes only catch the operation used to reduce the list in the form of the expected result.

This is not adapted to julia and I have to change these habits which is not a problem as julia is such a pleasure to use.

 

[DNF]

Actually, I see only a very marginal performance difference between your mapeBase_v4 (the first v4, don't know about the second v4) and the loop version, roughly 10%. Not sure why you're seeing a big difference.

As for the "crystal clear" thing: my background is almost 20 years of heavy Matlab use. Cryptic one-liners and contorted vectorizations are second nature to me. Quite often, reduce and generators make things easier and clearer, but not this time, even for me.

[Steven G. Johnson]

On Wednesday, October 26, 2016 at 3:29:39 PM UTC-4, DNF wrote:

Actually, I see only a very marginal performance difference between your mapeBase_v4 (the first v4, don't know about the second v4) and the loop version, roughly 10%. Not sure why you're seeing a big difference.

Perhaps he is using Julia 0.4 and you are using 0.5?  Higher-order functions are much faster in 0.5.

[DNF]

All higher-order functions? I thought it was mainly anonymous functions. Either way, that's a seriously big slowdown.

[Steven G. Johnson]

On Thursday, October 27, 2016 at 1:23:47 PM UTC-4, DNF wrote:

All higher-order functions? I thought it was mainly anonymous functions. Either way, that's a seriously big slowdown.

All higher-order functions can benefit in Julia 0.5, not just anonymous functions.  Because each function now has its own type, calling a higher-order function like reduce(...) can now compile a specialized version for each function you pass it, which allows it to do things like inlining. 

[Gregory Salvan]

Yes I'm using 0.4. Thank you for these infos.