Quadratic Programming in C#

I have manually translated and adapted the QuadProg solver for quadratic programming problems made by Berwin A. Turlach. His code was originally published under the GNU Library License, which has now been superseded by the GNU Lesser License. This adapted version honors the original work and is thus distributed under the same license.

Introduction

Despite the name, the terms linear or quadratic programming have little resemblance to the set of activities most people now know as programming. Those terms usually usually refers to a specific set of function optimization methods, i.e. methods which can be used to determine the maximum or minimum points of special kinds of functions under a given number of solution constraints. For example, suppose we would like to determine the minimum value of the function:

f(x, y) = 2x + y + 4

Under the constraints that x and y must be non-negative (i.e. either positive or zero). This may seem fairly simple and trivial, but remember that practical linear programming problems may have hundreds or even thousands of variables and possibly million constraints.

When the problem to be solved involves a quadratic function instead of a linear function, but still presents linear constraints, this problem can be cast as a quadratic programming problem. Quadratic functions are polynomial functions in each each term may have at most a total degree of 2. For example, consider the function

f(x, y, z) = 2x² + 5xy + y² - z² + x – 5.

Now let's check the sum of the degrees for each variable on the polynomial terms. We start by writing the missing terms of the polynomial

f(x, y, z) = 2 x2y0z0 + 5 x1y1z0 + 2 x0y2z0 - x0y0z2 + x1y0z0 - 5 x0y0z0

and then proceed to check the sum of the degrees at each term. In the first term, 2+0+0 = 2. For the second, 1+1+0 = 2, and so on. Those functions have a nice property that they can be expressed in a matrix form

f(x) = 1/2 xT Q x + cTx.

Here, x and c are vectors. The matrix Q is a symmetric matrix specifying how the variables combine in the quadratic terms of the function. If this matrix is positive definite, then the function is convex, and the optimization has a single, unique optimum (we say it has a global optimum). The coefficients c specify the linear terms of our function.

Source code

The available source code is based on a translation of the Fortran code written by Berwin A. Turlach. However, some modifications have been made. Fortran uses column-major ordering for matrices, meaning that matrices are stored in memory in sequential order of column elements. Almost all other languages use row-major ordering, including C, C++, Java and C#. In order to improve data locality, I have modified the code to use the transpose of the original matrices D and A. I have also modified the QP formulation adopted in the Goldfarb and Idnani paper to reflect the form presented in the introduction.

This code is part of the Accord.NET Framework. However, the version available in this blog post will most likely not be the most recently, updated, fixed and enhanced version of the code. For the latest version, be sure to download the latest version of the framework on the project site or through a NuGet package.

Using the code

The first step in solving a quadratic programming problem is, well, specifying the problem. To specify a quadratic programming problem, one would need two components: a matrix D describing the relationship between the quadratic terms, and a vector d describing the linear terms. Perhaps this would work better with an example.

Suppose we are trying to solve a minimization problem. Given a function, the goal in such problems is to find the correct set of function arguments which would result in the minimum possible value for the function. An example of a quadratic minimization problem is given below:

min f(x, y) = 2x² - xy + 4y² - 5x – 6y

subject to the constraints:

x - y = 5
x = 10

wolfram(generated with Wolfram Alpha)

However, note that this problem involves a set of constraints. The required solution for this minimization problem is required to lie in the interval specified by the constraints. More specifically, any x and y pair candidate for being a minimal of the function must respect the relations x - y = 5 and x >= 10. Thus, instead of lying in the unconstrained minimum of the function surface shown above, the solution lies slightly off the center of the surface. This is an obvious easy problem to solve manually, but it will fit for this demonstration.

wolfram2

As it can be seen (and also live demonstrated by asking Wolfram Alpha) the solution lies on the point (10,5), and the constrained minimum of the function is given by 170. So, now that we know what a quadratic programming problem looks like, how can we actually solve it?

Specifying the objective function

The first step in solving a quadratic programming problem is to specify the objective function. Using this code, there are three ways to specify it. Each of them has their own advantages and disadvantages.

1. Manually specifying the QP matrix.

This is the most common approach for numerical software, and probably the most cumbersome for the user. The problem matrix has to be specified manually. This matrix is sometimes denoted Q, D or H as it actually denotes the Hessian matrix for the problem.

The matrix Q is used to describe the quadratic terms of our problem. It is a n x n matrix, in which n corresponds to the number of variables in our problem, covering all possible combinations of variables. Recall our example given on the start of this section. We have 2 variables, x and y. Thus, our matrix Q is 2 x 2. The possible combinations for x and y are expressed in the table below.

  x y
x x*x x*y
y y*x y*y

To form our matrix Q, we can take all coefficients associated with each pair mentioned on the table above. The diagonal elements should also be multiplied by two (this is actually because the matrix is the Hessian matrix of the problem: it is the matrix of all second-order derivatives for the function. Since we have only at most quadratic terms, the elementary power rule of derivation “drops” the ² from the x² and y² terms – I think a mathematician would hit me with a stick for explaining it like this, but it serves well for a quick, non-technical explanation).

Remember our quadratic terms were 2x²  - 1xy + 4y². Writing the terms on their proper position and differentiating, we have:

 

As it can be seen, the matrix is also symmetric (and often, but not always, positive definite). The next step, more trivial, is to write a vector d containing the linear terms. The linear terms are –5x –6y, and thus our vector d can be given by:

Therefore our C# code can be created like this:

    double[,] Q =
    {   
        { +4,  -1 }, 
        { -1,  +8 },
    };
    
    double[] d = 
    {
        -5,
        -6 
    };

2. Using lambda expressions

This approach is a bit more intuitive and less error prone. However, it involves lambdas functions and some people find it hard to follow them. Another disadvantage is that we will lose the edit & continue debugging ability of visual studio. The advantage is that the compiler may catch some obvious syntax errors automatically.


    double x = 0, y = 0;

    QuadraticObjectiveFunction f =           // 2x² - xy + 4y² - 5x - 6y
        new QuadraticObjectiveFunction(() => 2*(x*x) - (x*y) + 4*(y*y) - 5*x - 6*y);

Note that the x and y variables could have been initialized to any value. They are only used as symbols, and not used in any computations.

3. Using text strings

This approach is more intuitive but a bit more error prone. The function can be specified using strings, as in a standard mathematical formula. However, since all we have are strings, there is no way to enforce static, compile time checking.


    QuadraticObjectiveFunction f = 
        new QuadraticObjectiveFunction("2x² - xy + 4y² - 5x - 6y");

Couldn’t be easier.

Specifying the constraints

The next step in specifying a quadratic programming problem is to specify the constraints. The constraints can be specified in almost the same way as the objective function.

1. Manually specifying the constraints matrix

The first option is to manually specify the constraints matrix A and vector b. The constraint matrix expresses the way the variables should be combined when compared to corresponding value on vector b. It is possible to specify either equality constraints or inequality constraints. The formulation used in this code is slightly different from the one used in Turlach’s original implementation. The constraints are supposed to be in the form:

A1 x = b1
A2 x = b2

This means that each line of matrix A expresses a possible combination of variables x which should be compared to the corresponding line of b. An integer variable m can be specified to denote how many of the first rows of matrix A should be treated as equalities rather than inequalities. Recall that in our example the constraints are given by 1x -1y = 5 and 1x = 10. Lets write this down in a tabular form:

# x y ? b
q1 1 -1 = 5
q2 1 0 = 10

Thus our matrix A and vector b can be specified as:

And not forgetting that m = 1, because the first constraint is actually an equality.

2. Using classes and objects

A more natural way to specify constraints is using the classes and objects of the Accord.NET Framework. The LinearConstraint class allows one to specify a single constraint using an object-oriented approach. It doesn’t have the most intuitive usage on earth, but has much more expressiveness. It can also be read aloud, it that adds anything! :-)


    List<LinearConstraint> list = new List<LinearConstraint>();

    list.Add(new LinearConstraint(numberOfVariables: 1)
    {
        VariablesAtIndices = new[] { 0 }, // index 0 (x)
        ShouldBe = ConstraintType.GreaterThanOrEqualTo,
        Value = 10
    });

    list.Add(new LinearConstraint(numberOfVariables: 2)
    {
        VariablesAtIndices = new int[] { 0, 1 }, // index 0 (x) and index 1 (y)
        CombinedAs     = new double[] { 1, -1 }, // when combined as 1x -1y
        ShouldBe       = ConstraintType.EqualTo,
        Value = 5
    });

The specification is centered around the notion that variables are numbered and have an associated index. For example, x is the zero-th variable of the problem. Thus x has an index of 0 and y has an index of 1. So for example, reading aloud the last constraint, it is possible to express how the variables at indices 0 and 1, when combined as 1x and –1y, should be equal to value 5.

2. Using lambda expressions

A more intuitive way to express constraints is again using lambda expressions. And again the problems are the same: some people find it hard to follow and we lose edit & continue.


    var constraints = new List<LinearConstraint>();
    constraints.Add(new LinearConstraint(f, () => x - y ==  5));
    constraints.Add(new LinearConstraint(f, () => x     >= 10));

3. Using text strings

Same as above, but with strings.


    var constraints = new List<LinearConstraint>();
    constraints.Add(new LinearConstraint(f, "x - y = 5"));
    constraints.Add(new LinearConstraint(f, "x >= 10"));

Finally, creating and solving the problem

Once we have specified what do we want, we can now ask the code for a solution. In case we have opted for manually specifying the matrix A, vector b and integer m, we can use:


    // Create the optimization problem
    var solver = new GoldfarbIdnaniQuadraticSolver(numberOfVariables:2, A, b, m);

In case we have opted for creating a list of constraints instead, we can use:


    // Create our optimization problem
    var solver = new GoldfarbIdnaniQuadraticSolver(numberOfVariables: 2, constraints: list);

After the solver object has been created, we can call Minimize() to solve the problem. In case we have opted for manually specifying Q and d, we can use:

    // Attempt to solve the problem
    double minimumValue = solver.Minimize(Q, d);

And in case we have opted for creating a QuadraticObjectiveFunction object, we can use:


    // Attempt to solve the problem
    double minimumValue = target.Minimize(f);

In either case, the solution will be available in the Solution property of the solver object, and will be given by:


    double value = solver.Value;    // f(x,y) = 170
    double x = solver.Solution[0];  // x = 10
    double y = solver.Solution[1];  // y = 5

Sample application

The Accord.NET Framework now includes a sample application demonstrating the use of the Goldfarb-Idnani Quadratic Programming Solver. It can be downloaded at the Accord.NET Framework site, and also comes together with recent versions of the framework (> 2.6).


Solver
Solver sample application included in the Accord.NET Framework.

Remarks

Because the code has been translated by hand (in contrast of using automatic translators such as f2c) there could be potential bugs in the code. I have tested the code behavior against R’s quadprog package and still didn’t find errors. But this does not mean the code is bug-free. As always, as is the case of everything else in this blog, this code is published in good faith, but I can not guarantee the correctness of everything. Please read the disclaimer for more information.

References

Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) for Unconstrained Optimization in C#

The Limited-memory Broyden-Fletcher-Goldfarb-Shanno method is an optimization method belonging to the family of quasi-Newton methods for unconstrained non-linear optimization. In short terms, it is an off-the-shelf optimizer for seeking either minimum or maximum points of a any differentiable and possibly non-linear function, requiring only an expression of the function and its gradient. The goal of this article is to provide and demonstrate a C# port of the L-BFGS method originally written by Jorge Nocedal in Fortran.

Introduction

Function optimization is a common problem found in many numerical applications. Suppose we have a differentiable function f : Rn → R and we would like to obtain its minimal or maximum value while traversing its space of possible input values. Those kind of problems arise often in parameter optimization of machine learning models and other statistic related applications (and in a myriad of other applications too – however, we will pay special attention to cases related to machine learning).

The problem of maximizing or minimizing a function can be simply stated as max f(x) or min f(x). When there aren’t any constraints limiting possible values for x, it is widely known from calculus that the maximum or a minimum of such a function would occur when the first derivatives f’(x) of the function f(x) in respect to x are zero. Newton’s method is a common method for finding the roots of a differentiable function and is indeed a suitable method for finding those values such that first derivatives of a function are zero.

 

Surface of the example function. Contour lines of the example function.

Example of a convex, but non-linear function f(x,y) = exp{-(x-1)²} + exp{-(y-2)²/2}. Images have been created using Wolfram Alpha.

 

However, while Newton’s method may work very well with functions of a single variable, the generalization to higher dimensions will require the replacement of the first derivative f’(x) with the function's gradient vector g of partial derivatives, and the second derivative f’’(x) with the inverse Hessian matrix H-1. The problem is that the Hessian is a square matrix with as many rows and columns as parameters of our function f, and, besides being very costly and challenging to compute, even its size may be infeasible to accommodate in main memory when dealing with very high dimensional problems.

