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