Planet Octave

November 28, 2015

Jordi Gutiérrez Hermoso

Octave code sprint 2015

So, let’s get this going!

#yop-poll-container-3_yp5659ba6ed7e75 { width:200px; background:#fff; padding:10px; color:#555; overflow:hidden; font-size:12px; } #yop-poll-name-3_yp5659ba6ed7e75 { font-size:14px; font-weight:bold; } #yop-poll-question-3_yp5659ba6ed7e75 { font-size:14px; margin:5px 0px; } #yop-poll-answers-3_yp5659ba6ed7e75 { } #yop-poll-answers-3_yp5659ba6ed7e75 ul { list-style: none outside none; margin: 0; padding: 0; } #yop-poll-answers-3_yp5659ba6ed7e75 ul li { font-style:normal; margin:0px 0px 10px 0px; padding:0px; font-size:12px; } #yop-poll-answers-3_yp5659ba6ed7e75 ul li input { margin:0px; float:none; } #yop-poll-answers-3_yp5659ba6ed7e75 ul li label { margin:0px; font-style:normal; font-weight:normal; font-size:12px; float:none; } .yop-poll-results-3_yp5659ba6ed7e75 { font-size: 12px; font-style: italic; font-weight: normal; margin-left: 15px; } #yop-poll-custom-3_yp5659ba6ed7e75 { } #yop-poll-custom-3_yp5659ba6ed7e75 ul { list-style: none outside none; margin: 0; padding: 0; } #yop-poll-custom-3_yp5659ba6ed7e75 ul li { padding:0px; margin:0px; font-size:14px; } #yop-poll-container-3_yp5659ba6ed7e75 input[type='text'] { margin:0px 0px 5px 0px; padding:2%; width:96%; text-indent:2%; font-size:12px; } #yop-poll-captcha-input-div-3_yp5659ba6ed7e75 { margin-top:5px; } #yop-poll-captcha-helpers-div-3_yp5659ba6ed7e75 { width:30px; float:left; margin-left:5px; height:0px; } #yop-poll-captcha-helpers-div-3_yp5659ba6ed7e75 img { margin-bottom:2px; } #yop-poll-captcha-image-div-3_yp5659ba6ed7e75 { margin-bottom:5px; } #yop_poll_captcha_image_3_yp5659ba6ed7e75 { float:left; } .yop_poll_clear { clear:both; } #yop-poll-vote-3_yp5659ba6ed7e75 { } .yop-poll-results-bar-3_yp5659ba6ed7e75 { background:#f5f5f5; height:10px; } .yop-poll-results-bar-3_yp5659ba6ed7e75 div { background:#555; height:10px; } #yop-poll-vote-3_yp5659ba6ed7e75 div#yop-poll-vote-3_yp5659ba6ed7e75 button { float:left; } #yop-poll-vote-3_yp5659ba6ed7e75 div#yop-poll-results-3_yp5659ba6ed7e75 { float: right; margin-bottom: 20px; margin-top: -20px; width: auto; } #yop-poll-vote-3_yp5659ba6ed7e75 div#yop-poll-results-3_yp5659ba6ed7e75 a { color:#555; text-decoration:underline; font-size:12px;} #yop-poll-vote-3_yp5659ba6ed7e75 div#yop-poll-back-3_yp5659ba6ed7e75 a { color:#555; text-decoration:underline; font-size:12px;} #yop-poll-vote-3_yp5659ba6ed7e75 div { float:left; width:100%; } #yop-poll-container-error-3_yp5659ba6ed7e75 { font-size:12px; font-style:italic; color:red; text-transform:lowercase; } #yop-poll-container-success-3_yp5659ba6ed7e75 { font-size:12px; font-style:italic; color:green; } .yop-poll-forms-display{}
When should the Octave 2015 code sprint be?

by Jordi at November 28, 2015 02:23 PM

August 31, 2015

Juan Pablo Carbajal

Archivando el trabajo

En esta clase veremos como crear y ejecutar un archivo m (mfile).

Todos los comandos que ejecutamos en GNU Octave pueden guardarse en un archivo de texto. Es importante que el archivo sea de texto plano (ascii), puesto que este es el formato que GNU Octave espera en los archivos de comandos. Esto es verdad en general para la mayoría de los lenguajes de programación.

Para empezar abrimos un archivo de texto con nuestro editor preferido (en mi caso, que trabajo en linux, uso Gedit). Enfatizo que no deben usar procesadores de texto como Office porque estos programas no guardan sus archivos en texto plano. Si estas en Windows puedes usar Notepad++.

El siguiente video muestra como trabajo en Ubuntu:

Crear y ejecutar mfile en octave-cli

El video muestra la creación de un mfile que corresponde a una función. Todavía no hemos aprendido a escribir funciones así que por el momento vamos a ver otro tipo de mfile: el archivo de comandos (script).

Abrí un archivo de texto y escribí los siguientes comandos:

disp ("Hola, este es un archivo de comandos! (script)");
a = 1:5;
txt = repmat (", %d",1,length(a)-1);
printf (["Los primeros %d números son: %d" txt "\n"],length(a), a);

Guarda este archivo notando el directorio donde quedará guardado (en el video utilicé el escritorio). Voy a suponer que el nombre que le diste al archivos es script_00.m
Ahora inicia una sesión de octave. Asegúrate que el archivo esta en el directorio actual: ejecuta pwd para ver el directorio actual y ls parta ver la lista de archivos.
Para ejecutar el script, simplemente lo llamamos por su nombre (sin extensión): 


¡Listo! Has creado tu primer script en Octave. La siguiente lista enumera los pasos para trabajar con scripts:
  • Abrir un archivo ASCII (texto plano) donde escribir los comandos.
  • Guardar el archivo con extensión .m (fijarse donde se guarda).
  • Llamar al archivo desde GNU Octave utilizando el nombre sin la extensión (asegurarse que pwd corresponde con el lugar done se guardó el archivo).
Es importante utilizar scripts (es la forma recomendada de trabajar), puesto que los comandos quedan archivados y podemos corregir errores y volver a ejecutarlos sin tener que escribir todo de nuevo.

Si estas usando una versión de Octave 3.8 o más nueva puedes utilizar la interfaz gráfica que incluye un editor de textos. El siguiente video muestra como crear un script que imprime los números del 1 al 10.

Crear y ejecutar mfile en octave-gui

La GUI de Octave puede configurarse en español si uds. lo desean. noten que esto solo cambia el idioma del menú, no el idioma del lenguaje de programación (funciones, ayuda y mensajes).

No duden en escribir sus preguntas y comentarios, aquí abajo o en el foro de discusión.

by Juan Pablo Carbajal ( at August 31, 2015 05:48 PM

August 27, 2015

Asma Afzal

Wrapping up

So, the GSoC period has come to an end now.

Project Goals  
My project was about creating Matlab compatible wrappers for the optim package.  Here is a brief list of my project goals.

1- lsqnonlin wrapping nonlin_residmin and residmin_stat
2- lsqcurvefit wrapping nonlin_curvefit and curvefit_stat
3- nlinfit wrapping nonlin_curvefit (it was initially decided to wrap leasqr but changed to avoid extra computations)
4- quadprog wrapping __qp__ instead of qp and returning lambda in the form of a structure as in Matlab
5- fmincon wrapping nonlin_min
6- Test and demos for the above functions
7- Stretch goals: previously decided to create other missing functions or perhaps additional backends but before midterm I decided instead to include optimoptions in my to do list.

The functions lsqnonlinlsqcurvefit and nlinfit are complete with tests and demos and integrated in the optim package. Since nlinfit is from the statistics package in Matlab, additional functions such as statset, statget were required for handling options. These functions are implemented with minor modifications in optimset, optimget and __all_opts__ as statset, statget and __all_stat_opts__ and are now a part of optim package.

The function quadprog required directly wrapping __qp__ instead of qp for the ordering of lambda. It is in the final stages of review and will soon be integrated.

fmincon has not been thoroughly reviewed yet. I will send it to Olaf after quadprog is committed to optim. 

Hiccups in the stretch goal
I couldn't create optimoptions in the GSoC time frame because it was a bit open ended and I had to come up with an object oriented design for the function. I was trying to understand how Matlab implements it for quite some time. Anyway, I didn't pursue it further and shifted my focus on the refinement of my almost complete functions to get them integrated in optim.  

Interesting Takeaways
This is my first experience of working with any open source organization and it's definitely a pleasant one. It's delightful to see people using my functions and possibly benefiting from my work [2-3]. :)

I think I have managed to meet all the goals I set before the start of GSoC. (Regrets? Well, I could have saved more time for optimoptions and it would've been better to discuss it way before than being stuck for a while.)

I'm extremely grateful to the Octave community, especially my mentors Olaf Till and Nir Krakauer for their unrelenting support. GSoC wouldn't have been possible without their constructive feedback. I have learned a lot from this experience.


by Asma Afzal ( at August 27, 2015 01:52 PM

August 20, 2015

Piotr Held

Summary of project

I have not written a post in a long while because I have had family issues. Although my progress was hampered, it was not completely halted. I was able to add some more features and clean some of the old code up. I have added two new sections with five new functions:
  • Testing for nonlinearity:
    • surrogates
    • endtoend
    • timerev
  • Spike trains:
    • spikeauto
    • spikespec
They are complete with documentation and some demos. The wiki has also been updated with a demo of surrogates. Thus the original plan has been completed apart from some of the functions in the final section Tutorial and randomize. There have been some additional functions: endtoend, spikeauto, spikespec and av_d2. These were added as a logical addition to the existing function.

What happened to randomize?

The first time I read the documentation and tried to plan the project I did not realize how big randomize is. Therefore, once I realized that it was a toolbox and not a function I needed to design how this toolbox should work with/in Octave. This took about a week of my time. I have not had any experience designing anything similar to randomize. I considered it an interesting challenge. I tried to create a hybrid that could use C++ and Octave functions to run simulated annealing. I wanted the main loop to be in C++ because this algorithm cannot be vectorized, so it will run slowly on Octave. 

My plan was to create some "runner" that can call some abstract methods of the base cost function class. I was able to accomplish this only partly, even though I introduced some polymorphism I only created one cost function (I did not finish the second one). This one cost function and the runner worked properly (as tested against the results from TISEAN). I was even able to do a type of pausing (like in other functions in the package), which allows the user to break the execution without having to kill Octave. All of this worked correctly with the one cost function. My plan was to modify the code to ensure the runner could manage any type of cost functions (that inherit from the base cost function class) once the second cost function was completed

Afterwards, I planned to create a base abstract class in Octave using classdef and similar keywords. Then I would have created an instance of the C++ class that would run the Octave methods using polymorphism. This was not possible however because classdef (and similar keywords) are yet not fully supported by Octave. They are not parsed properly by help, are not documented and not all of the Matlab functionality is present as of now.

Looking back at this attempt of porting the randomize toolbox I believe that with a good plan it could have taken more than half of my programming time this summer to create a complete toolbox with good tests and documentation. I do not regret trying to port it, but I think I can say that this part was not completed, even though much work was put into it.

All of the code created is placed in 'devel' folder. There are also some tests there that utilize the cgreen (cgreen++ actually) framework. It is almost completely undocumented, as I spent most of my time developing new code.

Other thoughts

I am very pleased with the support of the community. I expected they wouldn't have time to help, and at first I had trouble receiving the help I need, but once I understood what time (depending on who was online) to ask for help and which channel to use it became my biggest help in completing the project. My project needed help from people who knew things about how Octave was written in C++ and not just Octave developers. There were plenty of people who could assist me in just that.

One of the most exciting moments of the whole project was when I saw a problem on the mailing list that I encountered a few weeks earlier and could give some ideas on how to solve it.

If I were to do the project all over again and have to choose between the help of the community and the experience I have today after working on this project for over 3 months I would choose the community, because there are so many knowledgeable people there that know the inns and outs of Octave.

by Piotr Held ( at August 20, 2015 07:44 PM

August 19, 2015

Juan Pablo Carbajal

¿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 August 19, 2015 07:11 PM

August 16, 2015

Asma Afzal

Week 11 and 12: Integrating existing work in optim package.

A recap of the progress in two weeks:

  • I had to let go of optimoptions (for GSoC) mainly because of
    • time constraints 
    • and also because I don't have much experience with objected oriented programming. For optimoptions, I will have to come up with a design. I started with class implementation using classdef as in Matlab, but it is in its infancy in Octave and it could possibly be a limiting factor.
  • I am refining my existing functions and including tests and demos so they can be integrated in the optim package.
    • lsqnonlin and lsqcurvefit required additional options documentation and OutputFcn and Display setting. These two functions have now been successfully integrated in the optim package. [1]
    • Functions nlinfit and quadprog are under review.
    • I am working on fmincon now. Still have to discuss which backend should be used. lm_feasible can return Lagrange multipliers, gradient and hessian, but since it adheres to the constraints in all iterations, it behaves differently (from Matlab's algos) and sometimes less efficiently as octave_sqp, which only respects the constraints for the final result.

by Asma Afzal ( at August 16, 2015 08:28 AM

August 03, 2015

Asma Afzal

Week 10: Preliminary work on optimoptions

Thoroughly checking how optimoptions works in  Matlab.

options = optimoptions (SolverName)

Things to do:

  • Identify if Solver name is the right string or function handle

  • Cater for multiple Algorithms
     A subset of options for different algorithms.
  • Transfer relevant options of different solvers to modify/create option.
    oldoptions = optimoptions('fmincon','TolX',1e-10)
    newoptions = optimoptions('lsqnonlin',oldoptions)
  • Using dot notation or optimoptions to modify previously created options.
    (Second argument in optimoptions can be old options)
  • Display options:
    Set by user on top (for the current algo)
    Default options
    Options set for other algorithms.
Implementation ideas:

In Matlab, these two calls generate the same options object optim.options.Fmincon:

What would be more appropriate, IMO, will be to have a function optimoptions of the following format
opts = function optimoptions(SolverName,varargin)

    obj = optim.options.SolverName(varargin)

    opts = struct(obj)


This function will
  1. Instantiate the relevant class and request for default options from the solver.
  2. Compare the user provided options to add relevant options.
  3. Display options of the current algorithm.
  4. The output can also be returned in the form of a struct to be compatible with optimget.

by Asma Afzal ( at August 03, 2015 09:21 AM

July 24, 2015

Asma Afzal

Week 8 and 9: The ordering of lambda in quadprog

I was trying to dig through to figure out why the ordering of lambda in [1] does not match that of quadprog in Matlab.
The example shows how the values were different:

C = [0.9501    0.7620    0.6153    0.4057
    0.2311    0.4564    0.7919    0.9354
    0.6068    0.0185    0.9218    0.9169
    0.4859    0.8214    0.7382    0.4102
    0.8912    0.4447    0.1762    0.8936];
d = [0.0578
A =[0.2027    0.2721    0.7467    0.4659
    0.1987    0.1988    0.4450    0.4186
    0.6037    0.0152    0.9318    0.8462];
b =[0.5251
Aeq = [3 5 7 9];
beq = 4;
lb = -0.1*ones(4,1);
ub = 1*ones(4,1);
H = C ' * C;

f = -C ' * d;

[x, obj_qp, INFO, lambda] = qp ([],H,f,Aeq,beq,lb,ub,[],A,b);
lambda = 

Reordering lambda based on the length of constraints resulted in
lambda.eqlin = 
lambda.lower =

lambda.upper =

lambda.ineqlin =

Matlab however,  gave lambda.eqlin= -0.0165 for this example and lambda.lower =

There were two issues with the ordering:
  1.  the sign for Lagrange multipliers corresponding to linear equality constraint is always different from Matlab's
  2. The multipliers corresponding to the bound constraints as underlined above are swapped.
I tried several different examples to understand what is going on. For all the examples, the sign for lambda.eqlin was consistently different. Although, I still can't pinpoint why but for now I am just multiplying lambda.eqlin by -1.

For the swapping issue, I tried the same example with just the bound constraints:

[x, obj_qp, INFO, lambda] = qp ([],H,f,[],[],lb,ub)

but only considering lower bound constraints gave:
[x, obj_qp, INFO, lambda] = qp ([],H,f,[],[],lb,[])

which is how it is supposed to be. Tracing back, I found out that the ordering of lambda vector in qp.m [2] is not [equality constraints; lower bounds; upper bounds; other inequality constraints] like I previously assumed. From lines 287 and 288 in qp.m [2], the bound constraints are added to the inequality constraint matrix alternatively. So the issue wasn't swapping but understanding how the constraints are passed to __qp__.

In my code in [1], I had to make significant changes in the original qp.m code such as:
- Inequality constraint matrix has the order: [Linear inequality constraints; lower bounds; upper bounds]
- Check for too close bounds forming an equality constraint- This brings indexing issues as now the Lagrange multiplier value corresponding to bounds is in place of multipliers corresponding to linear equality constraints.
Also, Matlab only accepts too close bounds when using a medium scale algorithm and since the lower bound is approximately equal to the upper bound and it is considered as a single equality constraint, the single Lagrange multiplier is placed in the corresponding lambda.upper field while the corresponding lambda.lower value is zero.

Continuing the above example:

lb(4) = 0.3;
ub = 0.3*ones(4,1);

[x, obj_qp, INFO, lambda] = qp ([],H,f,[],[],lb,ub)
lambda =


Here, lb(4) = ub(4) and hence the constraint is treated as an equality constraint so the value for corresponding Lagrange multipliers is present on top (underlined)
I added checks for such cases and now my code in [1] gives same results as Matlab:


lambda =
  scalar structure containing the fields:
    lower =
    upper =


- qp.m strips off the -Inf constraints before passing them to __qp__. I am doing the same in quadprog. I have added further checks to make sure the multipliers are placed in the right positions in their respective fields. 

Plans for the next weeks:
- Get feedback from my mentor on the changes in quadprog.
- Begin intial work on optimoptions.


by Asma Afzal ( at July 24, 2015 05:40 AM

July 20, 2015

Piotr Held

Plans for randomize

This past week I have spent time trying to establish some testing framework for C++ methods and also trying to create a model for what the abstract base classes for cool, cost and permute should look like. I would like them to have the following methods/members:
  • Cost (Cost_fcn):
    • const Matrix *series-> pointer to the series for which a surrogate is generated
    • double cost-> current cost
    • void cost_transform (Matrix *)-> initial transformation which is used for better calculation of cost
    • Matrix cost_invert () const-> assigns to the input variable the inverse of the transformation performed above (to get the actual surrogate, not just a representation of it in a different form)
    • double cost_update (octave_idx_type nn1, octave_idx_type nn2, double cmax, bool &accept)-> perform quick update of cost (for a swap of elements under index nn1 and nn2) and decide if cost is smaller than maximum cost (cmax) if yes, accept new cost and return true otherwise reject new cost and return false
    • double cost_full()-> performs a full calculation of the cost, takes longer than cost_update()
    • getters/setters for cost
  • Cool (Cool_fcn):
    • double temp-> holds the current temperature
    • double cool (double cost, bool accept, bool &end)-> takes current cost, accept which holds whether the last cost_update() was accepted or not and returns the new temp and sets flag end to indicate if the simulated annealing is over
  • Permute (Permute_fcn):
    • Matrix *series-> holds the pointer to the series, and modifies the series only when exch() is called
    • void permute (octave_idx_type &n1, octave_idx_type &n2) const-> generates two indexes n1 and n2 that can be used to calculate Cost_fcn::cost_update()
    • void exch (octave_idx_type n1, octave_idx_type n2)-> exchange element under n1 with element under n2 in the series
Those are the methods I intend to have in the base/abstract classes which will be called by the Simulated Annealing runner code. I have not decided what that code should look like, but the current version seems to be working as well as the randomize program from TISEAN package.

I was also hoping to create a subclass of each of those abstract classes built to call GNU Octave code. This would allow the user to create their own functions without having to write anything in C++. However, this idea might not be practical for the following reasons:
  1. The example of Simulated Annealing provided in the TISEAN package (ver. 3.0.1), takes about 0.7 seconds to run using only C++ code and performs on average 900,000 calls to Cost_fcn::cost_update() and calling a simple function in Octave that many times (using the for loop) took 16 seconds
  2. I have trouble deciding how to neatly pass these functions/classes to randomize along with some parameters the user might want to include. I originally thought of using classdef - the new keyword introduced in Octave 4.0.0. I hoped to create an abstract Octave class and then let anyone subclass it to create their own cost, cool and permute classes. The problem is that classdef and all of the associated keywords are not documented, moreover according to Carnë Draug the help function will not recognize this new type of Octave class. So even if all of the needed functionality was available in Octave I might not be able to document it for the user
If obstacle 2. can be overcome it might still be beneficial for the package to create this type of functionality, regardless of how long the code will execute.

This week I plan to refine the design of the abstract classes as well as port more of the cost function options from TISEAN.

[Update]: I modified the design a bit and updated this post to fit the new code.

by Piotr Held ( at July 20, 2015 09:09 AM

July 14, 2015

Asma Afzal

Week 7: quadprog wrapping __qp__

I have wrapped quadprog on __qp__ instead of qp.m in [1].

Main differences between quadprog in [1] and qp.m.

- Input argument placement
  quadprog(H, f, A, b, Aeq, beq, lb, ub, x0, options)  =  qp (x0, H, f, Aeq, beq, lb, ub, [], A, b, options) 

- Check for empty inputs A and b
qp ([],H,f,Aeq,beq,lb,ub,[],A,[])

This works. qp simply ignores inequality constraints due to if checks
in lines  258, 266 and 275 of qp.m. Matlab gives an error if A is empty and b is not and vice versa.

quadprog (H, f, A, [], Aeq, beq, lb, ub)
Error: The number of rows in A must be the same as the length of b. I have added this check in line 181 in [1].

- Lambda output as a structure instead of a vector as in qp.m.

Ordering of lambda:
  • The order of lambda vector output (qp_lambda) from __qp__(in my code) is [equality constraints; inequality constraints; lower bounds; upper bounds]. 
  • The multipliers are present if the constraints are given as inputs so the size of qp_lambda depends on the size of constraints. 
  • Variables idx_ineq, idx_lb and idx_ub make sure I pick the right values. 

H = diag([1; 0]);
f = [3; 4];
A = [-1 -3; 2 5; 3 4];
b = [-15; 100; 80];
l = zeros(2,1);
[x, obj, info, qp_lambda] = qp ([], H, f, [], [],l,[],[], A, b)
[x,fval,exitflag,output,lambda] = quadprog (H, f, A, b,[],[],l,[])

qp_lambda =


lambda =
scalar structure containing the fields:

    lower =


    upper = [](0x0)
    eqlin = [](0x0)
    ineqlin =


Things to do:
  • Check the sign issue for lambda.eqlin (qp gives values -1* Matlab's)
  • Check if __qp__ changes the order of constraints. The values of lambda from qp.m in the last example in [2] are there but not coinciding with the respective constraints.
  • Move on to optimoptions.


by Asma Afzal ( at July 14, 2015 07:28 PM

July 09, 2015

Piotr Held

Progress report and plans

So far my progress has been as planned. Before the end of the midterm evaluation I was able to publish on my repository version 0.2.0 of the package, which included all of the functions from section Dimensions and entropies from the TISEAN documentation. As I mentioned in my previous post the functions that needed to be ported in this section are slightly different from what I wrote in my outline. The ported functions are:
  • d2
  • d1
  • boxcount
  • c2t
  • c2g
  • c2d
  • av_d2
I also wrote demos for most of those functions and updated the tutorial on the wiki page.

The first part of this week I spent improving on the build process. The function __c2g__ relies on C++ lambdas to work, therefore a configure script needed to be introduced to ensure the compiler has this capability. As was suggested by John Eaton, I tried to make the impact of not having that capability as small as possible. Currently if the compiler does not recognize C++ lambdas simply __c2g__ is not compiled and the function c2g does not work.

The plans

I was hoping to port all of the functions in the next section, Testing for Nonlinearity, by the end of the week. This might not be possible as randomize turned out to be a bigger function than I anticipated. It is actually not a function at all but, as the author of the TISEAN documentation puts it, "a toolbox". It generates surrogate data using simulated annealing. It needs to be supplied with three functions:
  1. the cost function 
  2. the cooling function -- how the temperature decreases
  3. the permutation function -- what to change every step
So currently if the user wants their own version of any of the functions above the user needs to write it in FORTRAN. My goal for this project would be to allow the user to write (use) their own octave function. The SA algorithm is an iterative method so using Octave code is not a good idea (as each line must be parsed when using for or while loops). As far as I understand the samin routine from the optim package will not suffice as it does not generate surrogate data, and has fewer options. Due to the size of this function it might take me some time to complete it.

I plan to tackle this problem as follows: I will rewrite in C++ the equivalent function to randomize_auto_exp_random and then try to refactor and modify the code to accept other functions. I plan to include all of the functions that are available in TISEAN in the Octave package, either through rewriting them or through linking to them. And I would like to make it easy for new functions to be added.

Further reading on randomize is available on the TISEAN documentation in the General Constrained RandomizationSurrogates paper Appendixrandomize description and randomize function extension description.

by Piotr Held ( at July 09, 2015 04:27 PM

July 02, 2015

Asma Afzal

Week 5 and 6: Refining fmincon

So my fmincon implementation is coming in shape [1]. 

[x,fval,exitflag,output,lambda,grad,hessian] = 

I came across a few issues which turned out to be bugs. Olaf pushed fixes in the central repository. Listing the issues for the record:
Setting gradc (the gradient of general equality/inequality functions)A bug in nonlin_min.m (and __nonlin_residmin__.m)
      objective_function = @ (p) p(1)^2 + p(2)^2;
 pin = [-2; 5];
 constraint_function = @ (p) p(1)^2 + 1 - p(2);
 gradc = @(p) [2*p(1);-1];
 [p, objf, cvg, outp] = nonlin_min (objective_function, pin, optimset
("equc", {constraint_function, gradc}))

error: function handle type invalid as index value
- Giving linear inequality/ equality constraints to lm_feasible. A bug in nonlin_min.m
      f = @(x) -x(1) * x(2) * x(3);
S = [1  -1;   2  -2;   2  -2]
b = [0;72];
x0 = [10;10;10];
[x,fval] = nonlin_min( f, x0, optimset ("inequc",{S,b}) )
      error: __lm_feasible__: operator -: nonconformant arguments (op1 is
3x1, op2 is 3x0)
-Any zero value in initial guess vector for nonin_residmin/nonlin_min gave an error. Required a minor change of sign(0)==0 in __dfdp__.m.
      k = 1:10;
func = @(x) 2 + 2 * k - exp (k * x(1)) - exp (k * x(2));
x0 = [0;0.5];
x = nonlin_residmin(func,x0)

warning: division by zero
warning: called from
    __dfdp__ at line 367 column 21
    __nonlin_residmin__> at line -1 column -1
    __lm_svd__ at line 191 column 9
    __nonlin_residmin__ at line 1125 column 21
    nonlin_residmin at line 98 column 21
    runlsqnonlin at line 9 column 3
error: svd: cannot take SVD of matrix containing Inf or NaN values

Functionality for Returning Hessian and Gradient
New options "ret_objf_grad" and "ret_hessian" to be introduced in nonlin_min (by Olaf). If anyone of these options is set to true, the 'outp' structure output argument of nonlin_min will contain additional fields .objf_grad and .hessian. My code currently checks this. 

Rearranging values of lambda in the fields of a structure.
- For lm_feasible, outp will contain an additional field lambda, a structure which contains Lagrange multipliers in fields separated by constraint type.

I added an additional feature in [1] to cater for the non linear constraint function set using deal();
      c = @(x) [x(1) ^ 2 / 9 + x(2) ^ 2 / 4 - 1;
        x(1) ^ 2 - x(2) - 1];
ceq = @(x) tanh (x(1)) - x(2);
nonlinfcn = @(x) deal (c(x), ceq(x));
obj = @(x) cosh (x(1)) + sinh (x(2));
z = fmincon (obj, [0;0], [], [], [], [], [], [], nonlinfcn)

z =

To do:
1- Write test cases/refined examples for for lsqnonlin, lsqcurvefit, nlinfit and fmincon.
2- Start wrapping quadprog to __qp__ instead of qp.m (because of the ordering of the lambda output).

by Asma Afzal ( at July 02, 2015 08:52 PM

June 26, 2015

Piotr Held

Progress report

The main goal of this post will be to create a progress report before the coming midterm assessment.
As I mentioned before I planned to complete the Dimensions and entropies section of the TISEAN documentation. This seems to be still a realistic goal.

Currently I have ported d2, av-d2, c2g, c2t along with writing documentation and demos for them. The current state of the tests needs improvement because they rely heavily on external files generated using the corresponding TISEAN programs. Because most of those functions/programs are closely linked I plan to improve on this feature once functions from the entire section are ported.

Currently I am working on c1 which already passes it's test. Once I complete it and write the documentation and demo, the only programs/functions that will need to be ported are boxcount and c2d. Once they are complete I plan to release version 0.2.0.

My elaborated proposal located on the octave wiki states that I planned to also port c2. Although the source code for such a program does exist in the TISEAN package (ver. 3.0.1), it does not seem to be mentioned in the documentation. Furthermore, installing the package on a computer does not give access to this program. Also, it seems to be redundant with other programs in the package. Therefore, I will not port it.

by Piotr Held ( at June 26, 2015 03:37 PM

June 24, 2015

Asma Afzal

Progress Update: Midterm Evaluation

Adding functions to the Optim package for Octave using existing back-ends.

Expected deliverables before midterm:
  • 'lsqnonlin' using 'nonlin_residmin'
Done in [1]. Differences in backends, nonlin_residmin uses "lm_svd" algorithm for optimization as currently the only backend. However, lsqnonlin in Matlab can choose between "trust-region-reflective" and "Levenberg-Marquardt" (LM) algorithms.
Another difference is in complex inputs. lm_svd does not support complex valued inputs whereas Matlab's LM algorithm can accept complex input arguments. One way of providing complex inputs to lsqnonlin in Octave is to split the real and imaginary parts into separate variables and running the optimization.  
  • 'lsqcurvefit' using 'nonlin_curvefit', 'nonlin_residmin', or 'lsqnonlin'
Done in [2] using nonlin_curvefit. lsqcurvefit is very similar to lsqnonlin with only a few minor interface differences. lsqcurvefit explicitly takes independent variables and the observations as inputs while these values can be wrapped inside the objective function while using lsqnonlin. Additional bounds for the optimized parameters can be specified. 
  • 'nlinfit' using 'leasqr',
I wrapped nlinfit on nonlin_curvefit and curvefit_stat as leasqr repeats the optimization to compute the additional statistics (Jacobian and Covariance matrix) while curvefit_stat saves this computation overhead. I have partially implemented nlinfit in [3] (It hasn't been thoroughly reviewed yet). Two missing features are: 1) Error models and Error parameters estimation, and 2) Robust Weight function. Meanwhile, no such functionality exists in the Octave's optimization backends for the missing features. My current implementation supports the input of scalar positive array of weights for robust regression.  
Since nlinfit is from the statistics toolbox of Matlab, it uses statset and statget to create and get options respectively. I created additional functions statset, statget and __all_stat_opts__ with minor changes to the code in optimset, optimget and __all_opts__.
  • 'fmincon' using 'nonlin_min',
In progress [4].

Future goals:
    1. Complete fmincon implementation.
    2. Create solver specific options using optimoptions and desirably still be able to use optimget. 
    3. Arranging lambda output for quadprog by wrapping it on __qp__ instead of qp.m
    4. Test cases for all the implemented functions.


by Asma Afzal ( at June 24, 2015 05:46 PM

Juan Pablo Carbajal

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 June 24, 2015 02:48 PM

June 22, 2015

Asma Afzal

Week4: fmincon wrapping nonlin_min

Time flies.. A third of the way through already..

fmincon/nonlin_min is the most elaborate function of all that I have previously 
implemented so before actual coding I would like to thoroughly check the the 
mapping of arguments and options. 

[x,fval,exitflag,output,lambda,grad,hessian] = fmincon (fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
[x,fval,exitflag,output] = nonlin_min (fun,x0,settings)

A,b Aeq, beq- Linear inequality and equality constraints
Inequality constraints:
Matlab standard: A * x - b <= 0 ,
Octave standard: m.' * x  + v >= 0. 
This implies: m=-A.' , v=b.
Set in Octave using 
optimset ("inequc", {m,v})
Similar for equality constraints:
optimset ("equc", {m,v})

lb, ub- Lower and upper bounds
Set in Octave using optimset
 optimset ("lbound", ..., "ubound",...)

nonlcon- Nonliner constraint function handle 
In Matlab, nonlinear constraints are given in a function with the following format

function [c,ceq,gradc,gradceq] = mycon(x)
c = ...     % Compute nonlinear inequalities at x.  
ceq = ...   % Compute nonlinear equalities at x.
% Optional output arguments:
gradc = ...   % Gradient of c(x).
gradceq = ...   % Gradient of ceq(x).
options = optimoptions('fmincon','GradConstr','on')

Alternative to nonlcon in nonlin_min
  optimset ("equc", {constraint_function})

Options- Options common to both fmincon and nonlin_min
User-supplied Gradient
In Matlab, the objective function must return the second output when the 
GradObj option is set:
User-supplied Hessian
In Matlab, the objective function must return the third output when the Hessian 
option is set:

Lambda- returned as a structure in Matlab but as a vector field of the "output" 
ouput argument in Octave.

Things to investigate:
1. Algorithm mapping. 
2. Exitflag mapping.
3. Setting gradc (the gradient of general equality/inequality functions). The second 
entry for the equc/inequc setting implements this feature in Octave as stated by 
optim_doc but I was unable to make it work properly.
4. Returning additional statistics (Hessian, Gradient)- I used 
residmin_stat/curvefit_stat previously for lsqnonlin and lsqcurvefit. No such function for nonlin_min.
5. Rearranging values of lambda in the fields of a structure.

P.s. this week I spent quite some time with Mercurial and how to possibly work in the optim package. There was a slight miscommunication/confusion but it is clear now and I will continue to publish my code on github..

by Asma Afzal ( at June 22, 2015 08:53 PM

June 18, 2015

Juan Pablo Carbajal

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 June 18, 2015 10:34 AM

June 14, 2015

Piotr Held

Improving the code

TISEAN was originally written as a command line set of programs. Because of this all the code is not very portable and many variables are global ones. So far this has been dealt with by creating local variables and extending the number of variables in function calls (in some cases up to 11). This is not optimal for code clarity, ease of maintaince and because many of the variables are passed as values (they are parameters) it also caused a slight slowdown in execution speed.

Due to all of these downsides I have contemplated possible solutions which I will attempt to describe.

Using structs

One idea that came to mind is to pack all of the global variables into a struct and pass the struct to all of the functions and obtain the global variables from this struct. This solution certainly solves the problem of passing so many parameters to functions. However, it does not improve code portability because every *.cc oct-file function needs its own struct. This solution is also problematic because all of the names of global variables have to be referenced now through the struct so center[i][j] would become[i][j] (obviously the name of the struct could be as short as p).

Using classes

Another quite simple solution would be to create a class. This class would have data members that were previously the global variables and function members that were the old functions called from old main(). As there are similarities between different TISEAN programs, it could be possible to even create a prototype class and inherit from it.

There are however downsides to this option as well. First of all, Octave code guidelines specifies that classes should be in separate files. This would mean creating 2 more files for each program that was ported using the C++ wrapper. Apart from that, the memory might have to allocated using new/delete, because the preferred method of using the macro OCTAVE_LOCAL_BUFFER might be difficult (or impossible) to apply to this case. This objection can be worked around in other ways, such as using Array classes to allocate the data and then get a pointer to them using fortran_vec().


Performing the aforementioned code improvement, although helpful, is not critical. Therefore any attempts to implement it will be deferred until after the functions outlined at the beginning of the project are complete.

Timeline update

So far I have been giving progress reports on the TISEAN porting project. This time, however, I would like to also compare the outlined schedule for the project with the actual progress made.

Since the last post I have additionally ported:
  • xzero
  • lyap_r
  • lyap_k
  • lyap_spec
In one of my first posts I stated that I would like to finish Dimensions and Entropies before the midterm assessment. Currently I have finished up Lyapunov Exponents and I plan to start working on Dimensions and Entropies this week. Since there are 2+ weeks to the Midterm Assessment I believe it is possible to complete all of the goals for this section of the project as planned.

by Piotr Held ( at June 14, 2015 08:58 PM

Asma Afzal

Week3: nlinfit, statset, stat_get and __all_stat_opts__

So this week I achieved the following milestones:

Wrapping nlinfit on nonlin_curvefit

[beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(X,Y,modelfun,beta0,options,Name,Value)

Implementation [1]:

I chose not to wrap nlinfit on lsqcurvefit because 
  1. We might end up wrapping lsqcurvefit on lsqnonlin eventually so it is undecided.
  2. The default options for lsqnonlin/lsqcurvefit are different from nlinfit. 
Missing features:
  1. RobustWgtFun - The field RobustWgtFun in options can be provided with a function handle which computes the residuals at every iteration. The backend optimization algorithm in Octave currently does not support this functionality.  
  2. Name-Value pairs. Currently the only implemented one is "weights",  which takes an array of weights for the weighted optimization. "ErrorModelInfo" and "ErrorParameters" are not implemented. The possible error models include, "constant", "proportional" and "combined". The error model also translates to a weight function which helps in reducing the effect of outliers.   
  3. ErrorModelInfo- output field which gives information about the error variance, and estimates the parameters of Error models. 
Setting options using statset

In Matlab, options for statistics toolbox are set using statset [2].  The functionality is almost identical to optimset but the different functions in Matlab are because of different toolboxes (statset for statistics and optimset for optimization toolbox).

Added functions:
statset.m, statget.m and __all_stat_opts__.m [3]-[5]
Creating these functions was pretty straight forward.

Still to do:
  1. Have to check if the weighted residual and weighted Jacobian output in octave is consistent with Matlab and further refine the functions with the feedback from my mentors.
  2. Move on to fmincon wrapping nonlin_min.

    by Asma Afzal ( at June 14, 2015 08:11 PM

    June 09, 2015

    Asma Afzal

    Week 2: lsqnonlin and lsqcurvefit

    A bit late blogging about week 2. 

    Almost completed functions lsqnonlin and lsqcurvefit. 
    Successfully mapped user-specified Jacobian. 
    In Matlab, if the Jacobian option is set to "on",  the model function must return 
    a second output which is the Jacobian function.
    In Octave, the Jacobian function handle is given to the dfdp option using optimset.

    Lsqnonlin function description:

    [x,resnorm,residual,exitflag,output,lambda,jacobian] = ...

    This function maps on:

    [x, residual, exitflag, output] = nonlin_residmin (fun, x0, options)

    Features of lsqnonlin in Octave:
    • Input arguments: Acceptable forms are lsqnonlin(fun,x0), lsqnonlin(fun,x0,lb,ub) and lsqnonlin(fun,x0,lb,ub,options)
    • Outputs
      • x, exitflag, residual and output currently same as nonlin_residmin.
      • resnorm=sum(residual.^2)
      • Lambda is computed using the complementary pivoting in __lm_svd__.
        It's values differ from Matlab's due to the difference in backends.
      • Jacobian is computed using the function residmin_stat ().

    Lsqcurvefit function description:

    [x,resnorm,residual,exitflag,output,lambda,jacobian] = ...

    This function maps on

    [x, fy, exitflag, output] = nonlin_curvefit (fun, x0, xdata, ydata, options)

    Features of lsqcurvefit in Octave:
    • Input arguments: Acceptable forms are lsqcurvefit (fun,x0,xdata,ydata), lsqcurvefit (fun,x0,xdata,ydata,lb,ub) and lsqcurvefit (fun,x0,lb,ub,xdata,ydata,options)
    • Outputs
      • x, exitflag, residual and output currently same as nonlin_curvefit.
      • residual = fy-ydata, resnorm = sum(residual.^2)
      • Lambda and Jacobian same as in lsqnonlin.
    There are only minor interface differences between lsqcurvefit and lsqnonlin.

    This week's plan:

    • Hopefully, with lsqnonlin and lsqcurvefit wrapped up, I'll move on to nlinfit
    • Three key challenges need to be addressed when wrapping nlinfit using nonlin_curvefit
      and curvefit_stat:
      • Weight functions: Currently, no such functionality exists in nonlin_curvefit,
         where a user can specify weight functions to perform Robust regression
        (weights computed using the specified function in every iteration).
      • Error Models and ErrorModelInfo
      • Setting options using statset instead of optimset or optimoptions. 

    by Asma Afzal ( at June 09, 2015 08:17 PM

    June 07, 2015

    Piotr Held

    Analyzing lfo-run

    I have written tests that compare lfo-run from TISEAN to the ported version lfo_run. The test that uses amplitude.dat works perfectly, but when I analyzed the results both programs/function gave for henon (Henon Maps) I ran into some problems. I will attempt to describe them.

    Input data

    The problems occur when analyzing a 1000 element Henon map (henon(1000)). For all of the implementations if I used a simple call with default parameters (m = 2, d= 1) the programs would quit due to a matrix singularity. The problems arose when (m = 4, d= 6) was used. With these parameters the program gave various results for various implementation methods.

    It is important to note that the prediction that I was testing tried to predict 1000 future elements (default for all implementations) on the basis of given 1000 elements.


    There are 3 implementations I used:
    1. The TISEAN implementation (uses lfo-run)
    2. The implementation similar to 1. but compiled as c++ and wrapped in enough code to run as m-file (uses __lfo_run__ and invert_matrix())
    3. The implementation that uses Matrix::solve() method
    I tried to find out if maybe method 1. and 2. do not differ due to a bug that was introduced while porting. I therefore ported it twice (the second time to a very rudimentary stage) and both times the same results were encountered. I do not understand why there is a discrepancy between between these two implementations.


    Since the goal of this project is to port TISEAN functions I have compare implementation 1. with 2. and 1. with 3, to see what differences I come across.

    A comparison between implementation 1. and 2. results in an error from implementation 2. The function generates about 700 elements (of the default 1000) and then throws an error that the forecast has failed.

    The comparison between implementation 1. and 3. is much more fruitful as the results are the same for about 150 elements and then they begin to differ (see Fig. 1.)
    Fig. 1 Comparison between the TISEAN implementation and using Matrix::solve() in TISEAN package from Octave
    These results can be achieved by cloning the repo, doing make run and running the script:

    cd tests/lfo_run/; test_lfo_run


    Before I give my suggestions for what is the cause for these discrepancies I would like to discuss another interesting discrepancy. This discrepancy is the maximum difference between the solution of the equation system obtained from implementation 2. and 3. When using implementation 3. for the forecast this value was 8.5e-14, but when using implementation 2. for the forecast this difference was 5e-13.

    I believe this is because the computational error is accumulated throughout the program. Each new forecast point is dependent on the previous ones. Moreover the Kahan algorithm (compensated summation) is never used in the TISEAN implementations. Even matrix multiplication (as seen e.g. in multiply_matrix())  uses the simple, but error accumulating for(...) sum += vec[i].

    As to why implementation 1. and 2. give different results I have two theories: either there still is a bug which I was unable to detect, or some compilation difference (e.g. a linked library) between the TISEAN program written in C (lfo-run) and the TISEAN package function written in C++(__lfo_run__).


    The question that poses itself is whether this warrants a rewriting of other TISEAN function that use the simple summation, or if this problem can be ignored. The authors of TISEAN said in their introduction that blindly using the programs they wrote may result in unintended or even wrong results. Trying to predict 1000 elements of a 1000 element Henon Map using first order local linear prediction might be considered a bad use-case.

    Progress report

    Since the last post I wrote a tutorial on I also ported:
    • lzo_gm
    • lzo_run
    • ikeda
    • lfo_run
    • lfo_ar
    • lfo_test
    • rbf
    • polynom
    During my work I discovered that polynom has similar functions (polynomppolyparpolyback) which provide extra options for performing polynomial fits. I will not include those functions in the project now, but they have a high priority once I finish all of the functions I outlined for this project.

    These newly ported functions aren't completely polished (some need demos and documentation) but they pass tests and don't have memory leaks. Once I clean these functions up and port xzero the last function in this section, I intend to create version 0.1.0 of the package. With this version I intend to branch out the repo to have a 'devel' and a 'stable' branch.

    Afterwards, I will add more information to the tutorials on the wiki page.

    by Piotr Held ( at June 07, 2015 03:16 PM

    June 02, 2015

    Asma Afzal


    Answering some questions to better understand the behavior of optimoptions.
    Although the answers have been discussed in [1], I am pasting some examples for clarity.
    *I was using the terms solver/algorithm in the wrong context before. Solver names are functions such as lsqnonlin, fmincon, etc. and one solver can have multiple algorithms such as interior-point, lev-mar, etc. *

    1) Is there an error or warning if optimoptions is used to assign an option not contained in the specified algorithm? 

    If the option belongs to a different algorithm of the solver, Matlab stores it as options "not" used by the current algorithm so if we change the algorithm to which the option belongs, we do not need to set it again. 

    For eg. The default algorithm of lsqnonlin is trust-region-reflective but when I try to set the option 'ScaleProblem' which is specific to 'levenberg-marquardt' algorithm, I get: 

    opts = optimoptions ('lsqnonlin', 'ScaleProblem', 'Jacobian')

    opts = 

      lsqnonlin options:

       Options used by current Algorithm ('trust-region-reflective'):
       (Other available algorithms: 'levenberg-marquardt')

       Set by user:
         No options set by user.


     Algorithm: 'trust-region-reflective'
        DerivativeCheck: 'off'
         .               .  
         .               .  
         .               .  
         .               .  
        TolX: 1.0000e-06

       Show options not used by current Algorithm ('trust-region-reflective')

    Set by user:
        ScaleProblem: 'jacobian'

    This gives the same result using dot notation:

    opts = optimoptions ('lsqnonlin')
    opts.ScaleProblem = 'Jacobian'

    2) ... or to assign an option that does not exist in the specified solver or any solver?

    This gives an error. Trying to set the option for SQP Algorithm, which is not used by lsqnonlin:

    opts=optimoptions('lsqnonlin', 'MaxSQPIter', 100)
    Error using optimoptions
    'MaxSQPIter' is not an option for LSQNONLIN

    using dot notation:

    opts.MaxSQPIter = 100;
    No public field MaxSQPIter exists for class optim.options.Lsqnonlin.
    3) If options are transfered to a different solver with optimoptions, are there errors or warnings if the new solver does not have some of these options?

    No errors or warnings. The options common to both solvers are copied. They could be options of different algorithms. For eg. 

    opts = optimoptions ('fmincon', 'Algorithm', 'sqp', 'TolX', 1e-10) 
    opts_lsq = optimoptions ('lsqnonlin', opts) 

    Options set in opts_lsq are: 
                    PrecondBandWidth: 0 
                    TolX: 1.0000e-10 

    The option PrecondBandwidth belongs to the trust-region algorithm of fmincon solver. 
    Another option copied from opts in opts_lsq belongs to the lev-mar algorithm of lsqnonlin. It is the stored option as mentioned in 1) 

      ScaleProblem: 'none' 
    4) Can options returned by optimoptions be used with optimset, and   vice versa? 

    This returns an error in both cases:

    opts = optimoptions ('lsqnonlin','MaxIter',100);
    opts = optimset (opts, 'TolX', 1e-6);
    Error using optimset
    Cannot use OPTIMSET to alter options created using OPTIMOPTIONS.
    Use OPTIMOPTIONS instead.

    opts = optimset ('TolX', 1e-6);
    opts = optimoptions ('lsqnonlin', opts);
    Error using optimoptions
    Invalid option name specified. Provide a string (such as 'Display').


    by Asma Afzal ( at June 02, 2015 07:51 PM

    May 31, 2015

    Asma Afzal

    lsqnonlin wrapping nonlin_residmin

    So the first week of GSoC is officially over.

    I was working on lsqnonlin. My code is accessible here:

    [x, resnorm, residual, exitflag, output, lambda, jacobian]...
    = lsqnonlin (fun, x0, lb, ub, options)

    [x, resid, cvg, outp] = nonlin_residmin (f, x0, settings)

    A recap of encountered problems:
    1. output lambda- it wasn't previously returned from the backend [1]
    2. options- In addition to optimset, Matlab uses optimoptions to set options specific to a certain optimization function. optimoptions
      • Creates and modifies only the options that apply to a solver
      • Shows your option choices and default values for a specific solver/algorithm
      • Displays links for more information on solver options and other available solver algorithms [2].
      octave currently does not have this functionality.  For more discussion on this, check [3]
    Things to do:
    1. check how user-specified jacobian is to be provided.
    2. check for matrix/complex inputs.
    3. come up with a plan for writing optimoptions in octave.

    by Asma Afzal ( at May 31, 2015 07:41 PM

    May 28, 2015

    Juan Pablo Carbajal

    El que busca encuentra

    El foro de discusión tiene una nueva entrada: una pregunta de Florencia que está resolviendo el desafío nro. 1.5.

    La pregunta de Floppy nos da pié para discutir sobre un tipo de funciones muy útiles. Funciones que nos permiten obtener la posición, dentro de un vector, de elementos que cumplen cierta condición lógica. Esto es muy similar a lo que hicimos en la clase ¿Verdadero o falso? cuando buscamos números pares.

    Find == Econtrar

    Las operaciones lógicas generan matrices con valores verdadero o falso. Así, la siguiente desigualdad 

    tf = [-1,5,-3] > 0

    Genera un vector de valores lógicos (verdadero o falso), en este caso [false,true,false]. Estos vectores lógicos son de muchísima utilidad: nos permiten filtrar matrices o buscar elementos que cumplan ciertas condiciones lógicas. Para realizar lo segundo usamos la función find.
    La función find (encontrar, en inglés) nos devuelve la posición de aquellos elementos que corresponden al valor "verdadero". Por ejemplo:

    find ([false, true, false])

    devuelve un 2, porque el segundo elemento es un valor verdadero. En Octave todo aquello que no es cero o vacío se considera verdadero. Es decir que

    find ([0,-3, 0])

    también devuelve 2, porque el segundo elemento es el único no cero. Cuando varios elementos son verdaderos, la función find devuelve todas las posiciones correspondientes, es decir

    find ([false, true, true])

    devuelve [2, 3].
    La función también acepta matrices como entradas y nos devuelve el índice lineal de los elementos verdaderos. En Octave las matrices están ordenadas por columnas (en Inglés se dice column-major order), es decir contamos de arriba para abajo y de izquierda a derecha. Por ejemplo en una matriz de 3x4 los indices lineales son

    Es decir que si ejecutamos

    find ([0 1 0; 0 0 1; 0 1 0])

    obtenemos [4,6,8]
    Para obtener las filas y columnas de los elementos verdaderos, tenemos que llamar a la función con dos argumentos de salida. Si ejecutamos

    [i, j] = find ([0 1 0; 0 0 1; 0 1 0]);
    ans =
       1   2
       3   2
       2   3

    obtenemos en i las filas y en j las columnas de los valores verdaderos.

    ¿Es miembro?

    Otra función muy útil es ismember (es miembro, en inglés). Esta función toma como entrada dos argumentos y nos dice si los elementos del primero están en el segundo. Solo voy a proveer un corto ejemplo y los invito a que lean la ayuda de la función (en Octave: help ismember) y pongan sus preguntas en el foro de discusión.
    Supongamos que queremos saber si la letra "a" está presente en un string

    x = "Hay letras a en esta frase?";

    podemos ejecutar

    tf = ismember ("a", x)
      tf =  1

    Claro que sí. La función también se puede usar para obtener la posición de las letras buscadas. Para esto intercambiamos los argumentos de entrada y buscamos el string x en la letra "a".

    pos = find (ismember (x,"a"))
      pos =
          2    9   12   20   24

    Pregunta 1. ¿Puedes explicar el porqué de este truquillo?

    La función ismember también puede darnos la posición de los elementos encontrados, pero cuando el elemento buscado se repite (como en nuestro caso), solo nos devuelve la última ocurrencia

    [tf, pos] = ismember ("a", x)
      tf  = 1
      pos = 24

    Esta funcionalidad es útil cuando no hay repeticiones y se ilustra mejor con el ejemplo en la ayuda de ismember

    a = [3, 10, 1];
    s = [0:9];
    [tf, pos] = ismember (a, s)
    tf = [1, 0, 1]
    s_idx = [4, 0, 2]

    Las funciones find e ismember son de mucha utilidad y tienen varios otros modos de uso. Lee la ayuda, explora y pregunta en el foro!

    by Juan Pablo Carbajal ( at May 28, 2015 07:26 AM

    May 27, 2015

    Piotr Held


    I am happy to announce that I completed the first two sections.

    Nonlinear Noise Reduction
    From this section I added:
    • ghkss
    Thus deprecating project.  It is important to note that data allocation in ghkss (only) is done via new/delete because of the way the function was originally written. In the future this will be replaced with OCTAVE_LOCAL_BUFFER (OLB) macro.

    Closer examination has revealed another interesting function in this section: nrlazy, which according to the documentation is similar to lazy. Because of this similarity porting it has low priority.

    Phase Space Representation
    From this section I added:
    • mutual
    • false_nearest
    • poincare
    The function poincare was discovered when further studying the documentation and was omitted in the initial plan, but as it seems important in the documentation it was ported.

    The function false_nearest is in a similar situation to ghkss, that is, new/delete is used instead of OLB. This will be improved on in the future.

    Since corr does not need to be ported, this section is complete.

    Nonlinear Prediction
    This next section is well on the way as the following have been ported:
    • upo
    • upoembed
    • lzo_test 
    The function predict turned out to be essentially the same with lzo_test. This had not been verified so far, but according to the documentation they do essentially the same. The additional option that predict has (flag -v) can be easily replaced using GNU Octave's std when calling lzo_test with parameter r. Therefore porting predict will be most likely unnecessary.

    The function upoembed is closely associated with upo. It takes the output of upo and creates a cell array of delay vectors for each orbit produced by upo. It was not mentioned in the original outline, but it is an important function for the package.

    The state of upo is not optimal. The original implementation only supported input up to 1e6 data points. This might not be a big problem as calculating upo on a 1e4 henon map takes about 8 seconds, so a 1e6 would take about 800 seconds ~ 15 min. Changing this might be problematic as the main data in the FORTAN program is in a common block (it is stored in a global variable), which cannot contain arrays of variable dimensions. The authors of TISEAN chose 1e6, but because of how the program is written making it unlimited would not be trivial. It might be beneficial to lift this number to 1e8 or 1e9 but one must keep in mind that because of how the FORTRAN program is written it will always allocate a local array of the maximum size possible (currently 1e6 elements). If the maximum input length is lifted to 1e8 or 1e9 the size of the data allocated by Octave and the local copy that the FORTRAN program uses can be a sizable amount of memory to allocate. Moreover it is important to note that each data point will be a real*8 (not a real*4). This brings me to the next point.

    FORTRAN data types
    This topic has been problematic for me from the very beginning. I had trouble realizing just how the dummy variables work. After some research I found that one can pass -freal-4-real-8 and promote all real*4 to real*8. This is beneficial as the input into a FOTRAN program is passed as doubles. This caused serious issues as I needed to ensure that when I call any TISEAN or SLATEC function/subroutine I needed to use real*4 instead of real*8. The solution I originally used was to copy the input variables into local variables. Apart from eliminating potential bugs and improving code complexity the previously mentioned flag also allows the FORTRAN programs to have the same type of precision expected from GNU Octave programs.

    Other things
    I spent my time also on other things. I significantly changed the Makefiles, removed all compilation warnings and everywhere except for ghkss and false_nearest moved from new/delete to OCTAVE_LOCAL_BUFFER or some type of Array.

    by Piotr Held ( at May 27, 2015 06:33 PM

    May 26, 2015

    Juan Pablo Carbajal

    ¿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 May 26, 2015 03:17 PM

    Antonio Pino

    On galleries and the beginning of summer

    After the community bonding period and before starting today the coding period, I will briefly list the transformation that has undergone my initial proposal: from just implementing new algorithms and then add them to GNU Octave, to various modifications of GNU Octave itself so that Higham's toolboxes run smoothly and in the end add the new algorithms. Sticking to what I said before, I expect to be doing the modifications (e.g. new bugs, patches, toolboxes) most of the first half of the coding period. From there we aim to go as far as we can about matrix functions, I will do so. These changes are noted in the new time line.
    During the community bonding period I have been setting up, you might have seen me at freenode. I have been becoming more and more acquainted with GNU Octave, and found out that the gallery function was broken, with unassigned variables and missing auxiliary functions. This function will prove useful to test matrix functions, because the eigenvalue decomposition strategy (if $A=VDV^{-1}$  then $f\left(A\right) = V f\left( D \right) V^{-1}$) yields a big error using the ill conditioned matrices gallery provides. Another example are the useful positive definite matrices that have a computable principal p-th root. gallery is indeed interesting for anyone looking for a matrix with a special characteristic to test an edge case of a function.
    Besides, I also have had the chance to see a plethora of Matlab-style short circuit operators, looping with infinite ranges, and even weird undocumented functions like superiorfloat (it returns either "single" or "double" strings depending on the input). What let Carnë to point me to Undocumented Matlab blog, where they document its unsupported hidden underbelly. More quirks (and their solutions) on the next post. 
    A final thank goes to The Project in general and my mentors (Carnë and Mario) in particular for this opportunity. I hope everyone pleasantly codes their summer away!

    by Antonio Pino Robles ( at May 26, 2015 12:44 AM

    May 25, 2015

    Antonio Pino

    GSoC 2015 - Matrix Functions in GNU Octave

    A brief intro

    First of all, let me introduce myself: I am Antonio Pino Robles—an Electronic Engineering student from the Basque Country—and I will be improving matrix functions in GNU Octave this summer, following Google Summer of Code program.
    The idea behind this is quite simple: given a square matrix $M\in \mathbb{C}^{n \times n}$ and a function $f$, GNU Octave will compute $f\left(M\right)$. You may think of them as an extension to scalar functions, i.e. starting from $f:\mathbb{C}\rightarrow \mathbb{C}$ compute $f:\mathbb{C}^{n \times n}\rightarrow \mathbb{C}^{n \times n}$. Their implementation is quite different, though. (Check Golub and van Loan's book[0] and the Short Course by Higham and Lin[1] for further info.)
    Let me note that matrix functions are already part of octave: expm, logmsqrtm in octave itself and funm,  trigonometric and hyperbolic matrix functions in the Linear-Algebra Octave-Forge package. There are also GPLed toolboxes by Nicholas J. Higham, namely the mctoolbox[2] and the mftoolbox[3]; furthermore, GPLed software from the NAMF group—led by N. J. Higham at The University of Manchester—is available as well.
    Hence, on a first part octave will be modified in order to run the toolboxes—as they are—smoothly , and then the existing implementations will be improved by means of updating their algorithms.
    Finally, for a more detailed description of the project please refer to my octave-wiki page:

    Agur bero bat!

    [0] G.H. Golub and C.F. Van Loan. Matrix Computations, 4th Edition. The Johns Hopkins University Press, Baltimore, USA, 2013.
    [1] Nicholas J. Higham and Lin Lijing, Matrix Functions: A Short Course, preprint, (2013).
    [2] N. J. Higham. The Matrix Computation Toolbox.
    [3] N. J. Higham. The Matrix Function Toolbox.

    by Antonio Pino Robles ( at May 25, 2015 11:05 PM

    May 23, 2015

    Asma Afzal

    Quadratic Programming quadprog, qp

    Just a post to review the things I've learned over the past week:

    Intuitive explanation of Lagrange multipliers [1]. Given a D-dimensional function to optimize and a set of a few equality constraints.
    • The vector normal to the function should be a scaled version of the vector normal to the constraint function (In other words, both normals are parallel to each other). This scaling factor is the Lagrange multiplier corresponding a particular constraint. 
    • $\nabla f=\lambda_1\nabla g+\lambda_2 \nabla h$, where $g$ and $h$ are two constraint functions.

    KKT Conditions: Generalized method of Lagrange multipliers to be applicable on Inequality constraints.  $g_i(x)-b_i\geq 0$

    • Feasibility- $g_i(x^*)-b_i\geq 0$
    • Stationarity- $\nabla f(x^*)-\sum\limits_i\lambda_i^*\nabla g_i(x^*)=0$
    • Complementary Slackness $\lambda_i^*(g_i(x^*)-b_i)=0$
    • Positive Lagrange multipliers $\lambda_i \geq 0, \forall i$

    In case we obtain negative Lagrange multiplier, the constraint corresponding to the most negative multiplier is removed and the optimization is performed again until all multipliers are positive.

    A bit about active set algorithm [2]:

    • Possibility that none of the constraints are active or may be some are active. We only need to solve for equality constraints that are active at the optimum (binding).
    • When we have an active set $S^*, x^* \in F$, where $F=\{x|A x^*\leq b\}$,$\lambda^* \geq0$, where $\lambda^*$ is the set of Lagrange multipliers for equality constraints $Ax=b$ 
    • Start at $x=x_0$ and initial active set
    • Calculate $x^*_{EQP}, \lambda_{EQP}^*$ which minimizes the EQP defined by the current active set. Two possible outcomes:
      1. $x^*_{EQP}$ is feasible  ($x_{EQP}^* \in F$). Set $x=x_{EQP}$ and check Lagrange mulitpliers $\lambda_{EQP}^*$. If positive, solution found! Otherwise, remove constraint with $\min(\lambda_{EQP}^*)$ and repeat.
      2. $x^*_{EQP}$is infeasible. We move as far as possible along the line segment from $x_0$ to $x^*_{EQP}$ while staying feasible. Add to $S$ the constraint we encounter that prevents further progress. This is the blocking constraint.

    Quadratic programming:

    $\min\limits_{x}\frac{1}{2} x^THx + x^Tq$
    $ A_{eq} x = b_{eq}$
    $lb \leq x \leq ub$
    $A_{lb} \leq A_{in}x \leq A_{ub}$

    What qp.m does:

     [xobjinfolambda] = qp (x0HqAeqbeqlbubA_lbA_inA_ub)     
    • Checks feasibility of initial guess $x_0$
    • Checks size of inputs and that they make sense.
    • Checks if bounds lb,ub too close or A_lb or A_ub too close. If they are very close then the inequality is treated as an equality constraint instead.
    • Checks if any bound is set to Inf or -Inf. qp simply strikes it off.
    • Calls backend solver __qp__ using null space active set algorithm. 
    The ordering of lambda

    • quadprog returns Lagrange multipliers in a structure (with fields upper, lower, eqlin, ineqlin) and the multipliers corresponding to the constraints not provided are left empty.  
    • In qp, lambda is a column vector with Lagrange multipliers associated to the constraints in the following order: [equality constraints; lower bounds; upper bounds; other inequality constraints]
    • The length of lambda vector output from qp depends on the number of different constraints provided as  input.
    • Two issues in wrapping qp.m
      1.  the order, i.e. the position of the bounds constraints within the inequality constraints) is not specified by qp.m. The code could change and the ordering too. 
      2. qp.m strips off the INF constraints before calling __qp__ but does not process the lambda (returned by __qp__) accordingly.
      • Solution:
        • If this order is "specified" then we could extract parts of lambda. Patch for Inf checks in qp output lambda will make things easier but is not critical.

    by Asma Afzal ( at May 23, 2015 05:40 AM

    May 20, 2015

    Mike Miller

    Birthday Resolutions - Review

    Last year on my birthday I decided to try setting some goals for self-improvement with a deadline of the following birthday. People typically set New Year’s resolutions for themselves, but I wanted to try something different. Partly because I’m a natural contrarian, but also because my birthday last year was unique and offered more than a few reasons for self-reflection. So with another birthday looming, it’s time now to review how this experiment worked out.

    First, because this was a particularly notable birthday I had decided to hold myself to 10 resolutions. So this was almost doomed to failure from the beginning, if success means hitting all 10, which I didn’t. If I want to do this again next year, I should definitely go with a smaller set of goals to better set myself up for success. Obvious.

    Some of my goals were broader than others, which made it harder to define a successful target to aim for. For example, my goal to attend more free software developer conferences (did) was a lot easier to define and complete than my goal to make more time for creative pursuits (didn’t).

    Despite these problems, I like how this experiment turned out. I was able to accomplish about half of my 10 goals (for some definition of “accomplish”), and I’m not one to dwell on the other half that didn’t get done. I also like pinning personal goals to my birthday, rather than arbitrarily to the start of the Western calendar year. It reminds me to not only celebrate my birthday but to keep trying to improve from year to year.

    Anyone else tried this? Any other non-traditional ideas for annual, or more frequent, resolutions and personal goals?

    May 20, 2015 11:33 PM

    May 15, 2015

    Piotr Held

    Progress report

    As I haven't had any significant roadblock or breakthroughs this week I wanted to give a little progress report on my work.
    1. Added functionality
    I have managed to add the following functions:

    • mutual
    • spectrum
    • lazy
    • delay
    • pca
    Along with their documentation, tests and a demo (for lazy). I was really happy that once I had produced some examples of how I want to port these functions the process of porting each one accelerated rapidly. 
    I am especially excited about the fact that I have now henon, delay, an equivalent of addnoise and project available as this allowed me to create a nice noise reduction demo for project (and for lazy, but the one for project is more impressive). Fig.1 sums those efforts up.
    Fig. 1 Noisy data and data cleaned up by project.
    2. Functions found to be non-equivalent
    I also spent a lot of my time (almost a week) researching which programs from TISEAN have a GNU Octave equivalent. Apart from my positive identifications, most of which I discussed in previous posts, I have made some negative ones. I found that both extrema and polynom have no GNU Octave equivalent. 
    There was a suggestion made that extrema might be similar to findpeaks from signal. The only problem is that findpeaks searches (and returns) all peaks, whereas extrema returns either minimums or maximums. It might be easier to implement it in Octave than to port, but this decision has not been made yet.
    The latter program polynom was compared to detrend, polyfit and wpolyfit. The results were disappointing. All of the GNU Octave functions attempt to fit a polynomial onto the data, whereas polynom tries to make a "polynomial ansatz" for the data. The results are vastly different as can be seen on Fig. 2.
    Fig. 2 Comparison of original data (green), polyfit fit (red), and polynom prediction (blue).
    Both programs were run to try to use a 4th order polynomial. 

    by Piotr Held ( at May 15, 2015 06:18 PM

    May 14, 2015

    Asma Afzal

    Nonlinear Regression and 'nlinfit'

    In MATLAB, all three fucntions 'lsqnonlin', 'lsqcurvefit' and 'nlinfit' are used to perform non-linear curve fitting.

    To better understand the differences and similarities in these functions, consider the model function:
    $y= \beta_1+\beta_2  \text{exp}(-\beta_3x)$

    We wish to estimate the $\beta=\{\beta_1,\beta_2,\beta_3\}$ for the set of independents {$x_i$} and observed values {$y_i$} such that the model fits the data.

    Both 'nlinfit' and 'lsqcurvefit' are very similar as we can pass the regression function to compute the parameters. 'lsqnonlin' on the other hand, solves optimization problems of the type $min_{\beta} \sum_k f_k(\beta)^2$, so we cannot directly specify the regression function and instead, an error function has to be provided.  This is shown in the code below:

    modelfun = @(b,x)(b(1)+b(2)*exp(-b(3)*x));
    b = [1;3;2]; %actual
    x = exprnd(2,100,1); %independents
    y = modelfun(b,x) + normrnd(0,0.1,100,1); %noisy observation
    beta0 = [2;2;2]; %guess
    beta = nlinfit(x,y,modelfun,beta0)
    beta = lsqcurvefit(modelfun,beta0,x,y)
    beta = lsqnonlin(@(b)err_fun(b,x,y),beta0) %err_fun = modelfun-y

    All three functions generate:

    beta =


    • lsqcurvefit is more superior in the sense that we can define the bounds for the design variable (unlike nlinfit) while inputting the observed values separately (unlike lsqnonlin). 
    • Nlinfit provides extra statistics such as covariance matrix of the fitted coefficients and information about error model.
    • As an alternative to defining weights for the observed values in 'nlinfit', 'RobustWgtFtn' option can choose from different pre-defined weight functions for robust regression (with robust regression, fitting criterion is not as vulnerable to unusual data as least squares weighting function.)


    by Asma Afzal ( at May 14, 2015 10:24 AM

    May 12, 2015

    Mike Miller

    Octave + Python: A New Hope

    As a fan of both Python and Octave for numerical computing, and an active Octave developer, I’m always excited to hear about projects in either environment that create new capabilities or open up new ways of looking for solutions to problems. So I am especially excited about a new project that has the potential to bring Octave and Python much closer together and to give users of either tool full use of the other.

    The broad goal of this project is to provide a two-way interface layer between Octave and Python. What does this mean specifically? Well, I expect a future version of Octave to have a function that will call Python functions, using an embedded Python runtime, with transparent conversion between native Octave types and Python / NumPy types. There will also be a Python module to do the inverse: allow Python code to call Octave functions, invoke an embedded Octave interpreter, and have automatic conversion between Python and Octave types.

    The way in which the seeds of this project came together very quickly is really interesting, and what I want to describe in this post. The first was in a mailing list side discussion in late March about the appropriateness of Octave and Matlab for teaching numerical programming. It was mentioned that recent versions of Matlab have a calling interface to Python. For years they had provided a similar interface to Java, but I had no idea that Python was now an option for Matlab users. I filed that away for later.

    Then there is the Octave symbolic package, which relies heavily on SymPy to do the actual symbolic computation, but interacts with Python and SymPy over a pipe. So that existing package would definitely benefit from having a Python interpreter embedded in Octave or in a loadable oct-file.

    And finally there was a post in early April from fellow Octave developer JordiGH, who wrote:

    I have a wild idea. I like Python, and I think Numpy and Scipy are a great tool. Interfacing Scipy with Octave is also a good thing. … I therefore propose to bring Pytave into Octave proper.

    Pytave is an already-existing project which provides a Python module that can call Octave functions. It worked with older versions of Octave years ago, but has not kept up with the Octave API. It did work, and it does have a lot of useful code for converting between Octave and Python types, lots of good groundwork to start building from.

    I’m not sure what led Jordi to think of this “wild idea” or share it with us, but it definitely inspired me to latch onto this project. The timing of his message, after the other previous uses and mentions of Python, and being just days before the start of my first PyCon experience, read to me like a call to action. This felt like a perfect confluence of events and ideas to bring Octave and Python together in a novel way.

    So, I have already put some effort into this, and am planning to do some more. I hope that I (and any other interested contributors) will be able to make some real progress on this Octave-Python interface during this summer. I will share some more specifics about the project in a followup post soon.

    Thoughts about this project? Interested in following our progress or contributing?

    May 12, 2015 03:41 AM

    May 07, 2015

    Piotr Held

    The problem with 'spectrum'

    The 'spectrum' function from TISEAN most likely needs to be rewritten in GNU Octave or there is no need for it. This is because linking to it does not seem like a good idea. This is because there is a suspicion it does not produce good results for some data inputs. 

    1. Where 'spectrum' works
    First it is important to note that 'spectrum' from TISEAN is basically a GNU Octave 'abs(fft(:))' with additional data manipulation/adjustment. This additional work is not an elegant one-line solution, which might warrant designating a separate function that would translate the Octave respective function into a form similar to the output of 'spectrum'. Although this might not be necessary since the data obtained from the Octave function is very similar to 'spectrum' (Fig. 1).
    Fig. 1 Unadjusted data from Octave
    After adjusting the data (which was done by analyzing the source code to determine what actions the TISEAN programs perform) it was possible to get a close fit with a small difference. An example of this type of adjustment is listed below (Fig. 2)
    Fig. 2 Adjusted data from Octave.

    As it is not a one-line fix to convert ' abs (fft (:))' into a similar format as 'spectrum' it will not be shown in the post. It is available in the 'tests/spectrum/test_spectrum.m' function located on the tisean package repo (here).

    2. Where the problem lies
    The problem is that when 'spectrum' is used to create a step response its results vary substantially from what is produced by Octave. The way the data looks suggests that there is something wrong with 'spectrum'. The adjusted version is situated below.
    For the most part, the data fits perfectly, but there seems to be a shadow on the bottom of the TISEAN data. If it is the case that there is a problem with 'spectrum' then its code should not be used in the future Octave package and should be rewritten or omitted (as similar results can be obtained from a simple Octave call).

    by Piotr Held ( at May 07, 2015 07:58 PM

    May 06, 2015

    Piotr Held

    Finding 'histogram' in GNU Octave

    Unlike 'corr' it is quite easy to find a representative for 'histogram' from TISEAN. It is 'hist' from GNU Octave. The data is almost the same with the exception that the TISEAN package normalizes by default so one needs to be careful when calling the respective functions. I will describe differences in the data and describe the differences in usage.

    1. Data comparison
    I have attached a comparison of the two data sets (from 'hist' and 'histogram' on one chart)
    Fig. 1 Comparison between 'hist' (Octave) and 'histogram' (TISEAN)
    When  one analyses the data close there is a slight discrepancy between the value on the 40th and 41st bar. But not only is it slight, it basically means that both programs assigned a certain value to two different bins, which should not be a major problem. All told we can say that both of those functions perform the same task.

    2. Usage comparison
    As mentioned before, usage varies on both functions.
     $ histogram amplitude.dat -b#n -o "amplitude_histogram.dat"
     [nn, xx] = hist (amplitude, #n, 1);  
    nn = transpose (nn); xx = transpose (xx)  
    amplitude_hist = [xx, nn];
    This way the data stored in 'amplitude_hist' is essentially the same with 'amplitude_histogram.dat'.

    by Piotr Held ( at May 06, 2015 10:46 AM

    Finding a 'corr' representative in Octave

    This article describes the methodology used to compare function from GNU Octave and the TISEAN package. To achieve the desired results the author assumes you have installed the TISEAN package (available here) and have downloaded amplitude.dat and have installed GNU Octave with the 'signal' package in version 1.3.0 or newer.

    1. Comparison

    Procedure taken to receive results:

      1. Generate amp_corr.dat using the TISEAN package 'corr' with the call:
     '$ corr amplitude.dat -D5000 -o "amp_corr.dat"'
      2. Generate similar autocorelation data using (in GNU Octave):
     'load amplitude.dat; [a,b] = xcorr(amplitude, 5000, 'coeff');'
         Then to save the data you can use:
     'idx = [rows(amplitude):2*rows(amplitude)-1]; 
      xcorr_res = a(idx);
      save "xcorr_res.dat" xcorr_res'
    There is a strong difference in the data. This might be because of the different methods used in both cases (as explained further in the methods used Section 2. Methods). Because of those differences the amplitudes of the data generated using 'xcorr' from 'signal' decreases linearly. Thus to compare the data the oscillation amplitude of the data generated by 'xcorr' must be amplified. This linear decrease was not proven but observed on the 'amplitude.dat' data.

    When a linear correction is applied:
     'mult = rows (amplitude) ./ (rows (amplitude) - [0:rows(amplitude)-1]);  
    xcorr_tisean_res = mult .* xcorr_res'
    Fig. 1 Difference between xcorr_tisean_res and amp_corr

    The resultant xcorr_tisean_res is close to the TISEAN 'corr' function, and the difference is smaller than 3% (see Fig. 1). The end of the data begins to change and this is most likely because there is no more data past 5000 and so the results vary. If a autocorrelation is calculated for less data (e.g.4500 instead of 5000) the difference is much less, as can be seen on the chart above.

    Even better results can be obtained for different data. We can generate a different set using the TISEAN package
     '$ar-model -s5000 amplitude.dat -p10 -o "amp_ar.dat"' 
    When the process described above is applied to this new data set ('amp_ar.dat') the resulting difference between 'xcorr' and 'corr' is shown on Fig. 2.
    Fig. 2 Difference between 'xcorr' and 'corr' on 'amp_ar.dat'

    Similarly to the previous case the data is the same for small ( < 4000) numbers but when they get close to the edge the difference becomes more pronounced.

    2. Methods

    The way TISEAN calculates autocorrelation in the 'corr' program is by using estimation method. It is described here:

    On the other hand the 'xcorr' function from the signal package uses the FFT (Fast Fourier Transform) method (it is described in the same Wikipedia article: here)

    This difference in methodology is the cause of the difference in the data results between both functions.

    3. Conclusions [edited]
    After more test we found 'corr' from TISEAN and 'xcorr' from 'signal' to perform the same autocorrelation and therefor it is not necessary to port it.

    It is important to note the different usage:
     $ corr amplitude.dat -Dn# -n -o "amplitude_tisean.dat" 
     [data, lags] xcorr (amplitude, n#, 'unbiased')   
    data = data(find (lags >0))
    Both of the usage noted above produce the same data.

    It is important to note the '-n' in the calling of the TISEAN program. It mean the data is not normalized. You can achieve similar data even when calling 'corr' with normalization, but it is more tricky:
     $ corr amplitude.dat -Dn# -o "amplitude_tisean.dat" 
     [data, lags] xcorr (center (amplitude), n#, 'coeff')   
    data = data(find (lags >0))  
    data = data .* (n# ./ (n# - (transpose ([0:n#-1]))))
    The results of this can be viewed in tests/corr with function test_corr.m (note: 'signal' package needed) available in from the tisean port package repo:

    by Piotr Held ( at May 06, 2015 09:58 AM

    May 05, 2015

    Asma Afzal

    Project Goals

    Here, I list down the project goals as stated in my Wiki page:

              Start of GSoC (May) 
    1. 'lsqnonlin' using 'nonlin_residmin'
    2. 'lsqcurvefit' using 'nonlin_curvefit', 'nonlin_residmin', or 'lsqnonlin',
    3. 'fmincon' using 'nonlin_min',
    4. 'nlinfit' using 'leasqr',
    5. Test cases for the above functions [10] .
    6. Instead of wrappers for top-level functions like qp, call back-end function (__qp__) to be able to extract lambda. See [11].
      Stretch Goals
    7.  Further missing functions in Optim package. See [12] Implement another back-end algorithm/add variant.

     6.  quadprog and lsqlin should call a private intermediate function instead of qp.m
     This private function should do the argument processing for calling __qp__. It could also be configured to call yet to be written alternatives to __qp__ .
    Among other things this should make ordering of the 'lambda' output feasible.

    ((I have yet to study __qp__ and how this will be done))

    by Asma Afzal ( at May 05, 2015 05:41 PM

    GSoC Acceptance :)

    I got selected to work on Octave GSoC! ( Yeaa! :) )

    Here is the link to my Wiki application:

    The project officially starts on 25th May, 2015. During this community bonding period, I will try to get a clear breakdown of the project goals with the help of my mentors.

    by Asma Afzal ( at May 05, 2015 05:04 PM

    April 29, 2015

    Piotr Held

    Estimated timeline

    Most of my discussion is based on my wiki post in which I delve into details how I want to approach the problem. This post is available here. I am aware that it might change as I go along but for now I do not have any new information which could help me estimate how much time each task takes.

    Taking under consideration the aforementioned wiki page I estimate that I will do as follow:
    1. For the next couple of weeks (May 6-17) I will try to test the Octave functions and compare them to the ones from the TISEAN package. Also I intend to do other tasks to facilitate my work later (e.g. fix the makefile)
    2. Next I will start porting according to the schedule outlined in the wiki. I would like to finish and clean up 'Dimensions and Entropies' before the midterm assessment.
    3. In the following 4 weeks or so I would like to finish the project up, obviously this might not be possible and that is why I leave myself a time 'cushion' if my estimates are mistaken.
    The biggest problem I foresee in the 1st stage is that FORTRAN programs use floats and Octave uses double, so the results might be different. This is a problem, especially in cases when data relies on some previously calculated data. I have tried different numbers and so far I saw that for 'henon' the difference between the result has been < 10^-3 (or <10^-2 later on) so I assumed they are the same. As of writing this I do not know how to elegantly compare the two data sets so as to compensate for the difference that results from floats and doubles.

    by Piotr Held ( at April 29, 2015 01:12 PM

    April 28, 2015

    Piotr Held

    First Post

    As I am not fully familiar with blogging I am posting this post only to educate myself in how it works.
    The next post I post will have a timeline.

    by Piotr Held ( at April 28, 2015 02:42 PM

    April 21, 2015

    Juan Pablo Carbajal

    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 April 21, 2015 02:47 PM

    April 11, 2015

    Juan Pablo Carbajal

    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 April 11, 2015 09:35 AM

    March 13, 2015

    Juan Pablo Carbajal

    ¡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 March 13, 2015 12:40 AM

    March 10, 2015

    Juan Pablo Carbajal

    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 March 10, 2015 09:34 PM

    March 09, 2015

    Juan Pablo Carbajal

    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 March 09, 2015 05:22 PM

    March 04, 2015

    Juan Pablo Carbajal

    ¡No soy cirujano pero necesito operar!

    Así es, necesitamos operar, pero no en el sentido médico. Necesitamos operar con números.
    En el ámbito de la programación llamamos operadores a un grupo de funciones que tienen una sintaxis particular (es decir que se escriben de una forma especial). Estas funciones generalmente coinciden con las bien conocidas operaciones de la aritmética:

     Operación Operador en Octave 
    Suma  +
    Resta  -
    Multiplicación  *
    División  / 

    Generalmente, cuando queremos escribir "Calcule la suma entre el número 3 y 4" anotamos 3 + 4 (3 más 4). Pero podríamos escribir +(3,4) ("sume 3 y 4") ó + 3 4 (eso se llama notación polaca) o incluso 3 4 + ("3 y 4, ahora sume", notación polaca inversa). Lo importante es que la operación requiere dos entradas, el 3 y el 4, y el símbolo que representa la operación es el operador. Hay operadores binarios, como el que acabamos de ver, pero también los hay unitarios (solo una entrada) y ternarios (tres entradas)... y así sucesivamente.
    La notación preferida, por ejemplo 3 + 4 para la suma, se llama notación de infijo (como prefijo, y sufijo) y el operador va entre medio de las entradas. Se darán cuenta que esta notación es muy clara para operadores binarios (con dos entradas).

    En la teoría de hoy discutiremos otros operadores ademas de los aritméticos, que en Octave pueden resultar muy útiles: especilamente cuando queremos hacer operaciones entre vectores o listas de números. No voy explicar todos los operadores! Les explicaré unos cuantos, al resto lo pueden investigar y aprender por su cuenta.


    Sumar dos vectores o listas de números

    Supongamos que hemos colectado la duración, en minutos, del 1er tiempo de varios eventos deportivos en un vector fila llamado primerT

    primerT = [20 45 30 30 46 22];

    En otro vector hemos colectado la duración del 2do tiempo de los mismos eventos

    segundoT = [27 49 35 30 45 20];

    Y ahora nos gustaría saber la duración total del partido. Tenemos que sumar la duración de los dos tiempos. en Octave hacemos

    TotalT = primerT + segundoT

    Como pueden ver el operador suma se aplica elemento a elemento.

    Multiplicar dos listas

    Supongamos que en otra ocasión hemos colectado los precios del pasaje de colectivo de diferentes años

    precio = [0.7, 1.2, 2.5, 4, 10];

    Y en otro vector hemos colectado la cantidad de veces que viajamos en colectivo y pagado el precio correspondiente

    veces = [10 8 5 5 2];

    Ahora queremos calcular la plata que gastamos en cada situación. Tenemos que multiplicar los dos vectores. En Octave, si ingresamos veces * precio recibimos un error. 

    Antes de seguir leyendo ¿Podés explicar el porqué?

    El problema es que Octave considera a las dos entradas vectores, y entre vectores hay una operación llamada producto escalar (tcc. producto punto, producto interno) que en Octave se calcula usando el operador * (asterisco). Para indicarle a Octave que lo que queremos es un producto elemento a elemento utilizamos el operador .* :

    veces .* precio

    y obtenemos el resultado deseado.


    Para elevar un número a una potencia, por ejemplo 3 al cuadrado, escribimos

    3 ^ 2

    el operador ^ representa la operación.

    Para calcular el cuadrado de una lista de números operamos con .^


    En general todo operador cuyo nombre empiece con un punto (.) realiza la operación elemento a elemento.

    Pregunta 1:
    ¿Qué operador utilizarías para dividir dos listas de números?


    Volvamos a las listas con duración de los tiempos de un deporte. Ahora vamos a organizar la información de otra manera que quizás es más fácil de visualizar.
    Primero ordenemos los deportes asociando a cada uno con un número:
    1. Football
    2. Handball
    3. Basketball
    Por cada deporte vamos a guardar la duración del 1er tiempo de diferentes eventos. A esta información la organizamos en forma de tabla:

       Football Handball  Basketball 
     1er Evento  45  20 15 
     2do Evento  49  20 18 
     3er Evento  46  25 18 
     4to Evento  45  22 15 

    Podemos ingresar esta tabla en Octave de la siguiente manera:

    primerT = [45 20 15; 49 20 18; 46 25 18; 45 22 15];

    Supongamos que hacemos lo mismo para el segundo tiempo

    segundoT = [49 21 15; 45 20 15; 46 20 15; 49 21 18];

    La duración total la seguimos calculando con primerT + segundoT. Noten que esta operación se ejecuta, nuevamente, elemento a elemento.

    Pregunta 2:
    ¿Cómo harías la multiplicación de estas dos tablas?

    Hay muchísimos más operadores para investigar ¡Disfruten!

    by Juan Pablo Carbajal ( at March 04, 2015 02:34 PM

    March 01, 2015

    Asma Afzal

    Hello World

    I am Asma Afzal, an Electrical Engineering PhD student at the University of Leeds, UK. I am applying for the GNU Octave GSoC project of creating MATLAB compatible wrappers for functions in the optim package. Additional stretch goals are still under discussion. These include, implementing additional optimization algorithms already present in C++/Fortran as they may be well-suited to some specific problems [1]. 


    by Asma Afzal ( at March 01, 2015 10:44 AM

    February 20, 2015

    Juan Pablo Carbajal

    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 February 20, 2015 07:07 PM

    December 24, 2014

    Juan Pablo Carbajal

    ¡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 December 24, 2014 04:16 PM

    December 21, 2014

    Juan Pablo Carbajal

    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 December 21, 2014 09:36 AM

    October 31, 2014

    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 October 31, 2014 02:26 AM

    October 21, 2014

    Jordi Gutiérrez Hermoso

    Customising Mercurial Like a Pro

    Consider the following Mercurial screenshot:


    This might look bewildering at first. It’s a custom hg wip command, a variation of hg log. You can see the DAG on the left, what looks like commit summaries in white, and what looks like usernames in mauve. This is not what stock Mercurial looks like, but in my opinion, this looks great. It has so much information packed into such a small space, all colour-coded for your convenience. Let’s take a look.

    Current commit


    The current commit is simply where you’re standing. It’s the commit that hg diff or hg status would use as reference for telling you how the working copy has changed. All that I did for hg wip was to highlight more clearly what this commit is. Note in the DAG that this commit is marked with an “@”.



    I have chosen to display bookmarks in green. Recall that bookmarks are simply movable markers that follow along with the commits that they’re pointing to. They’re very similar to git branches. Here I am using them to refer to each of my feature branches. There is a special “@” bookmark that indicates where upstream development is.



    In cyan we see the branch names. In this case there are only two commits on the stable branch. Branch names are permanent names tattooed on commits. This is useful for when you want to have a permanent record of where a commit was done.



    This part is one of my favourites. Mercurial phases indicate whether commits have been published yet or not, and whether they should be shared or not. Commits that have been published and therefore immutable are in orange. The commits I’ve been working on but haven’t been accepted into upstream yet are yellow drafts. Finally, a couple of local patches that should never be shared are secret commits in blue. I am using local revision numbers to identify commits, but it would also be possible to use global hashes.

    Selective about which commits to show

    Another cool thing about my hg wip command is that it only shows what I consider interesting information about my work-in-progress commits, hence the name of the command. To be precise, it only shows my draft and secret commits, current open heads of development, and the minimum common commits. No other intermediate commits are shown. Actually, usually I don’t even need the extra newlines, and I prefer this more compact format:


    So, how is this command built? It actually depends on a number of different Mercurial features, all useful on their own. All of the following snippets go into your ~/.hgrc config file.


    Revsets are a very useful language for selecting a subset of commits. Here you can see me talking about them. The revset I used for hg wip was

    wip = (parents(not public()) or not public() or . or head()) and (not obsolete() or unstable()^) and not closed()

    It looks like quite a mouthful, so let’s pick it apart:

    parents(not public())
    Take the parents of all the non-public commits (i.e. drafts and secrets).
    or not public() or . or head()
    Also take all non-public commits, the current commit (that’s “.”), and all heads.
    and (not obsolete() or unstable()^)
    But exclude obsolete and unstable commits (this only matters if you’re using the experimental Evolve extension).
    and not closed()
    And also exclude any closed branch heads

    The nice thing is that all of this complicated selection is also fast!

    Templating the output

    The next step is to show the output in the format that we want it. That’s what the templating engine is for. The template I used here was

    wip = '{label("log.branch", ifeq(branch, "default", "", branch))} {label("changeset.{phase}", rev)} {label("grep.user", author|user)}{label("log.tag", if(tags," {tags}"))} {bookmarks % "{ifeq(bookmark, currentbookmark, label('log.activebookmark', bookmark), label('log.bookmark', bookmark))} "}\n{label(ifcontains(rev, revset('parents()'), ''),desc|firstline)}'

    Whoa! That’s even worse than the revset. Unfortunately, due to the nature of the templating engine, inserting whitespace to make this more readable would also make this extra whitespace show up in the final output.

    If we are willing to define more intermediate templates and move them to an external template file, it can actually be made readable:

    wbrch = '{label("log.branch",
                     ifeq(branch, "default",
    wcset = '{label("changeset.{phase}",
    wuser = '{label("grep.user",
    wtags = '{label("log.tag",
                  if(tags," {tags}"))}'
    wbook = '{bookmarks % "{ifeq(bookmark, currentbookmark,
                                 label('log.activebookmark', bookmark),
                                 label('log.bookmark', bookmark))} "
    wdesc = '{label(ifcontains(rev, revset('parents()'), ''),
    changeset = '{wbrch} {wcset} {wuser} {wtags} {wbook}\n{wdesc}'

    Okay, this is a bit better, but it requires using an extra config file. The syntax is a bit tricky, but most of it is self-explanatory. Special keywords in {braces} get expanded to their corresponding values. There are a few functions like if() and ifcontains() to selectively output extra information if that information exists. For example, I use ifeq() to hide the branch name if this name is “default”.

    Once you’re in a function, keywords are already available. For example, rev gives you the revision number, and if you wanted global hashes instead, you could change that to node or shortest(node). You can even use a revset in the templating engine! In this case, for selecting the current commits and comparing it to the revision being printed. We use the parents() revset to select the current commit(s), because there may be more than one current commit if there’s an uncommitted merge state on the working directory. For template keywords that return a list of things, such as bookmark, we can iterate over elements of that list using “%” and then select different labels depending if the bookmark is the current bookmark or not.

    But what is this label function? Ah, well, labelling is how the colour extension decides how to colourise the output.

    Colourful Mercurial

    The finishing touches are now to pick colours for all of the labels we defined. For this we enable the colour extension and we configure it:



    #Custom colours

    #Colours for each label
    log.activebookmark = green bold underline

    changeset.public=orange bold
    changeset.secret=blue bold
    changeset.draft=brightyellow bold blue_background

    Here we set the mode to terminfo so that we can use 256 colours

    . Most modern terminals can support it, but sometimes you have to specify the $TERM enivronment variable to be xterm-256color or something similar. Then you can define your own colours from the 256 available, referring to the numbers on this chart. If you dislike my colour choices, this is where you can configure them.

    Define the command

    The final part is to just turn this on:

    wip = log --graph --rev=wip --template=wip

    This defines the hg wip command to be an alias to hg log together with the parameters for displaying the graph, using the wip revset, and the wip template.

    And voilà, fancy shmancy colourful hg command! Here is the complete addition to your ~/.hgrc all at once, for delicious copypasta:

    wip = (parents(not public()) or not public() or . or head()) and (not obsolete() or unstable()^) and not closed()

    wip = '{label("log.branch", ifeq(branch, "default", "", branch))} {label("changeset.{phase}", rev)} {label("grep.user", author|user)}{label("log.tag", if(tags," {tags}"))} {bookmarks % "{ifeq(bookmark, currentbookmark, label('log.activebookmark', bookmark), label('log.bookmark', bookmark))} "}\n{label(ifcontains(rev, revset('parents()'), ''),desc|firstline)}'



    #Custom colours

    #Colours for each label
    log.activebookmark = green bold underline

    changeset.public=orange bold
    changeset.secret=blue bold
    changeset.draft=brightyellow bold blue_background

    wip = log --graph --rev=wip --template=wip


    This command comes from ideas cobbled together from Steven Losh, Augie Fackler, and Sean Farley. They are all great contributors to Mercurial, and they have taught me so much! Thanks, guys!

    by Jordi at October 21, 2014 08:31 PM

    August 27, 2014

    Roberto Porcù

    Dense Output

    Dense Output

    Specifying specific output times for the solution, should not affect the internal time steps that the solver uses. The basic idea of the Dense Output concept is to provide the solution at a given time \(s \in [t, t+dt]\) with the same order of accuracy as the solutions computed at the internal time points by using suitable interpolation methods.
    Up to now only linear interpolation was performed and this significantly lowered the accuracy if a higher order solver was used.
    I then implemented a series of interpolation function:

    • linear_interpolation:
      x_out = linear_interpolation (t, x, t_out)
      Given the time span \(t=[t_0, t_1]\) and the function values \(x=[x_0, x_1]\), it returns the linear interpolation value \(x_{out}\) at the point \(t_{out}\).
    • quadratic_interpolation:
      x_out = quadratic_interpolation (t, x, der, t_out)
      Given the time span \(t=[t_0, t_1]\), the function values \(x=[x_0, x_1]\) and the derivative of the function at the point \(x_0\), it returns the quadratic interpolation value \(x_{out}\) at the point \(t_{out}\).
    • hermite_cubic_interpolation:
      x_out = hermite_cubic_interpolation (t, x, der, t_out)
      Given the time span \(t=[t_0, t_1]\), the function values \(x=[x_0, x_1]\) and the derivatives of the function at both points \(x_0\) and \(x_1\), it returns the 3rd order approximation \(x_{out}\) at the point \(t_{out}\) by performing Hermite interpolation.
    • hermite_quartic_interpolation:
      x_out = hermite_quartic_interpolation (t, x, der, t_out)
      Given the time span \(t=[t_0, t_1]\), the function values \(x=[x_0, x_{1/2}, x_1]\) (where \(x_{1/2}\) is the value of the function at the time \(t_0+dt/2\)) and the derivatives of the function at the extremes \(x0\) and \(x1\), it returns the 4th order approximation \(x_{out}\) at the point \(t_{out}\) by performing Hermite interpolation.
    • dorpri_interpolation:
      x_out = dorpri_interpolation (t, x, k, t_out)
      This interpolation method is specific for the Dormand-Prince Runge-Kutta scheme. Given the time span \(t=[t_0, t_1]\), the function value \(x=x_0\) and the vector \(k\) with the function evaluations required in the Dormand-Prince method, it returns the 4th order approximation \(x_{out}\) at the point \(t_{out}\). For more information on the method have a look at Hairer, Noersett, Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems (pag. 191-193).
    • hermite_quintic_interpolation:
      x_out = hermite_quintic_interpolation (t, x, der, t_out)
      Given the time span \(t=[t_0, t_1]\), the function values \(x=[x_0, x_{1/2}, x_1]\) and the derivatives of the function at each point, it returns the 5th order approximation \(x_{out}\) at the point \(t_{out}\) by performing Hermite interpolation.
    These methods are then used to perform the Dense Output according to the order of the solver chosen. This is the piece of code in integrate_adaptive.m that performs the interpolation:

    % if next tspan value is caught, update counter
    if( (z(end) == tspan(counter)) || (abs (z(end) - tspan(counter)) / ...
    (max (abs (z(end)), abs (tspan(counter)))) < 8*eps) )

    % if there is an element in time vector at which the solution is required
    % the program must compute this solution before going on with next steps
    elseif( vdirection*z(end) > vdirection*tspan(counter) )
    % initializing counter for the following cycle
    i = 2;
    while ( i <= length (z) )

    % if next tspan value is caught, update counter
    if( (counter <= k) && ...
    ( (z(i) == tspan(counter)) || (abs (z(i) - tspan(counter)) / ...
    (max (abs (z(i)), abs (tspan(counter)))) < 8*eps)) )

    % else, loop until there are requested values inside this subinterval
    while((counter <= k) && (vdirection*z(i) > vdirection*tspan(counter)))
    % choose interpolation scheme according to order of the solver
    switch order
    case 1
    u_interp = linear_interpolation ([z(i-1) z(i)], [u(:,i-1) u(:,i)], tspan(counter));
    case 2
    if (~isempty (k_vals))
    der = k_vals(1);
    der = feval (func, z(i-1) , u(:,i-1), options.vfunarguments{:});
    u_interp = quadratic_interpolation ([z(i-1) z(i)], [u(:,i-1) u(:,i)], der, tspan(counter));
    case 3
    % only ode23 - use k_vals
    u_interp = hermite_cubic_interpolation ([z(i-1) z(i)], [u(:,i-1) u(:,i)], [k_vals(:,1) k_vals(:,end)], tspan(counter));
    case 4
    % if ode45 used without local extrapolation this function doesn't require a new function evaluation.
    u_interp = dorpri_interpolation ([z(i-1) z(i)], [u(:,i-1) u(:,i)], k_vals, tspan(counter));
    case 5
    % ode45 with Dormand-Prince scheme:
    % 4th order approximation of y in t+dt/2 as proposed by Shampine in Lawrence, Shampine, "Some Practical Runge-Kutta Formulas", 1986.
    u_half = u(:,i-1) + 1/2*dt*((6025192743/30085553152)*k_vals(:,1) + (51252292925/65400821598)*k_vals(:,3) - (2691868925/45128329728)*k_vals(:,4) + (187940372067/1594534317056)*k_vals(:,5) - (1776094331/19743644256)*k_vals(:,6) + (11237099/235043384)*k_vals(:,7));
    u_interp = hermite_quartic_interpolation ([z(i-1) z(i)], [u(:,i-1) u_half u(:,i)], [k_vals(:,1) k_vals(:,end)], tspan(counter));

    % it is also possible to do a new function evaluation and the quintic hermite interpolator
    %f_half = feval (func, t+1/2*dt, u_half, options.vfunarguments{:});
    %u_interp = hermite_quintic_interpolation ([z(i-1) z(i)], [u(:,i-1) u_half u(:,i)], [k_vals(:,1) f_half k_vals(:,end)], tspan(counter));
    warning ('high order interpolation not yet implemented: using cubic iterpolation instead');
    der(:,1) = feval (func, z(i-1) , u(:,i-1), options.vfunarguments{:});
    der(:,2) = feval (func, z(i) , u(:,i), options.vfunarguments{:});
    u_interp = hermite_cubic_interpolation ([z(i-1) z(i)], [u(:,i-1) u(:,i)], der, tspan(counter));

    % add the interpolated value of the solution
    u = [u(:,1:i-1), u_interp, u(:,i:end)];

    % add the time requested
    z = [z(1:i-1);tspan(counter);z(i:end)];

    % update counters

    % if new time requested is not out of this interval
    if ((counter <= k) && (vdirection*z(end) > vdirection*tspan(counter)))
    % update the counter
    % else, stop the cycle and go on with the next iteration
    i = length (z) + 1;



    It is important to notice that:

    • The 1st order approximation doesn't require any additional function evaluation.
    • The 2nd order approximation may require the evaluation of the function at the current time. This can be avoided if the stepper already returns that value.
    • The only 3rd order solver implemented is ode23. The 3rd order approximation exploits the Runge-Kutta \(k\) values to avoid further function evaluations.
    • There are no 4th order schemes as yet implemented. However if ones were to use ode45 without local extrapolation then the dorpri_interpolation function can be used to obtain a 4th order approximation without any additional function evaluation. For any other 4th order scheme the hermite_quartic_interpolation function can be used.
    • For the 5th order method ode45, Shampine proposes to obtain a 4th order approximation at the middle point and to use quartic interpolation. It is however possible to directly do quintic interpolation but this require an additional function evaluation without (according to Shampine) a significant improvement.
    • For the higher order solvers (ode78), a suitable interpolator has not yet been implemented.
    Finally, since I wrote the interpolation functions in such a way that they are independent of the number of output points requested, a possible improvement would be to compute all the values of the dense output between \(t\) and \(t+dt\) all at once instead of one value at a time.

      by Jacopo Corno ( at August 27, 2014 03:00 AM

      August 18, 2014

      Eduardo Fernández

      Firm pencils down.

      During this week I have been reorganizing all the code, docs and tests in a better way for integrating into Octave. As Rik kindly suggested, I decided to organize things this way:

      • Inside libinterp/dldfcn directory I have created two files, and
      •  Within those files there are the dld functions that implements the each of the algorithms. They are ment to be built-in functions and follows the naming convention.

                    * contains __ilu0__() , __iluc__() and __ilutp__()

                    * contains __ichol0__() and __icholt__().

      • I have moved all the tests from .cc files to .m scripts so no tests are performed for built-in functions.

       The code is ready to be pulled from my repo to be reviewed :

      Practical examples: preconditioned conjugate gradient

      It is interesting to show how preconditioning techniques can improve the convergency of some iterative solvers. In this case I am running a Matlab example using the Poisson matrix (that is positive definite) obtained with gallery() function. The scritp:

      1. In this first case the convergency of pcg using ICHOL(0) algorithm, modified ICHOL(0) algorithm and without preconditioning are compared.

      N = 100;
      A = gallery ('Poisson', N);
      b = ones (size (A, 1), 1);
      tol = 1e-6; maxit = 100;
      [x0, fl0, rr0, it0, rv0] = pcg (A, b, tol, maxit);
      L1 = ichol (A);
      [x1, fl1, rr1, it1, rv1] = pcg(A, b, tol, maxit, L1, L1');
      opts.type = 'nofill'; opts.michol = 'on';
      L2 = ichol (A, opts);
      e = ones (size (A, 2), 1);
      norm(A * e - L2 * (L2' * e))
      [x2, fl2, rr2, it2, rv2] = pcg (A, b, tol, maxit, L2, L2');
      semilogy (0:maxit, rv0 ./ norm (b), 'b.');
      hold on;
      semilogy (0:it1, rv1 ./ norm(b), 'r.');
      semilogy (0:it2, rv2 ./ norm(b), 'k.');
      xlabel ('iterations');
      ylabel ('error');
      legend ('No Preconditioner', 'IC(0)', 'MIC(0)'); 



      2. In this second part of the script what is compared is the      convergency of ICHOLT algorithm with different values of droptol.
      L3 = ichol(A, struct('type', 'ict', 'droptol', 1e-1));
      [x3, fl3, rr3, it3, rv3] = pcg (A, b, tol, maxit, L3, L3');
      L4 = ichol (A, struct ('type', 'ict', 'droptol', 1e-2));
      [x4, fl4, rr4, it4, rv4] = pcg (A, b, tol, maxit, L4, L4');
      L5 = ichol (A, struct ('type', 'ict', 'droptol', 1e-3));
      [x5, fl5, rr5, it5, rv5] = pcg (A, b, tol, maxit, L5, L5');
      figure; semilogy (0:maxit, rv0 ./ norm (b), 'b-', 'linewidth', 2);
      hold on;
      semilogy (0:it3, rv3 ./ norm(b), 'b-.', 'linewidth', 2);
      semilogy (0:it4, rv4 ./ norm(b), 'b--', 'linewidth', 2);
      semilogy (0:it5, rv5 ./ norm(b), 'b:', 'linewidth', 2);
      ylabel ('error');
      xlabel ('iterations');
      legend ('No Preconditioner', 'ICT(1e-1)', 'ICT(1e-2)', ...
         'ICT(1e-3)', 'Location', 'SouthEast');


      As it can be seen Octave plots are the same as Matlab's ones.  Both lead to a decrease in the number of steps upt to convergence of the pcg method. ILU algorithms could also have been used here, but due to the simetry of the problem matrix ICHOL is faster.



      by Eduardo ( at August 18, 2014 08:19 PM

      Roberto Porcù

      FSAL - new stepper implementation

      As stated in the previous post, the implementation of the steppers as it was did not allow the possibility to exploit the FSAL (First Same As Last) property of the Bogacki-Shampine algorithm (ode23) and of the Dormand-Prince algorithm (ode45).
      The input and output arguments of the steppers have then be modified. As an example here is the runge_kutta_23 stepper:

      function varargout = runge_kutta_23 (f, t, x, dt, varargin)
      options = varargin{1};
      k = zeros (size (x, 1), 4);

      if (nargin == 5) % only the options are passed
      k(:,1) = feval (f, t , x, options.vfunarguments{:});
      elseif (nargin == 6) % both the options and the k values are passed
      k(:,1) = varargin{2}(:,end); % FSAL property
      k(:,2) = feval (f, t + (1/2)*dt, x + dt*(1/2)*k(:,1), options.vfunarguments{:});
      k(:,3) = feval (f, t + (3/4)*dt, x + dt*(3/4)*k(:,2), options.vfunarguments{:});

      %# computing new time and new values for the unkwnowns
      varargout{1} = t + dt; %t_next
      varargout{2} = x + dt.*((2/9)*k(:,1) + (1/3)*k(:,2) + (4/9)*k(:,3)); % return the 3rd order approximation x_next

      %# if the estimation of the error is required
      if (nargout >= 3)
      %# new solution to be compared with the previous one
      k(:,4) = feval (f, t + dt, varargout{2}, options.vfunarguments{:});
      varargout{3} = x + dt.*((7/24)*k(:,1) + (1/4)*k(:,2) + (1/3)*k(:,3) + (1/8)*k(:,4)); %x_est
      varargout{4} = k;


      And the call within the solver becomes:
      [s, y, y_est, k_vals] = stepper (func, z(end), u(:,end), dt, options, k_vals);

      where k_vals has to be initialized for the first iteration as f(t, x).
      This implementation will reduce the number of function evaluation for each step.

      Furthermore, after some tests in MATLAB, the return values for the solution and the estimate in the  runge_kutta_23 and runge_kutta_45 steppers have been swapped to automatically perform local extrapolation. The MATLAB functions are in fact of order 3 and 5 respectively.

      by Jacopo Corno ( at August 18, 2014 07:52 AM

      August 13, 2014

      Roberto Porcù

      Status of the code: bugfixes and new issues


      • odeset and odeget functions have been slightly modified to be compliant with MATLAB. Each MATLAB option is present and all the options are tested. The coding style has been adapted to the GNU-Octave standard.
      • ode_struct_value_check: this function has been introduced by Roberto in addition to odepkg_structue_check. The relation between the two functions has to be clarified: in particular it is necessary to understand if it is really necessary to have two different functions or one is sufficient.


      • The runge_kutta_78 stepper has been implemented.
      • Two 4th order steppers have been implemented: runge_kutta_45_dopri (Dormand-Prince coefficients) and runge_kutta_45_fehlberg (Fehlberg coefficients).


      • ode78 solver has been updated to the new structure. It now exploits the runge_kutta_78 stepper.
      • A series of tests has been added to each solver to check all the functionalities and the all options. This has made me possible to detect some bugs that have been corrected. In particular the adaptive timestep evaluation had some issues that lead to the use of too short timesteps. This has been corrected and now the algorithm proposed in [1] is used.
      • Furthermore the current implementation uses linear interpolation to evaluate the solution at the user specified times. This leads to a considerable loss in accuracy and is not consistent with MATLAB (which guarantees the same order of accuracy of the scheme used). In [1] different methods are proposed for the dense output: these will be used as a reference for the implementation of a better solution.
      • In the released version of odepkg some of the solvers perform local extrapolation, that is the higher-order estimate is chosen as the solution. With the new stepper structure, as it is now, this choice effects all the solvers. It have to be decided whether to perform it or not (MATLAB doesn't seem to do it, thus I suggest to avoid it).
      • MATLAB implementation of ode45 uses the Dormand-Prince (DP) coefficients. In the released version of odepkg there exits two solvers: ode45 that uses the Fehlberg coefficients and ode54 that uses the DP coefficients. To be consistent with MATLAB, ode45 now uses the DP method. This makes the runge_kutta_45_fehlberg stepper and the ode54 solver, as it is now, redundant. Either their elimination or a change of the solver might be considered. However one of the advantages of DP coefficients is the possibility to reuse the last function evaluation at a given step as the first evaluation of the subsequent one. This is not easily done with the stepper structure introduced by Roberto.


      • InitialStep option has been modified to be MATLAB compatible (it must be a positive scalar).
      • RelTol defalut value has been changed to 1e-3 instead of 1e-6 to be MATLAB compatible.
      • MaxStep option has been implemented.
      • NormControl option has been implemented.


      In addition to the general plan, a couple of issues need to be addressed:

      • Clarify the relation between ode_struct_value_check and odepkg_structue_check.
      • Decide if local extrapolation has to be used or not. My opinion (and the current implementation) is to avoid it to be compliant to what MATLAB seems to be doing.
      • Solve the dense output problem in a way that guarantees the consistency with MATLAB.
      • Consider if it's possible to reduce the number of function evaluation for the Dormand-Prince stepper (ode45) and the Bogacki-Shampine stepper (ode23) exploiting the FSAL property (first same as last).
      • Decide if in the future releases of odepkg ode54 has to be removed or maybe changed to become a Fehlberg solver.

      [1] E. Hairer, S.P. N{\o}rsett, G. Wanner, Solving Ordinary Differential Equations, 1993, Springer.

      by Jacopo Corno ( at August 13, 2014 02:28 AM

      August 12, 2014

      Eugenio Gianniti


      As said in my previous post, a missing feature in fem-fenics was the marking of subdomains. Indeed, I proposed an example that needed a file generated with a run of the corresponding Python code, which is not, honestly, the best approach. In order to address this issue, these days I have implemented a new class, subdomain, which can be used to mark mesh entities. In the following I will describe how to use this new functionality. Here is the code:

      pkg load fem-fenics msh

      ufl start Subdomains
      ufl fe = FiniteElement "(""CG"", triangle, 2)"
      ufl u = TrialFunction (fe)
      ufl v = TestFunction (fe)
      ufl a0 = Coefficient (fe)
      ufl a1 = Coefficient (fe)
      ufl g_L = Coefficient (fe)
      ufl g_R = Coefficient (fe)
      ufl f = Coefficient (fe)
      ufl a = "inner(a0*grad(u), grad(v))*dx(0) + inner(a1*grad(u), grad(v))*dx(1)"
      ufl L = g_L*v*ds(1) + g_R*v*ds(3) + f*v*dx(0) + f*v*dx(1)
      ufl end

      # Create mesh and define function space
      x = y = linspace (0, 1, 65);
      [msh, facets] = Mesh (msh2m_structured_mesh (x, y, 0, 4:-1:1));

      V = FunctionSpace ("Subdomains", msh);

      # Define boundary conditions
      bc1 = DirichletBC (V, @(x, y) 5.0, facets, 2);
      bc2 = DirichletBC (V, @(x, y) 0.0, facets, 4);

      # Define problem coefficients
      a0 = Constant ("a0", 1.0);
      a1 = Constant ("a1", 0.01);
      g_L = Expression ("g_L", @(x, y) - 10*exp(- (y - 0.5) ^ 2));
      g_R = Constant ("g_R", 1.0);
      f = Constant ("f", 1.0);

      # Define subdomains - Here are the edits #
      domains = MeshFunction ("dx", msh, 2, 0);
      obstacle = SubDomain (@(x,y) (y >= 0.5) && (y <= 0.7) && ...
                            (x >= 0.2) && (x <= 1.0), false);
      domains = mark (obstacle, domains, 1);

      # Define variational form
      a = BilinearForm ("Subdomains", V, V, a0, a1, domains);
      L = LinearForm ("Subdomains", V, g_L, g_R, f, facets, domains);

      # Assemble system
      [A, b] = assemble_system (a, L, bc1, bc2);
      sol = A \ b;
      u = Function ("u", V, sol);

      # Save solution in VTK format
      save (u, "subdomains");

      # Plot solution
      [X, Y] = meshgrid (x, y);
      U = u (X, Y);
      surf (X, Y, U);

      As you can see, it is basically the same as in the previous post, except the line used to import the meshfunction. I wrote in the corresponding comment where the edits are to be found. Now the workflow comprises these steps: first of all, a meshfunction needs to be created, then a subdomain, in the end we should mark cells. 

      The call to MeshFunction is something new, since it is now possible to instantiate a meshfunction given a mesh, the required topological dimension and the value to initialise it with. Moreover, the optional label "dx" means that it can be used in calls to BilinearForm or LinearForm to supply markers for subsets of the integration domains. In the example, this instruction returns a meshfunction of dimension 2, which means it holds values associated with each triangle in the mesh, initialised to be 0 in every entry.

      The subsequent instruction, instead, defines a subdomain, passing as arguments a function handle and a logical flag. The former will be the override of the dolfin::SubDomain::inside method, so it must return true for entities contained in the subset and false otherwise. In facts it checks whether the coordinates are inside the 2-interval defining the obstacle. The latter, instead, can be used to ask for a boundary subdomain, when set to true.

      At last, mark is called to set the entries corresponding to cells inside the subdomain to 1, so that the returned meshfunction now represents the obstacle: after these lines, the variable named domains assumes value 1 on cells inside the obstacle region and 0 outside. Thus, it is now possible to solve a problem whose formulation entails subdomains entirely using fem-fenics.

      by Eugenio Gianniti ( at August 12, 2014 10:54 PM

      August 11, 2014

      Eduardo Fernández

      Soft pencils down.

      It's been quite long since I posted here due to some personal situations. Anyway to sum up: I have finished ilu and ichol functions as I have planned in the beginning with great results.

      Things done after mid-term evaluation:
      • Implementing ICHOLT and ICHOL0 algorithms.
      • Fixing several bugs in ILU algorithms and introducing some enhancements for big sparse matrices with verly low densities.
      The files involved in ichol, within the repository, are:
      • src/
      • src/
      • ichol.m
      You can clone the code from the repo:

      Before going into the details of the algorithms' implementation, I want to point out some details about how ichol behave in MATLAB.

      1. In the real case the matrix must be symetric positive definite.  In the complex case the input matrix must be hermitian. That means: diagonal elements of the input and output matrix have to be non-zero, positive and real values. So that, at each iteration those conditions have to be fullfilled.
      2. If ichol is called just as L = ichol (A), Matlab ignores complex numbers and only work with their real part. Using L = ichol (A, setup) call, complex numbers are considered. Seriusly I do not understand why they do that and I have not followed that behaviour. Anyway if  to be 100% compatible I must change that, it would be only a line of code extra.

      Details of implementation


       In this file is located the implementation of ICHOL(0) algorithm. The zero-pattern of the output matrix is the same as the input one so it is known from the beginning how much  memory is needed to be allocated. The milu = ['on'|'off'] parameter indicates whether the dropped elements are added to the pivot or not (that keeps the colum sumation). 

      I will show two examples, one that corresponds to a big matrix with a very low density and the one that used Kai last year in his blog.

      Example 1:

        A = gallery ('poisson', 500);
        size (A)
        ans =

                250000   250000
        tic; L = ichol (A); toc;

        Elapsed time is 0.031718 seconds.
        density = nnz (A) / (size (A)(1))^2
        density =    1.9968e-05

        norm (A - L*L', 'fro') / norm (A, 'fro')
        ans =  0.0924207846384523

        ans =    2.28617974245061e-17

      It can be seen that the product L*L' is quite different from A, but the product L*L' will match A on its pattern (that is expected for the ICHOL(0) algorithm. The execution time is just given to give an idea of how fast the code is. It is executed in a i7 2.4GHz.

      Example 2:

       This example is taken from that post, written by Kai the past year. He faced problems with the michol option, obtaining different results from Matlab.

      A = [ 0.37, -0.05,  -0.05,  -0.07;
           -0.05,  0.116,  0.0,   -0.05;
           -0.05,  0.0,    0.116, -0.05;
           -0.07, -0.05,  -0.05,   0.202];

      A = sparse (A);
      opts.michol = 'on';
      L = ichol (A, opts);

      ans =

         0.60828   0.00000   0.00000   0.00000
        -0.08220   0.32014   0.00000   0.00000
        -0.08220   0.00000   0.32014   0.00000
        -0.11508  -0.18573  -0.18573   0.34607

      ans =

          0.6083         0         0         0
         -0.0822    0.3201         0         0
         -0.0822         0    0.3201         0
         -0.1151   -0.1857   -0.1857    0.3461

      Works fine.

      This file contains the implementation of ICHOLT algorithm. In this case the final structure of the output matrix is unknown. Therefore, a policy should be adopted for allocating memory. After trying different ways of doing that I end up using that one: 

        // max_len is the maximun length of ridx and data arrays for the output sparse matrix.
        max_len = sm.nnz ();
        max_len += (0.1 * max_len) > n ? 0.1 * max_len : n;

      What is done here is just to increment 10% of the actual size of the ridx and data internal arrays of the output sparse matrix. But only if that amount is larger than the dimension of the input matrix (n). In other case the increment in size is just n. That policy seems to work very well in every case I tested and do not slow down the process at all due to reallocations.

      Example 3:

      icholt accepts a parameter for controling the sparsity of the ouput matrix called droptol. If droptol = 0 then the complete factorization takes place. If we increase that value the output matrix will become more sparse as more elements will be dropped. Taking the same matrix than in example 1:

       A = gallery ('poisson', 500);
      opts.type= 'ict'

      % Complete factorization
      opts.droptol = 0;
      tic;L = ichol(A, opts);toc;
      Elapsed time is 46.0734 seconds. 
      norm (A - L*L', 'fro') / norm (A, 'fro')
      ans =    7.8595e-16

      % droptol = 1e-2
      tic;L = ichol(A, opts);toc;
      Elapsed time is 0.0650802 seconds.

      norm (A - L*L', 'fro') / norm (A, 'fro')
      ans =  0.016734

      % droptol = 1e-3
      tic;L = ichol(A, opts);toc;
      Elapsed time is 0.183416 seconds.

      norm (A - L*L', 'fro') / norm (A, 'fro')
      ans =  0.0021773

      % droptol = 1e-4
      tic;L = ichol(A, opts);toc;
      Elapsed time is 0.589693 seconds.
      norm (A - L*L', 'fro') / norm (A, 'fro')
      ans =    2.4820e-04

      As it can be seen, the higher the droptol parameter is, the sparser the matrix become. That lead to less execution times but on the other hand a higher error is obtained in the factorization. The complete factorization obviously have practically no error. Cool.

      Location of source files inside Octave core

      Now I've finished with the development of the algorithms, the final step is to integrate them into Octave core. For doing so I will create a subrepo of the default Octave repository and add the files. I have chosen the location for the functions looking at the last year repository Kai set.

      scripts/sparse: ilu.m ichol.m

      That is just a sugestion and should be revised and accepted by the maintainers.

      Future contributions

      There is a week left that I want to use it to start (and hopefully finish) the development of sprandsym extra parameters that Matlab have but Octave does not. As I submitted in the past a changeset for a similar functionality in sprand and sprandn, it will be much easier to implement for me.

      Also I am interested in developing some sparse linear solvers like minres and lsqr that Octave lacks. They are tightly related to the preconditioners I have been working on, and would be nice if they could be assigned to me for developing them.



      by Eduardo ( at August 11, 2014 06:47 AM

      August 09, 2014

      Eugenio Gianniti

      New features of meshfunction

      As you may recall from my last post, for DirichletBC to work in parallel runs I had to implement a new class, meshfunction. However it was still quite unfinished, with no way for the user to create one, except extracting it from a mesh produced by the msh package, no description to display, no way to save it. These days I have been tackling this issue: while at it I wondered what one could do with meshfunction and found out that they can come in handy when you are dealing with obstacles.

      At this link you can find a detailed explanation of the problem. It is a Poisson equation with variable diffusion coefficient on the unit square. Precisely, on [0.2, 1]x[0.5, 0.7] its value is 0.01, otherwise it is 1. The mentioned subset is the obstacle to diffusion, so we study its effect applying u = 0 on the y = 0 edge and u = 5 on y = 1. Here is the fem-fenics code:

      pkg load fem-fenics msh

      ufl start Subdomains
      ufl fe = FiniteElement "(""CG"", triangle, 2)"
      ufl u = TrialFunction (fe)
      ufl v = TestFunction (fe)
      ufl a0 = Coefficient (fe)
      ufl a1 = Coefficient (fe)
      ufl g_L = Coefficient (fe)
      ufl g_R = Coefficient (fe)
      ufl f = Coefficient (fe)
      ufl a = "inner(a0*grad (u), grad (v))*dx(0) + inner(a1*grad (u), grad (v))*dx(1)"
      ufl L = g_L*v*ds(1) + g_R*v*ds(3) + f*v*dx(0) + f*v*dx(1)
      ufl end

      # Create mesh and define function space
      x = y = linspace (0, 1, 65);
      [msh, facets] = Mesh (msh2m_structured_mesh (x, y, 0, 4:-1:1));

      V = FunctionSpace ("Subdomains", msh);

      # Define boundary conditions
      bc1 = DirichletBC (V, @(x, y) 5.0, facets, 2);
      bc2 = DirichletBC (V, @(x, y) 0.0, facets, 4);

      # Define problem coefficients
      a0 = Constant ("a0", 1.0);
      a1 = Constant ("a1", 0.01);
      g_L = Expression ("g_L", @(x, y) - 10*exp(- (y - 0.5) ^ 2));
      g_R = Constant ("g_R", 1.0);
      f = Constant ("f", 1.0);

      # Define subdomains
      domains = MeshFunction ("dx", msh, "cells.xdmf");

      # Define variational form
      a = BilinearForm ("Subdomains", V, V, a0, a1, domains);
      L = LinearForm ("Subdomains", V, g_L, g_R, f, facets, domains);

      # Assemble system
      [A, b] = assemble_system (a, L, bc1, bc2);
      sol = A \ b;
      u = Function ("u", V, sol);

      # Save solution in VTK format
      save (u, "subdomains");

      # Plot solution
      [X, Y] = meshgrid (x, y);
      U = u (X, Y);
      surf (X, Y, U);

      In the beginning there is the now familiar ufl block. As you might have noticed, subscripted measures appear in the definition of the bilinear form a and of the linear functional L. This is UFL notation for the integration on specific subsets of the computational domain. For instance, dx(1) is an integral over the subdomain marked with label 1, while ds(3) is an integral over the exterior edges marked with label 3. A third possibility, even if not used in this example, is to use dS for integrals on interior facets, which could be of use for interior penalty methods. Going back to the example, you can see that markers are used to enforce non-homogeneous Neumann conditions on the side edges and to assign the proper coefficient on the two subdomains.

      After defining the problem in UFL language, there are instructions to define the mesh, the function space, the essential boundary conditions and all the coefficients involved. All such lines come from fem-fenics before this summer or have been described in my previous posts, so I will not cover them in detail. The same applies for the assembly, solve and all the output in the end of the script. The only note is that the very last lines will error out in parallel runs: point-wise evaluations in DOLFIN can be performed only on local cells, but with meshgrid we are providing to every process the whole domain.
      The computed solution
      In between there are my latest efforts. At first, the brand new MeshFunction. With this, providing a mesh and a file name you can import a dolfin::MeshFunction. In this case it was saved in the XDMF format, here you can find the files needed to execute the script. DOLFIN uses this format for parallel input/output. It comprises a .h5 file storing data and a .xdmf with metadata useful to read the other one. The optional first argument is a string identifying the role of the returned meshfunction in the variational problem. In this case, with "dx" it will be searched for markers of the integrals on cells. All the previously mentioned measures are available, and "ds" is automatically attached to the meshfunction returned by Mesh. In the example this behaviour is exploited for the measure on edges.

      Afterwards, the mesh functions are passed as arguments to BilinearForm and LinearForm, so that the markers are available to assemble the system. In addition to the usual parameters, such as the name of the imported UFL problem, the function spaces and the coefficients, it is now possible to provide mesh functions properly labeled and they will be used.

      Currently fem-fenics allows for easily marking subdomains and exterior edges copying markers from the PDEtool representation returned by the functions of the msh package, which makes it quite tricky to properly identify the obstacle in the example. The approach used in the python interface to DOLFIN entails subclassing dolfin::Subdomain with the proper implementation of the inside method, then use an object of the derived class to mark a dolfin::MeshFunction. This could be an interesting feature to implement in the future also in fem-fenics.

      by Eugenio Gianniti ( at August 09, 2014 01:30 PM

      August 04, 2014

      Eugenio Gianniti

      MPI parallelism

      After quite a struggle, I have been able to obtain a working implementation of fem-fenics supporting MPI parallelism. Let's go through an example and highlight what has changed lately.

      pkg load fem-fenics msh

      ufl start Poisson
      ufl element = FiniteElement '("Lagrange", triangle, 1)'
      ufl u = TrialFunction (element)
      ufl v = TestFunction (element)
      ufl f = Coefficient (element)
      ufl g = Coefficient (element)
      ufl a = "inner (grad (u), grad (v))*dx"
      ufl L = f*v*dx + g*v*ds
      ufl end

      # Create mesh and define function space
      x = y = linspace (0, 1, 33);
      [mesh, facets] = Mesh (msh2m_structured_mesh (x, y, 1, 1:4));

      V = FunctionSpace ('Poisson', mesh);

      # Define boundary condition
      bc = DirichletBC (V, @(x, y) 0.0, facets, [2;4]);

      f = Expression ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
      g = Expression ('g', @(x,y) sin (5.0 * x));

      a = BilinearForm ('Poisson', V, V);
      L = LinearForm ('Poisson', V, f, g);

      # Compute solution
      [A, b] = assemble_system (a, L, bc);
      sol = A \ b;
      u = Function ('u', V, sol);

      # Save solution in VTK format
      save (u, 'poisson');

      The basic structure has remained the same. DOLFIN boasts the capability to be run both in serial and in parallel execution without intervening on the code, so I did my best to have the same behaviour from fem-fenics. The Poisson.m m-file above can be run either as you usually would do with any other m-file, or from the command line with an invocation such as:

      mpiexec -np 4 octave --eval Poisson

      Now, how is this possible? In the beginning, with the ufl block, the variational problem is defined in UFL language, written to an .ufl file and compiled via FFC. Since IO is performed, ufl.m ensures that only process zero will open and write to the file. Moreover, a MPI barrier makes sure that no process will proceed before the .ufl file is imported.

      As soon as the just-in-time compilation is over, there are two instructions to build the mesh, in this case on the unit square. For this, we rely on the msh package, which returns a PDE-tool-like representation of it. Mesh.oct must, then, convert it to DOLFIN internal representation and distribute it among processes. Here comes an issue: fem-fenics relies on markers present in the PDE-tool format to impose essential boundary conditions, and in serial runs dolfin::Mesh can store them, so that DirichletBC.oct needs just to know the boundary subset label. Unfortunately, this feature is not supported yet in parallel by the DOLFIN library, then Mesh.oct has been edited to return, if requested, also a meshfunction holding this information, in the example above facets. This way markers can be conveyed to DirichletBC.oct and boundary conditions can be applied on the correct edges.

      Further intervention was needed for the assembly and solve phase. In assemble_system.oct both the matrix and the vector are assembled locally on each portion of the mesh and, afterwards, gathered on process zero and joined, so that the system can be solved with the backslash instruction of Octave. In order to allow output in VTK format, in Function.oct the solution is split up and properly distributed among processes, so that each one holds the portion of degrees of freedom related to its subdomain and to the neighbouring vertices. After save.oct has written the solution to poisson.pvd and its auxiliary files, it can be visualised with ParaView.

      by Eugenio Gianniti ( at August 04, 2014 01:31 PM