Planet Octave

August 19, 2017

Michele Ginesi

Integral functions

Integral functionsDuring the last week I made few modifications to expint.m and wrote sinint.m and cosint.m from scratch. All the work done can be found on the bookmark expint of my repository.


As I mentioned here I rewrote expint.m from scratch before the GSoC. During the last week I moved the Lentz algorithm to a .cc function (in order to remain coherent with the implementations of gammainc and betainc) and added few tests.


The sinint function is present in the symbolic package, but is not present a numerical implementation in the core.
The sine integral is defined as $$ \text{Si} (z) = \int_0^z \frac{\sin(t)}{t}\,dt. $$ To compute it we use the series expansion $$ \text{Si}(z) = \sum_{n=0}^\infty \frac{(-1)^n z^{2n+1}}{(2n+1)(2n+1)!} $$ when the module of the argument is smaller than 2. For bigger values we use the following relation with the exponential integral $$ \text{Si} = \frac{1}{2i} (E_1(iz)-E_1(-iz)) + \frac{\pi}{2},\quad |\text{arg}(z)| \frac{\pi}{2}$$ and the following simmetry relations $$ \text{Si}(-z) = -\text{Si}(z), $$ $$ \text{Si}(\bar{z}) = \overline {\text{Si}(z)}. $$ The function is write as a single .m file.


As the sinint function, also cosint is present in the symbolic package, but there is not a numerical implementation in the core.
The cosine integral is defined as $$ \text{Ci} (z) = -\int_z^\infty \frac{\cos(t)}{t}\,dt. $$ An equivalent definition is $$ \text{Ci} (z) = \gamma + \log z + \int_0^z \frac{\cos t - 1}{t}\,dt. $$ To compute it we use the series expansion $$ \text{Ci}(z) = \gamma + \log z + \sum_{n=1}^\infty \frac{(-1)^n z^{2n}}{(2n)(2n)!} $$ when the module of the argument is smaller than 2. For bigger values we use the following relation with the exponential integral $$ \text{Ci} = -\frac{1}{2} (E_1(iz)+E_1(-iz)),\quad |\text{arg}(z)| \frac{\pi}{2}$$ and the following simmetry relations $$ \text{Ci}(-z) = \text{Ci}(z) -i\pi,\quad 0\text{arg}(z)\pi, $$ $$ \text{Ci}(\bar{z}) = \overline{\text{Ci}(z)} .$$ As for sinint, also cosint is written as a single .m file.

