5
$\begingroup$

The problem is relatively simple:

Given a 2D-binary array(grid), whereby every entry is either a 1 or -1, what is the absolute quickest way to take two coordinates points $o,p$ - which represent a line of distance $r_{op}$ - and move it through all possible valid grid points.

Let $c_1 = f(o) , c_2 = f(p)$, these return the value of the array at the corresponding points of $o,p$. After every move of the points, you want calculate the following:

$C = c_1\times c_2$

then take the average over the entire grid of $C$.

I have been using a nested Do-loop for my trial (see example code below), but it's way to slow for the size of grids i'm looking at and the number of different line lengths I am sampling across the grid, I thus am looking for a super speedy approach that I am probably unaware of.

    GetC[coords_, ArrayData_] := Module[{c1, c2},
      c1 = ArrayData[[coords[[1, 1]], coords[[1, 2]]]];
      c2 = ArrayData[[coords[[2, 1]], coords[[2, 2]]]];
     c1*c2];

    ftran[{{x1_, y1_}, {x2_, y2_}}, {xT_, yT_}] :=
    (
       Return[{{x1 + xT, y1 + yT}, {x2 + xT, y2 + yT}}];
    );


    grid = {1000, 1000};
    array = Table[(-1)^RandomInteger[{i, j}], {i, grid [[1]]}, {j,grid [[2]]}];

    o = {1, 1};
    p = {1, 2};

    line = {o, p}

    bounds = CoordinateBounds[line];
    xMax = grid[[1]] - bounds[[1, 2]];
    yMax = grid[[1]] - bounds[[2, 2]];
    cnt = 0
    c = 0;
    Do[
     Do[
        newLine = ftran[line, {i, j}];
        TempC = GetC[newLine, array];
        c = c + TempC;
        cnt++ 1;
       , {i, 0, xMax}];
      , {j, 0, yMax}]; // AbsoluteTiming

    averageC = c/cnt

For a 1000x1000 grid I get a timing of 17.5053s (of course this will change depending on your PC performance). To put things into perspective, my grids are roughly 2000x2000, and I'm looking at doing this procedure for many different lines ranging from $r_{op}=1$ to $r_{op} = 1000$. To be exact, this is 216341 possible different lengths. As you can see, if I cannot find a very quick way of rewriting my code, this task isn't really possible.

Using a time of 17s, this wold take ~ 42 days on my PC. Definitely not good enough.

Improve this question
$\endgroup$
4
  • $\begingroup$ I do not understand what you mean by the following: "what is the absolute quickest way to take two coordinates points o,p and move it through all possible valid grid points". What does "it" refer to here? What is the ultimate goal of this procedure, i.e. what are you going to do with the $c$ values? $\endgroup$
    – MarcoB
    Feb 9, 2017 at 15:35
  • $\begingroup$ @MarcoB Apologies, I probably should have been more explicit in my problem statement. I want to calculate the correlation between all possible pairs of points on a grid, and then see how the average correlation depends on the distance between the two points. Now, to do this, what you are essentially doing is moving a line(determined by points $(P,Q)$ through all possible pairs on the grid - hence the question as stated- and at each pair calculate the correlation between points, I.e return a value of 1 if $c_1=c_2$ or -1 if $c_1 \ne c_2$. I hope this clarifies things? $\endgroup$
    – Luca Pontiggia
    Feb 9, 2017 at 16:10
  • 1
    $\begingroup$ Have you looked at ListCorrelate? This is bound to be faster... $\endgroup$
    – bill s
    Feb 9, 2017 at 16:37
  • $\begingroup$ You can change the default functions in ListCorrelate -- the operations are defined by the g and h functions in the bottom two definitions in the help file. $\endgroup$
    – bill s
    Feb 9, 2017 at 17:24

3 Answers 3

Reset to default
10
$\begingroup$

Assuming your matrix consists of -1, 1 only, then the product of two entries, x and y, satisfies:

x y==2 Unitize[x+y]-1

So, we can use ListCorrelate and post process (using default parameters where possible should be faster than specifying g and h functions as mentioned in the comments) . Here is a function that does this:

corr[line_, array_]:=Module[{ker=Normal@SparseArray[line->1]},
    Mean @ Flatten @ If[Length[ker]==1,
        2Unitize@ListCorrelate[ker,array,{1,-1}]-1,
        2Unitize@ListCorrelate[ker,array,{1,-1},0,ListCorrelate,Plus,1]-1
    ]
]

Test (note that I simplified your random array generator to use RandomChoice):

SeedRandom[1];
grid={1000,1000};
array=RandomChoice[{-1,1},grid];
o={1,1};
p={3,1};

line={o,p}

bounds=CoordinateBounds[line];
xMax=grid[[1]]-bounds[[1,2]];
yMax=grid[[1]]-bounds[[2,2]];
cnt=0
c=0;
Do[Do[newLine=ftran[line,{i,j}];
TempC=GetC[newLine,array];
c=c+TempC;
cnt++ 1;,{i,0,xMax}];,{j,0,yMax}];//AbsoluteTiming

averageC=c/cnt
corr[line,array]//AbsoluteTiming

{{1,1},{3,1}}

0

{10.8666,Null}

-(699/499000)

{0.060691,-(699/499000)}

A speedup of almost 200 times.

Now that there are several answers, let's do a speed comparison. I rewrote the functions to use similar syntax:

corr[p_, q_, array_] := Module[{ker = Normal@SparseArray[{{1,1},{p,q}}->1]},
    Mean @ Flatten @ If[Length[ker]==1,
        2 Unitize @ ListCorrelate[ker, array, {1,-1}]-1,
        2 Unitize @ ListCorrelate[ker, array, {1,-1}, 0, ListCorrelate, Plus, 1]-1
    ]
]

george[p_, q_, array_]:= Mean @ Flatten @ (array[[ ;; -p, ;; -q ]] * array[[ p ;; , q ;; ]])

averageCompiled = Compile[{{x, _Integer}, {y, _Integer}, {grid, _Real, 2}}, 
    Module[{dimx, dimy},
        {dimx, dimy} = Dimensions[grid];
        Mean@Flatten@Table[
            grid[[i, j]] grid[[i + x - 1, j + y - 1]],
            {i, dimx - x + 1},
            {j, dimy - y + 1}
        ]
    ],
    CompilationTarget -> "C", RuntimeOptions -> "Speed"
];

Here is the data:

SeedRandom[1];
grid = {1000,1000};
array = Developer`ToPackedArray @ RandomChoice[{-1,1},grid];

And here is the comparision:

r1 = corr[3, 1, array]; //AbsoluteTiming
r2 = george[3, 1, array]; //AbsoluteTiming
r3 = averageCompiled[3, 1, array]; //AbsoluteTiming
N@r1 === N@r2 === r3

{0.065301,Null}

{0.007101,Null}

{0.01775,Null}

True

Note that the compiled version returns a real number, while the other two approaches return a rational.

The previous versions didn't work when the second point had a negative coordinate. Here's a variation that uses the difference in coordinates instead:

part[da_,db_,array_] := Mean @ Flatten @ (
    array[[ span[da], span[db]]] * array[[span[-da], span[-db]]]
)

span[x_]:=Switch[Sign[x],
    -1, 1 ;; x-1,
    0, 1 ;; -1,
    1, x+1 ;; -1
]

Compare the original corr and part:

corr[{{3,1}, {1,2}}, array] //AbsoluteTiming
part[-2, 1, array] //AbsoluteTiming

{0.076272,-(266/498501)}

{0.007252,-(266/498501)}

Improve this answer
$\endgroup$
4
  • $\begingroup$ I am gobsmacked!! This is amazing, I will try figure out what it is exactly you are doing. I only hope that one day I can utilize such methods in my code! $\endgroup$
    – Luca Pontiggia
    Feb 9, 2017 at 19:51
  • $\begingroup$ The second approach of @george2079 should be much faster than using ListCorrelate. $\endgroup$
    – Carl Woll
    Feb 9, 2017 at 20:06
  • $\begingroup$ The above methods don't seem to work for certain cases. Consider the line extended from point $o=(1,1)$ to $p=(-1,2)$. Since we don't want to deal with negative values in the array, we translate the line by $x' \rightarrow x+2$. This gives the new line from $o'=(3,1)$ to $p'=(1,2)$. Your original "corr[line_, array_]" could handle these cases, the latter (r1,r2,r3) cannot; I guess this is because the implicit assumption that all lines are given relative to $(1.1)$ is not withheld, but this would then exclude many possible cases. Eg. consider all $y>0$ integer solutions for $ x^2+y^2=625$. $\endgroup$
    – Luca Pontiggia
    Feb 10, 2017 at 5:59
  • $\begingroup$ That's some nifty code! One last comment to round all this off; can one extend this to multiple points? I expanded this in more detail in @george2079 answer below (excuse the taging mistake in that comment, I can't seem to edit it anymore). I am primarily interested in looking at three points, but would be curious if such code could be generalized to $n$ points. $\endgroup$
    – Luca Pontiggia
    Feb 10, 2017 at 8:06
4
$\begingroup$

just tightening up your code we get about a factor of two speed improvement..

grid = {100, 100};
o = {1, 1};
p = {1, 2};
line = {o, p}
bounds = CoordinateBounds[line];
xMax = grid[[1]] - bounds[[1, 2]];
yMax = grid[[1]] - bounds[[2, 2]];
c = Sum[
    (Times @@ (Extract[array, (# + {i, j}) & /@ line])), {i, 0, 
     xMax}, {j, 0, yMax}]; // AbsoluteTiming
averageC = c/((xMax + 1) (yMax + 1))

another approach

make two copies of the array zero padded to affect the shift and multiply the full matrices (this yields about a 100x speedup)

shiftx = p[[2]] - o[[2]];
shifty = p[[1]] - o[[1]];
(Flatten[(ArrayPad[array, {{shifty, 0}, {0, shiftx}}]) (ArrayPad[
        array, {{shifty, 0}, {shiftx, 0}}])] // 
    Total)/((xMax + 1) (yMax + 1))

To make this work you will need to work out the appropriate padding depending on the direction of the line.. (hope that makes sense).

another using Part (this should be robust in terms of handling any line direction)

lineproduct[array_, line_] := 
 Module[{p = If[line[[2]] < 0, -line, line]},
     If[p[[1]] >= 0,
      (array[[#[[2]] ;;, ;; -#[[1]]]] array[[;; -#[[2]], #[[1]] ;;]]) 
        &@{p[[2]] + 1, p[[1]] + 1},
      (array[[;; -#[[2]], ;; -#[[1]]]] array[[#[[2]] ;;, #[[1]] ;;]])
        &@{p[[2]] + 1, -p[[1]] + 1}
      ] // Flatten // Mean]
lineproduct[array, p - o]

or as per comment, lineproduct[Developer`ToPackedArray@array,p-o]

note if I understand the problem correctly the result should depend only on the difference p-o. this produces different results from the op if the first point 'o' is not {1,1}.

Improve this answer
$\endgroup$
4
  • $\begingroup$ You should clean this up. I would use Part instead of ArrayPad, something like: Mean@Flatten@(array[[;; -p, ;; -q]] * array[[p ;; , q ;;]]) This will be quite a bit faster than my ListCorrelate approach, as long as you convert the array into a packed array first. $\endgroup$
    – Carl Woll
    Feb 9, 2017 at 20:08
  • $\begingroup$ @CarlWoll good point, I worked that up. $\endgroup$
    – george2079
    Feb 9, 2017 at 21:25
  • $\begingroup$ Very clever approach! Would one be able to extend this to say three points? So let's say, given three points, $o,p,q$, move the now triangle around the grid and calculate the correlation given by $c_1*c_2*c_3$. I suppose you would do a padding in the direction of $o\rightarrow p$, then a separate padding from $o \rightarrow q$ and then I guess do a dot product on the resulting two flattened arrays? In terms of @george2079 cleanup suggestion, I would still need to think how I could implement the addition of a third point. $\endgroup$
    – Luca Pontiggia
    Feb 10, 2017 at 3:58
  • $\begingroup$ On second thought, this is not technically correct, as what I would be calculating is $(c_1\times c_2)(c_1\times c_3) \neq c_1\times c_2\times c_3$. $\endgroup$
    – Luca Pontiggia
    Feb 10, 2017 at 4:47
3
$\begingroup$

Some really good methods have already been discussed. Here is another way to write your algorithm, which is much shorter than what you wrote:

average[{x_, y_}, grid_] := Module[{sum, dimx, dimy},
  {dimx, dimy} = Dimensions[grid];
  N@Mean@Flatten@Table[
     grid[[i, j]] grid[[i + x - 1, j + y - 1]],
     {i, dimx - x + 1},
     {j, dimy - y + 1}
     ]
  ]

You can use it like this:

p = {3, 1};
SeedRandom[1];
array = RandomChoice[{-1, 1}, {1000, 1000}];
average[p, array] // AbsoluteTiming

{3.09915, -0.0014008}

I don't think o has to be defined at all. We can always assume that it is at position {1,1}.

Writing it in this way gives a speedup of about 8x. Now that we have it in this form we can easily see that it is compilable, and we compile it as C code:

averageCompiled = Compile[{{x, _Integer}, {y, _Integer}, {grid, _Real, 2}}, 
  Module[{dimx, dimy},
   {dimx, dimy} = Dimensions[grid];
   Mean@Flatten@Table[
      grid[[i, j]] grid[[i + x - 1, j + y - 1]],
      {i, dimx - x + 1},
      {j, dimy - y + 1}
      ]
   ], CompilationTarget -> "C", RuntimeOptions -> "Speed"]

If we test this, we find that we get a speedup of approximately 1000x compared with the original code:

averageCompiled[3, 1, array] // AbsoluteTiming

{0.025988, -0.0014008}

And we got the speed up without almost any effort at all. The only thing we did was to clean up your code and make it simpler, and then compile it.

It can also be written with Do in "C style", even more reminiscent of what you wrote:

averageCompiled = Compile[{{x, _Integer}, {y, _Integer}, {grid, _Real, 2}}, 
  Module[{sum, dimx, dimy},
   {dimx, dimy} = Dimensions[grid];
   sum = 0.;
   Do[
    sum += grid[[i, j]] grid[[i + x - 1, j + y - 1]],
    {i, dimx - x + 1},
    {j, dimy - y + 1}
    ];
   sum/((dimx - x + 1) (dimy - y + 1))
   ], CompilationTarget -> "C", RuntimeOptions -> "Speed"]

This is a tiny bit faster than the other compiled version, but the difference in negligible.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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