To overcome those limitations, several methods attempt to replace the Hessian matrix with an approximation extracted from the first partial derivatives alone. Those methods are called quasi-Newton methods, as they replace H in Newton’s method by an approximation instead of using the full Hessian. In case the function being optimized has any particular form or special characteristic which could be exploited, it is also possible to use more specialized methods such as the Gauss-Newton, Levenberg-Marquardt or even lower-bound approximation methods to avoid computing the full Hessian.

The quasi-Newton BFGS method for non-linear optimization builds an approximation for the Hessian matrix based on estimates extracted solely from first order information. However, even being cheaper to compute, the memory requirements are still high as the method still needs to accommodate a dense N x N matrix for the inverse Hessian approximation (where N is the number of variables for the function). To overcome this problem, the L-BFGS method, or the limited-memory version of the BFGS method, never forms or stores the approximation matrix, but only a few vectors which can be used to represent it implicitly.

Source code

The code presented here is based upon a direct translation of the original code by Jorge Nocedal. His original code was released under the permissive BSD license, and honoring the original license and the author's considerations, this code is released under the BSD as well.

L-BFGS Class diagrams

The L-BFGS method is implemented within the BroydenFletcherGoldfarbShanno class. Upon object creation, it is possible to specify a function and its gradient through either the use of Lambda expressions or by specifying handles for named functions. After the object has been initialized, a simple call to Minimize runs the standard L-BFGS method until convergence. The details for the Minimize method is given below.

    /// <summary>
    ///   Optimizes the defined function. 
    /// </summary>
    /// 
    /// <param name="values">The initial guess values for the parameters.</param>
    /// <returns>The values of the parameters which optimizes the function.</returns>
    /// 
    public unsafe double Minimize(double[] values)
    {
        if (values == null)
            throw new ArgumentNullException("values");

        if (values.Length != numberOfVariables)
            throw new DimensionMismatchException("values");

        if (Function == null) throw new InvalidOperationException(
            "The function to be minimized has not been defined.");

        if (Gradient == null) throw new InvalidOperationException(
            "The gradient function has not been defined.");


        // Initialization
        x = (double[])values.Clone();
        int n = numberOfVariables, m = corrections;


        // Make initial evaluation
        f = getFunction(x);
        g = getGradient(x);

        this.iterations = 0;
        this.evaluations = 1;


        // Obtain initial Hessian
        double[] diagonal = null;

        if (Diagonal != null)
        {
            diagonal = getDiagonal();
        }
        else
        {
            diagonal = new double[n];
            for (int i = 0; i < diagonal.Length; i++)
                diagonal[i] = 1.0;
        }


        fixed (double* w = work)
        {
            // The first N locations of the work vector are used to
            //  store the gradient and other temporary information.

            double* rho = &w[n];                   // Stores the scalars rho.
            double* alpha = &w[n + m];             // Stores the alphas in computation of H*g.
            double* steps = &w[n + 2 * m];         // Stores the last M search steps.
            double* delta = &w[n + 2 * m + n * m]; // Stores the last M gradient diferences.


            // Initialize work vector
            for (int i = 0; i < g.Length; i++)
                steps[i] = -g[i] * diagonal[i];

            // Initialize statistics
            double gnorm = Norm.Euclidean(g);
            double xnorm = Norm.Euclidean(x);
            double stp = 1.0 / gnorm;
            double stp1 = stp;

            // Initialize loop
            int nfev, point = 0;
            int npt = 0, cp = 0;
            bool finish = false;

            // Make initial progress report with initialization parameters
            if (Progress != null) Progress(this, new OptimizationProgressEventArgs
                (iterations, evaluations, g, gnorm, x, xnorm, f, stp, finish));


            // Start main
            while (!finish)
            {
                iterations++;
                double bound = iterations - 1;

                if (iterations != 1)
                {
                    if (iterations > m)
                        bound = m;

                    double ys = 0;
                    for (int i = 0; i < n; i++)
                        ys += delta[npt + i] * steps[npt + i];

                    // Compute the diagonal of the Hessian
                    // or use an approximation by the user.

                    if (Diagonal != null)
                    {
                        diagonal = getDiagonal();
                    }
                    else
                    {
                        double yy = 0;
                        for (int i = 0; i < n; i++)
                            yy += delta[npt + i] * delta[npt + i];
                        double d = ys / yy;

                        for (int i = 0; i < n; i++)
                            diagonal[i] = d;
                    }


                    // Compute -H*g using the formula given in:
                    //   Nocedal, J. 1980, "Updating quasi-Newton matrices with limited storage",
                    //   Mathematics of Computation, Vol.24, No.151, pp. 773-782.

                    cp = (point == 0) ? m : point;
                    rho[cp - 1] = 1.0 / ys;
                    for (int i = 0; i < n; i++)
                        w[i] = -g[i];

                    cp = point;
                    for (int i = 1; i <= bound; i += 1)
                    {
                        if (--cp == -1) cp = m - 1;

                        double sq = 0;
                        for (int j = 0; j < n; j++)
                            sq += steps[cp * n + j] * w[j];

                        double beta = alpha[cp] = rho[cp] * sq;
                        for (int j = 0; j < n; j++)
                            w[j] -= beta * delta[cp * n + j];
                    }

                    for (int i = 0; i < diagonal.Length; i++)
                        w[i] *= diagonal[i];

                    for (int i = 1; i <= bound; i += 1)
                    {
                        double yr = 0;
                        for (int j = 0; j < n; j++)
                            yr += delta[cp * n + j] * w[j];

                        double beta = alpha[cp] - rho[cp] * yr;
                        for (int j = 0; j < n; j++)
                            w[j] += beta * steps[cp * n + j];

                        if (++cp == m) cp = 0;
                    }

                    npt = point * n;

                    // Store the search direction
                    for (int i = 0; i < n; i++)
                        steps[npt + i] = w[i];

                    stp = 1;
                }

                // Save original gradient
                for (int i = 0; i < g.Length; i++)
                    w[i] = g[i];


                // Obtain the one-dimensional minimizer of f by computing a line search
                mcsrch(x, ref f, ref g, &steps[point * n], ref stp, out nfev, diagonal);

                // Register evaluations
                evaluations += nfev;

                // Compute the new step and
                // new gradient differences
                for (int i = 0; i < g.Length; i++)
                {
                    steps[npt + i] *= stp;
                    delta[npt + i] = g[i] - w[i];
                }

                if (++point == m) point = 0;


                // Check for termination
                gnorm = Norm.Euclidean(g);
                xnorm = Norm.Euclidean(x);
                xnorm = Math.Max(1.0, xnorm);

                if (gnorm / xnorm <= tolerance)
                    finish = true;

                if (Progress != null) Progress(this, new OptimizationProgressEventArgs
                    (iterations, evaluations, g, gnorm, x, xnorm, f, stp, finish));
            }
        }

        return f; // return the minimum value found (at solution x)
    }

The version of the code detailed here only supports Minimization problems. However, it is possible to note that any minimization problem can be converted into a maximization problem and vice-versa by taking the opposite of the function and its gradient.

Using the code

Code usage is rather simple. Suppose we would like to maximize the function g(x,y) = exp{-(x-1)²} + exp{-(y-2)²/2}:

Since we would like to perform a maximization, we first have to convert it to a minimization problem. The minimization version of the function is simply given by taking f(x,y) = –g(x,y):

Which is a convex function, as can be seen by plotting its surface. As it can be seen, the minimum of the function lies on the point (x,y) = (1,2). As expected, this point coincides with the roots of the partial derivative functions, as shown in the line plots below:

 

Example optimization cast as an minimization problem.

The minimization problem min f(x,y) = -(exp(-(x-1)²) + exp(-(y-2)²/2)), computed by Wolfram Alpha.

Partial derivative for x. Partial derivative for y.

The roots of the partial derivatives in respective to x (left) and y (right). It is possible to see that the roots occur at 1 and 2, respectively. Images computed using Wolfram Alpha. 

 

The first step in solving this optimization problem automatically is first to specify the target function. The target function f, in this case, can be specified through a lambda function in the form:

    Func<double[], double> f =
        (x) => -Math.Exp(-Math.Pow(x[0] - 1, 2)) - Math.Exp(-0.5 * Math.Pow(x[1] - 2, 2));

The next step is to specify the gradient g for the function f. In this particular case, we can manually compute the partial derivatives to be df/dx = -2 e^(-(x-1)^2) (x-1) and df/dy = -e^(-1/2 (y-2)^2) (y-2), in respect to x and y, respectively. Writing a lambda function to compute the gradient vector g, we have:

    Func<double[], double[]> g = (x) => new double[] 
    {
        2 * Math.Exp(-Math.Pow(x[0] - 1, 2)) * (x[0] - 1),
        Math.Exp(-0.5 * Math.Pow(x[1] - 2, 2)) * (x[1] - 2)
    };

We can also note that this example was rather simple, so the gradient vector was easy to calculate. However, the gradient could also have been computed automatically using Accord.NET's FiniteDifferences class. In either case, all we have to do now is to create our L-BFGS solver and call its Minimize() method to begin optimization:

    // Create the L-BFGS solver
    BroydenFletcherGoldfarbShanno lbfgs =
        new BroydenFletcherGoldfarbShanno(numberOfVariables: 2, function: f, gradient: g);

    // Minimize the function
    double minValue = lbfgs.Minimize(); // should be -2

After the optimization is complete, the solution parameter will be available in the Solution property of the lbfgs object, and will be equal to { 1.0, 2.0 }. And also in case one is interested in progress reports (such as in the case of optimizing very large functions), it is also possible to register a listener to the Progress event handler. The complete version of the sample application accompanying the source code is given below:

    // Suppose we would like to find the minimum of the function
    // 
    //   f(x,y)  =  -exp{-(x-1)²} - exp{-(y-2)²/2}
    //

    // First we need write down the function either as a named
    // method, an anonymous method or as a lambda function:

    Func<double[], double> f = (x) =>
        -Math.Exp(-Math.Pow(x[0] - 1, 2)) - Math.Exp(-0.5 * Math.Pow(x[1] - 2, 2));

    // Now, we need to write its gradient, which is just the
    // vector of first partial derivatives del_f / del_x, as:
    //
    //   g(x,y)  =  { del f / del x, del f / del y }
    // 

    Func<double[], double[]> g = (x) => new double[] 
    {
        // df/dx = {-2 e^(-    (x-1)^2) (x-1)}
        2 * Math.Exp(-Math.Pow(x[0] - 1, 2)) * (x[0] - 1),

        // df/dy = {-  e^(-1/2 (y-2)^2) (y-2)}
        Math.Exp(-0.5 * Math.Pow(x[1] - 2, 2)) * (x[1] - 2)
    };

    Console.WriteLine("Solving:");
    Console.WriteLine();
    Console.WriteLine("   min f(x,y)  =  -exp{-(x-1)²} - exp{-(y-2)²/2}");
    Console.WriteLine();

    // Finally, we can create the L-BFGS solver, passing the functions as arguments
    var lbfgs = new BroydenFletcherGoldfarbShanno(numberOfVariables: 2, function: f, gradient: g);

    // And then minimize the function:
    double minValue = lbfgs.Minimize();
    double[] solution = lbfgs.Solution;

    // The resultant minimum value should be -2, and the solution
    // vector should be { 1.0, 2.0 }. The answer can be checked on
    // Wolfram Alpha by clicking the following the link:

    // http://www.wolframalpha.com/input/?i=maximize+%28exp%28-%28x-1%29%C2%B2%29+%2B+exp%28-%28y-2%29%C2%B2%2F2%29%29

    Console.WriteLine("The minimum value of {0:0} occurs at the solution point ({1:0},{2:0})", 
        minValue, solution[0], solution[1]);

    

Conclusion

In this post, we showed how to use a reduced set of the Accord.NET Framework's to perform non-linear optimization. This routine is also used by the Conditional Random Fields and Hidden Conditional Random Fields trainers to optimize parameters for such models. The solver routines have been adapted from the original Fortran's source code from Nocedal, which, tracing back to a 2001 message from the author, also have been reported to be available under the public domain. Nevertheless, the reduced set of the framework available for download within this post is available under a BSD license, as an alternative to the version available within the Accord.NET Framework which is available only under a LGPL license.

As always, I hope someone will find it useful :-)

References


 
 
Copyright © 2010 Cesar Roberto de Souza - Some rights reserved. Otherwise specified, all works presented here are licensed under the Creative Commons Attribution-Noncommercial-Share Alike 3.0 Unported License.