by Michele Ginesi ( at August 19, 2017 02:32 AM

August 18, 2017

Joel Dahne

The Final Polish

We are preparing for releasing version 3.0.0 of the interval package and this last week have mainly been about fixing minor bugs related to the release. I mention two of the more interesting bugs here.

Compact Format

We (Oliver) recently added support for "format compact" when printing intervals. It turns out that the way to determine if compact format is enabled differs very much between different version of Octave. There are at least three different ways to get the information.

In the older releases (< 4.2.0 I believe) you use "get (0, "FormatSpacing")" but there appear to be a bug for version < 4.0.0 for which this always return "loose".

For the current tip of the development branch you can use "[~, spacing] = format ()" to get the spacing.

Finally in between these two version you use "__compactformat__ ()".

In the end Oliver, probably, found a way to handle this mess and compact format should now be fully supported for intervals. The function to do this is available here

Dot-product of Empty Matrices

When updating "dot" to support N-dimensional arrays I also modified it so that it behaves similar to Octaves standard implementation. The difference is in how it handles empty input. Previously we had

> x = infsupdec (ones (0, 2));
> dot (x, x)
ans = 0×2 interval matrix

but with the new version we get

> dot (x, x)
ans = 1×2 interval vector
   [0]_com   [0]_com

which is consistent with the standard implementation.

In the function we use "min" to compute the decoration for the result. Normally "min (x)" and "dot (x, x)" returns results of the same size (the dimension along which it is computed is set to 1), but they handle empty input differently. We have

> x = ones (0, 2);
> dot (x, x)
ans =
   0   0
> min (x)
ans = [](0x2)

This meant that the decoration would be incorrect since the implementation assumed they always had the same size. Fortunately the solution was very simple. If the dimension along which we are computing the dot-product is zero. the decoration should always be "com". So just adding a check for that was enough.

You could argue that "min (ones (0, 2))" should return "[inf, inf]" similarly to how many of the other reductions, like "sum" or "prod", return their unit for empty input. But this would most likely be very confusing for a lot of people. And it is not compatible with how Matlab does it either.

Updates on the Taylor Package

I have also had some time to work on the Taylor package this week. The basic utility functions are now completed and I have started to work on functions for actually computing with Taylor expansions. At the moment there are only a limited amount of functions implemented. For example we can calculate the Taylor expansion of order 4 for the functions $\frac{e^x + \log(x)}{1 + x}$ at $x = 5$.

## Create a variable of degree 4 and with value 5
> x = taylor (infsupdec (5), 4)
x = [5]_com + [1]_com X + [0]_com X^2 + [0]_com X^3 + [0]_com X^4

## Calculate the function
> (exp (x) + log (x))./(1 + x)
ans = [25.003, 25.004]_com + [20.601, 20.602]_com X + [8.9308, 8.9309]_com X^2 + [2.6345, 2.6346]_com X^3 + [0.59148, 0.59149]_com X^4

by Joel Dahne ( at August 18, 2017 05:01 PM

August 12, 2017

Michele Ginesi


betaincinvThe inverse of the incomplete beta function was present in Octave, but without the "upper" option (since it was missing in betainc itself). We decided to rewrite it from scratch using Newton method, as for gammaincinv (see my post on it if you are interested).
To make the code numerically more accurate, we decide which version ("lower" or "upper") invert depending on the inputs.
At first we compute the trivial values (0 and 1). Then the remaining terms are divided in two sets: those that will be inverted with the "lower" version, and those that will be inverted with the "upper" one. For both cases, we perform 10 iterations of bisection method and then we perform a Newton method.
The implementation (together with the new implementation of betainc) can be found on my repository, bookmark "betainc".

by Michele Ginesi ( at August 12, 2017 08:52 AM

August 11, 2017

Joel Dahne

Improving the Automatic Tests

Oliver and I have been working on improving the test framework used for the interval package. The package shares a large number of tests with other interval packages through an interval test framework that Oliver created. Here is the repository.

Creating the Tests

Previously these tests were separated from the rest of the package and you usually ran them with help of the Makefile. Now Oliver has moved them to the m-files and you can run them, together with the other tests for the function, with test @infsup/function in Octave. This makes it much easier to test the functions directly.

In addition to making the tests easier to use we also wanted to extend them to not only test scalar evaluation but also vector evaluation. The test data, input ad expected output, is stored in a cell array and when performing the scalar testing we simply loop over that cell and run the function for each element. The actual code looks like this (in this case for plus)

%! # Scalar evaluation
%! testcases = testdata.NoSignal.infsup.add;
%! for testcase = [testcases]'
%!   assert (isequaln (...
%!     plus ({1},{2}), ...
%!     testcase.out));
%! endfor

For testing the vector evaluation we simply concatenate the cell array into a vector and give that to the function. Here is what that code looks like

%! # Vector evaluation
%! testcases = testdata.NoSignal.infsup.add;
%! in1 = vertcat (vertcat ({:, 1});
%! in2 = vertcat (vertcat ({:, 2});
%! out = vertcat (testcases.out);
%! assert (isequaln (plus (in1, in2), out));

Lastly we also wanted to test evaluation of N-dimensional arrays. This is done by concatenating the data into a vector and then reshape that vector into an N-dimensional array. But what size should we use for the array? Well, we want to have at least three dimensions because otherwise we are not really testing N-dimensional arrays. My solution was to completely factor the length of the vector and use that as size, testsize = factor (length (in1)), and if the length of the vector has two or fewer factors we add a few elements to the end until we get at least three factors. This is the code for that

%! # N-dimensional array evaluation
%! testcases = testdata.NoSignal.infsup.add;
%! in1 = vertcat (vertcat ({:, 1});
%! in2 = vertcat (vertcat ({:, 2});
%! out = vertcat (testcases.out);
%! # Reshape data
%! i = -1;
%! do
%!   i = i + 1;
%!   testsize = factor (numel (in1) + i);
%! until (numel (testsize) > 2)
%! in1 = reshape ([in1; in1(1:i)], testsize);
%! in2 = reshape ([in2; in2(1:i)], testsize);
%! out = reshape ([out; out(1:i)], testsize);
%! assert (isequaln (plus (in1, in2), out));

This works very well, except when the number of test cases is to small. If the number of test is less than four this will fail. But there are only a handful of functions with that few tests so I fixed those independently.

Running the tests

Okay, so we have created a bunch of new tests for the package. Do we actually find any new bugs with them? Yes!

The function pow.m failed on the vector test. The problem? In one place $\&\&$ was used instead of $\&$. For scalar input I believe these behave the same but they differ for vector input.

Both the function nthroot.m and the function pownrev.m failed the vector test. Neither allowed vectorization of the integer parameter. For nthroot.m this is the same for standard Octave version so it should perhaps not be treated as a bug. The function pownrev.m uses nthroot.m internally so it also had the same limitation. This time I would however treat it as a bug because the function pown.m does allow vectorization of the integer parameter and if that supports it the reverse function should probably also do it. So I implemented support for vectorization of the integer parameter for both nthroot.m and pownrev.m and they now pass the test.

No problems were found with the N-dimensional tests that the vector tests did not find. This is a good indication that the support for N-dimensional arrays is at least partly correct. Always good to know!

by Joel Dahne ( at August 11, 2017 03:04 PM

August 09, 2017

Piyush Jain

Geometry Package (Octave)

Geometry package: Implement boolean operations on polygons

As part of GSoC 2017 , this project is intended to implement a set of boolean operations and supporting function for acting on polygons. These include the standard set of potential operations such as union/OR, intersection/AND, difference/subtraction, and exclusiveor/XOR. Other things to be implemented are the following functions: polybool, ispolycw, poly2ccw, poly2cw, poly2fv, polyjoin, and polysplit.

Bonding Period

After discussion with my mentor and keeping my proposal in mind, I tried to understand and enlist the tasks more specifically and in detail. I first tried to understand the conventions and other things about the organisation. I tried to understand the first basic thing which I will need throughout this project, i.e. how we create an oct-interface for the C++ codes to be executable in Octave. My first goal was to explore the possibility if the already implemented mex-interface geometry package can be improved in performance by replacing it with oct-interface. So, for understanding how these oct-files work and other things , I started to implement something in oct-files and getting familiar with it.

First Coding Phase

As stated, there is an already implemented Geometry3.0 package. This has mex interface for its functions. I tried to compare its performance with the oct interface version. For benchmarking, I first implemented my own function for polyUnion (using clipper library) with oct-interface (Find it here). Then, I compared its performance over a number of different sets of polgons (parametrized over number of vertices) and recorded the elapsed times with both the interfaces. On plotting a graph of Number of vertices V/s Elapsed time (for oct and mex), following obseravtions were made : </br>

  • The oct interface had a better performance over mex.
  • For 10000 vertices, the oct interface took about 0.008 seconds while the mex interface took about 0.014 seconds. This implies oct interface took 8X10e-3 seconds / 10e4 vertices i.e. 8X10e-7 seconds per vertex. For mex, it was 14X10e-7 seconds per vertex.
  • As it can be seen from the above data, the oct interface was not more than twice as good as mex interface.</br> From these observations, it was concluded that it is not worth to change the interface from mex to oct since there was not much improvement in performance. Thus, our next goal is now to incorporate new algorithms.

After spending decent time over checking out the new algorithm and understanding its implementation, I have now started to implement the polybool function. I also tried to compare its performance with the already implemented clipPolygon in the current geometry package. The new algorithm has a better performance than the old one.

The implementation of boolean operations on polygons, like DIFFERENCE , INTERSECTION , XOR and UNION with the new algorithm is almost done. A little work on tests, demos and documentation is still needed and I am working on that.

More about the algorithm by F. Martínez, A.J. Rueda, F.R. Feito

The algorithm is very easy to understand, among other things because it can be seen as an extension of the classical algorithm, based on the plane sweep, for computing the intersection points between a set of segments.When a new intersection between the edges of polygons is found, the algorithm subdivides the edges at the intersection point. This produces a plane sweep algorithm with only two kind of events: left and right endpoints, making the algorithm quite simple. Furthermore, the subdivision of edges provides a simple way of processing degeneracies. </br> Overall sketch of the approach for computing Boolean operations on polygons:

  • Subdivide the edges of the polygons at their intersection points.
  • Select those subdivided edges that lie inside the other polygon—or that do not lie depending on the operation.
  • Join the edges selected in step 2 to form the result polygon.</br>


    Let n be the total number of edges of all the polygons involved in the Boolean operation and k be the number of intersections of all the polygon edges. The whole algorithm runs in time O(n+k)log(n).</br>

After raw testing this new algorithm on several cases, I am now adding few tests in the m-script. Other than that, a demo has also been added in the m-script. The demo can be seen by the command demo clipPolygon_mrf. The tests have also been added. The test can be seen by test clipPolygon_mrf.

Second Coding Phase

After implementing the polybool function and checking it, we are planning to include it in the next release of geometry package. Now , to move forward, I am first importing some functions from last year GSoC’s repo and ensuring their MATLAB compatibility. Functions like poly2ccw, poly2cw, joinpolygons, splitPolygons have been created as aliases while ensuring their compatibility with their mathworks counterparts. Then after, there was some time invested on understanding the CGAL library and its implementation.

The further plan is to sync the matgeom package with the geometry package.

Third Coding Phase

Proceeding towards the next goal, the idea is to devise some way to automate the process somewhat, of syncing the matgeom and geometry package. The issue is that when a new release of geometry package is planned, there are some things which ahve been updated in matgeom but not in their geometry counterparts (if it exists). So, every time before releasing, so much time has to be invested in manually checking each edit and syncing it into the geometry.

To achieve this, first a workaround is implemented on a dummy repository - dummyMatGeom. Its master branch is matGeom (dummy) and there is another branch (named geometry) is created which contains geometry package (dummy). To test the entire procedure, go to the dummy repository dummyMatGeom , pull both branches in different folders, say “dummyMatGeom” for master branch and “dummyGeom” for geometry branch. Then follow the given steps as explained on the wiki page.


  • Clearly, the above procedure will only sync the script of the function, not it’s tests and demo, which are in separate folders in a Matlab package structure. Even if we try to concatenate their corresponding test/demo scripts with the function scripts (as it is in an octave package structure), there will be discrepancies because the notion or writing tests for octave and matlab packages are quite different. The way octave allows tests to work is unique to octave as explained here. SO, we can’t simply concatenate the Matlab test scripts with the functions.

  • Git doesn’t preserves the original version of geometry scripts and overwrites the whole file completely. For example :

  1. Original file at matGeom (upstream) ~~~~ % Bla bla % bla bla bla

function z = blabla (x,y) % Help of function for i=1:length(x) z(i) = x(i)*y(i); end

2. Ported to geometry

Copyright - Somebody

Bla bla

bla bla bla


Help of function

function z = blabla (x,y) z = x .* y; endfunction

3. Updated in matGeom

% Bla bla % bla bla bla

function z = blabla (x,y) % Help of function % updated to be more clear z = zeros (size(x)); for i=1:length(x) z(i) = x(i)*y(i); end

4. After syncing , the expected result is something like this :

Copyright - Somebody

Bla bla

bla bla bla


Help of function

updated to be more clear

function z = blabla (x,y) z = zeros (size(x)); z = x .* y; endfunction ~~~~

But, this doesn’t happen as expected. Git just finds the files which have been modified and overwrites those files completely. Considering the possibilities of the solutions, there are ways like git patch or git interactive which allows us to select the lines specifically which we want to be committed, but that would not serve our purpose as it would not be better than syncing it manually, file by file. Looking for a better solution to handle this !

Now, the further idea is to release geometry and I am getting involved into it to get a feel of how things are done.

August 09, 2017 07:30 PM

August 02, 2017

Michele Ginesi


table, th, td { border: 1px solid black; } betaincThe betainc function has two bugs reported: #34405 on the input validation and #51157 on inaccurate result. Moreover, it is missing the "upper" version, which is present in MATLAB.

The function

The incomplete beta function ratio is defined as $$I_x(a,b) = \dfrac{B_x(a,b)}{B(a,b)},\quad 0\le x \le 1,\,a>0,\,b>0,$$ where $B(a,b)$ is the classical beta function and $$B_x(a,b)=\int_0^x t^{a-1}(1-t)^{b-1}\,dt.$$ In the "upper" version the integral goes from $x$ to $1$. To compute this we will use the fact that $$\begin{array}{rcl} I_x(a,b) + I_x^U(a,b) &=& \dfrac{1}{B(a,b)}\left( \int_0^x t^{a-1}(1-t)^{b-1}\,dt + \int_x^1 t^{a-1}(1-t)^{b-1}\,dt\right)\\ &=&\dfrac{1}{B(a,b)}\int_0^1 t^{a-1}(1-t)^{b-1}\,dt\\ &=&\dfrac{B(a,b)}{B(a,b)}\\ &=&1 \end{array}$$ and the relation $$I_x(a,b) + I_{1-x}(b,a) = 1$$ so that $$I_x^U(a,b) = I_{1-x}(b,a).$$

The implementation

Even if it is possible to obtain a Taylor series representation of the incomplete beta function, it seems to not be used. Indeed the MATLAB help cite only the continuous fraction representation present in "Handbook of Mathematical Functions" by Abramowitz and Stegun: $$I_x(a,b) = \dfrac{x^a(1-x)^b}{aB(a,b)}\left(\dfrac{1}{1+} \dfrac{d_1}{1+} \dfrac{d_2}{1+}\ldots\right)$$ with $$d_{2m+1} = -\dfrac{(a+m)(a+b+m)}{(a+2m)(a+2m+1)}x$$ and $$d_{2m} = \dfrac{m(b-m)}{(a+2m-1)(a+2m)}x$$ which seems to be the same strategy used by GSL. To be more precise, this continued fraction is computed directly when $$x\dfrac{a-1}{a+b-2}$$ otherwise, the computed fraction is used to compute $I_{1-x}(b,a)$ and then it is used the fact that $$I_x(a,b) = 1-I_{1-x}(b,a).$$ In my implementation I use a continued fraction present in "Handboob of Continued Fractions for Special Functions" by Cuyt, Petersen, Verdonk, Waadeland and Jones, which is more complicated but converges in fewer steps: $$\dfrac{B(a,b)I_x(a,b)}{x^a(1-x)^b} = \mathop{\huge{\text{K}}}_{m=1}^\infty \left(\dfrac{\alpha_m(x)}{\beta_m(x)}\right),$$ where $$\begin{array}{rcl} \alpha_1(x) &=&1,\\ \alpha_{m+1}(x) &=&\dfrac{(a+m-1)(a+b+m-1)(b-m)m}{(a+2m-1)^2}x^2,\quad m\geq 1,\\ \beta_{m+1}(x) &=&a + 2m + \left( \dfrac{m(b-m)}{a+2m-1} - \dfrac{(a+m)(a+b+m)}{a+2m+1} \right)x,\quad m\geq 0. \end{array}$$ This is most useful when $$x\leq\dfrac{a}{a+b},$$ thus, the continued fraction is computed directly when this condition is satisfied, while it is used to evaluate $I_{1-x}(b,a)$ otherwise.
The function is now written as a .m file, which check the validity of the inputs and divide the same in the values which need to be rescaled and in those wo doesn't need it. Then the continued fraction is computed by an external .c function. Finally, the .m file explicit $I_x(a,b)$.


Next step will be to write the inverse. It was already present in Octave, but is missing the upper version, so it has to be rewritten.

by Michele Ginesi ( at August 02, 2017 04:12 AM

August 01, 2017

Michele Ginesi

Second period resume

Second period resumetd, th{ border: 1px solid } Here I present a brief resume of the work done in this second month.

Bessel functions

The topic of this month were Bessel functions. On the bug tracker it is reported only a problem regarding the "J" one (see bug 48316), but the same problem is present in every type of Bessel function (they return NaN + NaNi when the argument is too big).

Amos, Cephes, C++ and GSL

Actually, Bessel functions in Octave are computed via the Amos library, written in Fortran. Studying the implementation I discovered that the reported bug follows from the fact that if the input is too large in module, the function zbesj.f set IERR to 4 (IERR is a variable which describe how the algorithm has terminate) and set the output to zero, then return NaN when IERR=4. Obviusly, the same happen for other Bessel functions.
What I initially did was to "unlock" these .f files in such a way to give still IERR=4, but computing anyway the output, and modify in order to return the value even if IERR is 4. Then I tested the accuracy, together with other libraries.
On the bug report where suggested the Cephes library, so I tested also them, the C++ Special mathematical functions library, and the GSL (Gnu Scientific library). Unfortunately, both these alternatives work worse than Amos. I also tried to study and implement some asymptotic expansions by myself to use in the cases which give inaccurate results, unfortunately without success.
For completeness, in the following table there are some results of the tests (ERR in Cephes refers to cos total loss of precision error):
1e09 1e10 1e11 1e12 1e13
Amos 1.6257e-16 0.0000e+00 0.0000e+00 1.3379e-16 1.1905e-16
Cephes 2.825060e-08 ERR ERR ERR ERR
GSL 4.8770e-16 2.2068e-16 4.2553e-16 1.3379e-16 1.1905e-16
C++ 2.82506e-08 2.68591e-07 1.55655e-05 8.58396e-07 0.000389545

1e15 1e20 1e25 1e30
Amos 1.3522e-16 1.6256e-16 15.22810 2.04092
GSL 1.3522e-16 0 15.22810 2.04092

The problem with the double precision

As I explained in my last post the problem seems to be that there are no efficient algorithms in double precision for arguments so large. In fact, the error done by Amos is quite small if compared with the value computed with SageMath (which I'm using as reference one), but only if we use the double precision also in Sage: using more digits, one can see that even the first digit change. Here the tests:
sage: bessel_J(1,10^40).n()
sage: bessel_J(1,10^40).n(digits=16)
sage: bessel_J(1,10^40).n(digits=20)
sage: bessel_J(1,10^40).n(digits=25)
sage: bessel_J(1,10^40).n(digits=30)
sage: bessel_J(1,10^40).n(digits=35)
sage: bessel_J(1,10^40).n(digits=40)
The value stabilizes only when we use more than 30 digits (two times the digits used in double precision).

The decision

Even if we are aware of the fact that the result is not always accurate, for MATLAB compatibility we decided to unlock the Amos, since they are still the most accurate and, even more important, the type of the inputs is the same as in MATLAB (while, for example, GSL doesn't accept complex $x$ value). Moreover, in Octave is possible to obtain in output the value of IERR, thing that is not possible in MATLAB.
You can find the work on the bookmark "bessel" of my repository.


During these last days I also started to implement betainc from scratch, I think it will be ready for the first days of August. Then, it will be necessary to rewrite also betaincinv, since the actual version doesn't have the "upper" version. This should not be too difficult. I think we can use a simple Newton method (as for gammaincinv), the only problem will be, as for gammaincinv, to find good initial guesses.

by Michele Ginesi ( at August 01, 2017 09:30 AM

July 28, 2017

Joel Dahne

A Package for Taylor Arithmetic

In the last blog post I wrote about what was left to do with implementing support for N-dimensional arrays in the interval package. There are still some things to do but I have had, and most likely will have, some time to work on other things. Before the summer I started to work on a proof of concept implementation of Taylor arithmetic in Octave and this week I have continued to work on that. This blog post will be about that.

A Short Introduction to Taylor Arithmetic

Taylor arithmetic is a way to calculate with truncated Taylor expansions of functions. The main benefit is that it can be used to calculate derivatives of arbitrary order.

Taylor expansion or Taylor series (I will use these words interchangeably) are well known and from Wikipedia we have: The Taylor series of real or complex valued function $f(x)$ that is infinitely differentiable at a real or complex number $a$ is the power series
f(a) + \frac{f'(a)}{1!}(x-a) + \frac{f''(a)}{2!}(x-a)^2 + \frac{f'''(a)}{3!}(x-a)^3 + ....
From the definition it is clear that if we happen to know the coefficients of the Taylor series of $f$ at the point $a$ we can also calculate all derivatives of $f$ at that point by simply multiplying a coefficient with the corresponding factorial.

The simplest example of Taylor arithmetic is addition of two Taylor series. If $f$ has the Taylor series $\sum_{n=0}^\infty (f)_n (x-a)^n$ and $g$ the Taylor series $\sum_{n=0}^\infty (g)_n (x-a)^n$ then $f + g$ will have the Taylor series
\sum_{n=0}^\infty (f + g)_n (x-a)^n = \sum_{n=0}^\infty ((f)_n + (g)_n)(x-a)^n$
If we instead consider the product, $fg$, we get
\sum_{n=0}^\infty (fg)_n (x-a)^n = \sum_{n=0}^\infty \left(\sum_{i=0}^n (f)_n(g)_n\right)(x-a)^n.

With a bit of work you can find similar formulas for other standard functions. For example the coefficients, $(e^f)_n$, of the Taylor expansion of $\exp(f)$ is given by $(e^f)_0 = e^{(f)_0}$ and for $n > 0$
(e^f)_n = \frac{1}{n}\sum_{i=0}^{n-1}(k-j)(e^f)_i(f)_{n-i}.

When doing the computations on a computer we consider truncated Taylor series, we choose an order and keep only coefficients up to that order. There is also nothing that stops us from using intervals as coefficients, this allows us to get rigorous enclosures of derivatives of functions.

For a more complete introduction to Taylor arithmetic in conjunction with interval arithmetic see [1] which was my first encounter to it. For another implementation of it in code take a look at [2].

Current Implementation Status

As mentioned in the last post my repository can be found here

When I started to write on the package, before summer, my main goal was to get something working quickly. Thus I implemented the basic functions needed to do some kind of Taylor arithmetic, a constructor, some help functions and a few functions like $\exp$ and $\sin$.

This last week I have focused on implementing the basic utility functions, for example $size$, and rewriting the constructor. In the process I think I have broken the arithmetic functions, I will fix them later.

You can at least create and display Taylor expansions now. For example creating a variable $x$ with value 5 of order 3

> x = taylor (infsupdec (5), 3)
x = [5]_com + [1]_com X + [0]_com X^2 + [0]_com X^3

or a matrix with 4 variables of order 2

> X = taylor (infsupdec ([1, 2; 3, 4]), 2)
X = 2×2 Taylor matrix of order 2

ans(:,1) =

   [1]_com + [1]_com X + [0]_com X^2
   [3]_com + [1]_com X + [0]_com X^2

ans(:,2) =

   [2]_com + [1]_com X + [0]_com X^2
   [4]_com + [1]_com X + [0]_com X^2

If you want you can create a Taylor expansion with explicitly given coefficients you can do that as well

> f = taylor (infsupdec ([1; -2; 3, -4))
f = [1]_com + [-2]_com X + [3]_com X^2 + [-4]_com X^3

This would represent a function $f$ with $f(a) = 1$, $f'(a) = -2$, $f''(a) = 3 \cdot 2! = 6$ and $f'''(a) = -4 \cdot 3! = -24$.

Creating a Package

My goal is to create a full package for Taylor arithmetic along with some functions making use of it. The most important step is of course to create a working implementation but there are other things to consider as well. I have a few things I have not completely understood about it. Depending on how much time I have next week I will try to read a bit more about it probably ask some questions on the mailing list. Here are at least some of the things I have been thinking about

Mercurial vs Git?

I have understood that most of the Octave forge packages uses Mercurial for version control. I was not familiar with Mercurial before so the natural choice for me was to use Git. Now I feel I could switch to Mercurial if needed but I would like to know the potential benefits better, I'm still new to Mercurial so I don't have the full picture. One benefit is of course that it is easier if most  packages use the same system, but other than that?

How much work is it?

If I were to manage a package for Taylor arithmetic how much work is it? This summer I have been working full time with Octave so I have had lots of time but this will of course not always be the case. I know it takes time if I want to continue to improve the package, but how much, and what, continuous work is there?

What is needed besides the implementation?

From what I have understood there are a couple of things that should be included in a package besides the actual m-files. For example a Makefile for creating the release, an INDEX-file and a CITATION-file. I should probably also include some kind of documentation, especially since Taylor arithmetic is not that well known. Is there anything else I need to think about?

What is the process to get a package approved?

If I were to apply (whatever that means) for the package to go to Octave forge what is the process for that? What is required before it can be approved and what is required after it is approved?

[1] W. Tucker, Validated Numerics, Princeton University Press, 2011.
[2] F. Blomquist, W. Hofschuster, W. Krämer, Real and complex taylor arithmetic in C-XSC, Preprint 2005/4, Bergische Universität Wuppertal.

by Joel Dahne ( at July 28, 2017 09:56 PM

July 27, 2017

Enrico Bertino

Deep learning functions

Hi there,

the second part of the project is finishing. This period was quite interesting because I had to dive into the theory behind Neural Networks  In particular [1], [2], [3] and [4] were very useful and I will sum up some concepts here below. On the other hand, coding became more challenging and the focus was on the python layer and in particular the way to structure the class in order to make everything scalable and generalizable. Summarizing the situation, in the first period I implemented all the Octave classes for the user interface. Those are Matlab compatible and they call some Python function in a seamless way. On the Python side, the TensorFlow API is used to build the graph of the Neural Network and perform training, evaluation and prediction.

I implemented the three core functions: trainNetwork, SeriesNetwork and trainingOptions. To do this, I used a Python class in which I initialize an object with the graph of the network and I store this object as attribute of SeriesNetwork. Doing that, I call the methods of this class from trainNetwork to perform the training and from predict/classify to perform the predictions. Since it was quite hard to have a clear vision of the situation, I used a Python wrapper (Keras) that allowed me to focus on the integration, "unpack" the problem and go forth "module" by "module". Now I am removing the dependency on the Keras library using directly the Tensorflow API. The code in my repo [5].

Since I have already explained in last posts how I structured the package, in this post I would like to focus on the theoretical basis of the deep learning functions used in the package. In particular I will present the available layers and the parameters that are available for the training.  

Theoretical dive

I. Fundamentals

I want to start with a brief explanation about the perceptron and the back propagation, two key concepts in the artificial neural networks world. 


Let's start from the perceptron, that is the starting point for understanding neural networks and its components. A perceptron is simply a "node" that takes several binary inputs,  $ x_1, x_2, ... $, and produces a single binary output:

The neuron's output, 0 or 1, is determined by whether the linear combination of the inputs $  \omega \cdot x = \sum_j \omega_j x_j  $ is less than or greater than some threshold value. That is a simple mathematical model but is very versatile and powerful because we can combine many perceptrons and varying the weights and the threshold we can get different models. Moving the threshold to the other side of the inequality and replacing it by what's known as the perceptron's bias, b = −threshold, we can rewrite it as

$ out = \bigg \{ \begin{array}{rl} 0 & \omega \cdot x + b \leq 0 \\ 1 & \omega \cdot x + b > 0 \\ \end{array} $

Using the perceptrons like artificial neurons of a network, it turns out that we can devise learning algorithms which can automatically tune the weights and biases. This tuning happens in response to external stimuli, without direct intervention by a programmer and this enables us to have an "automatic" learning.

Speaking about learning algorithms, the proceedings are simple: we suppose we make a small change in some weight or bias and what see the corresponding change in the output from the network. If a small change in a weight or bias causes only a small change in output, then we could use this fact to modify the weights and biases to get our network to behave more in the manner we want. The problem is that this isn't what happens when our network contains perceptrons since a small change of any single perceptron can sometimes cause the output of that perceptron to completely flip, say from 0 to 1. We can overcome this problem by introducing an activation function. Instead of the binary output we use a function depending on weights and bias. The most common is the sigmoid function:
$ \sigma (\omega \cdot x + b ) = \dfrac{1}{1 + e^{-(\omega \cdot x + b ) } } $

Figure 1. Single neuron 

With the smoothness of the activation function $ \sigma $ we are able to analytically measure the output changes since $ \Delta out $ is a linear function of the changes $ \Delta \omega $ and $ \Delta b$ :
$ \Delta out \approx \sum_j \dfrac{\partial out}{\partial \omega_j} \Delta \omega_j + \dfrac{\partial out}{\partial b} \Delta b $

Loss function

Let x be a training input and y(x) the desired output. What we'd like is an algorithm which lets us find weights and biases so that the output from the network approximates y(x) for all x. Most used loss function is mean squared error (MSE) :
$ L( \omega, b) = \dfrac{1}{2n} \sum_x || Y(x) - out ||^2 $ ,
where n is the total number of training inputs, out is the vector of outputs from the network when x is input. 
To minimize the loss function, there are many optimizing algorithms. The one we will use is the gradient descend, of which every iteration of an epoch is defined as:

$ \omega_k \rightarrow \omega_k' = \omega_k - \dfrac{\eta}{m} \sum_j \dfrac{\partial L_{X_j}}{\partial \omega_k} $
$ b_k \rightarrow b_k' = b_k - \dfrac{\eta}{m} \sum_j \dfrac{\partial L_{X_j}}{\partial b_k} $

where m is the size of the batch of inputs with which we feed the network and $ \eta $ is the learning rate.


The last concept that I would like to emphasize is the backpropagation. Its goal is to compute the partial derivatives$ \partial L / \partial  \omega $ and $ \partial L / \partial b} $ of the loss function L with respect to any weight or bias in the network. The reason is that those partial derivatives are computationally heavy and the network training would be excessively slow.

Let be $ z^l $ the weighted input to the neurons in layer l, that can be viewed as a linear function of the activations of the previous layer: $ z^l = \omega^l a^{l-1} + b^l $ .
In the fundamental steps of backpropagation we compute:

1) the final error:
$ \delta ^L = \Delta_a L \odot \sigma' (z^L) $
The first term measures how fast the loss is changing as a function of every output activation and the second term measures how fast the activation function is changing at $ z_L $

2) the error of every layer l:
$ \delta^l = ((\omega^{l+1})^T \delta^{l+1} ) \odot \sigma' (z^l) $

3) the partial derivative of the loss function with respect to any bias in the net
$ \dfrac{\partial L}{\partial b^l_j} = \delta^l_j $

4) the partial derivative of the loss function with respect to any weight in the net
$ \dfrac{\partial L}{\partial \omega^l_{jk}} = a_k^{l-1} \delta^l_j $

We can therefore update the weights and the biases with the gradient descent and train the network. Since inputs can be too numerous, we can use only a random sample of the inputs. Stochastic Gradient Descent (SGD) simply does away with the expectation in the update and computes the gradient of the parameters using only a single or a few training examples. In particular, we will use the SGD with momentum, that is a method that helps accelerate SGD in the relevant direction and damping oscillations. It does this by adding a fraction γ of the update vector of the past time step to the current update vector.

II. Layers

Here a brief explanation of the functions that I am considering in the trainNetwork class


The convolution layer is the core building block of a convolutional neural network (CNN) and it does most of the computational heavy lifting. They derive their name from the “convolution” operator. The primary purpose of convolution is to extract features from the input image preserving the spatial relationship between pixels by learning image features using small squares of input data. 

Figure 2. Feature extraction with convolution (image taken form 
In the example in Fig.2, the 3×3 matrix is called a 'filter' or 'kernel' and the matrix formed by sliding the filter over the image and computing the dot product is called the 'Convolved Feature' or 'Activation Map' (or the 'Feature Map'). In practice, a CNN learns the values of these filters on its own during the training process (although we still need to specify parameters such as number of filters, filter size, architecture of the network etc. before the training process). The more number of filters we have, the more image features get extracted and the better our network becomes at recognizing patterns in unseen images.
The size of the Feature Map depends on three parameters: the depth (that corresponds to the number of filters we use for the convolution operation), the stride (that is the number of pixels by which we slide our filter matrix over the input matrix) and the padding (that consists in padding the input matrix with zeros around the border).


ReLU stands for Rectified Linear Unit and is a non-linear operation: $ f(x)=max(0,x) $. Usually this is applied element-wise to the output of some other function, such as a matrix-vector product. It replaces all negative pixel values in the feature map by zero with the purpose of introducing non-linearity in our network, since most of the real-world data we would want to learn would be non-linear.


Neurons in a fully connected layer have full connections to all activations in the previous layer, as seen in regular Neural Networks. Hence their activations can be computed with a matrix multiplication followed by a bias offset. In our case, the purpose of the fully-connected layer is to use these features for classifying the input image into various classes based on the training dataset. Apart from classification, adding a fully-connected layer is also a cheap way of learning non-linear combinations of the features. 


It is common to periodically insert a pooling layer in-between successive convolution layers. Spatial Pooling (also called subsampling or downsampling) reduces the dimensionality of each feature map but retains the most important information. In particular, pooling
makes the input representations (feature dimension) smaller and more manageable
reduces the number of parameters and computations in the network, therefore, controlling overfitting
makes the network invariant to small transformations, distortions and translations in the input image
helps us arrive at an almost scale invariant representation of our image 
Spatial Pooling can be of different types: Max, Average, Sum etc.

In case of Max Pooling, we define a spatial neighborhood (for example, a 2×2 window) and take the largest element from the rectified feature map within that window. 

Instead of taking the largest element we could also take the average.


Dropout in deep learning works as follows: one or more neural network nodes is switched off once in a while so that it will not interact with the network. With dropout, the learned weights of the nodes become somewhat more insensitive to the weights of the other nodes and learn to decide somewhat more by their own. In general, dropout helps the network to generalize better and increase accuracy since the influence of a single node is decreased.


The purpose of the softmax classification layer is simply to transform all the net activations in your final output layer to a series of values that can be interpreted as probabilities. To do this, the softmax function is applied onto the net intputs.
$ \phi_{softmax} (z^i) = \dfrac{e^{z^i}}{\sum_{j=0}^k e^{z_k^i}} $


Local Response Normalization (LRN) layer implements the lateral inhibition that in neurobiology refers to the capacity of an excited neuron to subdue its neighbors. This layer is useful when we are dealing with ReLU neurons because they have unbounded activations and we need LRN to normalize that. We want to detect high frequency features with a large response. If we normalize around the local neighborhood of the excited neuron, it becomes even more sensitive as compared to its neighbors. At the same time, it will dampen the responses that are uniformly large in any given local neighborhood. If all the values are large, then normalizing those values will diminish all of them. So basically we want to encourage some kind of inhibition and boost the neurons with relatively larger activations.

III. training options

The training function takes as input a trainingOptions object that contains the parameters for the training. A brief explanation:

Optimizer chosen for minimize the loss function. To guarantee the Matlab compatibility, only the Stochastic Gradient Descent with Momentum ('sgdm') is allowed

Parameter for the sgdm: it corresponds to the contribution of the previous step to the current iteration

Initial learning rate η for the optimizer

These are the settings for regulating the learning rate. It is a struct containing three values:
RateSchedule: if it is set to 'piecewise', the learning rate will drop of a RateDropFactor every a RateDropPeriod number of epochs. 

Regularizers allow to apply penalties on layer parameters or layer activity during optimization. This is the factor of the L2 regularization.

Number of epochs for training

Display the information of the training every VerboseFrequency iterations

Random shuffle of the data before training if set to 'once'

Path for saving the checkpoints

Chose of the hardware for the training: 'cpu', 'gpu', 'multi-gpu' or 'parallel'. the load is divided between workers of GPUs or CPUs according to the relative division set by WorkerLoad

Custom output functions to call during training after each iteration passing a struct containing:
Current epoch number, Current iteration number, TimeSinceStart, TrainingLoss, BaseLearnRate, TrainingAccuracy (or TrainingRMSE for regression), State


[1] [Ian Goodfellow, Yoshua Bengio, Aaron Courville] Deep Learning,  MIT Press,, 2016

by Enrico Bertino ( at July 27, 2017 09:34 PM

July 14, 2017

Joel Dahne

Ahead of the Timeline

One of my first posts on this blog was a timeline for my work during the project. Predicting the amount of time something takes is always hard to do. Often time you tend to underestimate the complexity of parts of the work. This time however I overestimated the time the work would take.

If my timeline would have been correct I would have just started to work on Folding Functions (or reductions as they are often called). Instead I have completed the work on them and also for functions regarding plotting. In addition I have started to work on the documentation for the package as well as checking everything an extra time.

In this blog post I will go through what I have done this week, what I think is left to do and a little bit about what I might do if I complete the work on N-dimensional arrays in good time.

This Week

The Dot Function

The $dot$-function was the last function left to implement support for N-dimensional arrays in. It is very similar to the $sum$-function so I already had an idea of how to do it. As with  $sum$ I moved most of the handling of the vectorization from the m-files to the oct-file, the main reason being improved performance.

The $dot$-functions for intervals is actually a bit different from the standard one. First of all it supports vectorization which the standard one does not

> dot ([1, 2, 3; 4, 5, 6], 5)
error: dot: size of X and Y must match
> dot (infsupdec ([1, 2, 3; 4, 5, 6], 5)
ans = 1x3 interval vector

  [25]_com   [35]_com   [45]_com

It also treats empty arrays a little different, see bug #51333,

> dot ([], [])
ans = [](1x0)
> dot (infsupdec ([]), [])
ans = [0]_com

Package Documentation

I have done the minimal required changes to the documentation. That is I moved support for N-dimensional arrays from Limitation to Features and added some simple examples on how to create N-dimensional arrays.

Searching for Misses

During the work I have tried to update the documentation for all functions to account for the support of N-dimensional arrays and I have also tried to update some of the comments for the code. But as always, especially when working with a lot of files, you miss things, both in the documentation and old comments.

I did a quick grep for the words "matrix" and "matrices" since they are candidates for being changed to "array". Doing this I found 35 files where I missed things. It was mainly minor things, comments using the "matrix" which I now changed to "array", but also some documentation which I had forgotten to update.

What is Left?

Package Documentation - Examples

As mentioned above I have done the minimal required changes to the documentation. It would be very nice to add some more interesting example using N-dimensional arrays of intervals in a useful way. Ironically I have not been able to come up with an interesting example but I will continue to think about it. If you have an example that you think would be interesting and want to share, please let me know!

Coding Style

As I mentioned in one of the first blog posts, the coding style for the interval package was not following the standard for Octave. During my work I have adapted all files I have worked with to the coding standard for Octave. A lot of the files I have not needed to do any changes to, so they are still using the old style. It would probably be a good idea to update them as well.

Testing - ITF1788

The interval framwork libary developed by Oliver is used to test the correctness of many of the functions in the package. At the moment it tests evaluation of scalars but in principle it should be no problem to use it for testing vectorization or even broadcasting. Oliver has already started to work on this.

After N-dimensional arrays?

If I continue at this pace I will finish the work on N-dimensional arrays before the time of the project is over. Of course the things that are left might take longer than expected, they usually do, but there is a chance that I will have time left after everything is done. So what should I do then? There are more thing that can be done on the interval package, for example adding more examples to the documentation, however I think I would like to start working on a new package for Taylor arithmetics.

Before GSoC I started to implement a proof of concept for Taylor arithmetics in Ocatve which can be found here. I would then start to work on implementing a proper version of it, where I would actually make use of N-dimensional interval arrays. If I want to create a package for this I would also need to learn a lot of other things, one of them being how to manage a package on octave forge.

At the moment I will try to finish my work on N-dimensional arrays. Then I can discuss it with Oliver and see what he thinks about it.

by Joel Dahne ( at July 14, 2017 04:59 PM

July 13, 2017

Joel Dahne

Set inversion with fsolve

This week my work have mainly been focused on the interval version of fsolve. I was not sure if and how this could make use of N-dimensional arrays and to find that out I had to understand the function. In the end it turned out that the only generalization that could be done were trivial and required very few changes. However I did find some other problems with the functions that I have been able to fix. Connected to fsolve are the functions ctc_intersect and ctc_union. The also needed only minor changes to allow for N-dimensional input. I will start by giving an introduction to fsolve, ctc_union and ctc_intersect and then I will mention the changes I have done to them.

Introduction to fsolve

The standard version of fsolve in Octave is used to solve systems of nonlinear equations. That is, given a functions $f$ and a starting point $x_0$ it returns a value $x$ such that $f(x)$ is close to zero. The interval version of fsolve does much more than this. It is used to enclose the preimage of a set $Y$ under $f$. Given a domain $X$, a set $Y$ and a function $f$ it returns an enclosure of the set
f^{-1}(Y) = \{x \in X: f(x) \in Y\}.
By letting $Y = 0$ we get similar functionality as the standard fsolve, with the difference that the output is an enclosure of all zeros to the function (compared to one point for which $f$ returns close to zero).

Example: The Unit Circle

Consider the function $f(x, y) = \sqrt{x^2 + y^2} - 1$ which is zero exactly on the unit circle. Plugging this in to the standard fsolve we get with $(0.5, 0.5)$ as a starting guess

> x = fsolve (@(x) f(x(1), x(2)), [0.5, 0.5])
x = 0.70711 0.70711

which indeed is close to a zero. But we get no information about other zeros.

Using the interval version of fsolve with $X = [-3, 3] \times [-3, 3]$ as starting domain we get

> [x paving] = fsolve (f, infsup ([-3, -3], [3, 3]));
> x
x ⊂ 2×1 interval vector

     [-1.002, +1.002]
   [-1.0079, +1.0079]

Plotting the paving we get the picture

which indeed is a good enclosure of the unit circle.

How it works

In its simplest form fsolve uses a simple bisection scheme to find the enclosure. Using interval methods we can find enclosure to images of sets. Given a set $X_0 \subset X$ we have three different possibilities
  • $f(X_0) \subset Y$ in which case we add $X_0$ to the paving
  • $f(X_0) \cap Y = \emptyset$ in which case we discard $X_0$
  • Otherwise we bisect $X_0$ and continue on the parts
By setting a tolerance of when to stop bisecting boxes we get the algorithm to terminate in a finite number of steps.


Using bisection is not always very efficient, especially when the domain has many dimensions. One way to speed up the convergence is with what's called contractors. In short a contractor is a function that can take the set $X_0$ and returns a set $X_0' \subset X_0$ with the property that $f(X_0 \setminus X_0') \cap Y = \emptyset$. It's a way of making $X_0$ smaller without having to bisect it that many times.

When you construct a contractor you use the reverse operations definer on intervals. I will not go into how this works, if you are interested you can find more information in the package documentation [1] and in these youtube videos about Set Inversion Via Interval Analysis (SIVIA) [2].

The functions ctc_union and ctc_intersect are used to combine contractors on sets into contractors on unions or intersections of these sets.

Generalization to N-dimensional arrays

How can fsolve be generalized to N-dimensional arrays? The only natural thing to do is to allow for the input and output of $f$ to be N-dimensional arrays. This also is no problem to do. While you mathematically probably would say that fsolve is used to do set inversion for functions $f: \mathbb{R}^n \to \mathbb{R}^m$ it can of course also be used for example on functions $f: \mathbb{R}^{n_1}\times \mathbb{R}^{n_2} \to \mathbb{R}^{m_1}\times \mathbb{R}^{m_2}$.

This is however a bit different when using vectorization. When not using vectorization (and not using contractions) fsolve expects that the functions takes one argument which is an array with each element corresponding to a variable. If vectorization is used it instead assumes that the functions takes one argument for each variable. Every argument is then given as a vector with each element corresponding to one value of the variable for which to compute the function. Here we have no use of N-dimensional arrays.


The only change in functionality that I have done to the functions is to allow for N-dimensional arrays as input and output when vectorization is not used. This required only minor changes, essentially changing expressions like
max (max (wid (interval)))

max (wid (interval)(:))

It was also enough to do these changes in ctc_union and ctc_intersect to have these support N-dimensional arrays.

I have made no functional changes when vectorization is used. I have however made an optimization in the construction of the arguments to the function. The arguments are stored in an array but before being given to the function they need to be split up into the different variables. This is done by creating a cell array with each element being a vector with the values of one of the variables. Previously the construction of this cell array was very inefficient, it split the interval into its lower and upper part and then called the constructor to create an interval again. Now it copies the intervals into the cell without having to call the constructor. This actually seems have been quite a big improvement, using the old version the example with the unit circle from above took around 0.129 seconds and with the new version it takes about 0.092 seconds. This is of course only one benchmark, but a speed up of about 40% for this test is promising!

Lastly I noticed a problem in the example used in the documentation of the function. The function used is

# Solve x1 ^ 2 + x2 ^ 2 = 1 for -3 ≤ x1, x2 ≤ 3 again,
# but now contractions speed up the algorithm.
function [fval, cx1, cx2] = f (y, x1, x2)
  # Forward evaluation
  x1_sqr = x1 .^ 2;
  x2_sqr = x2 .^ 2;
  fval = hypot (x1, x2);
  # Reverse evaluation and contraction
  y = intersect (y, fval);
  # Contract the squares
  x1_sqr = intersect (x1_sqr, y - x2_sqr);
  x2_sqr = intersect (x2_sqr, y - x1_sqr);
  # Contract the parameters
  cx1 = sqrrev (x1_sqr, x1);
  cx2 = sqrrev (x2_sqr, x2);

Do you see the problem? I think it took me more than a day to realize that the problems I was having was not because of a bug in fsolve but because this function computes the wrong thing. The function is supposed to be $f(x_1, x_2) = x_1^2 + x_2^2$ but when calculating the value it calls hypot which is given by $hypot(x_1, x_2) = \sqrt{x_1^2 + x_2^2}$. For $f(x_1, x_2) = 1$, which is used in the example, it gives the same result, but otherwise it will of course not work.


by Joel Dahne ( at July 13, 2017 10:24 AM

July 11, 2017

Michele Ginesi

Gnu Scientific Library

Bessel functions

During the second part of GSoC I have to work on Bessel functions. The bug is related in particular on the Bessel function of first kind and regard the fact that the function is not computed if the $x$ argument is too big.
During this week I studied the actual implementation i.e. the Amos library, written in Fortran. The problem is related to the fact that zbesj.f simply refuse to compute the result when the input argument is too large (the same problem happens for all bessel functions, both for real and complex arguments). What I did was to "unlock" the amos in such a way to compute the value in any case (which seems to be the strategy used by MATLAB). Then I compared the relative errors with the relative errors done by the Gnu Scientific Library (GSL).
At first, I would precise some limitation present in GSL:
  • The parameter $x$ must be real, while it can be complex in Amos and MATLAB. The same holds for the parameter $\alpha$.
  • The class of the output does not adapt to the class of the input, returning always a double.
Doing some tests, it seems that the amos work better in terms of accuracy (even if we are talking about errors which do not differ in the order). I concentrated on values in which Amos usually refuse to compute, since in every other zone of the $x-\alpha$ plane it is known the error is in the order of 1e-15. For $\alpha\in\{-1,0,1\}$ they return the same result, while for other values of $\alpha$, Amos are in general more accurate.
Anyway, I would remark the fact that there are not, as far as I know, accurate algorithms in double precision for very large magnitude. In fact, for such arguments, both Amos and GSL make a relative error of the order of the unity. This problem is evident when using SageMath to compute accurate values, e.g.
sage: bessel_J(1,10^40).n()
sage: bessel_J(1,10^40).n(digits=16)
sage: bessel_J(1,10^40).n(digits=20)
sage: bessel_J(1,10^40).n(digits=25)
sage: bessel_J(1,10^40).n(digits=30)
sage: bessel_J(1,10^40).n(digits=35)
sage: bessel_J(1,10^40).n(digits=40)
The values "stabilize" only when the number of digits is bigger than 30, far away from double precision.

Incomplete gamma function

I also gave a look on how to use GSL to eventually improve the incomplete gamma function. It is not possible, however, to use only GSL functions gsl_sf_gamma_inc_P and gsl_sf_gamma_inc_Q due to their limitations:
  • There is no the "scaled" option.
  • The value of $x$ must be real. This may be not a problem since even MATLAB does not accept complex value (while the version I worked on does).
  • The parameter $x$ must be positive. This is actually a problem of MATLAB compatibility.
  • The class of the output does not adapt to the class of the input, returning always a double.
I tested the accuracy of the GSL and they work better when $a1$ and $x\ll1$ so I think I will fix gammainc.m using the algorithm of gsl_sf_gamma_inc_P and gsl_sf_gamma_inc_Q for that values of the input arguments.

Incomplete Beta function

The actual incomplete beta function needs to be replaced. I've already studied the continued fraction exapnsion (which seems to be the best one to use). In GSL is implemented in a good way but still present two limitations:
  • Does not adapt to the input class (it always work in double precision).
  • There is no "upper" version. This is missing also in the current betainc function, but it is present in MATLAB. Unfortunately, it is not sufficient to compute it as $1 - I_x(a,b)$, it is necessary to find an accurate way to compute directly it.
So I will take inspiration to GSL version of the function but I think I will write beatinc as a single .m file.

by Michele Ginesi ( at July 11, 2017 09:01 AM

July 01, 2017

Michele Ginesi

First period resume

First period resume

First period resume

Here I present a brief resume of the work done in this first month.

1st week: Gammainc

The incomplete gamma function presented a series of inaccurate results in the previous implementation, so it has been decided to rewrite it from scratch. A big part of the work was done by Marco and Nir (here the discussion). I gave my contribution by fixing some problems in the implementation and by making a functioning commit (thanks to the suggestions given by Carne during the OctConf in Geneve).

2nd and 3rd week: Gammaincinv

The inverse of the incomplete gamma function was completely missing in Octave (and this was a problem of compatibility with Matlab). Now the function is present, written as a single .m file. The implementation consist in a simple, but efficient, Newton's method.

Last week: Betainc

When I first submitted my timeline, there was only a bug related to betainc, involving the input validation. During the first period of GSoC emerged a new bug of "incorrect result" type. From this, me and my mentor decided that it is necessary to rewrite the function, so I used this last week to study efficient algorithms to evaluate it. In particular, seems that the most efficient way is to use continude fractions, maybe after a rescale of the function. I'm already working on it, but I will complete it during the first week of august, after finishing the work on Bessel functions.

by Michele Ginesi ( at July 01, 2017 07:58 AM

June 30, 2017

Joel Dahne

One Month In

Now one month of GSoC has passed and so far everything has gone much better than I expected! According to my timeline this week would have been the first of two were I work on vectorization. Instead I have already mostly finished the vectorization and have started to work on other things. In this blog post I'll give a summary of what work I have completed and what I have left to do. I'll structure it according to where the functions are listed in the $INDEX$-file [1]. The number after the heading is the number of functions in that category.

Since this will mainly be a list of which files have been modified and which are left to do this might not be very interesting if you are not familiar with the structure of the interval package.

Interval constant (3)

All of these have been modified to support N-dimensional arrays.

Interval constructor (5)

All of these have been modified to support N-dimensional arrays.

Interval function (most with tightest accuracy) (63)

Almost all of these functions worked out of the box! At least after the API functions to the MPFR and crlibm libraries were fixed, they are further down in the list.

The only function that did not work immediately were $linspace$. Even though this function could be generalized to N-dimensinal arrays the standard Octave function only works for matrices (I think the Matlab version only allows scalars). This means that adding support for N-dimensional vectors for the interval version is not a priority. I might do it later on but it is not necessary.

Interval matrix operation (16)

Most of the matrix functions does not make sense for N-dimensional arrays. For example matrix multiplication and matrix inversion only makes sense for matrices. However all of the reduction functions are also here, they include $dot$, $prod$, $sum$, $sumabs$ and $sumsq$.

At the moment I have implemented support for N-dimensional arrays for $sum$, $sumabs$ and $prod$. The functions $dot$ and $sumsq$ are not ready, I'm waiting to see what happens with bug #51333 [2] before I continue with that work. Depending on the bug I might also have to modify the behaviour of $sum$, $sumabs$ and $prod$ slightly.

Interval comparison (19)

All of these have been modified to support N-dimensional arrays.

Set operation (7)

All of these functions have been modified to support N-dimensional arrays except one, $mince$. The function $mince$ is an interval version of $linspace$ and reasoning here is the same as that for $linspace$ above.

Interval reverse operation (12)

Like the interval functions above, all of the functions worked out of the box!

Interval numeric function (11)

Also these functions worked out of the box, with some small modifications to the documentation for some of them.

Interval input and output (9)

Here there are some functions which require some comments, the ones I do not comment about have all gotten support for N-dimensional arrays.

I think that this function does not make sense to generalize to N-dimensions. It could perhaps take an N-dimensional arrays as input, but it will always return a row vector.  I have left it as it is for now at least.

$disp$ and $display$
These are functions that might be subject to change later on. At the moment it prints N-dimensional arrays of intervals in the same way Octave does for normal arrays. It's however not clear how to handle the $\subset$ symbol and we might decide to change it.

Interval solver or optimizer (5)

The functions $gauss$ and $polyval$ are not generalizable to N-dimensional vectors. I don't think that $fzero$ can be generalized either, for it to work the functions must be real-valued.

The function $fsolve$ can perhaps be modified to support N-dimensional vectors. It uses the SIVIA algorithm [3] and I have to dive deeper into how it works to see if it can be done.

For $fminsearch$ nothing needed to be done, it worked for N-dimensional arrays directly.

Interval contractor arithmetic (2)

Both of these functions are used together with $fsolve$ so they also depend on if SIVIA can be generalized or not.

Verified solver or optimizer (6)

All of these functions work on matrices and cannot be generalized.

Utility function (29)

All of these for which it made sense have been modified to support N-dimensional arrays. Some of them only works for matrices, these are $ctranspose$, $diag$, $transpose$, $tril$ and $triu$. I have left them as they were, though I fixed a bug in $diag$.

API function to the MPFR and crlibm libraries (8)

These are the functions that in general required most work. The ones I have added full support for N-dimensional arrays in are $crlibm\_function$, $mpfr\_function\_d$ and $mpfr\_vector\_sum\_d$. Some of them cannot be generalized, these are $mpfr\_matrix\_mul_d$, $mpfr\_matrix\_sqr\_d$ and $mpfr\_to\_string\_d$. The functions $mpfr\_linspace\_d$ and $mpfr\_vector\_dot\_d$ are related to what I mentioned above for $linspace$ and $dot$.


So summing up the functions that still require some work are
  • Functions related to $fsolve$
  • The functions $dot$ and $sumsq$
  • The functions $linspace$ and $mince$
Especially the functions related to $fsolve$ might take some time to handle. My goal is to dive deeper into this next week.

Apart from this there are also some more things that needs to be considered. The documentation for the package will need to be updated. This includes adding some examples which make use of the new functionality.

The interval package also did not follow the coding style for Octave. All the functions which I have made changes to have been updated with the correct coding style, but many of the functions that worked out of the box still use the old style. It might be that we want to unify the coding standard for all files before the next release.

[1] The $INDEX$ file
[2] Bug #51333
[3] The SIVIA algorithm

by Joel Dahne ( at June 30, 2017 04:38 PM

Enrico Bertino

End of the first work period

Hi all,

this is the end of this first period of GSoC! It was a challenging but very interesting period and I am very excited about the next two months. It is the first time where I have the opportunity to make something that has a real value, albeit small, on someone else's work. Even when I have to do small tasks, like structure the package or learn Tensorflow APIs, I do it with enthusiasm because I have very clear in mind the ultimate goal and the value that my efforts would bring. It may seem trivial, but for a student this is not the daily bread :) I really have fun coding for this project and I hope this will last until the end!

Speaking about the project, I've spent some time wondering which was the best solution to test the correct installation of the Tensorflow Python APIs on the machine. The last solution was putting the test in a function __nnet_init__ and calling it in the PKG_ADD (code in inst/__nnet_init__ in my repo [1]).

Regarding the code, in this last days I tried to connect the dots, calling a Tensorflow network from Octave in a "Matlab compatible" way. In particular, I use the classes that I made two weeks ago in order to implement a basic version of trainNetwork, that is the core function of this package. As I explained in my post of June 12, trainNetwork takes as input the data and two objects: the layers and the options. I had some some difficulty during the implementation of the Layer class due to the inheritance and the overloading. Eventually, I decided to store the layers in a cell array as attribute of the Layer class. Overloading the subsref, I let the user call a specific layer with the '()' access, like a classic array. With this kind of overloading I managed to solve the main problem of this structure, that is the possibility to get a property of a layer doing for example layers(1).Name

classdef Layer < handle
  properties (Access = private)
    layers = {};

  methods (Hidden, Access = {?Layers})
    function this = Layer (varargin)
        nargin = numel(varargin);
        this.layers = cell(1, nargin);
      for i = 1:nargin
        this.layers{i} = varargin{i};

  methods (Hidden)
    function obj = subsref(this, idx)
      switch idx(1).type
          case '()'
              idx(1).type = '{}';
              obj = builtin('subsref',this.layers,idx);
          case '{}'
              error('{} indexing not supported');
          case '.'
            obj = builtin('subsref',this,idx);

    function obj = numel(this)
      obj = builtin('numel',this.layers);

    function obj = size(this)
      obj = builtin('size',this.layers);

Therefore, I implemented the same example of my last post in a proper way. In tests/script you can find the function cnn_linear_model which consists simply in:

1. Loading the datasets

2. Defining layers and options
  layers = [ ...
    imageInputLayer([28 28 1])
  options = trainingOptions('sgdm', 'MaxEpochs', 1);

3. Training
net = trainNetwork(trainImages, trainAngles, layers, options);
4. Prediction
acc = net.predict(testImages, testAngles) 

TrainNetwork is a draft and I have not yet implemented the classe seriesNetwork, but I think it's a good start :) In next weeks I will focus on the Tensorflow backend of the above mentioned functions, with the goal of having a working version at the end of the second period!



by Enrico Bertino ( at June 30, 2017 09:32 AM

June 24, 2017

Enrico Bertino

Package structure

Hello! During the first period of GSoC I have worked mostly on analyzing the Matlab structure of the net package in order to guarantee the compatibility throughout the whole project. The focus of the project is on the convolutional neural networks, about which I will write the next post.

Regarding the package structure, the core will be composed by three parts:

  1. Layers: there are 11 types of layers that I defined as Octave classes, using classdef. These layers can be concatenated in order to create a Layer object defining the architecture of the network. This will be the input for the training function. 
  2. Training: the core of the project is the training function, which takes as input the data, the layers and some options and returns the network as output. 
  3. Network: the network object has three methods (activations, classify and predict) that let the user compute the final classification and prediction. 

Figure 1: conv nnet flowchart 

I have already implemented a draft for the first point, the layers classes [1]. Every layer type inherits some attributes and methods from the parent class Layers. This is useful for creating the Layer object: the concatenation of different layers is always a Layer object that will be used as input for the training function. For this purpose, I overloaded the cat, horzcat and vertcat operators for Layers and subsref for Layer. I need to finalize some details for the disp methods of these classes.

Figure 2: Layers classes definitions 

 is used in every class for the parameter management and the attributes setting.

The objects of these classes can be instantiated with a corresponding function, implemented in the directory inst/. Here an example for creating a Layer object 

> a = imageInputLayer([2,2,3]); # first layer
> b = convolution2dLayer(1,1); # second layer
> c = dropoutLayer(1); # third layer
> layers = [a b c]; # Layer object from layers concat
> drop = layers(3); # Layer element access
> drop.Probability # Access layer attribute
ans = 0.50000

All functions can be tested with the make check of the package.

The next step is a focus on the Tensorflow integration, with Pytave, writing a complete test for a regression of images angles and comparing the precision and the computational time with Matlab.


p.p1 {margin: 0.0px 0.0px 0.0px 0.0px; font: 14.0px Monaco; color: #f4f4f4; background-color: #000000; background-color: rgba(0, 0, 0, 0.95)} span.s1 {font-variant-ligatures: no-common-ligatures} p.p1 {margin: 0.0px 0.0px 0.0px 0.0px; font: 14.0px Monaco; color: #f4f4f4; background-color: #000000; background-color: rgba(0, 0, 0, 0.95)} span.s1 {font-variant-ligatures: no-common-ligatures}

by Enrico Bertino ( at June 24, 2017 08:55 PM

Train a Convolutional Neural Network for Regression


I spent the last period working mostly on Tensorflow, studying the APIs and writing some examples in order to explore the possible implementations of neural networks. For this goal, I chose an interesting example proposed in the Matlab examples at [1]. The dataset is composed by 5000 images, rotated by an angle α, and a corresponding integer label (the rotation angle α). The goal is to make a regression to predict the angle of a rotated image and straighten it up.
All files can be found in tests/examples/cnn_linear_model in my repo [2].

I have kept the structure as in the Matlab example, but I generated a new dataset starting from LeCun's MNIST digits (datasets at [3]). Each image was rotated by a random angle between 0° and 70°, in order to keep the right orientation of the digits (code in dataset_generation.m). In Fig. 1 some rotated digits with the corresponding original digits.

Figure 1. rotated images in columns 1,3,5 and originals in columns 2,4,6

The implemented linear model is:
$ \hat{Y} = \omega X + b $,
where the weights $\omega$ and the bias $b$ will be optimized during the training minimizing a loss function. As loss function, I used the mean square error (MSE):
$  \dfrac{1}{n} \sum_{i=1}^n (\hat{Y_i} - Y_i)^2 $,
where the $Y_i$ are the training labels. 

In order to show the effective improvement given by a Neural Network, I started to make a simple regression feeding the X variable of the model directly with the 28x28 images. Even if for the MSE minimization a close form exists, I implemented an iterative method for discovering some Tensorflow features (code in For evaluate the accuracy of the regression, I consider a correct regression if the difference between angles is less than 20°. After 20 epochs, the convergence was almost reached, giving an accuracy of $0.6146$.

Figure 2. rotated images in columns 1,3,5 and after the regression in columns 2,4,6

I want to analyze now the improvement given by a feature extraction performed with a convolutional neural network (CNN). As in the Matlab example, I used a basic CNN since the input images are quite simple (only numbers with monochromatic background) and consequently the features to extract are few.
  • INPUT [28x28x1] will hold the raw pixel values of the image, in this case an image of width 28, height 28
  • CONV layer will compute the output of neurons that are connected to local regions in the input, each computing a dot product between their weights and a small region they are connected to in the input volume. This results in volume such as [12x12x25]: 25 filters of size 12x12
  • RELU layer will apply an element-wise activation function, such as the $max(0,x)$ thresholding at zero. This leaves the size of the volume unchanged ([12x12x25]).
  • FC (i.e. fully-connected) layer will compute the class scores, resulting in volume of size [1x1x1], which corresponds to the rotation angle. As with ordinary Neural Networks, each neuron in this layer will be connected to all the numbers in the previous volume.
Figure 3. CNN linear model architecture

We can visualize the architecture with Tensorboard where the graph of the model is represented.

Figure 4. Model graph generated with Tensorboard

With the implementation in, the results are quite satisfying: after 15 epochs, it reached an accuracy of $0.75$ (205 seconds overall). One can see in Fig. 4 the marked improvement of the regression.

Figure 5. rotated images in columns 1,3,5 and after the CNN regression in columns 2,4,6

With the same parameters, Matlab reached an accuracy of $0.76$ in 370 seconds (code in regression_Matlab_nnet.m), so the performances are quite promising

In the next post (in few days), I will integrate the work done up to now, calling the Python class within Octave and making a function that simulates the behavior of Matlab. Leveraging the layers classes that I made 2 weeks ago, I will implement a draft of the functions trainNetwork and predict making the Matlab script callable also in Octave.

I will also care about the dependencies of the package: I will add the dependency from Pytave in the package description and write a test as PKG_ADD in order to verify the version of Tensorflow during the installation of the package.


by Enrico Bertino ( at June 24, 2017 02:36 PM

June 22, 2017

Joel Dahne

Vectorization and broadcasting

At the moment I'm actually ahead of my schedule and this week I started to work on support for vectorization on N-dimensional arrays. The by far biggest challenge was to implement proper broadcasting and most of this post will be devoted to going through that. At the end I also mention some of the other things I have done during the week.

Broadcasting arrays

At the moment I have implement support for broadcasting on all binary functions. Since all binary functions behave similarly in respect to broadcasting I will use $+$ in all my example below, but this could in principle be any binary function working on intervals.

When adding to arrays, $A, B$, of the same size the result is just an arrays of the same size with each entry containing the sum of the corresponding entries in $A$ and $B$. If $A$ and $B$ does not have the same size then we try to perform broadcasting. The simplest form of broadcasting is when $A$ is an arrays and $B$ is a scalar. Then we just take the value of $B$ and add to every element in $A$. For example

> A = infsupdec ([1, 2; 3, 4])
A = 2×2 interval matrix
   [1]_com   [2]_com
   [3]_com   [4]_com
> B = infsupdec (5)
B = [5]_com
> A + B
ans = 2×2 interval matrix
   [6]_com   [7]_com
   [8]_com   [9]_com

However it is not only when one of the inputs is a scalar that broadcasting can be performed. Broadcasting is performed separately for each dimension of the input. We require either that the dimensions are equal, and no broadcasting is performed, or that one of the inputs have that dimension equal to $1$, we then concatenate this input along that dimension until they are of equal size. If for example $A$ has dimension $4\times4\times4$ and $B$ dimension $4\times4\times1$ we concatenate $B$ with itself along the third dimension four times to get two arrays of the same size. Since a scalar has all dimensions equal to 1 we see that it can be broadcasted to any size. Both $A$ and $B$ can also be broadcasted at the same time, along different dimensions, for example

> A = infsupdec (ones (1, 5, 2))
A = 1×5×2 interval array
ans(:,:,1) =
   [1]_com   [1]_com   [1]_com   [1]_com   [1]_com
ans(:,:,2) =
   [1]_com   [1]_com   [1]_com   [1]_com   [1]_com
> B = infsupdec ([1, 2, 3, 4, 5; 6, 7, 8, 9, 10])
B = 2×5 interval matrix
   [1]_com   [2]_com   [3]_com   [4]_com    [5]_com
   [6]_com   [7]_com   [8]_com   [9]_com   [10]_com
> A + B
ans = 2×5×2 interval array
ans(:,:,1) =
   [2]_com   [3]_com   [4]_com    [5]_com    [6]_com
   [7]_com   [8]_com   [9]_com   [10]_com   [11]_com
ans(:,:,2) =
   [2]_com   [3]_com   [4]_com    [5]_com    [6]_com
   [7]_com   [8]_com   [9]_com   [10]_com   [11]_com

The implementation

I'll go through a little bit about my implementation. I warn you that I'm not that familiar with the internals of Octave so some things I say might be wrong, or at least not totally correct.

Internally all, numerical, arrays are stored as a linear vector and the dimensions are only metadata. This means that the most efficient way to walk through an array is with a linearly increasing index. When $A$ and $B$ have the same size the most efficient way to sum them is to linearly go through the arrays. In pseudo code

// Calculate C = A + B
for (int i = 0; i < numel (A); i++) {
  C(i) = A(i) + B(i);

This works fine, and apart from unrolling the loop or doing optimizations like that it is probably the most efficient way to do it.

If $A$ and $B$ are not of the same size then one way to do it would be to simply extend $A$ or/and $B$ along the needed dimensions. This would however require coping a lot of data, something we want to avoid (memory access is expensive). Instead we try to be smart with our indexing to access the right data from both $A$ and $B$.

After asking on the IRC-channel I got pointed to this Octave function which performs broadcasting. My implementation, which can be found here, is heavily inspired by that function.


Here I compare the performance of the new implementation with the old one. Since the old one could only handle matrices we are limited by that. We can measure the time it takes to add two matrices $A$, $B$ with the code

tic; A + B; toc;

We do 10 runs for each test and all times are in seconds.

Addition of large matrices

Case 1: A = B = infsupdec (ones (1000, 1000));
       Old         New
       0.324722    0.277179
       0.320914    0.276116
       0.322018    0.276075
       0.318713    0.279258
       0.332041    0.279593
       0.318429    0.279987
       0.323752    0.279089
       0.317823    0.276036
       0.320509    0.280964
       0.320610    0.281123
Mean:  0.32195     0.27854
Case 2: A = B = infsupdec (ones (10, 100000));
        Old         New
        0.299321    0.272691
        0.297020    0.282591
        0.296460    0.274298
        0.294541    0.279661
        0.298306    0.277274
        0.301532    0.275531
        0.298163    0.278576
        0.298954    0.279868
        0.302849    0.275991
        0.297765    0.278806
Mean:   0.29849    0.27753

Case 3: A = B = infsupdec (ones (100000, 10));
        Old         New
        0.286433    0.279107
        0.289503    0.278251
        0.297562    0.279579
        0.292759    0.283311
        0.292983    0.281306
        0.290947    0.282310
        0.293025    0.286172
        0.294153    0.278886
        0.293457    0.278625
        0.296661    0.280804
Mean:   0.29275     0.28084

Broadcasting scalars

Case 4: A = infsupdec (ones (1000, 1000));
             B = infsupdec (1);
        Old         New
        0.298695    0.292419
        0.298158    0.292274
        0.305242    0.296036
        0.295867    0.291311
        0.296971    0.297255
        0.304297    0.292871
        0.298172    0.300329
        0.297251    0.291668
        0.299236    0.294128
        0.300457    0.298005
Mean;   0.29943     0.29463

Case 5: A = infsupdec (1);
             B = infsupdec (ones (1000, 1000));
         Old         New
        0.317276    0.291100
        0.316858    0.296519
        0.316617    0.292958
        0.316159    0.299662
        0.317939    0.301558
        0.322162    0.295338
        0.321277    0.293561
        0.314640    0.291500
        0.317211    0.295487
        0.317177    0.294376
Mean:   0.31773     0.29521

Broadcasting vectors

Case 6: A = infsupdec (ones (1000, 1000));
             B = infsupdec (ones (1000, 1));
        Old         New
        0.299546    0.284229
        0.301177    0.284458
        0.300725    0.276269
        0.299368    0.276957
        0.303953    0.278034
        0.300894    0.275058
        0.301776    0.276692
        0.302462    0.282946
        0.304010    0.275573
        0.301196    0.273109
Mean:   0.30151     0.27833

Case 7: A = infsupdec (ones (1000, 1000));
             B = infsupdec (ones (1, 1000));
         Old         New
        0.300554    0.295892
        0.301361    0.294287
        0.302575    0.299116
        0.304808    0.294184
        0.306700    0.291606
        0.301233    0.298059
        0.301591    0.292777
        0.302998    0.290288
        0.300452    0.291975
        0.305531    0.290178
Mean:   0.30278     0.29384

We see that in all cases the new version is faster or at least equally fast as the old version. In the old version the order of the input made a slight difference in performance (case 4 vs case 5). In the new version both inputs are treated in exactly the same way so we no longer see that difference.

Possible improvements

In theory the cases when we broadcast a scalar could be the fastest ones. If $B$ is a scalar we could, in pseudo code, do something similar to
// Calculate C = A + B with B scalar
for (int i = 0; i < numel (A); i++) {
  C(i) = A(i) + B;

This is however not implemented at the moment. Instead we use the ordinary routine to calculate the index for $B$ (since it is a scalar it will always evaluate to $1$). If we would like to optimize more for this case we could add a check for if $A$ or $B$ are scalars and then optimize for that. Of course this would also make the code more complicated, something to watch out for. At the moment I leave it like this but if we later want to optimize for that case it could be done.

Other work

Apart from the work to fix the broadcasting for binary functions there were very little to do for many of the functions. All binary functions that use this code, and all unary functions using an even simpler code, worked directly after fixing the oct-files. Some of them required small changes to the documentation but other than that the octave-scripts were fine. So mainly it has been a matter of actually going through all files and check that they actually did work.

Bug #51283

When going through all the functions I noticed a bug in the interval version of $\sin$,

 > sin (infsupdec (0))
ans = [0]_com
> sin (infsupdec ([0, 0]))
ans = 1×2 interval vector
   [0, -0]_com   [0, -0]_com

The second version here is wrong, $-0$ should never be allowed as a value for the supremum of an interval. I was able to track this down to how Octaves $\max$ function works, see bug #51283. As Oliver writes there the exact behaviour of the $\max$-function is not specified in IEEE Std 754-2008 so we cannot rely on that. To solve this I have added a line to manually set all $-0$ to $+0$ in the supremum of the interval.


by Joel Dahne ( at June 22, 2017 05:12 PM

June 20, 2017

Michele Ginesi

Timetable: modification

Timetable: modification

Timetable: modification

According to my timetable (that you can find here), during this last week of June, I should've work on the input validation of betainc. Since a new bug related to this function has been found and, moreover, the actual implementation doesn't accept the "lower" or "upper" tail (as MATLAB do), me and my mentor decided to use this week to start studying how to rewrite betainc (main references will be [1] and [2]) and to use the last part fo the GSoC to actually implement it. In this way, my timetable remain almost identical (I will use July to work on Bessel functions) and I will be able to fix also this problem.

[1] Abramowitz, Stegun "Handbook of Mathematical Functions"
[2] Cuyt, Brevik Petersen, Vendonk, Waadeland "Handbook of Continued Fractions for Special Functions"

by Michele Ginesi ( at June 20, 2017 07:01 AM

June 16, 2017

Joel Dahne

Construction and Printing

This week I have started to work on methods for constructing and printing N-dimensional arrays of intervals. In my timeline I estimated that this work would take 2 weeks. However in this first week I have managed to complete most of the work. I will give some comments on how I have worked with the Mercurial repository, how the work went and different things I encountered along the path.

Working with Mercurial

This is essentially the first time I'm using Mercurial for revision control, though I have used git before. However I quickly learned how to use it for the basic tasks that I need, committing, comparing files and checking the history. As mentioned in a previous post you can find my repository here [1].

Coding style

When I started to work with the files I realized that they did not follow Octaves coding standard [2]. After a short discussion on the mailing list we decided that I will update the files I change to follow the standard coding style. Usually it is not a good idea to change coding style and add functionality in the same commit. However most of the changes to coding style are only white space changes so they can be ignored using the -w flag in Mercurial. Thus we decided that as long as the coding style changes are only such that it is ignored with -w I will do it in the same commit as the added functionality. If there are some coding style changes that's not only white space, the most common example is to long lines, I do a commit with only changes to the coding style first. So if you want to take a look at the functionality I have added you will probably want to use the -w flag. Note however that I have not updated the coding style for any files I have not changed otherwise.


Normally I do one commit for each file, though in many cases the bare intervals and the decorated intervals have almost identical functions and in that case I commit changes to them both at the same time. Of course it also happens that I have to go back and do more changes to a files, in that case I just do another commit.

The actual work

The work went much faster than I expected. The main reason for this is that Octave has very good support for indexing. For example expressions like

isnai(x.inf <= x.sup) = false;

works just as well for matrices as for N-dimensional arrays. In fact the constructor for bare intervals even worked for N-dimensional arrays from the beginning, there I only had to do slight modification to the documentation and add some tests!

Not all functions were that easy though. Some functions that have not been updated in a while clearly assumed the input was a matrix, for example in $hull$

sizes1 = cellfun ("size", l, 1);
sizes2 = cellfun ("size", l, 2);

In most cases I only needed to add more general indexing, often times even making the code clearer.

In some functions all I had to do was to remove the check on the input data so that it would accept N-dimensional arrays. This was true in for example $cat$ were all I had to do was to remove the check and do some minor modifications to the documentation.

I can conclude with saying that Octave has great support for working with N-dimensional arrays. Since internally the data for intervals are stored only as arrays I could make good use of it!

Noteworthy things

While most functions were straight forward to modify some required some thought. How should they even work for N-dimensional input?


When modifying the $disp$-function I chose to mimic how Octave handles displaying N-dimensional arrays. I noticed that this is different from how Matlab handles it. In Matlab we have

> x = zeros (2, 2, 2)

x(:,:,1) =

     0     0
     0     0

x(:,:,2) =

     0     0
     0     0

while in Octave it's

> x = zeros (2, 2, 2)
x =

ans(:,:,1) =

   0   0
   0   0

ans(:,:,2) =

   0   0
   0   0

I don't know the choice behind Octaves version. At least at first glance I think I prefer the way Matlab does it. But since I'm working in Octave I chose that style.

The next question was how to handle the subset symbol, $\subset$. The interval package uses $=$ or $\subset$ depending on if the string representation is exact or not. For example

> x = infsup (1/2048, 1 + 1/2048);
> format short; x
x ⊂ [0.00048828, 1.0005]
> format long; x
x = [0.00048828125, 1.00048828125]

How should this be handled for N-dimensional arrays? One way would be to switch all $=$ to $\subset$ is the representation is not exact. Another to use $\subset$ on all submatrices that does not have an exact string representation. The third way, and how it is implemented now, is to only change the first $=$ to $\subset$, the one after the variable name. Like this

> x(1,1,1:2) = infsup (1/2048, 1 + 1/2048)
x ⊂ 1×1×2 interval array

ans(:,:,1) =   [0.00048828, 1.0005]
ans(:,:,2) =   [0.00048828, 1.0005]

This might be a bit odd when you first look at it, on some places we use $=$ and on some $\subset$. Though I think it somehow makes sense, we are saying that $x$ is a subset of the $1\times1\times2$ interval array given by

ans(:,:,1) =   [0.00048828, 1.0005]
ans(:,:,2) =   [0.00048828, 1.0005]

which actually is true. Anyway I will leave like this for now and then we might decide to switch it up later.

linspace and mince

The standard implementation of $linspace$ only supports scalar or vector input. It could be generalized to N-dimensional arrays by for example returning a N+1-dimensional array were the last dimension corresponds to the linearly spaced elements. But since this has not been done in the standard implementation I will at least wait with adding for intervals.

The function $mince$ can be seen as a interval generalization of $linspace$. It  takes an interval and returns an array of intervals whose union cover it. This could similarly be expanded to N dimensions by creating the array along the N+1 dimension. But again we choose to at least wait with adding this.

meshgrid and ndgrid

The interval package already has an implementation of $meshgrid$. But since it previously did not support 3-dimensional arrays it had to fit 3-d data in a 2-d matrix. Now that it supports 3-d data it can output that instead.

Currently the interval package does not implement $ndgrid$. When I looked into it I realized that the standard implementation of $ndgrid$ actually works for interval arrays as well. I have not looked into the internals but in principle it should only need the $cat$ function, which is implemented for intervals. Further I noticed that the standard $meshgrid$ also works for intervals. However the interval implementation differs in that it converts all input to intervals, were as the standard implementation allows for non-uniform output. Using the interval implementation of $meshgrid$ we have

> [X Y] = meshgrid (infsup (1:3), 4:6)
X = 3×3 interval matrix

   [1]   [2]   [3]
   [1]   [2]   [3]
   [1]   [2]   [3]

Y = 3×3 interval matrix

   [4]   [4]   [4]
   [5]   [5]   [5]
   [6]   [6]   [6]

but if we fall back to the standard implementation (by removing the interval implementation) we get

> [X Y] = meshgrid (infsup (1:3), 4:6)
X = 3×3 interval matrix

   [1]   [2]   [3]
   [1]   [2]   [3]
   [1]   [2]   [3]

Y =

   4   4   4
   5   5   5
   6   6   6

Note that the last matrix is not an interval matrix.

So the question is, should we implement a version of $ndgrid$ that converts everything to intervals or should we remove the implementation of $meshgrid$? It's at least most likely not a good idea that the functions are different. I think that removing the implementation of $meshgrid$ makes most sense. First of all it's less code to maintain, which is always nice. Secondly you can manually convert all input to the function to intervals if you want uniform output. If you do not want uniform output then the standard implementation works were as the interval implementation does not, so the standard implementation is more general in a sense.

We have to choose what to do, but for now I leave it as it is.

Non-generalizable functions

From what I have found there is no way to create a 3-dimensional array in Octave in the same way you can create a 2-dimensional one with for example

M = [1, 2; 3, 4];

Instead higher dimensional arrays have to be created using other functions, for example $reshape$ or $zeros$, or by specifying the submatrices directly

M(:,:,1) = [1, 2; 3, 4];
M(:,:,2) = [5, 6; 7, 8];

This means that the functions $\_\_split\_interval\_literals\_\_$, which is used to split a string like $"[1, 2; 3, 4]"$ into its separate components, cannot really be generalized to N dimensions.


by Joel Dahne ( at June 16, 2017 02:07 PM

Michele Ginesi


The gammaincinv function

The inverse of the incomplete gamma function


Given two real values $y\in[0,1]$ and $a>0$, gammaincinv(y,a,tail) gives the value $x$ such that \[P(x,a) = y,\quad\text{ or }\quad Q(x,a) = y\] depending on the value of tail ("lower" or "upper").


The computation is carried out via standard Newton method, solving, for $x$ \[ y - P(x,a) = 0 \quad\text{ or }\quad y - Q(x,a) = 0. \] The Newton method uses as exit criterion a control on the residual and on the maximum number of iterations.
For numerical stability, the initial guess and which gammainc function is inverted depends on the values of $y$ and $a$. We use the following property: \[\text{gammaincinv}(y,a,\text{`upper'})=\text{gammaincinv}(1-y,a,\text{`lower'})\] which follows immediately by the property \[P(x,a) + Q(x,a) = 1.\] We start computing the trivial values, i.e. $a=1,y=0$ and $y=1$. Then we rename the variables as follows (in order to obtain a more classical notation): if tail = "lower", then $p=y,q=1-y$; while if tail = "upper", $q=y,p=1-y$. We will invert with $p$ as $y$ in the cases in which we invert the lower incomplete gamma function, and with $q$ as $y$ in the cases in which we invert the upper one.
Now, the initial guesses (and the relative version of the incomplete gamma function to invert) are the following:
  1. $p\dfrac{(0.2(1+a))^a}{\Gamma(1+a)}$. We define \[r = (p\,\Gamma(1+a))^{\frac{1}{a}}\] and the parameters \begin{align*} c_2 &= \dfrac{1}{a+1}\\ c_3 &= \dfrac{3a+5}{2(a+1)^2(a+2)}\\ c_4 &= \dfrac{8a^2+33a+31}{3(a+1)^3(a+2)(a+3)}\\ c_5 &= \dfrac{125a^4+1179a^3+3971a^2+5661a+2888}{24(1+a)^4(a+2)^2(a+3)(a+4)}\\ \end{align*} Then the initial guess is \[x_0 = r + c_2\,r^2+c_3\,r^3+c_4\,r^4+c_5\,r^5\] and we invert the lower incomplete gamma function.
  2. $q\dfrac{e^{-a/2}}{\Gamma(a+1)}, a\in(0,10)$. The initial guess is \[x_0 = - \log (q) - \log(\Gamma(a))\] and we invert the upper incomplete gamma function.
  3. $a\in(0,1)$. The initial guess is \[x_0 = (p\Gamma(a+1))^{\frac{1}{a}}\] and we invert the lower incomplete gamma function.
  4. Other cases. The initial guess is \[x_0 = a\,t^3\] where \[t = 1-d-\sqrt{d}\,\text{norminv}(q),\quad d=\dfrac{1}{9a}\] (where norminv is the inverse of the normal distribution function) and we invert the upper incomplete gamma function.
The first three initial guesses can be found in the article "Efficient and accurate algorithms for the computation and inversion of the incomplete gamma function rations" by Gil, Segura ad Temme; while the last one can be found in the documentation of the function igami.c of the cephes library.
You can find my work on my repository on bitbucket here, bookmark "gammainc".


Working on this function I found some problems in gammainc. To solve them, I just changed the flag for the condition 5 and added separate functions to treat Inf cases.

by Michele Ginesi ( at June 16, 2017 10:39 AM

June 06, 2017

Michele Ginesi


The gammainc function

The incomplete gamma function


There are actually four types of incomplete gamma functions:
  • lower (which is the default one) $$P(a,x) = \frac{1}{\Gamma(a)}\int_{0}^{x}e^{-t}t^{a-1}$$
  • upper$$Q(a,x) = \frac{1}{\Gamma(a)}\int_{x}^{+\infty}t^{a-1}e^{-t}dt$$
  • scaled lower and scaled upper which return the corresponding incomplete gamma function scaled by $$\Gamma(a+1)\frac{e^x}{x^a}$$


To compute the incomplete gamma function, there is a total of seven strategies depending on the position of the arguments in the $x-a$ plane.
  1. $x=0,a=0$. In this case, the output is given as follows:
    • for "lower" and "scaledlower" is $1$
    • for "upper" and "scaledupper" is $0$
  2. $x=0,a\neq0$. In this case:
    • for "lower" is $0$
    • for "upper" and "scalelower" is $1$
    • for "scaledupper" is $+\infty$
  3. $x\neq0,a=0$. In this case:
    • for "lower" is $1$
    • for "scalelower" is $e^x$
    • for "upper" and "scaledupper" is $0$
  4. $x\neq0,a=1$. In this case:
    • for "lower" is $1-e^x$
    • for "upper" is $e^{-x}$
    • for "scaledlower" is $\frac{e^x-1}{x}$
    • for "scaledupper" is $\frac{1}{x}$
  5. $0 x\leq36,a\in\{2,3,4,\ldots,17,18\}$. In this case we used the following expansion (which holds true only for $n\in\mathbb{N}$): denote $\gamma(a,x):=\Gamma(a)P(a,x)$, then $$\gamma(n,x) = (n-1)!\left( 1-e^{-x}\sum_{k=0}^{n-1} \frac{x^k}{k!}\right).$$
  6. $x+0.25 a $ or $x0$, or $x$ not real or $|x|1$ or $a5$. Here we used the expansion $$\gamma(a,x) = e^{-x} x^a \sum_{k=0}^{+\infty} \frac{\Gamma(a)}{\Gamma(a+1+k)}x^k.$$
  7. for all the remaining cases was used a continued fraction expansion. Call $\Gamma(a,x) = \Gamma(a) Q(a,x)$, the expansion is $$\Gamma(a,x) = e^{-x}x^a\left( \frac{1}{x+1-a-}\frac{1(1-a)}{x+3-a-}\frac{2(2-a)}{x+5-a-}\cdots \right).$$ This case is handled by an external .cc function.

My contribution

Actually I dind't participate from the start in the (re)making of this function: most of the work done on it is due to Nir and Marco (see the discussion on the bug tracker).
My contribution was, before the GSoC, to adapt the codes in such a way to make the type of the output (and the tolerances inside the algorithms) coherent on the type of the input (single, int or double). Then I corrected few small bugs and made a patch helped by Carne during the OctConf in Geneve.
During this first week of GSoC I fixed the input validation and added tests to this function, finding some problem in the implementation. Most of them depended on the fact that the continued fractions were used in cases in which they are not optimal. To solve these problems I changed the conditions to use the series expansion instead of the continued fractions.
Now gammainc works properly also with complex argument (for the $x$ variable), while Matlab doesn't accept non-real input.
You can find my work on my repository on bitbucket here, bookmark "gammainc". I will work on this same bookmark the next two weeks while I will work on gammaincinv.

by Michele Ginesi ( at June 06, 2017 09:09 AM

Carnë Draug

Function help text in scripts and command line

I guess I'm one of the few [[citation needed]] that writes Octave programs instead of only things meant to run from the Octave console. You end up with a few hundred lines of Octave and tens of functions in a single file. For debugging you just source the whole thing.

I always thought that Octave's help function only handled functions in their own m files, the help text being the first or second block of comments that does not start with Author or Copyright. But I was wrong. You can document functions in the Octave command line and you can also document functions in a script. Octave seems to pick a block of comments before a function definition as that function help text. Like so:

octave-cli:1> ## This is the help text of the foo function.
octave-cli:1> function x = foo (y, z), x = y+z; endfunction
octave-cli:2> ## This is the help text of the bar function.
octave-cli:2> ##
octave-cli:2> ## Multi-line help text are not really a problem,
octave-cli:2> ## just carry on.
octave-cli:2> function x = bar (y, z), x = y+z; endfunction
octave-cli:3> man foo
 'foo' is a command-line function

 This is the help text of the foo function.

octave-cli:4> man bar
 'bar' is a command-line function

 This is the help text of the bar function.

 Multi-line help text are not really a problem,
 just carry on.

It also works in script file with one exception: the first function of the file which picks the first comment which will be the shebang line:

$ cat > an-octave-script << END
> #!/usr/bin/env octave
> ## This is the help text of the foo function.
> function x = foo (y, z)
>    x = y+z;
> endfunction
> ## This is the help text of the bar function.
> ##
> ## Multi-line help text are not really a problem,
> ## just carry on.
> function x = bar (y, z)
>   x = y+z;
> endfunction
$ octave-cli
octave-cli:1> source ("an-octave-script")
octave-cli:2> man bar
'bar' is a command-line function

This is the help text of the bar function.

Multi-line help text are not really a problem,
just carry on.

octave-cli:3> man foo
'foo' is a command-line function

!/usr/bin/env octave

I reported a bug about it (bug #51191) but I was pleasantly surprised that functions help text were being identified at all.

June 06, 2017 12:00 AM

June 01, 2017

Carnë Draug

My Free software activities in May 2017


Bugs and patches:

  • Fixed (bug #47115) where label2rgb in the Image package would not work with labelled images of type uint32 and uint64. Reason was that behind the scenes ind2rgb from the Octave core was being used, and Octave imitates Matlab in that labelled images can only be up to uint16. That sounds like a silly limitation so I dropped it from core in ind2rgb and ind2gray. However, there's other places where that limitation is in place and needs to be lifted such as cmunique, cmpermute, and the display functions.
  • Fixed hist3 handling of special cases with NaN values on the rows that define the histogram edges and centring values in an histogram when there's only one value. Started as side effect of fixing (bug #51059) .
  • (bug #44799) - document imrotate does not work with the spline method.
  • Fix input check for fspecial e65762ee9414
  • Review wiener2, function for adaptive noise reduction, from Hartmut and extended it for N dimensions.
  • Fixed mean2 Matlab incompatibilities for some weird cases (bug #51144)
  • Make imcast support class logical.
  • Review otsuthresh, function to compute threshold value of an histogram via Otsu's method, by Avinoam Kalma
  • Review affine2 and affine3, classes to represent 2d and 3d affine transforms, submitted by Motherboard and Avinoam . Started generalising code to N dimensions to avoid code duplication. Only one method missing (maybe you wanna help?).

As sysadmin:

  • Finished updating the wiki for Octave, due to issues with pygments.


  • Started packaging libsystem-info-perl, got stuck due to a licensing issue

June 01, 2017 12:00 AM

May 31, 2017

Mike Miller

Blog Relaunch, More Pythonic With Pelican

I just finished converting this blog to Pelican, a static site generator built entirely on Python. The conversion process was so easy, and Pelican is such a pleasure to work with, that I wish I had done this much sooner instead of letting this blog languish for over a year and a half.

To reboot the blog and kick things off again, I'll describe just how easy it was to convert to Pelican, compare it to what I was working with before, and lay out some ideas for the future.

Starting with Pelican

So Pelican is a static site generator, which means it builds a complete website from a template and a bunch of files representing content. The content is usually written in a plain text markup language like Markdown. Pelican has explicit support for common blog features like articles, pages, categories, comments, and syndication. The blog was already written in Markdown and built using a static site generator, and I wanted to try redoing it with Pelican.

Creating a Pelican site is simple, just run pelican-quickstart, fill in some content, and run pelican to compile it. I only had to copy in the existing Markdown content and massage the metadata headers into Pelican's format. I also wanted to preserve the URL scheme I had already been using, and found that it was trivial to tweak the default config files to do exactly that. The result was a new site that mostly replicated the content and URLs of the existing blog.

I noticed a small difference in the rules that Pelican uses for transforming a page title into a URL. There is one post that contains a "+" plus sign. My previous blog software converted this into the word "plus", while Pelican drops it by default. Fortunately, it's simple to add custom URL transformation rules, and I got this post where it needed to be. If generic transformation rules don't work, Pelican also allows an article to specify exactly what its URL should be in its metadata.

I also wanted to be able to add a delimiter in articles to mark the beginning of the article as a summary. By default, Pelican pulls out the first 50 words as the summary, probably cutting off in mid-sentence. I wanted more control, without having to duplicate the summary in a header field. Fortunately, there's a plugin for that. The pelican-plugins repository has dozens of plugins, including one that recognizes a delimiter comment marking the end of the page summary.

At this point, setting aside appearance, the old blog was converted perfectly, preserving all of the content, URLs, summaries, feeds, tags, and comments. I had wanted to know if I could use Pelican to publish my blog, and in about an hour I had answered that easily, and it was basically ready to use. Of course the next question was, what else can I do?

Ditching Octopress

I'll take a small detour here and try to explain why I did this conversion in the first place.

When I started this blog three years ago, I wanted to use a static site generator. I settled on Octopress at the time because it seemed simple and included all of the features I wanted. Octopress is based on Jekyll, which was and is probably the most popular static site generator. Octopress looked like Jekyll with a bunch of integrated plugins and theming problems already solved. It looked like a quick way to get started with Jekyll without having to learn Ruby and mess around with CSS, just write and publish. And yes, it was pretty easy.

One problem that was quickly apparent is that Octopress was basically a preconfigured Jekyll project. To customize Octopress for your own blog, you essentially have to fork Octopress and make changes to it. There was no way to run upstream Octopress against your own configuration and content.

Another problem that's a direct corollary, Octopress is based on a now-ancient version of Jekyll. It's impossible to take advantage of any Jekyll improvements over the years because Octopress used Jekyll 0.12 from 2012. If Octopress updates to a newer version, the only way to take advantage would be to rebase your local fork onto Octopress upstream changes and hope for the best.

Not only are configuration changes made in the Octopress repository itself, but content is put there by default. So a typical Octopress blog is a fork of the Octopress project, with local config files, and with all articles and pages mixed in there too. I actually created a subrepo for just the content, to try to maintain the content separately from the Octopress repo. This seemed better at first, but I still had config changes in one repo and content in the other, and they're not exactly independent. Add in a customized fork of a theme, and that's three separate repositories to maintain one simple blog.

I saw that Octopress had been working on a redesign to fix these problems and work with the latest Jekyll, but it seems like development has stalled. Maybe Jekyll is good enough now that Octopress isn't really needed anymore? Maybe the authors moved on and the Octopress community dissolved? Whatever the state of the project, I had a blog built with a tool that was cumbersome to maintain, wasn't being updated, and I didn't really understand. I stopped updating my blog due to a confluence of many life events, and being difficult to maintain wasn't helping.

I was at PyCon 2017 last week, and one of the strongest messages that I took with me was to try to find a Pythonic way to do everything. Since Python can do so many things and make them so easy, why not do more things in Python? I made a note to find and try out a best-practices Pythonic static site generator to reboot my blog, and I found Pelican. And here we are.

Customizing Pelican

Pelican feels extremely Pythonic to me, because it made it extremely easy to build a replacement blog from scratch using sane defaults and working with my existing Markdown. The default behavior is easy and intuitive and makes something that just works. Now, beyond using the defaults, customizing Pelican is done with Python modules, and there is an active community building reusable plugins and themes for Pelican.

To add plugins and choose a theme, I only needed to point at the pelican-plugins and pelican-themes repositories and add a little configuration. Most of the themes are listed with screenshots and live demos at I was very happy to find the pelican-bootstrap3 theme because it is clean and minimal, it's based on the widely used Bootstrap framework, and it is extremely customizable. With just a little exploration of the config space, I now have something that behaves and looks more like I wanted the original blog to and less like the Pelican default theme.

I had to stop myself from customizing too much, write this post, and update this blog after a long hiatus. It's functional and pretty much does what I want for now, but there is still more customization to be done. I would like to install a favicon. I want to get the built-in site search feature working. I would like to look at the Bootswatch themes for Bootstrap that can easily be dropped into this theme. I might add a tag cloud once I have some more content.

For now the most important thing is that the blog is up and running again. I hope the switch to Pelican will lower the barrier to updating a bit and help keep things moving. Sometimes the right tools make the difference. I'm also curious what else Pelican can be used for. Are there any downsides that I haven't run into yet? Other Pythonic tools for compiling content into web sites, documents, presentations?

by Mike Miller at May 31, 2017 10:40 PM

Joel Dahne

The first day and my repository

Yesterday I handed in the last exam for the semester(!) and today I started to work on the project full time. I begun by setting up a repository for my code. At the moment I have only cloned the interval package, I'm still to make my own changes to it.

I started to read the constructor for intervals to see how it works. While doing that I realized that there is much less to do than I thought. The constructor can actually handle creating intervals from numerical arrays of any dimension. The function displaying the interval can only handle up to two dimension but if you look at the internal state it actually has more dimensions. It is not perfect though, it does not work for decorated intervals and it cannot handle strings and cells as input. But still, it makes it much easier for me.

I will continue to work and hopefully I have actually committed something by the end of the week. Next week I will be away but after that I will code for the rest of the summer.

by Joel Dahne ( at May 31, 2017 05:38 PM

May 30, 2017

Michele Ginesi




During this first period I searched for the bugs related to special functions in the bug tracker. I found these four bugs that need to be fixed:
  1. gammainc: the upper version rounds down to zero if the output is small. For this bug there is a lot of work already done by Marco and Nir. I will just fix the last few things (like the input validation).
  2. gammaincinv: after finishing gammainc, I should start working on the inverse (which is missing in Octave).
  3. betainc: it is necessary to fix the input validation.
  4. besselj: it gives NaN if the input is too large. Actually, also other versions of the Bessel function have the same problem.
My idea of the timetable is the following
  • First part (May 30 - June 26: 4 weeks):
    • 1st week: finish to fix gammainc (input validation, documentation, add tests)
    • 2nd and 3rd weeks: write gammaincinv (via Newton's method)
    • 4th week: fix the input validation for betainc (and, if needed, add tests and fix the documentation)
  • Second part (July 1 - July 24: 3 weeks): fix the Bessel functions. Some ideas comprehend to test other libraries in addition to the amos (which is the library currently in use) and try to implement new algorithms. Add tests and revise the documentation.
  • Final part (July 29 - August 21: 3 weeks): add tests to the other special functions to make sure they work properly, trying to fix the problems that, eventually, will be found and revise the documentation if needed.

by Michele Ginesi ( at May 30, 2017 03:54 AM

May 28, 2017

Joel Dahne


This post will be about specifying a timeline for my work this summer. As I mentioned in the introductory post the work will be about implementing support for higher dimensional arrays in the interval package. To begin with I have divided the work into 5 different parts:

1. Construction and Printing

This part will be about implementing functions for creating and printing arrays. It will mainly consist of modifying the standard constructor and all the different functions used for printing intervals so they can handle higher dimensional arrays.

2. Vectorized Functions

Here I will work on generalizing all functions supporting vectorization to also support arrays of higher dimensions.

3. Folding Functions

Here I will work on generalizing all functions implementing some sort of folding to support higher dimensions. By folding I mean taking a multidimensional array as input and returning an array of lower dimension. Examples of these functions are $sum$, $prod$ and $dot$.

4. Plotting

I'm not sure what support Octave has for plotting higher dimensional arrays, but if there are some functions which could also be useful for intervals I will try to implement them here.

5. Documentation

I'll write the documentation alongside the rest of the work. In the end I will try to add some usage examples and integrate it better with the standard documentation.

So these are the parts in which I have divided my work in and the timeline will be

  • Phase 1 (30/5 - 30/6)
    • Week 1: Setting up
    • Week 2-3: Construction and Printing
    • Week 4: Vectorized Functions
  • Phase 2 (3/7 - 28/7)
    • Week 5: Continue on Vectorized Functions 
    • Week 6-7: Folding Functions
    • Week 8: Plotting
  • Phase 3 (31/7 - 25/8)
    • Week 9: Continue on Plotting
    • Week 10-11: Documentation
    • Week 12: Finishing up!

My first week will be rather short since I have an exam due in the middle of the week. After that I will also be away for a week not working on the project, I have no counted that week in the timeline above.

by Joel Dahne ( at May 28, 2017 08:30 PM

May 21, 2017

Joel Dahne


Hi! I'm Joel Dahne and I will be using this blog to share the progress of my work for Google Summer of Code 2017 under GNU Octave where I will be working on improving the interval package.

About me

Currently I'm a master student in mathematics at Uppsala University, Sweden, I'm just about to complete my first year and have one more to go. I first encountered interval numerics with my bachelor thesis, during which I tried to generalize the method described in this paper to higher dimensions. The code from the project is available on GitHub, Working on that project ignited my interest for interval numerics and I want to work to increase its availability to the common user.

The project

In short my project is about adding support for N-dimensional arrays in the interval package. At the moment the package only has support for up to 2-dimensional arrays. Most of the time this is enough but sometimes it is limiting. The goal is to have the syntax identical to the one for normal arrays. Switching from floats to intervals should be as easy as adding $infsup$ in front of the variable.

During my time preparing I have also noticed that some of the folding functions (e.g. $sum$, $prod$ and $dot$) handle some edge cases differently than the standard implementation. Hopefully this can be resolved at the same time.

This blog

During my project I will try to update this blog with my current progress. It's my first time writing a blog so we will see how it goes. The next step for me is to create more detailed plan of how I should structure my work during the summer, this will most likely be the subject of my next post.

by Joel Dahne ( at May 21, 2017 10:36 PM

May 16, 2017

Enrico Bertino

Introductory Post

Hello, I'm Enrico Bertino, and this is the blog I'll be using to track my work for Google Summer of Code 2017. The project will consist in rifare the neural network package with a focus on convolutional neural networks. Here my proposal.

Something about myself: I’m an Mathematical Engineering student and I live in Milan. After a BSc in Politecnico di Milano I did a master double degree both in Ecole Centrale de Nantes, France, and Politecnico di Milano. In the two years in France I attended engineering classes with spec in virtual reality and in Italy I am enrolled in the last year of MSc major Statistics and Big Data.

During my project, I will leverage on the Pytave project in order to exploit Python code for the implementation of the package. This is almost necessary because nowadays there is a strong tendency to open-researching in the deep learning field, where even the biggest companies release their frameworks and open their code. The level of this resources is very high, mostly for frameworks and APIs on Python, such as Tensorflow. This integration will allow us to guarantee a continuous updating on neural networks improvement and leading to a greater focus on Octave interface. In this sense, letting Octave coexist with Tensorflow, we will exploit their synergy and reach fast consistent goals.

As of now, I tested the integration of Tensorflow via Pytave. I have come across some problems of variable initialization in Ubuntu 16.04, whereas everything works fine on macOS Sierra. I will perform different tests in order to be be sure it works before starting the implementation of the package. Here my repository.

by Enrico Bertino ( at May 16, 2017 01:53 PM

May 13, 2017

Michele Ginesi


The expint function

The exponential integral

This is the first function I worked on in Octave: you can find the discussion here. Even if it will not be part of my GSoC project, I think it would be interesting to show how the function has been rewritten.

Definition and main properties

The canonical definition for the exponential integral is $$E_i (z) = -\int_{-z}^{+\infty} \frac{e^{-t}}{t}dt$$ but, for Matlab compatibility, expint compute the function $$E_1(z) =\int_{z}^{+\infty}\frac{e^{-t}}{t}dt.$$ The two definitions are related, for positive real values $z$, by $$E_1(-z) = -E_i(z) - \mathbb{i}\pi.$$ More in general, the function $E_n$ is defined as $$E_n(z) = \int_1^{+\infty}\frac{e^{-zt}}{t^n}$$ and the following recurrence relation holds true: $$E_{n+1}(z) = \frac{1}{n}\left(e^{-z}-zE_n(z)\right).$$ Moreover $$E_n(\bar{z}) = \overline{E_n(z)}.$$


To compute $E_1$, I implemented three different strategies:
  • Taylor series (Abramowitz, Stegun, "Handbook of Mathematical Functions", formula 5.1.11, p 229): $$E_1(z) = -\gamma -\log z -\sum_{n=1}^{\infty}\frac{(-1)^nz^n}{nn!}$$ where $\gamma\approx 0.57721$ is the Euler's constant.
  • Continued fractions (Abramowitz, Stegun, "Handbook of Mathematical Functions", formula 5.1.22, p 229): $$E_1(z) = e^{-z}\left(\frac{1}{z+}\frac{1}{1+}\frac{1}{z+}\frac{2}{1+}\frac{2}{z+}\cdots\right)$$ or, in a more explicit notation $$E_1(z)=e^{-z}\frac{1}{z+\frac{1}{1+\frac{1}{z+\frac{2}{1+\frac{2}{z+\ddots}}}}}$$ This formulation has been implemented using the modified Lentz algorithm ("Numerical recipes in Fortran 77" p.165).
  • Asymptotic series ("Selected Asymptotic Methods with Application to Magnetics and Antennas" formula A.10 p 161): $$E_1(z)\approx \frac{e^{-z}}{z}\sum_{n=0}^{N}\frac{(-1)^n n!}{z^n}.$$ A difficulty about this approximation is that the series is divergent, and that is the reason why the sum goes up to a certain $N$ "big" and not up to infinity.
Then I tested them in the square $[-100,100]\times[-100,100]$ of the complex plane, comparing the result with the ones given by the Octave symbolic package. With these informations, I identified three zones of the complex plane: in each zone one strategy is better than the others.
Then the final function divides the input into three parts, accordingly to the position of the same in the complex plane, and compute the function using the opportune strategy.

by Michele Ginesi ( at May 13, 2017 07:16 AM

May 10, 2017

Michele Ginesi


Introduction to the project


About me

My name is Michele Ginesi, I obtained a Bachelor's degree in Applied Mathematics in Verona, Italy. Now I am getting a Master degree in Mathematics in the same university.
I was selected to partecipate to the Google Summer of Code under GNU Octave for the project Make specfuns special again.

About the project

Special functions are an interesting and important topic in Mathematics, and so it is fundamental to have a way to compute them in an accurate way.
Some examples of special functions (whith some important application of the same) are:
  • Gamma function $\Gamma$. This is one of the most important: it can be viewed as the extension of the factorial function ($\Gamma(n)=(n-1)!$, for $n\in\mathbb{N}$) and it is a component in various distribution functions in probability theory.
  • Beta function $B$. This was the first known Scattering amplitude in String theory.
  • Bessel functions. These are the canonical solutions $u(r)$ of the Bessel's differential equation $$r^2\frac{d^2u}{dr^2}+r\frac{du}{dr}+(r^2-\alpha^2)u=0 $$ for any complex number $\alpha$. At a first approach may seems that they are an end in themselves, but actually the Bessel equation describes the radial component of the two dimensional wave equation ($u_{tt}-c^2\Delta u = 0$) applied on a disc.
The most common strategies used to approximate special functions are Taylor series and continued fractions, sometimes asymptotic series. Some particular functions can be implemented via recurrence formula or other type relations with other functions (for example, $B(z,w) = \Gamma(z)\Gamma(w)/\Gamma(z+w)$).
This project will be divided into three main parts:
  1. Fix the already known bugs related to special functions (e.g. #48316, #47800, #48036).
  2. When the known bugs will be fixed, I will proceed to add new tests to make sure that all the functions are accurate.
  3. Fix the new problems/bugs that, eventually, will be found during the second phase.
The main reference for this work will be Handbook of Mathematical Functions by Irene Stegun and Milton Abramowitz which contains all the functions I will work on completed with (almost) all the expansions that are needed to implement them; and Handbook of Continued Fractions for Special Functions by Annie Cuyt, Vigdis Brevik Petersen, Brigitte Verdonk, Haakon Waadeland and William B. Jones.
To test the functions I will use, in addition to the Handbook of Mathematical Functions, Sage and the Octave symbolic package to get a reference value.
Here you can find my repository on which I will work.

by Michele Ginesi ( at May 10, 2017 05:24 AM

March 14, 2017

John W. Eaton

Timeline of Octave Releases

Someone asked me in private email why randi isn’t available in Octave. They were using 3.2.4 on Ubuntu 12.something and quoted a 6-year old email from the mailing list that said it was added in Octave 3.4. Maybe they are a time traveler from the past? Anyway, after fumbling around for the n-th time trying to find out when 3.2.4 was released (the answer is, a long time ago) I created this page to make those kinds of searches easier.

by jwe at March 14, 2017 07:25 PM

Silent Running, or, Back to the X-Term of the Future

I bought one of these:

put 8GB of memory and a 120GB SSD in it and mounted it on the wall.  So now I have (more or less) a modern-day quad-core X-Terminal with a gigabit network connection running Debian and Gnome.

Next, I moved my former noisy, heat-generating deskside system to the basement and, voila, sweet silence!  All I hear now is the residual ringing in my ears that is probably the result of attending too many loud rock-n-roll shows in my youth.  But that’s a story for a different time…

My home directory remains on the system that used to sit beside the desk.  It has a 6-core AMD CPU and a pair of disks in a RAID-1 configuration.  Accessing the files over the gigabit network seems fine, at least for normal usage.  And even though the new ZBOX has pretty modest specs, it is more than sufficient for most of the things I do: email, web browsing, and running Octave for testing.  I don’t use it  for building Octave.  For that I have a couple of other systems also located down in the basement where I can’t hear them.

Total cost was about $250.  I figured that wasn’t too much to pay for silence.  And if it turned out that the performance wasn’t great, I could always save up a few more pennies to buy a different model that has a faster processor and/or more memory and then reuse this box as a firewall – it does have 2 separate ethernet ports, as if it were intended to be used for that purpose.  Or maybe as an HTPC since it also has an HDMI output and Intel graphics that are supposed to support up to 4k graphics (I don’t know, my monitor works great with this system but it is just 1980 x 1200).

Setup was pretty easy and Debian installed without any problems from a USB stick.

It takes about 35 seconds to boot up to a gdm3 login screen.  That includes a few seconds each for POST and grub delays.  I could eliminate those and save a few seconds but then it would be harder to enter the BIOS or change boot parameters.  And how often do I reboot anyway?

From the login screen to a terminal window is another 15 seconds or so.  And that’s all with a 1.6GHz Celeron processor.  If I’d splurged and spent a few hundred more for the Core i5 version, I guess I’d be up and running in about 30 seconds?

by jwe at March 14, 2017 03:16 PM

August 30, 2016

Cristiano Dorigo

CGS, TFQMR and final revision of all the codes

Dear all,

During these weeks I worked on improving CGS, which I rewrite from scratch. Then I implemented TFQMR, which is not available in Octave. Finally I made a revision on the codes that I improved the past weeks.
All the file revisioned (gmres, pcg, bicg, bicgstab) and the new ones (cgs and tfqmr) are available here:

in the folder definitive_codes


As written in my second post, the actual Octave's CGS implementation is not very clear. Indeed I didn't understand most of the code (I was not able to recognize how the is was adapted). Moreover I tested it and the results are not correct (in my second post there is an example of this).
Then I decided to rewrite if from scratch following closely the algorithm given by Saad in the reference book suggested.
I wrote it following the pattern that I used in the other codes (i.e. to wrote three subscripts, each for the different cases: all matrices, all functions handle or mixed). Iit has also the same characteristics of the other methods improved, for example: the approximation returned is the one with the minimum residual, there are five types of flag (0 convergence, 1 maximum number of iterations reached, 2 preconditioner singular, 3 stagnation, 4 division by zero during the process)... I don't go in deep about this details to make less boring this post, since they are already explained in my previous posts. 


In my last post I wrote that the mentors told me to implement a new algorithm not available in Octave. Since I had no suggestion about which method write, and the community had no particular needs, I chose the transpose-free qmr algorithm because it is the only methods treated in the Saad book.
Also with this algorithm, I use the same pattern of the others (i.e. the three subscripts) and synchronized in such a way it has the same type of input and output as the others.

One mention to the number of iterations effectively performed. In the Saad book, the method has different behaviours if the iteration number is even or if it is odd, but at every iteration the approximation is computed. Making some test in Matlab with its tfmqr, I noticed that if the method converges, the vector of the residuals in long two times the number of iteration that Matlab tells that are necessary to convergence. Then I think that Matlab counts as one iteration the block (odd iteration + even iteration), i.e. I think that Matlab_it = odd_Saad_it + even_Saad_it. In this way the iterations performed by Matlab is two times the effective iterations performed. I decided to count the iterations in this way to be compatible.


Finally I checked (one more time) the codes that I improved in these weeks (i.e. pcg, gmres, bicg, bicgstab).
Indeed I updated the documentation (for example: making all similar, giving more details about the preconditioning, adding the reference book,...).
Then I checked the codes and I tried to make it more readable and to adapt it at the Octave conventions.
Then I added some tests in all the codes, mostly to check if all the subscripts works (i.e. checked every combination of A matrix/function, M1 matrix/function, M2 matrix/function), but also to check that the flags works well (I checked only the flag 0, 1 and 2, since the others are difficult to reproduce) and also to check that these methods solves also complex linear systems.

The mentors suggest me to make a patch with the methods improved (pcg, gmres, bicg, bicgstab and cgs), in such a way it can be easily the review of them and if it is in time (maybe) they can be included in Octave 4.2.
This patch is available here

and it is the file SoCiS16_Improve_iterative_methods.diff

I made this patch from the Octave changeset 22406:dc4e2203cd15 of the Octave 4.1.0+ version.


Since the official deadline for the SoCiS  is tomorrow (the 31-th of August), this is my last post in this blog.

To make a small summary, during this summer of codes I make an improving on
most of the iterative methods to solve linear systems. In particular I improved pcg, gmres, bicg, bicgstab and cgs.
I wrote for all of these methods three subscripts to make them more efficient, dividing the particular cases "matrix", "function handle" and "mixed", in such a way that, according to the input data, the method uses the appropriate subscripts to make computations.
I tried to make these methods more compatible with Matlab (for example noticing that Matlab gives as output the approximation with the minimum residual and not the last performed, as used in Octave, or adding the possibility of "flag =2", i.e. singular preconditioner).
I tried to fix some "non-clear" situation that I noticed, for example pcg does not noticed if the matrix was not Hermitian positive definite if it was a complex matrix, or to face the problem of a division by zero in bicg, bicgstab and cgs.
Then I implemented tfqmr which is available in Matlab, but not in Octave.

Unfortunately, there are other two methods that needs an improve (pcr and qmr), but due to unexpected problems with the above methods I had not enough time to improve also these. When my scholastic commitments are not so big, I want to try to improve also these two methods.

I want to thank all the Octave community that gives me this fantastic opportunity to face a big project as this and then to learn  much theoretically and practical things about these iterative solvers for linear systems.
I want to thank specially the mentors, who were always available to give suggestions, advices and ideas of how to proceed in the work, mostly when came some unexpected problems.
I hope that this improving will be useful to Octave, and to contact me if there is something not clear or wrong.
I also hope to continue to help the Octave community in the following times.

One more thanks to all of you.

by Cristiano Dorigo ( at August 30, 2016 10:34 AM

August 23, 2016

Francesco Faccio

Summary of the work - GSoC

MathJax TeX Test PageHello!

In this final blog post I want to summarize the work done so far for the GSoC. First of all I want to thank Octave community for their support:
when I had some problems I have always found someone in the mailing list or in the IRC channel willing to help me.

In my public repository I pushed three final commits:

The first commit I pushed regards the work we did in the first part of the GSoC and it is a patch for Octave which allows to include IDA module of Sundials library as a dependence of Octave.

I wrote some configuration tests, checking for sundials_ida and sundials_nvecserial headers and functions and a macro to test whether sundials_ida has been configured with double precision realtype. I included a test for klu library in suitesparse, because it is required for sundials' sparse solver.
If all these conditions are satisfied, you can use ode15i and ode15s in Octave (you can see in my previous post how to configure sundials and octave correctly and how to run benchmark problems).

The second commit regards my work on solvers ode23 and ode45 which were already in Octave.
Since there were a lot of duplicate code, we decided to refactor the solvers, exporting the check of the options into two private functions:

$[defaults, classes, attributes] = odedefaults (n, t0, tf)$

is a function which returns a struct defaults, containing the default values of all the solvers, and structs classes and attributes, containing the arguments for validateattributes.

$options = odemergeopts (useroptions, options, classes, attributes, fun_name);$

is a function which loops over default fields and, if a user provides a new value, uses validateattributer or validatestring to test the value and overwrites the default one.

odeset.m and odeget.m have been modified in such a way that the user can add new special fields to an ODE option structure. They don't depend on private functions ode_struct_value_check and known_option_names, but use inputParser to perform input checking. This version of odeset doesn't include the option to catch a misspelled option name yet, so that 'Rel' doesn't become automatically 'RelTol'.

We wrote function decic, which can compute consistent initial conditions for ode15i, close to initial guesses y0 and yp0.

The third commit contains files ode15i.m, ode15s.m and a common interface


The main functionalities of ode15i have been achieved. Here's a list of the options accepted by the solver:

RelTol: You can specify a relative error tolerance as a positive scalar. The default value is 1e-3.

AbsTol: You can specify an absolute error tolerance as a positive scalar or vector. The default value is 1e-6.

NormControl: Sundials controls the relative error using the norm of the solution at each step, so option NormControl is always "on".

OutputFcn, OutputSel: You can specify them to plot the components of the solution you want after each successful integration step. The default output function is $@odeplot$, but you can provide a custom output function.

Refine: You can increase the number of output points using the refine factor. uses the IDAGetDky function to reach the highest possible order of interpolation.

Stats: Only standard statistics are printed (number of successful steps, number of failed attempts and the number of function evaluations).

Events: If you provide an event function, it must be of the form $ [value,isterminal,direction] = EventFcn(t,y,yp)$. The solver localizes a change of sign on your event function and uses linear interpolation to extimate the time when the event occurs. A way to improve this functionality will be to use the rootfinding function of Sundials in order to have an higher order of interpolation.

MaxStep: Use this option to specify an upper bound for the integration step. The default value is $0.1*abs(tf - t0)$.

InitialStep: Use this option to specify the first integration step of the solver. Its default value depends on the tspan and the initial slope.

Jacobian: You can specify the Jacobian as a function which returns the modified Jacobian $[\frac{\partial F}{\partial y}, \frac{\partial F}{\partial y'}] = jac (t, y, y')$ or as a cell array containing the two constant Jacobian matrices $ \{[\frac{\partial F}{\partial y}], [\frac{\partial F}{\partial y'}] \}$.

In both cases, if the two parts of the Jacobian are sparse, the Jacobian will be stored in CCS format.
If you don't provide the Jacobian, it will be extimated internally by Sundials, using difference-quotient approximation.
In this implementation you can't provide only one part of the Jacobian yet, but you have to specify both.

JPattern, Vectorized: These options haven't been developed yet.

MaxOrder: You can specify an upper bound for the order of the solver. The default value is 5.

This solver can be extended in order to increase the compatibility with Matlab:
- completing the missing options
- allowing the user to solve equations with complex number
- allowing the user to use single precision variables.


This solver uses the same dld-function of ode15i. Inside the m-file, the semi-explicit DAE is written as an implicit one and solve as ode15i.
The solver uses the same options of ode15i, but in addition, it has:

NonNegative: this option has not been implemented yet. It can be implemented using IDASetConstraint function.

Events: Same as ode15i, but it must have the form $[value,isterminal,direction] = EventsFcn(t,y)$.

Jacobian: You can provide a function which returns the Jacobian of the form $[\frac{\partial F}{\partial y}] = jac(t,y)$ or a constant Jacobian matrix.

Mass: You can provide a constant Mass matrix or a function handle of the form $M(t,y)$ or $M(t)$. If the Mass matrix depends on the state, the Jacobian Matrix will be computed internally from Sundials. If you provide a sparse Jacobian and a sparse Mass matrix, they will be stored in CCS format.

MStateDependence: Specify MStateDependence to "strong" or "weak" if the Mass matrix has the form $M(t,y)$. Use the value "none" if the Mass matrix has the form $M(t)$ or it is constant.

MvPattern: This option has not been developed yet.

MassSingular: This option has not been developed yet. The solver solves the problem both with and without a singular mass matrix in the same way.

InitialSlope: You can specify the initial slope of the problem. The default value is a vector of zeros.

BDF: this implementation of ode15s uses always BDF method.

This solver can be extended in order to increase the compatibility with Matlab:
- completing the missing options
- allowing the user to solve equations with complex number
- allowing the user to use single precision variables.

This code has not been merged yet, since Octave is concluding a release process.

Concluding, I really enjoied working on this project and there is still a lot of work to do. I will work in the next weeks to extend the functionalities of ode15i and ode15s. I plan to continue to contribute to Octave development, maybe both in ode and in some statistics project.

by Francesco Faccio ( at August 23, 2016 09:04 AM

Chiara Segala

Final Evaluation

Dear all,

GSoC16 is at the end, I want to thank you for giving me the opportunity to participate in this unique experience.

My goal was to add two solvers to the odepkg package, two exponential integrators, that are a class of numerical methods for the solution of partial and ordinary differential equations.
  • Here is a summary of the work done in this months:

For the midterm evaluation (repository of the midterm work): 

After a first period of study, I started with implementations.
I start with the implementation of the φ functions (phi1m.m, phi2m.m, phi3m.m, phi4m.m), necessary to calculate the matrix functions in the two schemes (exprk.m, exprb.m) for a general exponential Runge-Kutta and Rosenbrock integrator. These schemes have been very useful to verify the correctness of the official methods. I applied them to four different exponential methods (exprk2.m, exprk3.m, exprb3.m, exprb4.m).
Then I wrote a code that implements the augmented matrix  Ã (augmat.m) described in theorem 2.1 of [HAM 11]. This matrix will be given as input to expmv and will allow us to implement the exponential integrators by evaluating a single exponential of this augmented matrix avoiding the need to compute any φ functions.

For the final evaluation (REPOSITORY OF THE FINAL WORK):

After the midterm evaluation, I implemented the official exponential method that will be in odepkg, it is the advanced third order method with adaptive time stepping from Runge-Kutta family. Unlike the general schemes, for this code I used the function expmv so that I could evaluate a single exponential of an augmented matrix instead of computing any φ functions. I wrote a code for the stepper (exp_runge_kutta_23) and for the solver (odexprk23), according to odepkg standards.
Then I changed the file augmat.m to calculate a "double" augmented matrix, and in this way, I have been able to calculate, with a single call to expmv, both the solution to the next step, and the lower order solution for the estimation of the error.
In the same way I implemented the official exponential fourth order method with adaptive time stepping from Rosenbrock family (exp_rosenbrock_34, odexprb34)
I finally tested the order of my methods and I translate in C++ a for loop ( of the Nick Higham function expmv to improve the speed of an algorithm.

  • In the repository of the final work you will find a patch for odepkg.

I added to the package the following files

/inst/odexprk23.m : Solve a set of stiff Ordinary Differential Equations (stiff ODEs) with a third-order exponential Runge-Kutta method.
/inst/odexprb34.m : Solve a set of stiff Ordinary Differential Equations (stiff ODEs) with a fourth-order exponential Rosenbrock method.
/inst/steppers/exp_runge_kutta_23.m : Stepper for the odexprk23 solver.
/inst/steppers/exp_rosenbrock_34.m : Stepper for the odexprb34 solver.
/inst/utilities/expmv.m : Compute the matrix exponential times vector or matrix (2010 Awad H. Al-Mohy and Nicholas J. Higham).
/inst/utilities/augmat.m, estimates.m, expmv_tspan.m, select_taylor_degree_aug.m, NormAm.m, select_taylor_degree.m, theta_taylor, theta_taylor_half, theta_taylor_single : Functions related to expmv.
/src/ : for loop of expmv. 

  • I export the patch from
odepkg changeset: 542:3f74a6e911e1,
octave version: 4.1.0+, changeset: 22348:9deb86bb5632.
My patch exponential_integrators.diff requires a changeset >= 22345.
To apply the patch follow the commands

hg import exponential_integrators.diff
hg archive ../odepkg-0.9.1.tar.gz

run Octave 4.1.0+

pkg install odepkg-0.9.1.tar.gz
pkg load odepkg 

to test the patch:

demo odexprk23
demo odexprb34 (slightly slower)

if everything goes well, the demo gives a plot on convergence of methods.

  •  Demos give the following warnings
warning: odexprk23: unknown field 'LinOp' in ODE options
warning: odexprk23: unknown field 'GFun' in ODE options

warning: odexprb34: unknown field 'DFdt' in ODE options  
I haven't fixed the warnings because another GSoC16 student is working on the package and because outside the scope of my project, but the problem can quickly be solved collaborating at the official end of this experience.

by Chiara Segala ( at August 23, 2016 03:09 AM

August 22, 2016

Amr Mohamed


GNU Octave – Geometry Package – GSoC 2016

This Repository ( ) is a clone(fork) of the Geometry Package which is a part of the Octave Forge Project (

This fork adds new functions to the official Geometry Package as part of GSoC (Google Summer of Code) 2016.

The official Geometry Package can be found here:

Added Functions

  • polyjoin: Convert cell arrays of polygons into NaN delimited coulmn vectors.
  • polysplit: Convert NaN separated polygon vectors to cell arrays of polygons.
  • ispolycw: Check whether the polygon(s) has(have) clockwise orientation or counterclockwise orientation.
  • poly2cw: Convert Polygons to clockwise contours(polygons).
  • poly2ccw: Convert Polygons to counterclockwise contours(polygons).
  • polybool: Perform Boolean Operations on polygons.
    Valid Operations are: (intersection,and,&), (union,or,|,+,plus), (subtraction,minus,-), (exclusive-or,xor).
  • poly2fv: Convert Polygons into Faces and Vertices Representation.


Community Bonding:

As mentioned in my proposal , I started working on the project 3 weeks before the official start of the GSoC program to compensate the three weeks that overlap with my final examinations.
I started my implementing the polysplit & polyjoin functions as they don’t depend on other functions.
These functions change the representation of the polygons from one form to another.
Since the polysplit & polyjoin functions don’t depend on Boost Geometry, they are implemented as .m files and placed in the inst/polygons2d/ directory.

First Coding Period:

During the first two weeks of this coding period, I continued the implementation of the polysplit/polyjoin functions adding help messages and regression tests.
Then, i worked on the first function that is implemented as a cc file which is the ispolycw function.

As all the package functions were implemented as .m files, i created a new directory src/ so that i can add .cc functions to the package.

The ispolycw function depends on calcualting the signed area of the polygon using the Boost Geometry function area(polygon).
The orientation of Self-intersecting polygons is defined as the respective orientation of the leftmost top point.If these points are co-linear ,  the signed area of the polygon is used to determine its orientation.

To compile the ispolycw function and generate .oct file that can be run by Octave, A new Makefile was added to the src/ directory.
The Makefile initially used the socket’s package implementation of the Makefile.

Two weeks after the start of the GSoC and we were already ahead of schedule. So we decided to start working on the polybool function as it is the main function of the whole project.
We decided to start by implementing the polybool function to support doing  boolean operation between only two polygons.

Second Coding Period:

By the start of the second coding period, we changed the implementation of the polybool  function to use the Boost Geometry multipolygon concept so that it can support applying boolean operations on multiple polygons with inner holes.

During the first three weeks the  second coding period, simple benchmarks were used to test the polybool  function.
And consequently multiple implementation bugs were fixed. One of these bugs was that a single inner ring was initially matched to multiple outer rings that surrounds this inner ring.
So a boolean array is currently used to ensure that each inner ring is matched to one and only one outer ring.

Then, i wrote two simple functions poly2cw poly2ccw which are used to toggle the orientation of polygons and thus convert inner rings (holes) into outer ones and vice versa.

Then in the fourth week of the second coding period, i decided to make better understanding of Makefiles , how to use them to build only changed files, how to organize dependencies between targets and how to use variables.
Also, i tried to use the dissolve extension of the Boost Geometry to handle self intersecting polygons.
It took me a few days to manage to use it as the Boost Geometry version on my system (Ubuntu 14.04) was 1.54 while the dissolve extension is developed based on Boost Geometry libray of version 1.60
After that, i adapted the ClipperLib benchmarks to Octave so that we can compare my implementation of these functions to the ClipperLib / other programs performance.

By the end of July, i  started working on the poly2fv function which is used to triangulate polygons so that they can be better visualized using the patch command.
This function returns a 2D matrix consisting of the triangles’ vertices and another 2D matrix representing the indexes of each triangle respectively.

For the first two weeks of August, i reviewed all the newly implemented package functions, added missing regression test , improved the help messages , added more input format checks , added warning messages when some boolean operations fail to be done mainly because of self-intersections that can’t be resolved.

Then during the last two weeks of the program, i had to learn how to use autotools so that i can use a configure script to test whether the Boost library is installed on the system, and check whether the dissolve function can be used to resolve self-intersecting polygons or not.
It was a great experience to manage to use the autotools to manage the generation of the final Makefile that is required to generate .oct files from the implemented .cc files.

Week Task Working hours / week
22-4 / 29-4 Understanding the organisation conventions – Prepare a complete plan for each function (Input/Output requirements – Optimal implementation based on memory/time metrics) 20 hour/week
30-4 / 6-5 Implementing the polysplit function and testing it. 20 hour/week
7-5 / 13-5 Implementing the polyjoin function and testing it. 14 hour/week
14-5 / 10-6 * Finalize the implementation of the conversion functions (polysplit-polyjoin) 7 hours /week
11-6 / 19-6 Implementing and testing the ispolycw function. * 40 hours/week
20-6 / 27-6 Midterm Evaluation week 40 hours/week
25-6 / 1-7 Working on the poly2cw and ispolycw scripts 40 hours/week
2-7 / 8-7 Performing severe testing to all the implemented functions 30 hours/week
9-7 / 15-7 Working on the polybool script. 30 hours/week
16-7 / 22-7 40 hours/week
23-7 / 29-7 Implementing and testing the poly2fv function. 40 hours/week
30-7 / 5-8 A 2 week buffer for finalizing the scripts , debugging them and writing the documentation. 40 hours/week
6-8 / 12-8 40 hours/week
13-8 / 23-8 Tidying the code, writing tests, improving the  documentation and submitting the code sample. 40 hours/week

* Final Examinations – Schedule isn’t announced yet – will inform the mentor whenever it is announced.


All the newly added functions has the following set of conventions:

  1. Outer Rings of the polygon has a clockwise orientation, while Inner Rings has a counterclockwise orientation.
  2. Description of the polygons is represented by two arguments X and Y. Each argument can be represented using Cell Arrays or nan Delimited Vectors.
  3. The functions don’t require that the last point is the same as the first one.

Example: To represent a square as shown in figure:

— NaN Delimited Vectors:

x=[0 0 3 3 NaN 1 2 2 1];

y=[0 3 3 0 NaN 1 1 2 2];

— Cell Arrays:

x={[0 0 3 3],[1 2 2 1]};

y={[0 3 3 0],[1 1 2 2]};


polybool & poly2fv functions use the Boost Libraries.

The recommended version is 1.60

Getting Started

To test the updates to the Geometry package:

  • Make sure you have Boost Libraries installed.
  • Clone the package
hg clone
  • Run ./bootstrap in the root directory of the package.
  • Create a tar file
tar -zcf geometry.tar.gz octave-geometry/
  • Run Octave.
  • In Octave Command Prompt, write
pkg install geometry.tar.gz
  • Load the package
pkg load geometry


All the benchmarks where run on a Ubuntu 16.04 computer with Boost Library 1.60.

All the scripts can be found here:

As my blog doesn’t support tables, Kindly check the Tests/Notes which can be found as a single pdf file here:


-poly2fv depends on polytri library which isn’t developed to throw exceptions when it can’t triangulate polygons so sometimes poly2fv calls may cause Octave to crash.

by amrkeleg at August 22, 2016 04:07 PM

Abhinav Tripathi


During Google Summer of Code 2016, I worked on the Symbolic Project under the organisation Octave. The project was directed towards adding a new IPC mechanism into symbolic – Pytave. During the course of GSoC work was done in 2 repositories (octsympy on github & pytave on bitbucket).

Code Changes:

All of the changes done by me can be viewed at the following links:

Summary of changes:

As we were to use pytave in Symbolic we worked on it to make it more mature. Add more functionalities to provide a better interface for Symbolic to deal with. Also many new features for conversion to and from python objects were introduced. The most important addition would be @pyobject and its functionalities.
On symbolic side, we used the features of pytave as much to our benefits as possible. We currently only have extra conversion for @sym objects which are specific to Symbolic. For rest of the octave types we use pytave‘s internal conversion.

Detailed PR wise changes:

To know about changes done before August, please read my other posts on this blog in the sequence: post-1, post-2, post-3 and post-4. Here I will summarise PR wise work done in August only.

In pytave, the PRs #33, #31 and #32 were merged which were under review in july. The PR 33 did some cleanup and added tests. PR 31 adapted existing tests since after inclusion of @pyobject we were now converting  python – lists, dicts and tuples; to @pyobject in octave. PR 32 fixed some doctests which were either not working in Python 3 or gave errors due to the latest @pyobject changes.
PR #24 added a new functionality to specify a namespace to be passed in py* functions to ensure that any changes in symbols/objects happens only in that namespace. This required some refactoring which was done in the same PR in a separate commit.
Then PR #40 added a method to allow comparison of 2 pyobjects. This would be useful in writing tests that need a desired output.

On the Symbolic side, PR #512 is under review which takes care of updating out existing methods to use the new pytave features. We have had a good amount of changes in pytave which now allow us to do most of the things without the need to execute and python code on our own (using pyexec). We now do it using py.* calls.

Goals not met:

Almost all the goals were met during the project. Only 1 major goal was not met which was building pytave on windows. Efforts were made towards that but later on it was more fruitfull to add new functionalities that would make the projects more mature. I would definetly try to work on this afterwards and hopefully in near future we will have fully working pytave on windows too.


by genuinelucifer at August 22, 2016 03:50 PM

Barbara Lócsi


GSoC 2016 is about to end so in this post I will sum up the work we had done during this period. I will write about 3 major topics: eig, gsvd and chol2inv. The original project was about eig, and when it was done we moved on to the other ones.


The goals of the project were:
Certain calling forms of the eig function are currently missing, including:
  •  preliminary balancing
  • computing left eigenvectors as a third output
  •  choosing among generalized eigenvalue algorithms
  •  choosing among return value formats of the eigenvalues (vector or matrix)

Calling forms:

We finished these goals before the end of GSoC, and most of the work is already merged to the main octave repository:
In my repository there are some other to be reviewed commits that reduces code duplication with the use of templates.
You can find some information about the individual tasks on this blog in previous posts, and also here:


The function gsvd was imported from Octave-Forge linear-algebra package. However it uses deprecated LAPACK functions: *ggsvd. Moreover it is not compatible with Matlab’s gsvd. For more information check these discussions:
  1. Deprecated LAPACK routines (GSVD)
  2. new gsvd function incompatible w/Matlab
  3. bug #48807: new gsvd function is incompatible w/Matlab
Work already merged to the main Octave repository:
Other commits that are not yet merged:


chol2inv couldn’t handle sparse matrixes (bug # 36437), we submitted a fix for this bug and tests for chol2inv.

GSoC 2016 with GNU Octave

This was an exeptional opportunity to gain experience and get into open source development and learn about Octave. I am really happy for this project and the time we spent working on it. I would like to thank for my mentors Carnë Draug and Nir Krakauer, without their help and guidance non of this would have been possible.

by Barbara Lócsi ( at August 22, 2016 02:02 PM

August 16, 2016

Francesco Faccio

How to replicate test


Here's a little guide to replicate the test I made for ode15i


Download last version of Sundials (only IDA module is required (IDA version > 2.8.0), but SundialsTB, idas and kinsol will be used for some test). Follow the install guide of Sundials to build IDA module with KLU (you will have to specify the lib and include directories where you have previously installed KLU). Use SUNDIALS_PRECISION as double (default) in configuration.


Clone my octave repo you can find here:
and build Octave from source. KLU is already a Octave dependence, but you can specify other versions you have installed at configuration time.


Here you can find the scripts and C functions I used to test ode15i.
You can test my oct-file running the scripts in folder Octave and the same scripts in Matlab folder for testing Matlab's ode15i.

To test the C examples of Sundials compile the C file you can find in folder C, linking against libklu, libsundials_ida and libsundials_nvecserial. If you installed these libraries in default path, type:

$ gcc -c idaHeat2D_klu.c idaRoberts_dns.c
$ gcc -o idaHeat2D_klu idaHeat2D_klu.o -lsundials_ida
-lsundials_nvecserial -lklu -lm
$ gcc -o idaRoberts_dns idaRoberts_dns.o -lsundials_ida
-lsundials_nvecserial -lklu -lm
$ ./idaRoberts_dns
$ ./idaHeat2D_klu

Sundials provides some mex interfaces for solving the benchmark problems. In order to test that, you have to edit file SundialsTB/install_STB.m changing mex compiler (near row 20):

mexcompiler = 'mkoctfile --mex';

Install Octave and run from Octave install_STB.m script. Choose to compile IDAS and kinsol interface and install the Toolbox. A startup file is generated where you chose to install SundialsTB: run it.

Now you could be able to run the script in MEX folder, which uses Sundials mex interface for solving the dense problem.

For problem ida_klu in Octave and Matlab folder, you can choose to test the solver with a sparse jacobian function, a dense jacobian function, a cell jacobian dense or a sparse jacobian dense.

Since I added the cell jacobian option in ode15i recently and I changed some parameters in the benchmark problems, I have some new results:

Problem Robertson eq (tspan from zero to 4e10, abstol = [1e-8, 1e-14,1e-6], reltol = 1e-4, dense methods):

Matlab: 0.41s
Octave: 0.06s
C impl: 0.002s
Mex impl: 0.12s

Without jacobian:

Matlab: 0.45s
Octave: 0.07s
Mex impl: 0.15s

Problem Heat2D (tspan from zero to 10.24, reltol = 1e-5, abstol = 1e-8, sparse methods):

With a sparse jacobian function:

Matlab: 0.34s
Octave: 0.13s
C impl: 0.005s

With a dense jacobian function:

Matlab: 0.44s
Octave: 0.15s

With a dense jacobian cell:

Matlab: 0.45s
Octave: 0.08s

With a sparse jacobian cell:

Matlab: 0.34s
Octave: 0.07s

Without a jacobian:

Matlab: 0.56s
Octave: 0.47s

It seems a quite goode result, considering that C and Mex implementations are specificly designed for solving these problems. Since I use a m-interface before calling my oct-file, there is a little overhead, but it will be negligible for bigger more complex problems.

For these test I configured Octave with -O2 flag and all the results refer to the code till this commit

by Francesco Faccio ( at August 16, 2016 01:43 PM

August 03, 2016

Amr Mohamed

Blog Post

Hi all,
I am so happy that our project is on the right track.
We have managed to implement the final function (poly2fv) on time.

This function perform Constrained Delaunay triangulation that i useful for plotting polygons with inner holes.
The function is implemented as a cc source file ,and it uses the polytri library which can be found here:

We have also adapted some of the ClipperLib Benchmarks to Octave and here is a screenshot for one of them:



by amrkeleg at August 03, 2016 04:16 PM

August 02, 2016

Cristiano Dorigo


Dear all,

In this last period I worked on the BICG and BICGSTAB algorithm.
You can find the codes in the folder bicg_bicgstab at

Since the main loop of both of them have more or less 30-40 lines of code, I wrote them from scratch, because in my opinion was easier instead to modify and to adapt the original code at the many situations that occurred.

The strategy that I used is the same of PCG and GMRES: for each algorithm I wrote three sub-scripts to subdivide the cases which: (a) the matrix and the preconditioners are all matrices, (b) all functions handle or (c) some are matrices and some functions handle.

As for PCG and GMRES, I tried to make these codes compatible with Matlab.

In my BICG and BICGSTAB codes there are not actually the documentation, but I'll add them as soon as possible

Main modifications on BICGSTAB and BICG

Many of the changes that I applied at the BICGSTAB and BICG were quite similar to those that I made in PCG and GMRES. I make a little list to summarize them: 

  • Now the approximation of the solution X returned is the one that has the minimum residual and not the last one computed by the algorithm. Also the output variable ITER contains the iteration in which X was computed, not the number of iterations performed.
  • In the original Octave BICGSTAB code the stagnation criterion is
    resvec(end) == resvec(end - 1).
    Instead in the BICG code it is: res0 <= res1 (where res1 is the actual residual and res0 the previous)
    I changed this criterion with:
    norm(x - x_pr) <= norm(x)*eps
    (where x_pr is the previous iteration approximation)
    This criterion is more efficient, since it checks the variation of the approximations and not the variation of the residual, and it seems more similar to the "stagnation" described in the Matlab help. Moreover in BICGSTAB (and in BICG) there are many fluctuations of the residuals (they are not decrescent in general), so it's better to check directly the approximations.
    Due to these fluctuactions, the BICG criterion is wrong, since in this manner the Octave BICG usually stops also if it is not necessary.
  • I added the "flag = 2" criterion, the one that Matlab describes as "preconditioner matrix ill-conditioned", that in Octave is skipped. The check that I implemented is for the singularity of the preconditioner matrix M, since seems that in Matlab this flag appears only in this situation.
Now I will focus on some situations regarding BICGSTAB:
  •  I added the "half iterations" as in Matlab: studying the Saad implementation, I noticed that the approximation X is made in two steps at every iteration, indeed it is defined as:
    x = x + alpha*p + omega*s
    Then to be compatible with Matlab (since p/alpha and s/omega are not related and computed in different situations), I set x = x + alpha*p (and its residual) as the half iteration, and x = x + alpha*p + omega*s as the complete iteration.
    In this way  the "real" total number of iteration performed (recalling that ITER contains the iteration which X was computed) are
    (length(RESVEC) - 1) / 2.   
  • There was a problem when the matrices A, M1, M2 and/or the right-hand-side b are complex, indeed many times the residual blow up and the method does not converge. In the Saad reference book (where it was taken the Octave BICGSTAB, I think), the method is formulated only for the real case. The only difference is that using this formulation, some inner product are inverted (e.g. it uses <v, w>, instead of <w, v>). In the real case there is no problem since the inner real product is symmetric, but in the complex not, indeed the result of <v, w> is the complex conjugate of <w, v>.
  • In the Matlab help I noticed that there is the possibility of "flag = 4". I tried to reproduce it, then I used the same input data in the Octave BICGSTAB: the approximation X was all composed by NaN, and also for the residuals. Then I study in deep the Octave code and the reference book. I discovered that ih this situation there was some diviosn by zero probably caused by the Non-Symmetric Lanzcos Biorthogonalization (that is used in both BICG and BICGSTAB).
    With this Lanczos algorithm is possible the situation of a "serious breakdown" (not going into deep details, but it happens when a certain inner product is zero, instead of a non-zero number), and there is not easy examples to reproduce it. The refernce book suggest two possibility to avoid this situation: to restart the Lanczos biortho, or to use a difficult variant called "Look-Ahead Lanczos".
    I tried to study the "Look-Ahead" but in the Saad book is only mentioned the general idea, and it is not simple to implement and to understand quickly. Then I tried the second strategy it is not efficient (an easy experiment to try this strategy is: run BICGSTAB to find the last iteration without NaN, then run another time the BICGSTAB with the output X as initial guess), indeed the method does not converge and the residual does not decrease.
    Then (since in Matlab is not mentioned if it uses or not the Look Ahead) in the BICGSTAB implementation that I wrote I don't use any strategy and I set "flag = 4" when the guilty inner product is zero.
    To see an example of this situation in the Octave BICGSTAB see the file example_bicgistab.m
Instead, during the implementation of BICG there was not (for now) strange situations. For both BICG and BICGSTAB I tried the same tests that are given in their Octave implementations and I add some test to try the "flag =2" and "flag = 4" situations.


I know that to fix PCG, GMRES, BICG and BICGSTAB I used a lot of time. This because in their implementation emerged a lot of strange situations that needed a deep revision of these algorithms and many times also a discussion (with the mentors or the community) of how to proceed to solve particular problems. Also the mentors told me that they expected less time/work to improve and make a revision of the codes already implemented.
They also told me to spend the remaning time of the SoCiS to implement some methods not present in Octave. For this, I want to ask you if there is one of the suggested algorithms (minres, symmlq, symmetr, tfqmr, bicgstabl or lsqr) that you prefer to implement.

As usual, if there is something wrong/not clear, please contact me via e-mail or write a comment.


by Cristiano Dorigo ( at August 02, 2016 06:52 AM

August 01, 2016

Abhinav Tripathi


Work in the past month has been distributed over symbolic and pytave side.

The work that was merged to symbolic in this time was mostly to convert objects from python to octave types. PR #474 and #491 worked in this regard.  PR #490 helped to make the tests compatible with Pytave IPC. Then we also added Pytave IPC to the build system (travis) with the PR #489.
Also, we worked to improve the passing of variables to python by using pytave. PR #465 has mostly achieved the goal. But it is failing on travis. The cause was finally detected in PR #502 to be the passing of struct to pytave. I couldn’t get much time to work on this as I was working on pytave. And hopefully when pytave matures then this error will be resolved in pytave itself.

In pytave, there were some bug fixes which caused octave to segfault. In the PR #12, the issue was addressed when an octave value could not be converted to python which threw an exception and crashed octave. In the PR #15, the issue where the exception text from Python could not be properly extracted was handled. This was later improved with the PR #30 for the cases of crash with Python 3. In PR #29, the issue was handled when a large number (integer) was being converted from python to octave then boost threw positive_overflow exception.

Also, there were some feature improvements too. PR #13, added the feature that checks for a function called via pycall into main namespace and builtins namespace too for existence. After pyobject was merged into pytave, most of the work was related to pyobject. PR #17 completed the functionality of returning pyobject from py* functions by adding the same mechanism to pycall which was implemented in pyeval. PR #16 added the feature to check length of a pyobject. PR #25 dropped the existing conversion mechanism which converted python list, tuples and dicts to octave matrices and structs. Now they return a pyobject instead.

There are a few PRs currently under review. PR #33 removes some of the unused code from pycall and adds a few BIST tests which showcase changes due to PR #25. PR #32 fixes some doctests which are also affected by the changes in the same PR. PR #31 aims at adapting all the methods of pyobject to the changes which dropped the conversion of lists, dicts and tuples. In PR #24, we are trying to add a new feature to allow users to specify the namespace in which pyexec/pyeval would run their code. Hence allowing isolation of different codes to have the same names/functions by running in different namespaces. Currently it is buggy and doesn’t actually work as expected. The symbols in a local namespace are also accessible from the global namespace. Moreover any 2 local namespace are sharing the same symbols. Currently it just works in the case where 1 global and 1 local namespaces are required. In such case, they are having separate symbol tables (dictionary actually as it is python).

During this time we had one setback. We tried a lot to build pytave on windows but failed. It has been stopped currently but we will resume it after GSoC ends. Tatsuro helped a lot in this. Thanks to him we have gotten really close to successfully building pytave on windows. The complete details about what was tried and which errors occurred can be found on my wiki page where I described Building Pytave on windows.

The project has been going smoothly because of a great participation from the mentors and octave community. Hopefully we will be able to have a working release of symbolic with pytave IPC (atleast for linux) before the end of GSoC.

by genuinelucifer at August 01, 2016 06:04 PM

Francesco Faccio

Performance test


As klu library is now linked to ode15i, I could test the performance of the code both in the dense and sparse case. I tested my oct-file against Matlab ode15i m-file, the C implementation of sundials examples and the m-file relying on the mex interface of sundials. Although the code still contains unefficient use of memory, I think that the results are pretty interesting.

Before going into details, I have to specify that C and mex implementations are written for solving these specific problems, so their code is highly optimized (specific Jacobian matrices are constructed directly inside functions required by IDA, while in my oct-file I receive a function that when evaluated it returns the Jacobian matrices I have to pass to IDA). In other words, solving a problem requires more time than solving the problem.

The first test regards the solution of Robertson chemical kinetics problem, in which the Jacobian is a 3x3 dense matrix. I found a solution for RelTol = 1e-4, AbsTol = [1e-8, 1e-14, 1e-6], tspan = [0, 0.4] and I got (in mean):

0.015s for the oct-file
0.13s for the mex-file
0.002s for the C-implementation
0.19s for Matlab ode15i.

The second test regards the solution of a 2D heat equation. In this problem the Jacobian was a 100x100 sparse matrix with 330 non zero elements. I tested with RelTol = 1e-5, AbsTol = 1e-8 and tspan =[0, 0.01, 0.02, 0.04, 0.08, 0.16, 0.32, 0.64, 1.28, 2.56, 5.12, 10.24] and got:

0.42s for the oct-file
0.005s for the C-implementation
0.45s for Matlab ode15i.

Unfortunately Sundials doesn't provide a mex interface for klu, so I couldn't test a sparse mex file version for this problem.

We can now make some considerations:
-Oct-file is 10 times faster than mex-file and this confirms that mex files cannot reach the performance of an equivalent (in this case, of a more general) otc-file
-Our solver is already faster than Matlab m-file
-The C implementation is the fastest, but code is designed to solve only a specific problem. We can try to get closer to this speed.

In order to improve the performance of the solver, the next step is to try to reduce looping and the use of memory and to write a DAE class, using functors and lambda functions to avoid the use of global pointers.

by Francesco Faccio ( at August 01, 2016 03:54 AM

July 26, 2016

Chiara Segala

odexprb34 solver

As I did with the odexprk23 solver from the Runge-Kutta family, I implemented the official exponential fourth order method with adaptive time stepping from Rosenbrock family described in my postAlso in this case I used the function expmv and I wrote a code for the stepper (exp_rosenbrock_34), in which there are three calls to expmv function. The first two for the calculation of the second and third stage
U_{n2} = u_n + \frac{1}{2} h_n \varphi_1( \frac{1}{2} h_n J_n) F(t_n,u_n) + \frac{1}{2} h_n \varphi_2(\frac{1}{2} h_n J_n) \frac{1}{2} h_n v_n ,
U_{n3} = u_n +  h_n \varphi_1( h_n J_n) [F(t_n,u_n) + D_{n2}] + h_n \varphi_2 ( h_n J_n) h_n v_n ,
[Atilde, eta] = (1/2*dt, Jn, [Fn, 1/2*dt*vn]);
X = myexpmv (1/2*dt, Atilde, [zeros(size (Atilde)-1, 1); 1/eta])(1:d, :);

U(:, 2) = x(:) + X;

[Atilde, eta] = augmat (dt, Jn, [Fn+D(:, 2), dt*vn]);
X = myexpmv (dt, Atilde, [zeros(size (Atilde)-1, 1); 1/eta])(1:d, :);

U(:,3) = x(:) + X;

The third call to expmv is for both the solution to the next step, and the lower order solution for the estimation of the error

u_{n+1} = u_n + h_n \varphi_1(h_n J_n) F(t_n,u_n) + h_n \varphi_2(h_n J_n) h_n v_n \dots
\dots + h_n \varphi_3(h_n J_n) [ 16 D_{n2} - 2 D_{n3}] + h_n \varphi_4(h_n J_n) [ -48 D_{n2} + 12 D_{n3}] ,
\hat{u}_{n+1} =  u_n + h_n \varphi_1(h_n J_n) F(t_n,u_n) + h_n \varphi_2(h_n J_n) h_n v_n + h_n \varphi_3(h_n J_n) [ 16 D_{n2} - 2 D_{n3}] ,
v = [Fn, dt*vn, 16*D(:, 2)-2*D(:, 3), -48*D(:, 2)+12*D(:, 3)];
w = [Fn, dt*vn, 16*D(:, 2)-2*D(:, 3)];
[Atilde, eta] = augmat (dt, Jn, v, w);
p = size (v, 2);
q = size (w, 2);
ep = [[zeros(p-1, 1); 1/eta; zeros(q, 1)], [zeros(p+q-1, 1); 1/eta]];
X = [eye(d, d), zeros(d, p+q)]*myexpmv(dt, Atilde, [sparse(d, 1), sparse(d, 1); ep]);
X1 = X(:, 1);
X2 = X(:, 2);
x_next = x(:) + X1;
x_est = x(:) + X2;

I wrote also a code for the solver (odexprb34), according to odepkg standards and  I tested the order of the method with the usual example described in my post.

by Chiara Segala ( at July 26, 2016 07:40 AM

July 18, 2016

Juan Pablo Carbajal


¿Qué es una variable?

En lenguajes de programación, una variable es un símbolo (e.g. una letra:"x" o una palabra:"valor1" ) que se utiliza para guardar información.
En GNU Octave es muy fácil crear variables, simplemente
asignamos la información deseada a la variable. 

Por ejemplo:

x = 1

guarda el valor "1" en la variable llamada x. El signo = se utiliza para asignar valores, uno podría leer el código de arriba como "equis guarda el valor uno".

Pregunta 1: 
¿Cúal es el valor de la variable y en el siguiente código?

x = 1;
y = x + 1

Las variables en GNU Octave no solo pueden guardar un valor, sino que también pueden guardar listas de valores (vectores y matrices o tablas). 


x = [1,2,3]
y = "Hola Salta la linda!"

La variable x guarda un
vector fila
con componentes 1, 2 y 3. La variable y guarda un arreglo de caracteres, en otras palabras: un vector de letras! 
A estos objetos se los llama "strings" en Inglés y yo prefiero esta palabra que la palabra "arreglo", porque con esta última parece que estamos haciendo negocios sucios (chiste ;D).

Pregunta 2:
¿Cúal es la diferencia entre u y v en el siguiente código?

u = [0.1 -1 5];
v = [0.1; -1; 5];

Bueno, con estos introdujimos rápidamente variables en GNU Octave. Hagan sus preguntas! (y no se olviden de experimentar! También hay mucha ayuda en la web).

Hasta la próxima!

PS: No se olviden que todo el código de las clases se puede ejecutar en una sesión de GNU Octave.

by Juan Pablo Carbajal ( at July 18, 2016 07:57 PM

Vectores y matrices

En esta clase introducimos dos elementos fundamentales del calculo númerico, el vector y la matriz.

En la clase de Variables vimos que una variable puede guardar varios valores al mismo tiempo. Por ejemplo
u = [0,1,2,3,4,5]
Nota: los decimales en Octave se escriben con punto, no con coma. Es decir 1/2 se escribe 0.5 y no 0,5.

En el ejemplo, u es una variable que tiene 6 elementos o componentes. A estas variables, que organizan los valores de forma horizontal, las llamamos vectores fila. Las llamamos vectores fila porque los valores están organizados en una fila. Si ejecutamos el siguiente código:
size (u)

Octave nos dice:
ans =
1 6

OK, vamos por partes... Primero el comando size(u) nos devuelve el tamaño de u (size quiere decir tamaño en Inglés). Y la respuesta fue 1 6 ¿Como leemos la respuesta? 
El tamaño de un vector se indica de la forma "filas, columnas". En este caso 1 6 nos dice que u tiene 1 fila y 6 columnas.

Pregunta 1:
¿Cuantas filas y cuantas columnas tienen los siguientes vectores?
a = [0,1,2]
b = [0.1 0.5 -0.1]
c = [-3,2]
d = [83]
e = [0,1] 
tiene 1 fila y 2 columnas.

¿Cómo armamos un vector columna? Es decir, ¿Cómo organizamos valores de forma vertical?
v = [0;1;2;3]
Este es un vector columna, donde los valores se organizan verticalmente. Quiere decir que para crear un vector columna tenemos que separar las componentes con punto y coma.

Pregunta 2:
¿Cuantas filas y cuantas columnas tienen los siguientes vectores?
a = [0;1;2]
b = [0.1; 0.5; -0.1]
c = [-3;2]

Para acceder a la 2da componente del vector a = [0 -1 2] ejecutamos lo siguiente
y Octave responde
ans =

Supongamos que a(2) == -1 fue un error y deseamos cambiar ese valor por 1. Simplemente hacemos
a(2) = 1

Con toda esta información ya puedes armar tablas, o en términos matemáticos, una matriz. Una matriz no es otra cosa que la concatenación de filas o columnas ¿A alguien se le ocurre cómo armar una matriz? Por ejemplo ¿Cómo crearías la siguiente matriz en Octave?
0 1
1 0
0 1

Noten que uds. ya usan las matrices es su vida diaria, por ejemplo una tabla como esta:
 1. Literatura4
 2. Química 5
Puede entenderse como una matriz que tiene 2 columnas (Materia, Nota) y dos filas (cada materia).

¿Preguntas? Escriban en la sección discusión o en los comentarios más abajo.

En la próxima clase veremos como representar la posición de un objeto como un vector y su movimiento en el tiempo como una matriz.

by Juan Pablo Carbajal ( at July 18, 2016 07:57 PM

¿Para qué me sirven los vectores?

En la clase de vectores comentaba como se puede crear un vector en GNU Octave, pero quedó en el tintero la motivación para el uso de vectores ¿Para qué me pueden servir esas listas de numeritos que llamamos vectores?

Cartas de amor...

Supongo que alguna vez han mandado una carta por correo postal. Para que la carta llegue a destino necesitamos proveer cierta información además del nombre del destinatario. Supongamos que la carta la enviamos dentro de Argentina. Necesitamos indicar la provincia a donde va la carta, la cuidad, la calle y el número de casa. En principio podríamos organizar la información en de esta forma

carta = [ "Buenos Aires","Lomas de Zamora", "Ceferino Namuncura", "150" ];

Ok, esto ya se parece a un vector pero las componentes son strings en vez de números. No hay drama, simplemente podemos indicar la provincia con un número del 1 al 23 (Argentina tiene 23 provincias) y hacer algo parecido con las ciudades y las calles. El número de casa lo podemos usar directamente. Es decir que el destino de nuestra carta podría representarse con el vector

carta = [ 23, 1352, 12345, 150 ];

donde cada componente fue reemplazada por un número según una tabla especificada. La primera componente del vector nos indica la provincia, la segunda nos muestra las ciudad, la cuarta nos indica la calle y la quinta componente la altura de la calle ¡Estas cartas viven en un espacio de 5 dimensiones!
Un conjunto de estas cartas se puede organizar en una lista de vectores, una matriz:

cartas = [ 23, 1352, 12345, 150; ...
          5, 130, 4, 756; ...
          12, 7, 2341, 29 ];  

Cada vector fila de esa matriz representa una carta. La primera columna de la matriz nos muestra las provincias, la segunda columna nos muestra las ciudades, la cuarta nos indica la calle y la quinta columna la altura de la calle. 

Panadero, quiero pan!

Ok, quizás es ejemplo fue muy abstracto. El siguiente ejemplo me lo pasaron por Facebook, también es abstracto, pero quizás es más fácil de entender. 

Pensemos por un momento en recetas de cocina, en particular en la lista de ingredientes de cada receta. En la tabla puse dos ejemplos de ingredientes para hacer pan

 Pan de cazuela  Pan fuerte
 1.Harina (gr)
 250 300
 2.Sal (gr) 5 20
 3.Levadura (gr) 5 20
 4.Aceite (ml) 0 50
 5.Agua (ml) 175 50
 6.Miel (ml) 0 100
 7.Huevos 0 0

Esta tabla ya nos muestra una representación de estas recetas como vectores. En este caso son dos vectores columna que en GNU Octave podríamos anotar como

recetas = [ 250 300; 5 20; 5 20; 0 50; 175 50; 0 100; 0 0 ];

Los ingredientes del pan de cazuela son la primera columna de esta matriz y los del pan fuerte son la segunda. Noten que estos vectores son de 7 dimensiones (7 ingredientes), pero vale preguntarse si la dimensión Huevos tiene sentido en el espacio de recetas de pan... ¿Hay recetas de pan que utilizan huevos?

Los dos ejemplos anteriores son bastante abstractos y si los has entendido me puedo imaginar que ya tienes muchas ideas de como representar otras cosas como vectores o matrices. 

Trayectorias espaciales

El siguiente ejemplo es paradigmático en el campo del cálculo numérico, y es muy relevante para este curso: representar la posición de objectos en el espacio!

Imagina una hormiga caminando sobre una hoja de papel A4 (ancho: 21cm, largo: 29.7cm). Para nuestra conveniencia hemos marcado una de las esquinas de la hoja como punto de referencia. En cualquier instante podemos leer la distancia entre la hormiga y la esquina de la hoja en las dos dimensiones del papel: ancho y largo. Ponemos la hormiga en un punto del papel y una vez por segundo miramos donde se encuentra. Podríamos observar algo así:

trayectoria hormiga

Figura 1: Posiciones de la hormiga en el papel. Cada punto representa la posición de la hormiga en las coordenadas del papel observadas cada 1 segundo. La línea punteada es la trayectoria de la hormiga observada con más frecuencia.

La posición luego del primer segundo la guardamos en una matriz con una fila y con dos columnas (un vector fila!)

p(1,1:2) = [ 10.7 14.8 ];

Luego de otro segundo obtenemos la nueva posición y la guardamos en la segunda fila de la matriz:

p(2,1:2) = [ 11.1 14.8 ];

Y así sucesivamente. Lego de 10 segundos, la matriz tiene 10 filas y 2 columnas. Esta matriz representa la trayectoria de la hormiga a intervalos de 1 segundo. Cada vector fila nos indica la posición en un dado instante de tiempo.

p = [ 10.7 14.8; ...
      11.1 14.8; ...
      11.4 15.1; ...
      11.2 15.0; ...
      11.0 15.3; ...
      11.8 15.8; ...
      11.6 16.8; ...
      11.8 17.0; ...
      11.0 17.3; ...
      10.8 17.2];

Creo que estos tres ejemplos deberían entusiasmarte para pensar qué cosas pueden representarse utilizando vectores y matrices y cómo hacerlo. Vale la pena preguntarse como podemos utilizar estas representaciones. Por ejemplo, usando la matriz de la hormiga ¿Cuál es la distancia entre la primera y la ultima posición de la hormiga?

¿Se te ocurren otros ejemplos? ¡Anótalos en los comentarios!

Nota: El archivo adjunto es un script de GNU Octave que genera la trayectoria de la hormiga y crea la figura que se muestra en esta clase.

by Juan Pablo Carbajal ( at July 18, 2016 07:57 PM

¡Animarse a animar!

Todos disfrutamos de una buena animación, aquí veremos como podemos animar nuestros datos en GNU Octave. Les presento algunos ejemplos a modo de motivación:

Para lograrlo vamos a conocer la función comet y, para hacer nuestras propias animaciones, vamos mirar un poco en detalle los objetos generados por la función plot.

La función comet

En su forma más intuitiva, esta función toma como entradas tres argumentos: 2 vectores (o listas de valores) y una valor escalar (un número). Los dos vectores son los pares de valores (x,y) que se irán dibujando uno a uno. El argumento escalar define el tiempo de espera antes de mostrar el siguiente punto. Probemos lo siguiente:

T   = 10;
fps = 25;
t   = linspace (0, T, fps*T).';

l = 0.5;
w = 2*pi*0.5;

R = exp (-l * t);
x = R .* cos (w * t);
y = R .* sin (w * t);

fig = figure (1);
set (fig, "Name", "Espiral");

comet (x, y, 1/fps);

Si ejecutas el código, deberías ver algo como lo que se muestran en el siguiente video:
Pregunta 1:
¿Entiendes todo el código?

Anatomía de una animación
Para realizar cualquier tipo de gráfico en Octave necesitamos una figura (la ventana donde vemos el gráfico, creada con la función figure). Esta figura contiene un par de ejes coordenados (un "hijo", que se guarda en el campo de la figura llamado, en inglés, "children"); que en Octave se llaman axes y se que crean con la función homónima. Dentro de estos ejes podemos tener objetos gráficos, como líneas, puntos y polígonos.

Como se darán cuenta, las figuras pueden tener mas que un hijo, es decir una figura puede contener varios ejes permitiendo poner varios gráficos en la misma figura.

Resumiendo: una figura contiene un par (o más) de ejes coordenados, los cuales a su vez contienen objetos gráficos
Veamos un ejemplo

close all
fig = figure (1, "Name", "Figura 1");  # Figura con nombre
ax  = axes ();                         # fig adopta ejes automáticamente.

x = randn (10,1);                      # Graficamos algo
h = plot (ax, x, '.');                 # y lo ponemos en los ejes creados

Por supuesto, todo esto (excepto el nombre de la figura) se hace automáticamente simplemente ejecutando plot (x,'.') ¡Verifícalo!
En la última línea del código, le pedimos a la función plot que nos devuelva el objeto gráfico. Podemos ver las propiedades de este objeto ejecutando get(h) y obtendremos una lista enorme de propiedades del objeto gráficos y sus valores actuales. Fíjate que en esa lista hay campos con los nombres ydata y xdata. Estos campos contienen los puntos que se muestran en el plot. Si modificamos estos campos, modificaremos el gráfico. Por ejemplo podemos desplazar los valores un lugar hacia la derecha haciendo

t = get (h, "xdata");
set (h, "xdata", shift (t, -1));

Pregunta 2
¿Puedes modificar el código para realizar la animación que se muestra abajo (click para ver animación)?
gif animado

Hemos visto los principios de animación: prepara el gráfico y luego modifica su contenido.
En general es mala idea rehacer el gráfico a cada instante, pues es muy lento y las animaciones quedan mal.
Puedes animar cualquier propiedad de los objetos gráficos, ejes y figuras.

Si tienes alguna idea y necesitas ayuda no dudes en preguntar en el foro

by Juan Pablo Carbajal ( at July 18, 2016 07:57 PM

Ecuaciones lineales? ¡No problem!

Uno de los problemas más recurrentes en el cálculo numérico es el de resolver sistemas de ecuaciones lineales. ¿Qué es una ecuación lineal? Si no conocen este nombre todavía, no se preocupen. Vamos a ilustrar todo desde el principio y de manera básica. También veremos una aplicación fundamental de las matrices.

La ecuación lineal

Los problemas lineales son muy comunes en la vida cotidiana. Son tan comunes que en el colegio nos enseñan a resolverlos con una regla: la famosa, y a veces odiada, regla de tres simple.
Esta regla resume la facultad de resolver proporciones. Ejemplo:

Problema 1:
Si un paquete de harina pesa 10 kg cuando está lleno ¿Cuánto pesa cuando solo queda un tercio del paquete?

En cada generación y en cada escuela se enseñan distintas formas de plantear este problema, pero la idea fundamental es que el peso del paquete es proporcional a la cantidad de harina. Esto quiere decir que mientras más(menos) harina, más(menos) pesa el paquete. Pero no solo implica que las dos magnitudes aumentan juntas, sino también que aumentan juntas de una manera particular: cualquier aumento de harina implica un aumento de igual proporción en el peso del paquete. Esto es lo que se llama una relación lineal

Donde P es el peso del paquete y V es la fracción de llenado del paquete (1=lleno, 0=vacío). El factor k es la proporcionalidad entre P y V. En nuestro problema la proporción es:

Otros ejemplos de proporcionalidades son la relación entre el volumen y la masa de un líquido (la proporcionalidad se llama densidad del líquido); entre el volumen y la altura de un cilindro (la proporcionalidad está dada por el área de la base del cilindro); entre el la altura de un fajo de billetes y la cantidad de billetes (la proporcionalidad es el espesor de cada billete), etc... seguro pueden imaginarse en un montón de relaciones lineales.

Una relación lineal que aparece en física del colegio secundario es la relación entre distancia recorrida en linea recta y rapidez.

Donde x es la posición final, x0 la posición inicial, v es la rapidez y Δt el tiempo transcurrido. Como ven, la distancia recorrida es proporcional al tiempo transcurrido; también podemos decir que la distancia recorrida es proporcional a la rapidez.

Sistemas de ecuaciones lineales

Consideren el siguiente problema.

Problema 2:
Encuentra los números tales que: si sumamos el doble del primero con la mitad del segundo, obtenemos 10 y si restamos el primero al segundo obtenemos 5.

Escribamos este problema en ecuaciones

Los números que estamos buscando son x,y  que al momento representan incógnitas. Podemos organizar las ecuaciones de la siguiente manera:

donde el punto indica producto escalar entre el vector fila y el vector columna (ver clase sobre el producto escalar). Ok, si miran estas ecuaciones con fuerza, se darán cuenta que se pueden simplificar un poco mas:

Es decir, los sistemas de ecuaciones lineales se pueden escribir como la multiplicación entre una matriz (que tiene los factores de las ecuaciones lineales) y el vector de incógnitas, igualada al vector de datos. Es decir, en general, los sistemas lineales se pueden escribir de la forma siguiente


Donde A es la matriz de factores, x el vector de incógnitas y b el vector con los datos. A los vectores los marqué con negritas, pero en el futuro no lo voy a hacer, excepto para enfatizar o para evitar confusión entre escalares y vectores.

Resolviendo sistemas de ecuaciones lineales

Cuando trabajamos con escalares (con números) sabemos como despejar una incógnita: se aplican las operaciones inversas hasta que la incógnita queda despejada. Si tenemos una ecuación lineal, podemos despejar de la siguiente forma

Como ven, multiplicamos ambos lados de la igualdad por la inversa multiplicativa de a. Pero para resolver el sistema de ecuaciones (1) ¿Cuál es la inversa multiplicativa de una matriz? La respuesta es un poco larga y no la vamos a discutir hoy (pueden investigar en la internet sobre "matriz inversa"). Lo bueno es que GNU Octave sabe como calcular la inversa de una matriz, así que van a poder resolver sistemas lineales antes de aprender la teoría. Con suerte esto los va a motivar para estudiar y entender la teoría.

Lo que si es necesario saber es que una matriz, a diferencia de un escalar, no siempre tiene inversa multiplicativa. En otras palabras ¡Hay problemas lineales que no tienen solución!

Pregunta 1:
¿Podés construir un sistema de dos ecuaciones lineales que no tenga solución?

Resolviendo el problema calculando la inversa

Si estamos seguros que el problema tiene solución, entonces la inversa de la matriz de factores tiene que existir y en Octave la calculamos con la función inv(). En nuestro ejemplo

A  = [2 0.5; -1 1];
b  = [10; 5];
Ai = inv (A);
x  = Ai * b

x =


Pregunta 2:
Si ahora tratás de calcular la inversa del problema imposible de la Pregunta 1 ¿Qué obtenés?

Resolviendo el problema sin calcular la inversa

Otra forma de resolver el problema en Octave es utilizar el operador \ (división izquierda). Este es un operador inteligente y elige la mejor forma de resolver el problema, desde un punto de vista numérico (ya veremos que significa esto). En la mayoría de los casos este operador no calcula la inversa explícitamente, puesto que dicha operación es un poco problemática. El código para resolver el sistema (1) con este operador es:

A = [2 0.5; -1 1];
b = [10; 5];
x = A \ b

x =


Pregunta 3:
Si ahora tratás de resolver el problema imposible de la Pregunta 1 ¿Qué obtenés?

¿Se chocan o no se chocan? Problemas de encuentro

Advertencia: el siguiente problema puede traerte recuerdos del colegio secundario :D.

Problema 2:
Un móvil A sale de un punto con una rapidez de 60 Km/h y luego de 30 minutos parte otro B del mismo punto a una rapidez de 80 Km/h ¿Cuándo alcanzará el segundo móvil al primero y a qué distancia del punto de partida se encuentran?.

Primero planteamos las ecuaciones de cada móvil:

Donde el subíndice indica a que móvil nos referimos. Estas ecuaciones nos dicen donde estarán los móviles A y B luego de un tiempo t. Vamos a reescribir estas ecuaciones para poner el problema en forma matricial (e.d. usando matrices). Primero notemos que nuestro objetivo es conocer en que tiempo t, las posiciones son iguales, e.d. xA = x= x
Reemplazando en las ecuaciones de arriba y poniendo las incógnitas del mismo lado obtenemos:

Ahora la matriz de factores es evidente y el problema se puede escribir como

El código de Octave que resuelve el Problema 2 es:

vA  = 60; t0A = 0; 
vB  = 80; t0B = 0.5;   % 30 minutos son media hora
A   = [ vA -1; vB -1 ];
b   = [ vA*t0A; vB*t0B ];
sol = A \ b

sol =

Los móviles se encuentran 2 horas luego de la partida de A (t=2) y se encuentran a 120 kilómetros del punto de partida (x = 120). Para verificar que la solución es correcta podemos reemplazar en las ecuaciones iniciales (o multiplicar A*sol y verificar que obtenemos b).

¡A resolver problemas!
Si tienen preguntas participen en el foro de discusión.

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

Funciones para generar valores

GNU Octave provee un gran número de funciones que son útiles para resolver nuestros problemas numéricos.
No tiene sentido enumerar aquí al enormidad de funciones de Octave, para eso está el manual y la ayuda interactiva que describimos en una clase pasada.  Los invito a que usen el foro de discusión para hacer preguntas del tipo: ¿Con qué función puedo hacer XXXX?

En nuestras clases irán apareciendo funciones, cuya funcionalidad la explicaremos en su debido momento (en realidad una vez que saben el nombre de la función pueden ir al manual o utilizar la ayuda interactiva).
En esta clase, introduciré dos funciones que nos permiten generar valores que incrementan de forma lineal o logarítmica.

Valores igualmente espaciados

La función linspace toma tres argumentos de entrada:

x = linspace (inicio, final, numero);

El argumento inicio indica el primer valor y el argumento final el último. El argumento numero indica cuantos valores generar entre inicio y final, todos ellos igualmente espaciados.

valores igualmente espaciados

Pregunta 1:
Si quiero valores que van desde 0 a 16 con un espacio entre ellos de 2, puedo ejecutar el código:

x = 0:2:16;

¿Cómo hay que llamar a linspace para generar los mismos valores?  ¿Puedes explicar la solución?

Valores logarítmicos

La función logspace funciona igual que la función linspace pero en vez de distribuir valores igualmente espaciados, distribuye los exponentes de los valores. Por ejemplo

x = logpspace (-1, 3, 5);

genera los valores 

x =

   1.0e-01   1.0e+00   1.0e+01   1.0e+02   1.0e+03

Si observan bien, pueden notar que es el exponente el que está igualmente espaciado (-1,0,1,2,3). Es decir, x(1) es 1x10-1 mientras que x(end) es 1x103.

valores espaciados logaritmicamente
Pregunta 2:
¿Qué valores obtengo si ejecuto el siguiente código?

x = log (logspace(-3,3,10))

Pregunta 3:
¿Cúal es el comando que genera valores logarítmicos entre 0.5 y 125?

No duden en dejar sus preguntas, sugerencias y comentarios en el foro!

¡Hasta la próxima!

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

Repitiendo instrucciones: lazos "for"

Seguimos con el tema solicitado por Matías: control de flujo. Ahora vamos a ver como repetir un grupo de instrucciones utilizado el lazo for (yo les digo "loops", no lazos).
La idea es sencilla: uno tiene que repetir una acción muchas veces y no quiere escribir el código correspondiente a cada una de ellas.
Por ejemplo, supongamos que queremos sumar todos los valores en un vector:

x = rand (100,1)

El vector x contiene 100 números entre 1 y 0 elegidos al azar (ver ayuda de la función rand y preguntar en el foro si algo no se entiende). El objetivo es sumar todos estos valores. La forma ingenua de hacer esto sería escribir las 100 líneas:

y = x(1);
y += x(2);
y += x(3);

Y así hasta sumar todo los números. Está claro que esto no nos lleva muy lejos. Los lazos nos permiten repetir un grupo de comandos sin tener que escribir todas estas líneas. Utilizando el loop for podemos hacer la suma:

x = rand (100,1);
y = 0;
for i = 1:100
    y += x(i);

El loop define una variable de iteración i e indica que esta variable va desde 1 a 100. Esta variable se utiliza luego para indexar el vector x y acceder a cada uno de sus valores. Para vectores de otros tamaños el programa es casi el mismo

x = rand (500,1);
y = 0;
for i =1:500
    y += x(i);

Pregunta 1:
¿Cómo podemos escribir el loop para que funcione con cualquier tamaño del vector x?

El precio de la iteración

GNU Octave se basa en un paradigma de programación vectorial. Esto quiere decir que operaciones entre vectores y matrices son veloces y las iteraciones explicitas, como el loop for, son lentas. La idea es que en vez de hacer una iteración, podemos hacer operaciones muy rápidas utilizando más memoria. Es decir, cambiamos  iteraciones por espacio en memoria.

El ejemplo anterior se puede escribir de manera vectorial de la siguiente manera:

x = rand (500,1);
y = sum (x);

Este ejemplo en particular no cuesta más memoria que el loop porque en ambos casos construimos el vector x. Sin embargo, el loop podría ahorrar memoria si lo escribiésemos así:

y = 0;
for i =1:500
    y += rand ();

En este caso nunca construimos el vector x y ahorramos esa memoria. Esto podría ser útil en el caso que necesitas una cantidad tan grande de números aleatorios que no entran en la memoria de tu PC (si tienes 1Gb de memoria puedes tener alrededor de 1 millón de estos números). Sin embargo esto tiene un precio en tiempo y lo podemos verificar utilizando las funciones tic y toc. El siguiente programa compara los tiempos de ejecución de las dos formas de hacer esta suma:

N = 1000;
disp ('Loop')
y = 0;
for i =1:1000
    y += rand ();
disp ('Vectorizado')
x = rand (1000,1);
y = sum (x);

En mi PC esto da 

Elapsed time is 0.0164838 seconds.
Elapsed time is 7.60555e-05 seconds.

¿Cuánto da en la tuya? Puedes descargar el script al pié de la página.

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

¡Ayuda! ¡Socorro!

GNU Octave tiene un sistema de ayuda integrado. En esta clase veremos como usar el sistema de ayuda y otras fuentes de información. La mayoría de esta ayuda está en Inglés, y es importante saber leer este idioma a la hora de programar. Sin embargo muchas de las personas que trabajan en GNU Octave hablan español y están dispuestos a ayudar!

En la terminal

Mientras trabajamos en Octave hay varias formas de pedir ayuda sobre una función. Cada forma se aplica en diferentes casos.

Si sabes el nombre de la función tu mejor amigo es el comando help. Supongamos que sabemos que la función que queremos usar es la arco tangente (la inversa de la función tangente). En Octave, hay dos funciones que calculan la arco tangente: atan y atan2. Veamos la ayuda de cada función (en Octave 3.8.2):

> help atan
'atan' is a built-in function from the file libinterp/corefcn/

 -- Mapping Function: atan (X)
     Compute the inverse tangent in radians for each element of X.

     See also: tan, atand.

> help atan2
'atan2' is a built-in function from the file libinterp/corefcn/

 -- Mapping Function: atan2 (Y, X)
     Compute atan (Y / X) for corresponding elements of Y and X.  Signal
     an error if Y and X do not match in size and orientation.

     See also: tan, tand, tanh, atanh.

¿Entiendes la diferencia entre estas dos funciones?

A medida que vayas acumulando más experiencia con Octave, podrás recordar los nombres de la funciones que te resultan más útiles. Sin embargo puede ser que a veces no te acuerdes el nombre completo de la función pero solo el principio. En este caso podemos utilizar el autocompletar para obtener el nombre de las funciones. Supongamos que nos acordamos que una función empezaba con "sum" pero no recordamos el resto. en la terminal de Octave ingresamos este texto y luego presionamos dos veces seguidas la tecla TAB. El resultado es el siguiente

> sum
sum     sumblk  summer  sumsq

Octave nos da una lista de todas las funciones que empiezan con "sum" ¡Esto es de gran ayuda para encontrar nombres de funciones!

Otra situación bastante común es cuando no sabemos el nombre de la función pero sabes que el texto de ayuda contiene una frase/palabra clave que si conocemos. En este caso una herramienta muy útil es lookfor (significa "buscar" en Inglés). Consideremos un ejemplo. Supongamos que queremos buscar si Octave tiene una función para calcular la suma de una lista de números. Como frase a buscar (en Inglés) elegimos "sum of" (suma de) y ejecutamos:

> lookfor "sum of"
trace    Compute the trace of A, the sum of the elements along the main diagonal.
hypot    Compute the element-by-element square root of the sum of the squares of X and Y.
cumsum   Cumulative sum of elements along dimension DIM.
sum      Sum of elements along dimension DIM.
sumsq    Sum of squares of elements along dimension DIM.
md5sum   Calculate the MD5 sum of the file FILE.

Octave nos responde con una lista de nombres de funciones y la primera linea de la ayuda correspondiente. Una vez que tenemos los nombres de las funciones podemos utilizar el comando help para investigar en mas detalle que es lo que hace cada función.

Finalmente, puedes acceder al manual de Octave en la terminal utilizando el comando doc.

En internet

Como siempre estar conectado hace la vida más fácil. En internet podemos utilizar el manual online de octave para buscar ayuda sobre las diferentes funciones y funcionalidades. Es un buen lugar para consultar y puedes descargar la versión en PDF si acaso no tienes conexión todo el tiempo. He puesto algunas referencias en castellano en la sección de lecturas recomendadas. También puedes visitar el Wiki de octave.

Otro recurso excepcional es el chat room #octave en Aquí puedes charlas con otros usuarios y desarrolladores de todas partes del mundo. No es raro ver conversaciones en castellano, italiano y alemán. Puedes accerder al chat room utilizando un sitio web o cualquier cliente de IRC (por ejemplo Pidgin).

Finalmente, si prefieres conversaciones por e-mail puede utilizar la lista de mails

¡Utilizando todos estos recursos no puedes sentirte solo a la hora de aprender GNU Octave!

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

¡Por favor Sr. Kirchhof, resuelva mi circuito!

En esta entrega vamos a utilizar lo que discutimos sobre ecuaciones lineales para resolver problemas con circuitos eléctricos.
En la figura siguiente vemos un circuito compuesto de 4 resistores y una fuente de voltaje.
El dibujo de la izquierda muestra el circuito con la fuente conectada al nodo 1 y el nodo 4 conectado a tierra. El dibujo de la derecha muestra un circuito equivalente. Nosotros podemos controlar la fuente y por lo tanto conocemos el voltaje Vi entre los nodos 1 y 4. Todos las otras 4 caídas de voltaje V12, V13, ... son incógnitas que tenemos que averiguar. También son desconocidas las 5 corrientes: la corriente I entregada por la fuente y  la corriente circulando por los 4 resistores I12, I13, .... Planteemos las ecuaciones de este circuito y resolvamos con Octave.

Ley Kirchhof del voltaje

La ley de Kirchof del voltaje (que es una ley que se puede derivar de las ecuaciones del electromagnetismo), dice que en ausencia de campos magnéticos, la suma de las caídas de voltaje en un circuito cerrado (en un lazo, en un loop) debe ser 0.
En la figura de arriba, a la derecha, he marcado dos lazos con flechas rojas. La dirección de las flechas es arbitraria y simplemente indica en que dirección el voltaje cae. Lo único importante es que una vez que hemos dado un sentido a las flechas no lo podemos cambiar (si lo cambiamos mientras escribimos las ecuaciones, estas van a estar mal). Es decir que si V12 es positivo, quiere decir que el nodo 1 tiene un voltaje más alto que el nodo 2. Si este voltaje es negativo, entonces es al revés, el nodo 2 tiene un voltaje mas alto que el nodo 1.
En el dibujo de la izquierda, marque otro lazo que no involucra a la fuente.

Escribamos las ecuaciones:

Las dos primeras corresponden a los lazos en la figura de la derecha. El primer lazo empieza en el nodo 1, va al 2 luego al 4 y vuelve al 1.  Cuando circulo en la misma dirección que la flecha el voltaje aparece en la ecuación con signo positivo, y cuando circulo en contra del sentido de la flecha, el voltaje aparece con signo negativo (¡voy en contra!). En el caso del segundo lazo (la segunda ecuación) me paseo por los nodos 4,3,1 y de regreso  4. al tercera ecuación es dar una vuelta pasando por todos los resistores. 
¿Se dan cuenta que la suma de las dos primeras ecuaciones genera la tercera? Esto quiere decir que la tercera ecuación no nos brinda nueva información: sabiendo las dos primeras ecuaciones, ya sabemos lo que dice la tercera. Por lo tanto la tercera ecuación no la vamos a utilizar (podrían no utilizar cualquier otra).
Aplicando nuestra habilidad para transformar sistemas de ecuaciones lineales en productos de matrices podemos obtener:

En este sistema parecen las 4 incógnitas de voltaje. En Octave podemos escribir la matriz de factores de la ley de Kirchhof de voltajes (LKV), de la siguiente manera:

LKV = [1 0 1 0; 0 1 0 1];

Ya tenemos las leyes de Kirchhof de voltaje, ahora veamos las corrientes.

Ley Kirchhof de las corrientes

La idea detrás de esta ley es que la corriente se conserva. Es decir que cuando corrientes se unen o se separan en un nodo, la suma se tiene que conservar. Para escribir estas ecuaciones nos paramos en un nodo y vemos que corrientes llegan al nodo (positivas) y que corrientes se van del nodo (negativas). En este caso también obedecemos el sentido de las flechas que definimos al principio. Repito, es importante mantener el sentido de las flechas hasta el final, pero el sentido es arbitrario. Una vez que elijen un sentido para la flecha, no lo cambien hasta que terminen.
Hagamos un ejemplo: la corriente que llega al nodo 1 viene directo de la fuente (corriente I); la corriente I13 también llega al nodo, mientras que la corriente I12 se va del nodo. Todo esto nos da la primera ecuación del siguiente sistema:

Las demás ecuaciones las pueden escribir siguiendo el mismo método para cada nodo. La ultima ecuación, nodo 4, también se puede deducir de las tres primeras (¿Pueden explicar cómo?), por lo tanto no la vamos a utilizar. El sistema escrito en forma matricial es:

Noten que tenemos 3 ecuaciones y 5 incógnitas de corriente.

En Octave podemos cargar la matriz de factores para la ley de Kirchhof de las corrientes (LKC) con el código:

LKC = [-1 1 0 0 1; 1 0 -1 0 0; 0 -1 0 1 0];

Hasta aquí hemos coleccionado 5 ecuaciones: 2 de LKV y 3 de LKC. Sin embargo tenemos 9 incógnitas, 4 voltajes y 5 corrientes. Para poder encontrar tantas ecuaciones como incógnitas necesitamos recurrir a la relación entre voltaje y corriente determinada por los resistores.

Relación entre corriente y voltaje

En el k-ésimo resistor, la relación entre la caída de voltaje en el resistor y la corriente fluyendo a través de el, viene dada por:

Como ven es también una relación lineal: la caída de voltaje es proporcional a la corriente. Escribiendo esta ecuación para cada resistor, obtenemos el siguiente sistema de ecuaciones:

Para construir esta matriz en Octave, es necesario dar el valor de resistencia de los resistores. Aquí elijo cualquier valor, uds. pueden usar otros:

R = [1 1 3 3];
VIR = [eye(4,4) diag(-R) zeros(4,1)];

En la última linea de comando utilizo 3 funciones que no hemos discutido aún. Para entender que hacen tenemos que mirar un poco la forma de las matrices en el sistema de ecuaciones. Si se fijan en la primera parte de la matriz podrán ver que hay unos (1) formando una diagonal. En la segunda parte también tenemos una diagonal formada por el negativo de las resistencias. La última columna de la matriz está llena de ceros. Ok, a las funciones:

eye(n,n): Esta función genera una matriz con n filas y n columnas que tiene la diagonal llena con unos.

Pregunta 1:
¿Qué matriz se genera si llamamos a esta función con dos argumentos distintos? Por ejemplo eye(2,3).

diag(v): Esta función genera una matriz cuya diagonal se llena con los elementos del vector v.

Pregunta 2:
¿Podés generar la matriz eye(4) con la función diag()?

zeros(n,m): Esta función genera una matriz con n filas y m columnas llenas de ceros.

Ensamblando todo el sistema

Quizás ya se han dado cuenta que ahora tenemos tantas ecuaciones como incógnitas. Ahora podemos armar el sistema de ecuaciones completo

Con color resalté cada uno de los bloques de la matriz y su procedencia. Los bloques marcados con un cero grande, indican que allí solo hay ceros.
En Octave podemos ensamblar esta matriz y resolver el sistema (¡Ojoi! antes de resolver necesitan darle un valor al voltaje Vde la fuente) con el siguiente código:

A   = [LKV zeros(2,5); zeros(3,4) LKC; VIR];
neq = size (A,1); # Numero de ecuaciones
b   = [-Vi; Vi; zeros(neq-2,1)];

# Resolver el sistema
x  = A \ b;
V  = x(1:4);
I  = x(5:8);
Ii = x(9);

disp ("Voltages"); disp (V); disp ("Corrientes"); disp (I); disp ("Corriente de la fuente"); disp (Ii)

Noten como he vuelto a utilizar la función zeros() para rellenar los bloques con ceros. La última linea de este código muestra los resultados en pantalla ¿Coinciden con tus cálculos manuales?

El archivo adjunto contiene todo el código aquí utilizado.

¿Encontraste un error en este post? ¿Tienes preguntas? Utiliza el foro de discusión.

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

Tic,tac, tic, toc! ¿Cuánto tarda mi programa?

Matías Fajardo pidió que discutamos control de flujo de programas en GNU Octave. Es un pedido razonable así que vamos a ir en esa dirección. 
Con el tema de control de flujo tocamos algo importante para lenguajes como GNU Octave. Para entender la importancia del tema necesitamos primero poder medir el tiempo que tarde en correr un programa. Por esta razón vamos a empezar nuestro camino hacia control de flujo con unas funciones muy sencillas: tic y toc.


El objetivo es el siguiente: queremos medir cuanto tarda la ejecución de un programa. Por ejemplo ¿Cuánto tarda la creación de una matriz en Octave? Es decir ¿Cuánto tarda la ejecución del comando
rand (n)
para distintos valores de n?

Una forma de medir el tiempo de ejecución es usar las funciones tic y toc:
n = 100;
rand (n);
La función tic inicia un cronómetro y la función toc muestra le valor actual del cronómetro. 
En mi computadora los resultados para n igual a 1e2, 1e3, y 1e4 son
Elapsed time is 0.000303984 seconds.
Elapsed time is 0.024765 seconds.
Elapsed time is 1.57575 seconds.

Pregunta 1:
¿Qué tiempos obtenés vos?

La respuesta de toc nos dice el tiempo segundos que transcurrió desde la última vez que llamamos a la función tic. Este tiempo se llama tiempo de reloj de pared (wall-clock time, en Inglés) y corresponde más o menos al tiempo real que transcurrió en tu vida mientras se ejecutaba el código.

Pregunta 2:
¿Qué resultados obtenés si corres el siguiente código?
n = 5e2;
rand (n);

n = 5e2;
clear all
rand (n);
¿Notás algo raro?

¡Preguntas y comentarios en el foro!

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

Si, si en cambio, sino: IF-ELSEIF-ELSE

En la clase sobre lógica discutimos un poco la forma de generar valores verdadero/false, usualmente llamados valores booleanos (se lee bu-leanos y viene del nombre de un inglés llamado George Boole).

Casi todos los lenguajes de programación ofrecen una forma de controlar la ejecución de código utilizando como criterio valores booleanos. Es decir, cierta parte del código se ejecuta si una condición lógica es verdadera. Una forma ejemplar de realizar esta ejecución selectiva es utilizando IF-ELSE (si, sino) o IF-ELSEIF (si, si en cambio) o una combinación de ambos.


En Octave, todos aquellos comandos dentro de un bloque if se ejecutan solo si a condición dada es verdadera. 

if (x == 0)
  disp ("Equis es igual a cero");

El bloque if se empieza en la linea donde llamamos a esta función y termina con endif (o simplemente end, es lo mismo). En este caso, se imprime a la pantalla la oración "Equis es igual a cero" cuando la comparación x == 0 es verdadera. La condición puede ser tan compleja como sea necesario (aunque esto puede hacer el código difícil de entender), por ejemplo

if (x == 0 || x == 5)
  disp ("Equis es igual a cero o igual a cinco");

Si la condición fuese falsa, el código dentro del bloque nunca se ejecuta. En cualquier caso, una vez evaluado el bloque (con o sin ejecución del código correspondiente), el programa continua. El siguiente ejemplo imprime "Equis es igual a cinco" dependiendo del valor de x, pero siempre imprime "Terminé".

if (x == 5)
  disp ("Equis es igual a cinco");
disp ("Terminé.");

Pregunta 1:
¿En que situación se imprime a pantalla el mensaje "¡Muy loco!"?

if (x == 0 && x == 5)
  disp ("¡Muy loco!");

¿Veremos alguna vez el mensaje en pantalla?


El bloque if se puede extender para que cierto código se ejecute solo cuando la condición es falsa. Considera el siguiente código

disp ("Eqis es: ")

if (x > 0)
    disp ("positivo.");

if (x <= 0)
    disp ("negativo o cero.");

Dado que la segunda condición es lo opuesto a la primera podemos simplificar usando else

disp ("Eqis es: ")
if (x > 0)
    disp ("positivo.");
    disp ("negativo o cero.");

Pregunta 2:
¿Cómo harías para separar en "positivo", "cero" o "negativo" el ejemplo anterior?


Una respuesta a la pregunta anterior viene de la mano de elseif (si en cambio), que nos permite evaluar una condición adicional.

disp ("Eqis es: ")
if (x > 0)
    disp ("positivo.");
elseif (x == 0)
    disp ("cero.");
    disp ("negativo.");

Pregunta 3:
¿Podés dar otros ejemplos que produzcan lo mismo que el ejemplo anterior?

Debe notarse que el código dentro de un elseif se ejecutara solo y solo si la condición del if es falsa. En el siguiente ejemplo jamás veremos el mensaje "¡Ahá!" en pantalla aunque la condición del elseif es siempre verdadera.

if (true)
 disp ("Sí, es verdad.");
elseif (true)
 disp ("¡Ahá!");

Esto quiere decir que las condiciones en un bloque IF-ELSEIF deberían que ser mutuamente exclusivas (nunca ambas pueden ser verdaderas al mismo tiempo). Esto no es necesario, como se muestra en el siguiente ejemplo, pero la lógica resultante es confusa

if (x > 0 && x <=1)
    disp ("Entre cero y uno (incluido).");
elseif (x > 0.5)
    disp ("Más que un medio.");

Pregunta 4:
¿Podés simplificar la lógica del ejemplo anterior?

Espero recibir sus preguntas y dudas en el foro!

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

¿Verdadero o falso?

En el cole te enseñaron una cosa que llaman lógica formal (o solamente lógica). Si te gustó tenés suerte, probablemente tu profesor/a era muy buen@ y entendía el tema (¡yo tuve mucha suerte!). Lamentablemente la mayoría de los estudiantes con los que charlo odian la lógica. En esta clase voy a tratar de darles un ejemplo de porqué la lógica es super importante y que no es tan aburrida como muchos creen.

Esta clase puede ayudarte con los desafíos del nivel 1, si todavía no los resolviste.

"La lógica no sirve para nada"

¡Cuantas veces habrás escuchado esa frase o alguna otra similar! Como la mayoría de las herramientas matemáticas, la lógica no es más que una formalización de algo que hacemos naturalmente. Voy a tomar un ejemplo mundano: seleccionar fruta o verdura cuando vamos de compras.

Hace unos días, en uno de mis viajes de aprovisionamiento al mercadito de la esquina, observaba con cierta admiración a una mujer mayor seleccionar las manzanas que ponía en su bolsa de las compras. Su meticulosidad era sorprendente y me preguntaba si esa mujer no habría sido una matemática en sus años de actividad laboral.
Cada vez que la señora juntaba una manzana, la miraba de muchos ángulos diferentes, la golpeaba suavemente con un dedo, la olía y finalmente decidía si la manzana volvía a la góndola o si se convertía en parte de su compra.
Podemos hacer un modelito sencillo del proceso de selección y pensar que la señora ejecutaba una función que tomaba como entrada los aspectos de la manzana bajo observación. Esta función evalúa tres condiciones y decide si la manzana se compra o se devuelve:
  1. Tiene buena forma?
  2. Suena bien?
  3. Huele bien?
La respuesta a estas preguntas es si o no. Podemos re-escribir las preguntas como afirmaciones (o condiciones) y decidir si estas afirmaciones son verdaderas (en inglés: true) o falsa (en inglés: false):
  1. La manzana tiene buena forma.
  2. La manzana suena bien al ser golpeada suavemente.
  3. La manzana huele bien.
Si la manzana en cuestión hace que todas estas afirmaciones sean verdaderas, entonces la compramos. De lo contrario, si cualquiera de estas afirmaciones es falsa, devolvemos la manzana a la góndola.

Símbolos y operadores

En el ejemplo introduje los dos símbolos fundamentales de la lógica: true y false. Estos son los nombres que se utilizan en GNU Octave para estos valores, pero vos podés usar cualquier par de valores que te guste: (V,F); (0,1); (-1,1); (0,5); (blanco, negro); etc.. Para poder ejercer la lógica necesitamos saber como operar con estos valores.

Tomemos dos afirmaciones, una falsa y otra verdadera, para ejemplificar. 
  1. La afirmación verdadera es: la palabra "hacha" tiene 5 letras
  2. La afirmación falsa:  La palabra "hacha" no tiene ninguna letra "h".
Verifiquemos los valores lógicos de estas afirmaciones en GNU Octave. El siguiente código guarda el string "hacha" en una variable llamada p (de palabra) y evalúa la 1ra afirmación.

p = "hacha";
# Cuenta la cantidad de letras
nLetras = length (p);
# Evalua la 1ra afirmacion 
nLetras == 5

La primera linea de código guarda la palabra en la variable. La linea luego del primer comentario devuelve la longitud del string p, es decir la cantidad de letras, que guardamos en a variable nLetras (recuerda que las lineas que empiezan con el caracter # o %, son ignoradas por GNU Octave, son lineas de comentarios, y tu también puedes ignorarlas. Son solo ayuda). La última linea de código pregunta si nLetras es igual a 5. Notar que el operador == no asigna el valor a nLetras, sino que compara los valores y responde verdadero si estos son iguales (experimenta un poco con este operador!).
Si ejecutas esas líneas de código en una sesión de Octave obtendrás la respuesta 

ans = 1

En Octave el valor true se representa con un 1 (ejecuta true en tu sesión de Octave para verificar) por lo tanto el valor de la 1ra afirmación es verdadero. Veamos el código de la 2da afirmación:

all (p != "h")

¿Simple, eh? Dejo como ejercicio entender que esta ocurriendo aquí. Lo que necesitas saber es que all en inglés significa "todo/todos" y que el operador != es lo opuesto a == (verifica desigualdad en vez de igualdad) ¿Qué ocurre cuando comparamos el string p con el caracter "h"?

La respuesta de Octave a esta linea de código es 0. Indicando que, efectivamente, la segunda afirmación es falsa.

Ok, tenemos una afirmación verdadera y una falsa, y vamos a guardar estos valores de verdad en un vector:

afirmacion = [true, false]

De esta forma afirmacion(1) es verdadera (true, 1) y afirmacion(2) es falsa (false, 0).

¿Qué dirías de la nueva afirmación más abajo?
  • La palabra "hacha" tiene 5 letras y la palabra "hacha" no tiene ninguna letra "h"
¿Es verdadera o falsa?

El operador que hemos usado para construir esta nueva afirmación es el "y", en inglés "and", y que en lógica llamán conjunción. En Octave podemos utilizar la función and o el operador &, veamos: 

# operador
afirmacion(1) & afirmacion(2)
# funcion
and (afirmacion(1), afirmacion(2))

Ejecuta este código ¿Estas de acuerdo con Octave sobre el valor de verdad de la nueva afirmación?
¿Cuál sería el resultado si las dos afirmaciones fuesen verdaderas?

Veamos otra afirmación
  • La palabra "hacha" tiene 5 letras o la palabra "hacha" no tiene ninguna letra "h"
Para construir esta nueva afirmación utilizamos la disyunción inclusiva, el "o". Pero no el "o" exclusivo (que algunos escriben "ó" para que quede claro), sino un "o" que está contento si esto es verdad, si lo otro, o si ambos son verdad. Este operador en Octave se escribe | y la función es or:

# operador
afirmacion(1) | afirmacion(2)
# funcion
or (afirmacion(1), afirmacion(2))

Experimenta con este operador para entender como funciona. 

Los operadores & y |, al igual que las respectivas funciones toman como entrada dos afirmaciones, son operadores binarios (en el sentido que toman 2 entradas). Otros operadores binarios que ya conocés son el + y el *, ambos toman dos números como entrada y devuelven el resultado.

Ejercicio 1: ¿Puedes construir el resultado de los operadores & y | para todas las entradas posibles?

Ejercicio 2: ¿Puedes escribir la afirmación que la señora de la historia evaluaba para seleccionar las manzanas?

Ejemplo numérico

Ok, toda esta cháchara sobre lógica, pero ¿Cómo nos puede ser útil para los desafíos?

En los desafíos de nivel 1 es necesario filtrar o seleccionar ciertos números según sus propiedades. Por ejemplo ¿Cómo podemos seleccionar números pares?
Lo que sabemos es que un número par se escribe como

K = 2*N

donde N puede ser cualquier numero entero. Todo bien con esto, pero para cada número K que me dan, tendría que buscar un N, entre todos los enteros, que cumpla con esta condición...posible, pero medio pesado. De manera equivalente podemos decir que un número es par si al dividirlo en 2 el resto es cero. Es decir, un número par es divisible por 2. GNU Octave tiene una función que nos devuelve el resto de una división entera, la función rem (del inglés "remainder", que quiere decir "resto" o "sobrante") . Como ejemplo de uso tomemos la division de 7 en 2 partes. 7 se puede repartir entre 2 en 3 partes iguales, pero nos sobra 1. Veamos que dice a la función:

rem (7, 2)
ans = 1

El primer argumento es el numero a dividir (el dividendo) y el segundo argumento es el divisor.

Ahora vamos a utilizar esta función y nuestro conocimiento de lógica para seleccionar los números pares de una lista de números:

l = 0:3:15;
tf = rem (l,2) == 0
tf =
   1   0   1   0   1   0

La lista l contiene números enteros desde 0 hasta 15 en incrementos de a 3 y Octave nos dice que la afirmación "el número es divisble por 2" es verdadera para el primer número, falsa para el segundo, verdadera para el tercero y así sucesivamente. 
¿Cuales son estos números? Claramente 0, 6 y 12. En Octave podemos obtener estos números utilizando el vector lógico al que hemos llamado tf. Prueba lo siguiente


Sigue experimentando con este tipo de ejercicios hasta que entiendas como funciona la cuestión. Si tienes preguntas o sugerencias no dudes en escribirlas en el foro de discusión o los comentarios.

Adjunto un archivo con todos los comandos que hemos utilizado en esta clase.

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

Producto escalar

En la clase sobre operadores vimos que el operador asterisco * parece no funcionar con vectores. Hoy veremos que este operador realiza el famoso producto escalar entre vectores o el producto entre matrices.

Longitud de un vector (módulo o norma)

Para indicar la posición de un punto en el plano (digamos una bolita en una cancha) necesitamos dar al menos dos magnitudes. Por ejemplo, tomando un punto como referencia (una esquina de la cancha), podemos indicar las distancias al punto desde ese origen en dirección horizontal y en dirección vertical. Estos dos valores los podemos acomodar en un vector [x y].

¿Cuál es la distancia en linea recta desde el origen hasta la bolita? La respuesta viene de la mano del teorema de Pitágoras

Este valor es la longitud del vector (t.c.c. módulo del vector o norma del vector). Dentro de la raíz cuadrada tenemos el cuadrado del módulo del vector, la suma del producto de cada componente por si misma. 
Con toda libertad podemos preguntarnos como se llama la misma operación entre dos vectores diferentes. Digamos que tenemos nuestro vector [x y] y un otro vector [u w]. Haciendo las cuentas obtenemos

Este valor es claramente un escalar (un número, no un vector), como lo es la distancia. A esta operación se la denomina producto escalar entre vectores. El producto escalar es una función que toma dos vectores como entrada y devuelve un valor escalar, de allí su nombre.

¿Matriz o vector?

Los matemáticos son tipos muy ocurrentes y se inventaron una forma de escribir el producto escalar de manera muy concisa. Primero decidieron diferenciar entre vectores fila y vectores columna (como ya vimos en otra clase) y luego decidieron que cuando uno multiplica un vector fila por un vector columna, se debe calcular el producto escalar. Es decir, multiplicamos el primer elemento del vector fila por el primer elemento del vector columna y lo sumamos al producto de los segundos elementos.

aimación producto escalar
Cambiar de vector filas a columnas en Octave es muy sencillo.

x  = [1 2];
xT = x'

xT =

y la operación se llama
transponer. Es decir la comillas simples ' es un operador unitario (toma una sola entrada) que devuelve el vector transpuesto: cambia filas por columnas o viceversa. Ahora podemos calcular el módulo cuadrado y la longitud del vector x definido más arriba:

mod2 = x * x';
long = sqrt (mod2)

Pregunta 1:
¿Qué hace la función sqrt()? ¿Por qué se llama así?

Ya nos vamos dando cuenta que no hay mucha diferencia entre un vector y una matriz. En efecto, el producto entre matrices es el producto escalar de los vectores fila de la primera matriz, por los vectores columna de la segunda.

A = [1 2; 1 -1];
# vectores filas
a1 = A(1,:); # primera fila
a2 = A(2,:); # segunda fila

B = [1 0; 0 1];
# vectores columnas
b1 = B(:,1); # primera columna
b2 = B(:,2); # segunda columna

C  = A * B;
Cp = [a1*b1 a1*b2; a2*b1 a2*b2];

assert (Cp == C)

La matriz Cp se construye usando el producto escalar de las filas y las columnas. La función assert produce un error si la condición que se pasa por argumento no es verdadera. Al no recibir ningún error, entendemos que las dos matrices son iguales.

¡Un vector es tan solo una matriz con una sola fila o columna!

¿Preguntas? Participa en el foro.

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

Crea tus propias funciones

A esta altura están en posición de escribir sus propias funciones. Pero ¿Qué es una función?

Una función es una entidad que recibe una serie de argumentos como entradas, opera sobre ellos y luego devuelve ciertos resultados. Las funciones cumplen un rol fundamental: encapsulan tareas. Al encapsular tareas, en especial aquellas que se repiten seguido, nuestros programas son mas fáciles de leer y a veces hasta más eficientes.

Considera la situación de saludar a personas que encuentras durante tu día. Supongamos que tu frase elegida como saludo es "Hola, <nombre>",  donde <nombre> es el nombre de pila de la persona a al que saludamos. Dado que a esto lo vamos a hacer muy seguido, podríamos definir una función (o rutina) que haga lo siguiente:

1. Toma el nombre de la persona a saludar.
2. Adjunta "Hola, ".
3. Imprima el saludo
4. Devuelva un resultado con el saludo.

Esta función en octave sería

function frase = saludar (nombre)
  frase = ["Hola, " nombre];
  disp (frase);

En una sesión de octave utilizaríamos esta función de la siguiente manera

x = "María";
y = saludar (x);

En la pantalla veríamos "Hola, María" y el contenido de la variable y sería exactamente esa frase.

Pregunta 1
¿Puedes definir una función que tome dos entradas y devuelva la suma de las mismas?

Estructura de una función

El siguiente esquema muestra la estructura de una función

function <resultados> = <nombre de función> ( <entradas> )

Las palabras y símbolos en negrita son necesarios y obligatorios.

Las palabra clave function es necesaria (y obligatoria) para definir una función. Un archivo que contiene una función debe siempre tener como primer comando ejecutable esta palabra clave. El archivo donde guardamos la función siempre tiene que llamarse  <nombre de funcion>.m , es decir el nombre del archivo debe coincidir con le nombre de a función que allí se define.
Las variables que se pasarán al ámbito desde donde se ha llamado la función (p. ej. la sesión de Octave) se devuelven como un vector en <resultados>. En el caso de la función saludar.m de más arriba, devolvíamos solo una variable. Pero una función puede devolver un montón de cosas o nada. Ejemplos

function [x, y] = ceroyuno ()
 x = 0;
 y = 1;

function diceHola ()
  disp ("Hola");

De manera similar, las variables de entrada  se pasan en una lista en <entradas>. Las entradas pueden ser muchas o ninguna (como en los ejemplos anteriores). Ejemplo

function z = resta (x,y)
  z = x - y;

Pregunta 2:
¿Cual es el nombre del archivo .m donde debemos guardar la función resta definida en el ejemplo anterior?

Pregunta 3:
¿Como se define una función que devuelva la cantidad de letras "a" en un string?


Cuando discutimos como guardar nuestro trabajo, vimos que desde Octave podemos ejecutar el contenido de un archivo de texto. Los comandos en ese archivo tienen acceso a las variables que existen en el ámbito (en inglés scope) de la sesión de Octave; esto quiere decir que los scripts pueden "ver" y "editar" las variables que creamos directamente en la sesión, y que si estos scripts crean nuevas variables estas quedaran en nuestra sesión luego de que el script haya terminado de ejecutarse.

Una función, en cambio, crea su propio ámbito de variables. Las variables que son creadas dentro de una función no son accesibles desde la sesión de Octave, o por lo menos no directamente. Excepto en casos particulares, las variables creadas dentro de una función son eliminadas cuando la función termina.

Con estas ayuditas ya puedes empezar crear tus propias funciones. El mundo de las funciones es mucho más amplio de lo que he descrito aquí ¡No te olvides de explorar!. Espero tus dudas y preguntas en el foro.

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM