Message Boards Message Boards

0
|
4685 Views
|
2 Replies
|
1 Total Likes
View groups...
Share
Share this post:

Basic question on Functional Programming equivalent to For[..] loop

Posted 11 years ago

This is probably a simple question. I have a for loop that I use to sum the effects of charges at a test point. The test points are in one list and the charge locations are in another list. The code below illustrates a basic case.

Each item in chargeList contains the coordinate of the charge and the charge. chargeList[[j]] = {x,y,z,q}. Each item in testPointList contains the coordinates of the test point. testPointList[[i]]={x,y,z}.

For each test point, I need to sum chargeList[[j]][[4]] divided by the distance of each charge to the test point. While the procedural nested For[...] loop approach works, it can be very slow for large lists.

chargeList = Table[{-3, yPos, 1.54, 1}, {yPos, -3, 3, 0.5}]
testPointList = 
 Table[ {1.34 Cos[\[Theta]], 1.34 Sin[\[Theta]], 1.54}, {\[Theta], 0, 
   2 Pi, Pi/10}]

Clear[zField];
pointCount = Length[testPointList];
chargeCount = Length[chargeList];
zField = {};

For[i = 1, i <= pointCount, i++,

  zPot = 0;

  For[j = 1, j <= chargeCount, j++,
        zPot += 
     chargeList[[j]][[4]]/
      Sqrt[(testPointList[[i]][[1]] - 
           chargeList[[j]][[1]])^2 + (testPointList[[i]][[2]] - 
           chargeList[[j]][[2]])^2 + (testPointList[[i]][[3]] - 
           chargeList[[j]][[3]])^2]]

   AppendTo[
    zField, {testPointList[[i]][[1]], testPointList[[i]][[2]], 
     zPot}]];
POSTED BY: Doug Kimzey
2 Replies

We can look at this in pieces. Let's first examine your code:

For[i = 1, i <= pointCount, i++, zPot = 0;
  For[j = 1, j <= chargeCount, j++, 
    zPot += chargeList[[j]][[4]]/
      Sqrt[(testPointList[[i]][[1]] - 
           chargeList[[j]][[1]])^2 + (testPointList[[i]][[2]] - 
           chargeList[[j]][[2]])^2 + (testPointList[[i]][[3]] - 
           chargeList[[j]][[3]])^2]] AppendTo[
    zField, {testPointList[[i]][[1]], testPointList[[i]][[2]], 
     zPot}]];

The outer loop can be changed in a Table command:

zField1 = Table[zPot = 0;
   For[j = 1, j <= chargeCount, j++, 
    zPot += chargeList[[j]][[4]]/
      Sqrt[(testPointList[[i]][[1]] - 
           chargeList[[j]][[1]])^2 + (testPointList[[i]][[2]] - 
           chargeList[[j]][[2]])^2 + (testPointList[[i]][[3]] - 
           chargeList[[j]][[3]])^2]] ;
   {testPointList[[i]][[1]], testPointList[[i]][[2]], zPot}
   , {i, 1, pointCount}];

As you can see, zField1 is now being assigned as a total result rather than appended to. Then let's take one more step and remove the inner For loop and replace it with the use of Total and Map:

zField2 = Table[
   zPot = 
    Total[Map[#[[4]]/
        Sqrt[(testPointList[[i]][[1]] - #[[1]])^2 + \
(testPointList[[i]][[2]] - #[[2]])^2 + (testPointList[[i]][[3]] - \
#[[3]])^2] &, chargeList]];
   {testPointList[[i]][[1]], testPointList[[i]][[2]], zPot}
   , {i, 1, pointCount}];

Here, Map applies the operation you were doing item by item and Total sums the individual elements. Next, Mathematica has the Norm function built in and this can replace your Sqrt of the sum of the squares:

zField3 = Table[
   zPot = 
    Total[Map[#[[4]]/Norm[testPointList[[i]] - #[[1 ;; 3]]] &, 
      chargeList]];
   {testPointList[[i]][[1]], testPointList[[i]][[2]], zPot}
   , {i, 1, pointCount}];

I prefer not to index anything if I don't have to. Here, I've created a function, midFunc, which does some of the work, which I will use in a Map call:

zFieldPiece[testPoint_, charges_] := {testPoint[[1]], testPoint[[2]], 
  Total[Map[#[[4]]/Norm[testPoint - #[[1 ;; 3]]] &, chargeList]] }

zField4 = Map[zFieldPiece[#, chargeList] &, testPointList];

The speed up advantage has largely been gained by removing the For loops. The trouble with a For loop is that the code creates iterators that must be interpreted. Some of the changes I made, however, were strictly for readability. Norm is also optimized when you get to larger lists though I doubt you'd see much of a difference in this case. One other item you may notice is that I avoided assignment statements when at all possible. Assignment statements are extra steps and while this is again, not likely to be a major difference, can add up in very large data sets. You may not that I only used Set exactly once and SetDelayed exactly once in the final result.

POSTED BY: Jason Grigsby
Posted 11 years ago

Thank you for a working example program with sample data. That makes it so much easier to understand the question and propose and test a solution.

Since you want the result of "doing the same thing" to a list of elements, the first thought is to use Map with a function operating on an element.

Since part of your operation is totaling up a calculation over all items in a list, the second thought is to use another Map in Total[ Map[ ] ]

Since your calculation is roughly finding a distance between two points, the third thought is Norm of a difference.

Can you take those hints, study the documentation for each of those functions and see if you can find your solution?

I don't know if it is going to speed it up enough, but it removes a lot of your subscripting, which is often slow with big lists, and should reduce your code down to four lines, two creating your lists, one for a "helper function" to keep the two Map's from confusing their naming of # in the anonymous functions and one for the final Map. That does seem to reproduce exactly the same list as your zField.

If even this isn't fast enough then some very bright algorithm guys came up with a method that is far faster than n^2 for large clouds of interacting particles. That was perhaps twenty years ago and unfortunately I don't have any good keywords to search for to find that. But that would probably be a much more complicated and subtle method than just doing a summation over n^2 pairs of points.

POSTED BY: Bill Simpson
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract