Posts Tagged ‘programming’

Using parfor in Matlab

February 27, 2018 Leave a comment

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 speed-up 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:

  1. 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.
  2. 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.
  3. 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

you can do the following:

parfor i=1:N
   array(i) = func(i);

function res = func(i)

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');

If you recieve an error, then run

c = parcluster('local');

and add to path just the right folders for your code to work.


Metapost – or how to code an image

November 9, 2017 Leave a comment

What software do you use when you need to draw a nice image? I often need to draw things related to research and I never managed to efficiently use a software which uses the mouse or touchpad to draw and modify things. Moreover, if you need to add some mathematical text to the figure things get even more complicated. For me it is more natural to use code to generate graphics, if this is possible. Using something like Matlab to draw is possible. The advantage is that once you have a working code which produces what you want, you can easily modify it. If you have an image where you need to repeat things, using loops can facilitate the job (programmers will understand…).

This is were Metapost comes into play. I found this a long time ago while searching for a tool to build nice graphics for my master thesis. Being used to LaTeX for typesetting math, it was a natural way to draw using code. I do not claim it is the best or the simplest way, but for me it works. Most importantly, it allows me to build high quality, vectorized graphics, which are memory efficient. The advantage of vector graphics is that you can zoom as much as you want and you’ll never see the pixels.

There are a bunch of places where you can learn Metapost. I’ll put here some references I use a lot. The easiest way to start drawing in Metapost is to take a look at some existing codes and modify them. You can find lots of examples in this tutorial. If you have a linux system, you can compile metapost .mp files to obtain high quality pdfs using the command mptopdf. If not, then you can use the online Metapost editor found here. In any case, here are some of my codes for some drawings in Metapost.

Here is the code for my website logo: logo If you want to see the lossless vectorized pdf click to see the following file: logo-0    You can see that it has a border which changes from one color to another. This was done using a loop and a parameter to vary the color as we go along the boundary.

% copy from here to use the online previewer
u:=25; % 25 = 25bp = 25 PostScript points = 25/72 in
wi:=10; % width in units u
he:=7; % height in units u
hoehe:=he*u; % height
%for i=0 upto he:
%draw (0, i*u)--(breite, i*u) withcolor .7white;
%for j=0 upto wi:
%draw (j*u, 0)--(j*u, hoehe) withcolor .7white;
path p,q;
pickup pensquare scaled 20;
draw p;
numeric c,d,detail;
color a,b,co;
for i=1 upto (detail/2):
q := subpath (h*(i-1)/detail,h*i/detail) of p;
draw q withcolor (2*i/detail)[a,b];
for i=detail downto (1+detail/2):
q := subpath (h*(i-1)/detail,h*i/detail) of p;
draw q withcolor (2*(detail-i)/detail)[a,b];
label ("B",(1/15)*(2.6u,2.7u)) scaled 15 withcolor co ;
label ("2" infont defaultfont scaled 8,(4.9u,4.1u))
withcolor co;
% end copy here for online previewer

Here’s another example of a figure where I needed a grid of disks aligned on top of another figure. Using loops in Metapost allowed me to get what I wanted: cm

For the high quality pdf click here: cioranescu-murat-0  and see the code below

u:=25; % 25 = 25bp = 25 PostScript points = 25/72 in
wi:=10; % width in units u
he:=7; % height in units u
hoehe:=he*u; % height

path p,pa;
draw p;
fill p withcolor .8white;
pair a,b;
numeric c,d;
a:= point c of p;
b:= point d of p;
pickup pencircle scaled 1.5;
%draw a;
%draw b;
path q,r,s;
q:= subpath (c,d) of p;
%draw q withcolor red;
%r:= b..((a+b)/2 shifted (.5(a-b) rotated -90))..a;
%draw r;
%fill buildcycle(r,q) withcolor .5red;
%fill pa withcolor .9999blue;
%label.rt(btex $\Omega$ etex,.5(u,3u))  scaled 2;

draw (-u,9u)-- (10u,9u)--(10u,-2u)--(-u,-2u)--cycle;
for i=0 upto 9:
for j=-1 upto 8:
    draw fullcircle scaled 0.3u shifted (i*u, j*u);
    fill fullcircle scaled 0.3u shifted (i*u, j*u) withcolor white;

There are many ways today to draw what you want. Metapost is one of the tools available, if you like coding.

Project Euler – Problem 264

July 28, 2017 Leave a comment

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 x^2+y^2=n, where n is a positive integer, not necessarily a square… It is possible to enumerate all solutions to the following equation for fixed n, simply by looping over x and y. This helps you find all lattice points on the circle of radius \sqrt{n}.

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

June 11, 2017 7 comments

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:



FreeFem to Matlab – fast mesh import

October 14, 2016 Leave a comment

I recently wrote a brief introduction to FreeFem++ in this post. FreeFem is a software designed for the numerical study of partial differential equations. It has the advantage of being able to easily define the geometry of the domain, construct and modify meshes, finite element spaces and solve problems on these meshes.

I use Matlab very often for numerical computations. Most of the numerical stuff I’ve done (take a look here if you want) was based on finite differences methods, fundamental solutions and other classical techniques different from finite elements methods. Once I started using finite elements I quickly realized that Matlab is not that easy to work with if you want some automated quality meshing. PDEtool is good, but defining the geometry is not easy. There is also a simple tool: distmesh which performs a simple mesh construction for simple to state geometries. Nevertheless, once you have to define finite element spaces and solve problems things are not easy…

This brings us to the topic of this post: is it possible to interface Matlab and FreeFem? First, why would someone like to do this? Matlab is easier to code and use than FreeFem (for one who’s not a C expert…), but FreeFem deals better with meshes and solving PDE with finite elements. Since FreeFem can be called using system commands, it is possible to call a static program from Matlab. FreeFem can save meshes and variables to files. Let’s see how can we recover them in Matlab.

There is a tool called “FreeFem to Matlab” developed by Julien Dambrine (link on Mathworks). There’s also a longer explanation in this blog post. I recently tried to use the tool and I quickly found that it is not appropriate for large meshes. It probably scans the mesh file line by line which makes the loading process lengthy for high quality meshes. Fortunately there’s a way to speed up things and I present it below. I will not cover the import of the data (other than meshes) since the function importdata from the FreeFem to Matlab tool is fast enough for this.

Read more…

ICPC 2015 World Final Problem B

May 21, 2015 Leave a comment

This the the solution to Problem B from the International Collegiate Programming Contest. The list of problems can be found here. The idea is to find the maximal area of the intersection of two moving polygons. The inputs give the initial positions of two convex polygons and the vectors giving their speed in the plane.

One idea of a solution is as follows: take a maximal time {T} and a discretization of the interval {[0,T]} by dividing it into {N} parts. Iterate over the translations at times given by this discretization, compute at each step the area of the intersection of the two polygons (if there is any intersection at all), and in the end find the time {t_0} for which this area is maximized. Now, even if the discretization step {T/N} is large (greater than the demanded precision of {10^{-3}}), we can conclude that {t_0} is an approximation of the final time with an error smaller than {T/N}. This is due to the fact that the function representing the area of the intersection has two monotonicity intervals, as a consequence of the fact that the polygons are convex. On the first interval, the intersection area is increasing, on the second one it is decreasing. Thus, once we have a discrete maximum, we are close to the real maximum.

Now, all we are left to do in order to achieve any desired precision is to refine the search near this discrete maximum, find a new, better approximation, and refine again, until we have enough precision. Of course, one first problem with this algorithm is the initial search using a maximal time {T}. It is possible that if {T} or {N} are not large enough, then we do not detect any intersection of the two polygons. Thus, an initial guess, based on a mathematical argument is needed in order to reduce the search of the optimal time to an interval which is small enough to have enough initial precision to detect the discrete maximum.

The algorithm presented below, uses Matlab’s predefined function polybool which can compute the intersection of two polygons. Of course, this is an overkill, since dealing with intersections of convex polygons is not that complicated (but still, I didn’t have enough time to play with the problem, in order to provide a more optimized version). I do not treat the search for the initial time interval. As I think about it, I guess a n argument based on finding some line intersections should give us a narrow enough time interval (with some care for the case when the two speed directions are collinear). The algorithm presented below solves the sample cases, but could fail in more general situations.

