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[1] * x[1]^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
```

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[1] * x[1]^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