Using parfor in Matlab
We all know that loops don’t behave well in Matlab. Whenever it is possible to vectorize the code (i.e. use vectors and matrices to do simultaneous operations, instead of one at a time) significant speedup is possible. However, there are complex tasks which cannot be vectorized and loops cannot be avoided. I recently needed to compute eigenvalues for some 10 million domains. Since the computations are independent, they could be run in parallel. Fortunately Matlab offers a simple way to do this, using parfor.
There are some basic rules one need to respect to use parfor efficiently:
 Don’t use parfor if vectorization is possible. If the task is not vectorizable and computations are independent, then parfor is the way to go.
 Variables used in each computation should not overlap between processors. This is for obvious reasons: if two processors try to change the same variable using different values, the computations will be meaningless in the end.
 You can use an array or cell to store the results given by each processor, with the restriction that processors should work on disjoint parts of the array, so there is no overlap.
The most restrictive requirement is the fact that one cannot use the same variables in the computations for different processors. In order to do this, the simplest way I found was to use a function for the body of the loop. When using a matlab function, all variables are local, so when running the same function in parallel, the variables won’t overlap, since they are local to each function.
So instead of doing something like
parfor i = 1:N commands ... array(i) = result end
you can do the following:
parfor i=1:N array(i) = func(i); end function res = func(i) commands...
This should work very well and no conflict between variables will appear. Make sure to initialize the array before running the parfor, a classical Matlab speedup trick: array = zeros(1,N). Of course, you could have multiple outputs and the output array could be a matrix.
There is another trick to remember if the parpool cannot initialize. It seems that the parallel cluster doesn’t like all the things present in the path sometimes. Before running parfor try the commands
c = parcluster('local'); c.parpool
If you recieve an error, then run
restoredefaultpath c = parcluster('local'); c.parpool
and add to path just the right folders for your code to work.
Plotting 3D level sets in Paraview
Surfaces can be represented as certain levelsets for some 3D functions. Given a set of points with values attached to them a level set associated to a certain number will separate the points into two sets: with values higher or lower than the number considered. It is nice to have good looking plots when working with the levelset method in shape optimization. Paraview is a nice, open source framework, which has the right tools in order to produce high quality plots.
I’ll briefly describe below how to use Paraview to make some nice pictures of levelsets. First of all, you’ll need your data in some format that Paraview can understand. I use vtk file format for which there is a nice automated interface in the software I use for the optimization (FreeFem++). In the vtk file you need to have a set of points and a scalar value attached to them.
If you want to create levelsets associated to certain values, follow the steps below:

 Load your file (containing points with at least a scalar value) and click Apply.
 Next, go to Filters/Alphabetical and select “Cell Data to Point Data” (if you forget to do this, you’ll get a rough surface where you see the discretization). Click Apply.
 Then apply the Contour filter (by clicking the button or going in Filters/Alphabetical). You’ll have to select the field for which the contours will be drawn and then put in the values of the levelsets you want to see. Click Apply. An example can be seen below.
If your levelset cuts the boundary of the domain, Paraview will draw a hole there. If you want to have a closed region instead, you need to use the IsoVolume filter instead of the Contour one. The difference is that you need to specify two values and Paraview will draw the surface enclosing the points corresponding to these values. Many other features are directly available: you can color the level set following another scalar value, you can set the lighting, etc. You can also symmetrize your geometry using the Reflect filter. Below you can see a result built from my work.
You can also create animations in a pretty straightforward way. Just go to View and select the Animation box. Then you’ll see the animation options. Add a Camera object with Orbit field selected. You’ll be presented with multiple options, like the center of rotation, direction of the vertical and initial position. Once everything is set, click the play button to see the animation. Then go to File/Save Animation to save it to a file.
I heard that Paraview could to many things when dealing with visualization aspects, but I hesitated to use it until now since the interface is not straightforward. The use of Filters is not clear in the beginning, but after playing with some examples, everything becomes really easy to use. The next step is to automatize all this using scripts.
Happy New Year!
A hint for Project Euler Pb 613
The text for Problem 613 can be found here. The hint is the following picture 🙂
Project Euler – Problem 264
Today I managed to solve problem 264 from Project Euler. This is my highest rating problem until now: 85%. You can click the link for the full text of the problem. The main idea is to find all triangles ABC with vertices having integer coordinates such that
 the circumcenter O of each of the triangles is the origin
 the orthocenter H (the intersection of the heights) is the point of coordinates (0,5)
 the perimeter is lower than a certain bound
