# Unable to Efficiently Compute the Weighted Average, with weight $\overline{r}$, for points in a closed region?

In a previous post, I wanted to compute the weighted average, with weight $\overline{r}(x,y)$, of the whole cartesian plane. This was not possible since the integrals did not converge.

Instead I'm computing the weighted average of points inside Region $\mathcal{C}$.

The definition of $\overline{r}(x,y)$ is the average ray function

Suppose we have curve $\mathcal{C}$ and a ray, from $\mathbf{p}\in\mathcal{C}$, rotating full circle at $\Delta \theta$ increments while intersecting $\mathcal{C}$. We will denote the intersections as $\left\{\mathbf{q}_i\right\}_{i=1}^{n}$ (for $n$ intersections) and the length of the ray as $r_i = \|{\mathbf{q}_i}\mathbf - \mathbf p\|$.

Then as $\Delta \theta\to0$, the average radius function is $\overline{r}(\mathbf{p})=\left(\sum\limits_{i=1}^{n}r_i\right)/{n}$

For star shaped curves, one can use an alternate and simpler definition

$$\overline{r}(\mathbf{p})=\frac1{2\pi}\oint_{\mathbf {q}_1\in\mathcal C}\|\mathbf {q}_1-\mathbf p\|\,\mathrm d\theta.$$

Thanks to @Rahul, in this post, I have a code that can implement $\overline{r}(x,y)$. For example, if $\mathcal{C}$ is defined by $x^2+y^2+\sin(4x)+\sin(4y)=4$

curve = DiscretizeRegion[
ImplicitRegion[
x^2+y^2+Sin[4*x]+Sin[4*y] == 4, {{x, -3, 3}, {y, -4, 4}}], {{-3, 3},
{-4, 4}}, AccuracyGoal -> 8]
q = MeshCoordinates[curve];
edges = First /@ MeshCells[curve, 1];
signedAngle[a_, b_] := Arg[(Complex @@ a)/(Complex @@ b)]
1/(2 \[Pi]) Abs[Sum[Module[{q1, q2, r, d\[Theta]}, q1 = q[[First@e]];
q2 = q[[Last@e]];
r = EuclideanDistance[p, (q1 + q2)/2];(*midpoint approximation*)
d\[Theta] = signedAngle[q1 - p, q2 - p];
r d\[Theta]], {e, edges}]]


Edit: A faster code can be implemented using @HenrySchumacer's answer in the this post.

avgRadiusFunction[curve_MeshRegion] :=
Module[{q, e1, e2}, q = MeshCoordinates[curve];
{e1, e2} = Transpose[DeveloperToPackedArray[MeshCells[curve, 1][[All, 1]]]];
With[{q11 = q[[e1, 1]], q12 = q[[e1, 2]], q21 = q[[e2, 1]], q22 = q[[e2, 2]]},
Compile[{{p, _Real, 1}},
Module[{u11, u12, u21, u22, r, dθ, r2squared, r2, p1, p2},
p1 = Part[p, 1];
p2 = Part[p, 2];
u11 = q11 - p1;
u12 = q12 - p2;
u21 = q21 - p1;
u22 = q22 - p2;
r2squared = (u21^2 + u22^2);
dθ = ArcTan[(u11 u21 + u12 u22)/r2squared, (u12 u21 - u11 u22)/r2squared];
r2 = Sqrt[r2squared];
(RotateRight[r2] + r2).dθ/(4. Pi)
],
RuntimeAttributes -> {Listable},
Parallelization -> True]
]
];
curve = DiscretizeRegion[ ImplicitRegion[x^2+y^2+Sin[4*x]+Sin[4*y] == 4, {{x, -3, 3}, {y, -4, 4}}], {{-3, 3},{-4, 4}}, AccuracyGoal -> 8]


We can plot the contour of $\overline{r}(x,y)$ on $\mathcal{C}$.

k = ContourPlot[f[{x, y}], {x, -3, 3}, {y, -3, 3},
Exclusions -> None, Contours -> 100, PlotLegends -> Automatic]


Hence the weighted average in $\mathcal{C}$ is

$$\begin{array}{cc} {x_w=\frac{\int\int_{\mathcal{C}}x\overline{r}(x,y) \ dA }{\int\int_{\mathcal{C}}\overline{r}(x,y) \ dA}} & { y_w=\frac{\int\int_{\mathcal{C}} y\overline{r}(x,y) \ dA}{\int\int_{\mathcal{C}}\overline{r}(x,y) \ dA} }\end{array}$$

However, even for the simplest double integrals, I'm unable to approximate using NIntegrate In less than an hour.

x_w=(NIntegrate[x*avgRadius[{x,y}],{x,-10,10},



How do I do accurately and efficiently approximate the weighted average of the $x$ and $y$-coordinates over curve $\mathcal{C}$ with weights $\overline{r}(x,y)$?

"Accurately" means the error should be less than .00001. "Efficently" means computing the integral under $15$ minutes? I am unable to do so using NIntegrate

• Your main issue is that your function is too slow. Also you shouldn't keep calling avgRadius[{x, y}] over-and-over. Cache that result. I'd suggest replacing the EuclideanDistance with its appropriate symbolic definition, calculating symRad = avgRadius[{x, y}] and trying to Compile that result (using Compile[ ... , Evaluate@symRad, ...]). This will hopefully be fast enough to actually compute this integral. Commented Jan 13, 2018 at 18:11
• You should also point to this old post of yours. Commented Jan 13, 2018 at 18:25
• @b3m2a1 Why is there a singularity? Commented Jan 13, 2018 at 18:56
• @Arbuja never mind it's just a strong peak (which makes more sense). If you just do a Table and ListPlot3D Mathematica clips it. Then there's a different issue in the NIntegrate. All I know is it's unhappy with your function for one reason or another. Commented Jan 13, 2018 at 19:02
• Also your AccuracyGoal is probably way overkill. Try dropping it to 3 or 4. Commented Jan 13, 2018 at 19:04

With the function avgRadiusFunction from this post (hardened version!) and using the mass matrix from the finite element method (and coarsening the mesh considerably, this can be done as follows:

implicitcurve =  ImplicitRegion[ x^2 + y^2 + Sin[4*x] + Sin[4*y] == 4, {{x, -3, 3}, {y, -4, 4}}];
ρ = 0.01;
ρfine = 0.000001;
curvecoarse = DiscretizeRegion[implicitcurve, {{-3, 3}, {-4, 4}},  MaxCellMeasure -> ρ]
curvefine = DiscretizeRegion[implicitcurve, {{-3, 3}, {-4, 4}},  MaxCellMeasure -> ρfine];


As integration domain, I use simply the bounded domain enclosed by curvecoarse.

Needs["NDSolveFEM"]
domain = ToElementMesh[curvecoarse, "MeshOrder" -> 2, MaxCellMeasure -> ρ];
domain["Wireframe"]


This computes the mass matrix:

mass = Module[{vd, sd, cdata, x, y, u},
vd = NDSolveVariableData[{"DependentVariables","Space"} -> {{u}, {x, y}}];
sd = NDSolveSolutionData[{"Space"} -> {domain}];
cdata = InitializePDECoefficients[vd, sd, "MassCoefficients" -> {{1}}];
DiscretizePDE[cdata, InitializePDEMethodData[vd, sd], sd]["MassMatrix"]
];


And this computes the integral $\int_{\operatorname{domain}} \bar r(x,y) \, \operatorname{d} x \, \operatorname{d} y$

{x,y}=Transpose[domain["Coordinates"]];


21.4127

The $x_w$ and $y_w$ can be obtained like this:

xw = (x.mass.avgradii)/totalavgradius


-0.0689768

-0.0690877

Of course, you are free to decrease with the edge length parameter ρ which controls the accuracy. For example, with ρ = 0.001;ρfine = 0.000001; I get

xw = (x.mass.avgradii)/totalavgradius


-0.0690308

-0.0690308

• Nice. I clearly need to get comfortable with FEMs. Commented Jan 13, 2018 at 20:19
• How do I use this to solve $x_w$ and $y_w$? Do I take s=x*avgRadiusFunction[curve] and apply Total[mass.s[curvfine][domain["Coordinates"]]] Commented Jan 13, 2018 at 20:56

I mentioned this in the comments, but we'll tackle this via compilation and a more manageable AccuracyGoal:

curve =
DiscretizeRegion[
ImplicitRegion[
x^2 + y^2 + Sin[4*x] + Sin[4*y] ==
4, {{x, -3, 3}, {y, -4, 4}}], {{-3, 3}, {-4, 4}},
AccuracyGoal -> 4];
q = MeshCoordinates[curve];
edges = First /@ MeshCells[curve, 1];

1/(2 π) Abs[
Sum[Module[{q1, q2, d1, d2, r, dθ}, q1 = q[[First@e]];
q2 = q[[Last@e]];
r = Norm[p - (q1 + q2)/2];(*midpoint approximation*)

d1 = q1 - p;
d2 = q2 - p;
dθ =
Arg[(d1[[1]] + I*d1[[2]])/(d2[[1]] + I*d2[[2]])];
r dθ], {e, edges}]];
compSym =
Compile[{{x, _Real}, {y, _Real}},
RuntimeOptions -> {"EvaluateSymbolically" -> False},
CompilationTarget -> "C"
];


Then we can see how long it takes on a reasonably small grid:

AbsoluteTiming[
heights =
Table[
compSym[x, y],
{x, -10, 10, .1},
{y, -10, 10, .1}
];
][[1]]

3.49649


Which is still mad slow. So the NIntegrate will probably take a looong time.

And so first we check to see if we messed up and introduced a "MainEvaluate":

<< CompiledFunctionTools
pstr = CompilePrint[compSym];
StringContainsQ[pstr, "MainEvaluate"]

False


And it turns out we didn't.

Now plotting those heights:

ListPlot3D[heights, PlotRange -> All]


This seems like the kind of NIntegrate should be able to handle. So:

NIntegrate[
compSym[x, y],
{x, -10, 10},
{y, -10, 10}
]

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

142.029


And so there's something in your domain making NIntegrate unhappy.

The other integrals can also be handled via compilation, e.g.:

compXSym =
Compile[{{x, _Real}, {y, _Real}},
RuntimeOptions -> {"EvaluateSymbolically" -> False},
CompilationTarget -> "C"
];

compYSym =
Compile[{{x, _Real}, {y, _Real}},

But fundamentally NIntegrate is unhappy with something in your problem and it's worth trying to figure out what that is.