# A Functional way to write map and then reduce?

#1

I am trying to implement GDQ method in Julia. Writing code in an imperative schema is a straightforward practice but I feel stuck in rewriting it with functional idioms below is the part of code which I am hanged in. Is it possible to rewrite it in functional schema in an efficient way?

Imperative code:

``````function coeff_1(points)
m1_ = m1(points)
a_ij(i, j) = m1_[i]/(points[i] - points[j])/m1_[j]

count = length(points)
result = zeros((count, count))

for i = 1:count
for j = 1:count
if (i != j)
result[i, j] = a_ij(i, j)
result[i, i] -= result[i, j]
end
end
end
result
end
``````

Incomplete functional functional implementation:

``````function coeff_1_functional(points)
m1_ = m1(points)
a_ij(i, j) = m1_[i]/(points[i] - points[j])/m1_[j]

count = length(points)

result = [a_ij(i, j) for i = 1:count, j = 1:count]
# the below part is where I stuck, It is both imperative and inefficient
for i = 1:count
result[i, i] = -sum(result[i, :])
end
result
end
``````

With Regards

#2

Or maybe this is what you are looking for.

``````result = [ i!=j?a_ij(i,j):2*a_ij(i,j) for i = 1:count, j = 1:count ] - diagm( [sum([a_ij(i,j) for j =1:count ]) for i = 1:count] )
``````

#3

This is just a straightforward expression and may not help with the efficiency issue. I won’t really comment on the efficiency issue.

``````function coeff_f(points)
m = m1(points)
ix = 1:length(points)

a(i,j) = m[i]/(points[i] - points[j])/m[j]
d(i) = -sum(a(i,j) for j = ix if i != j)

[ ((i == j) ? d(i) : a(i,j)) for i = ix, j = ix ]
end
``````

But might it help to define the auxiliary function for the diagonal entries? Put the condition only there and then compute the diagonal separately.

#4

By putting the diagonal condition in the auxiliary, I mean this: instead of NaN, make `a(i,j)` zero on the diagonal, so it can be used more freely.

``````function coeff_g(points)
m = m1(points)
ix = 1:length(points)

a(i,j) = (i == j) ? 0.0 : m[i]/(points[i] - points[j])/m[j]
d(i) = -sum(a(i,j) for j = ix)

[ a(i,j) for i = ix, j = ix ] + diagm([ d(i) for i = ix ])
end
``````

#5

Out of curiosity, what is the problem with the first version of the code?

#6

If you refer to my two versions, I think the second version is just a nicer expression. Making `a(i,j)` defined, as `0.0`, on the diagonal made it usable without the conditions.

#7

You can make transposed versions of m1_ and points, then broadcast:

``````function set_diagonal!(matrix, vector)
set_diagonal_at_index!(i) =
matrix[i, i] = vector[i]
set_diagonal_at_index!.(1:length(vector) )
matrix
end

function coeff_1_functional(points)
m1_ = m1(points)
result = m1 ./ ( points .- points') / m1'
column_sums = sum(a, 2)
set_diagonal!(result, column_sums)
end
``````

I haven’t tested it out, but I think it’s about right. Best in 0.6 which has lazy transposes

#8

I think there is something with diagonal elements.
Maybe the else part of ternary operator should be zero.

#9

It’s completely correct and efficient, just for a simple(functional) and efficient implementation.

#10

Nice implementation full of tricks.
Thank you.

#11

Yes, that’s so nice.
I just tested it, but it seems to be memory hungry.
Is it the trade of being functional, or there is a better way?

#12

The following should also do the work. It is both short and hopefully quick. Unfortunately, because `m_1` isn’t detailed in the question, it is hard to make sure everything is 100% correct. So, I’ve googled GDQ to look for the value of `m_1` and added a function to calculate it which might be cleaner than the current implementation you have.

This also defines the useful `setdiag!` function, but makes it both mutating and returning the mutated object, which plays nicely in composition with further processing.

The use of `inv.(m1)` is important because doing 1/x is slow on most computers and is better done only once per element of `m1`.

The code:

``````setdiag!(m,v) = begin foreach(i->m[i,i]=v[i],1:length(v)) ; m end

m_1(points) = prod(setdiag!(points.-points',ones(length(points))),2)

function coeff_1_functional(points)
m1 = m_1(points)
result = setdiag!(m1.*inv.(points.-points').*inv.(m1)',zeros(length(points)))
return setdiag!(result,-sum(result,2))
end
``````

#13

I had implemented `m_1` as this:

``````function m1(points)
function m1_map(row_point)
row_reduce(prod_, point) = ifelse(point != row_point, (row_point - point) * prod_, prod_)
reduce(row_reduce, 1, points)
end

map(m1_map, points)
end
``````

The most interesting thing about your vectorized implementation is that parts are putted together in a nice way, some code reuse but again it is not as efficient as imperative one. Maybe here the best way is to use C like code.

Thank you.

#14

Now with the `m1` implementation, I could check both functions give the same results (Yep!). Benchmarking shows my `m_1` is 10x faster. As for the `coeff_1`, the loop code in the original question 2x faster than my `coeff_1_functional` in Julia 0.5 and about the same on Julia 0.6 (which enjoys better optimization of vectorized operations).

PS benchmark used `points` vector 10 elements long.