7
$\begingroup$

Suppose I have a table with some data that are 'lumped' around some central location, such as the dummy example

i0 = 22.3;
j0 = 34.1;
table = Table[
  Exp[-((i - i0)^2 + (j - j0)^2)/10^2]
  , {i, 1, 50}, {j, 1, 50}
  ]

which kind of looks like

ListPointPlot3D[table, PlotRange -> Full]

Mathematica graphics

I would like to recover the centroid of the "mass" determined by the signal in my table, i.e. in the case above, the tuplet {i0,j0}, in the cleanest way possible. I've managed to do it in the obvious but expensive way with a bunch of explicit sums and products with explicit Ranges, but it feels like there should be a built-in function that will do this ─ and I've as yet been unable to find it.

Can this be done? If so, how?

$\endgroup$
11
$\begingroup$
wd = WeightedData[Tuples @ Range @ Dimensions @ table, Join @@ table]
Mean @ wd

{22.3232, 33.9072}

Also

Total[MapIndexed[#2 # &, table / Total[table, 2], {2}], 2] (* and  *)
Dot[Join @@ table , Tuples[Range @ Dimensions @ table]] / Total[table, 2]

{22.3232, 33.9072}

$\endgroup$
8
$\begingroup$
i0 = 22.3; j0 = 34.1;
table = Table[
   Exp[-((i - i0)^2 + (j - j0)^2)/10^2], {i, 1, 50}, {j, 1, 50}];

x = Array[#1 &, Dimensions[table]];
y = Array[#2 &, Dimensions[table]];

pt = {Total@Flatten[x table], Total@Flatten[y table]}/
  Total[table, 2]
(* {22.3232, 33.9072} *)

ListDensityPlot[table, Epilog -> {Red, Point@Reverse[pt]}, 
 PlotRange -> All]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Thanks both. I accepted the 'proper' built-in as that's what the OP asks for, but this is faster for my instance and it is clean enough that it makes for readable code. $\endgroup$ – Emilio Pisanty Jul 11 '18 at 12:19
6
$\begingroup$

I can provide some improvement if it is about speed.

Let's generate a larger data set (I use Compile merely to speed it up a bit).

i0 = 22.3; j0 = 34.1;
m = 2500;
n = 1500;
x = Subdivide[1., 50, m - 1];
y = Subdivide[1., 50, n - 1];
table2 = Partition[#, n] &@
    Compile[{{X, _Real, 1}, {i0, _Real}, {j0, _Real}},
      Exp[-((X[[1]] - i0)^2 + (X[[2]] - j0)^2)/10^2],
      RuntimeAttributes -> {Listable}
      ][Tuples[{x, y}], i0, j0];

Szabolcs' approach

AbsoluteTiming[
 pt = {
    Total@Flatten[Transpose[ConstantArray[x, n]] table],
    Total@Flatten[ConstantArray[y, m] table]
    }/Total[table, 2]
 ]

{0.733425, {22.3288, 33.8733}}

The problem is that before summation, some large arrays have to be constructed and multiplied. But the summations and matrix-matrix products can also be expressed by cheaper matrix-vector products:

AbsoluteTiming[
 pt2 = With[{buffer = ConstantArray[1., m].table},
   {x.(table.ConstantArray[1., n]), 
     buffer.y
   }/(ConstantArray[1., n].buffer)
   ]
 ]

{0.044209, {22.3288, 33.8733}}

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.