function Prob1_2015(p1,p2)

% choose initial polygons; see the text of the problem
p1 = [6 3 2 2 4 3 6 6 6 7 4 6 2 2 2];
p2 = [4 18 5 22 9 26 5 22 1 -2 1];

%p1 = [4 0 0 0 2 2 2 2 0 1 1];
%p2 = [4 10 0 10 2 12 2 12 0 -1 1];

np1 = p1(1);
np2 = p2(1);
cp1 = p1(2:1+2*np1);
cp2 = p2(2:1+2*np2);

vp1 = p1(end-1:end);
vp2 = p2(end-1:end);

cp1 = cp1(:);
cp1 = reshape(cp1,2,np1);

xp1 = cp1(1,:);
yp1 = cp1(2,:);
cp2 = cp2(:);
cp2 = reshape(cp2,2,np2);

xp2 = cp2(1,:);
yp2 = cp2(2,:);

%set precision
prec = 1;

%here there should be a clever initial choice of 
%the starting time and stopping time
%this choice works well in this case
start = 0;
stop  = 5;

while prec>1e-6
n = 100;

timex = linspace(start,stop,n);
areas = zeros(size(timex));

for i = 1:n
   t = timex(i);
   xt1 = xp1+t*vp1(1);
   yt1 = yp1+t*vp1(2);
   xt2 = xp2+t*vp2(1);
   yt2 = yp2+t*vp2(2);
   [x,y] = polybool('intersection',xt1,yt1,xt2,yt2);
   areas(i) = polyarea(x,y);
[m,I] = max(areas);
if m<1e-6
   % if no polygonal intersection is detected
   % there is nothing further to do

tapp = timex(I);
fprintf('Precision %f | Time %f\n',prec,tapp);
start = tapp-prec;
stop  = tapp+prec;
prec  = prec/10;
%this draws the final position of the polygon
%as well as their intersection
axis equal
hold on
axis equal
hold off

Here’s what you get for the first sample test provided in the questions:polygons3

 >> Prob1_2015
    res = 4.193548

ICPC 2015 World Final Problem A

May 21, 2015 Leave a comment

This is the solution to Problem A from the International Collegiate Programming Contest. The list of problems can be found here. This first problem consists simply of reading the parameters of the function defined below, and computing its values on the set {\{1,2,...,n\}}. Then, you need to find the maximum decrease in the function values.

The inputs are parameters {p,a,b,c,d,n}, the function is

\displaystyle \text{price}(k) =p(\sin(ak+b)+\cos(ck+d)+2).

and the values to be considered are {p(1),p(2),...,p(n)}. Below is a Matlab code which works, at least for the given Sample Inputs. This should be optimized by removing the for loop. It is instant for n=1000 fg=000000 but it should be modified to work until n=10^6 fg=000000.

function res = Prob1_2015(p,a,b,c,d,n)
dis = 1:n;

vals = p*(sin(a*dis+b)+cos(c*dis+d)+2);

mat = zeros(n,1);
for i = 1:n
mat(i) = vals(i)-min(vals(i:end));

res = max(mat)

The sample outputs are:

 >> Prob1_2015(42,1,23,4,8,10)  
 >> Prob1_2015(100,7,615,998,801,3)
 >> Prob1_2015(100,432,406,867,60,1000)

New version which works for large n. It suffices to eliminate the computation of the min at each of the phases of the iteration. This solves the problem in 0.2 seconds for n= 10^6

function res = Prob1_2015(p,a,b,c,d,n)
dis = 1:n;
vals = p*(sin(a*dis+b)+cos(c*dis+d)+2);
mat       = zeros(n,1);
lastmax   = vals(1);
for i = 1:n
  lastmax   = max(lastmax,vals(i));
  mat(i)    = lastmax-vals(i);
res = max(mat)
Categories: matlab, Olympiad Tags: , ,
%d bloggers like this: