I’m in need of a good mixing rule function. A mixing rule is a type of of weighed mean. if we have a property p and some weights x, a mixing rule is like this:

`p_mixing = sum(op(pi,pj)*xi*xj for i,j = 1:N)`

the op(pi,pj) is a type of mean operation, where op(pi,pi)= pi (for example, an aritmethic mean)

using the properties of symmetry and that fact, i made a function that does such mixing rule

``````function mixing_rule(op,x,p)
N = length(x)
res1 = p[N]*x[N]^2
for i = 1:N-1
res1 += p[i]*x[i]^2
for j=(i+1):N
res1 += 2*x[i]*x[j]*op(p[i],p[j])
end
end
return res1
end
``````

can it be faster?, i’m not convinced of my methods, or, there exists some function in some package that does this? sounds like it can be implemented by a mapreduce or somethig

This takes less than 1/6th the time on my machine:

``````function mixing_rule_2(op,x,p)
N = length(x)
@boundscheck length(p) == N || throw(DimensionMismatch("Dimensions of x and p don't match"))
@inbounds begin
res1 = zero(p * x^2)
for i = 1 : N
res1 += p[i] * x[i]^2
@simd for j = 1 : i - 1
res1 += 2 * x[i] * x[j] * op(p[i], p[j])
end
end
end
return res1
end
``````
• use `@inbounds` and `@simd` to elide bounds checks and get the inner loop to vectorize

• iterate over `j = 1 : i - 1` instead of `j=(i+1):N` in the inner loop for better cache locality

You can probably do quite a bit better still (e.g. with threading), but this is a nice quick win.

Full code including benchmark setup.
``````function mixing_rule(op,x,p)
N = length(x)
res1 = p[N]*x[N]^2
for i = 1:N-1
res1 += p[i]*x[i]^2
for j=(i+1):N
res1 += 2*x[i]*x[j]*op(p[i],p[j])
end
end
return res1
end

function mixing_rule_2(op,x,p)
N = length(x)
@boundscheck length(p) == N || throw(DimensionMismatch("Dimensions of x and p don't match"))
@inbounds begin
res1 = zero(p * x^2)
for i = 1 : N
res1 += p[i] * x[i]^2
@simd for j = 1 : i - 1
res1 += 2 * x[i] * x[j] * op(p[i], p[j])
end
end
end
return res1
end

using BenchmarkTools
N = 10_000
x = rand(N)
p = rand(N)
op = (x, y) -> (x + y) / 2

@assert mixing_rule(op, x, p) ≈ mixing_rule_2(op, x, p) atol=1e-5
@btime mixing_rule_2(\$op, \$x, \$p)

@btime mixing_rule(\$op, \$x, \$p)
``````
4 Likes

Thanks!, it worked flawlessly, to optimize further, i found this:

``````@boundscheck checkbounds(p,N)
``````

and it shaves some aditional time, but nothing compared to your original function 