I will not give detailed advice or codes. You can already find a program online for this problem (I won’t tell you where) and it can serve to verify the final code, before going for the final result. Anyway, following the hints below may help you get to a solution.
The initial idea has to do with a geometric relation linking the points A, B, C, O and H. Anyone who did some problems with vectors and triangles should have come across the needed relation at some time. If not, just search for important lines in triangles, especially the line passing through O and H (and another important point).
Once you find this vectorial relation, it is possible to translate it in coordinates. The fact that points A, B, C are on a circle centered in O shows that their coordinates satisfy an equation of the form , where is a positive integer, not necessarily a square… It is possible to enumerate all solutions to the following equation for fixed , simply by looping over and . This helps you find all lattice points on the circle of radius .
Once these lattice points are found one needs to check the orthocenter condition. The relations are pretty simple and in the end we have two conditions to check for the sum of the x and y coordinates. The testing procedure is a triple loop. We initially have a list of points on a circle, from the previous step. We loop over them such that we dont count triangles twice: i from 1 to m, j from i+1 to m, k from j+1 to m, etc. Once a suitable solution is found, we compute the perimeter using the classical distance formula between two points given in coordinates. Once the perimeter is computed we add it to the total.
Since the triple loop has cubic complexity, one could turn it in a double loop. Loop over pairs and construct the third point using the orthocenter condition. Then just check if the point is also on the circle. I didn’t manage to make this double loop without overcounting things, so I use it as a test: use double loops to check every family of points on a given circle. If you find something then use a triple loop to count it properly. It turns out that cases where the triple loop is needed are quite rare.
So now you have the ingredients to check if on a circle of given radius there are triangles with the desired properties. Now we just iterate over the square of the radius. The problem is to find the proper upper bound for this radius in order to get all the triangles with perimeter below the bound. It turns out that a simple observation can get you close to a near optimal bound. Since in the end the radii get really large and the size of the triangles gets really large, the segment OH becomes small, being of fixed length 5. When OH is very small, the triangle is almost equilateral. Just use the upper bound for the radius for an equilateral triangle of perimeter equal to the upper bound of 100000 given in the problem.
Using these ideas you can build a bruteforce algorithm. Plotting the values of the radii which give valid triangles will help you find that you only need to loop over a small part of the radii values. Factoring these values will help you reduce even more the search space. I managed to solve the problem in about 5 hours in Pari GP. This means things could be improved. However, having an algorithm which can give the result in “reasonable” time is fine by me.
I hope this will help you get towards the result.
Project Euler 607
If you like solving Project Euler problems you should try Problem number 607. It’s not very hard, as it can be reduced to a small optimization problem. The idea is to find a path which minimizes time, knowing that certain regions correspond to different speeds. A precise statement of the result can be found on the official page. Here’s an image of the path which realizes the shortest time:
Project Euler tips
A few years back I started working on Project Euler problems mainly because it was fun from a mathematical point of view, but also to improve my programming skills. After solving about 120 problems I seem to have hit a wall, because the numbers involved in some of the problems were just too big for my simple bruteforce algorithms.
Recently, I decided to try and see if I can do some more of these problems. I cannot say that I’ve acquired some new techniques between 20122016 concerning the mathematics involved in these problems. My research topics are usually quite different and my day to day programming routines are more about constructing new stuff which works fast enough than optimizing actual code. Nevertheless, I have some experience coding in Matlab, and I realized that nested loops are to be avoided. Vectorizing the code can speed up things 100 fold.
So the point of Project Euler tasks is making things go well for large numbers. Normally all problems are tested and should run within a minute on a regular machine. This brings us to choosing the right algorithms, the right simplifications and finding the complexity of the algorithms involved.