ComputSimu

Computing the Goodwin-Staton integral

2016-11-07T09:50:00.000+01:00

In this post we will explore different approaches for the evaluation of the Goodwin-Staton integral, which is one of the generalizations of the complementary error function. The Goodwin-Staton integral is defined as $$\begin{equation}G(z) = \int_0^{\infty} \frac{e^{-t^2}}{t+z} \; dt \quad |\mathrm{arg}\, z| < \pi \end{equation}$$ From now on, we consider the case where $z$ is real positive number. Note that this special functions can be expressed in terms of other elementary special functions; the error function and the exponential integral. Thus, $$\begin{equation}G(x) = -\frac{e^{-x^2}}{2}(\pi i \mathrm{erf}(ix) + E_i(x^2)) \end{equation}$$ Additionally, the previous representation can be slightly simplified by considering the following connection formulas: $E_i(x) = -E_1(-x)$ and $\mathrm{erfi}(x) = -i\mathrm{erf}(ix)$, where $\mathrm{erfi}(x)$ is the imaginary error function, both valid when $x > 0$. In particular, $E_1(x)$ for real negative argument is not currently supported in Scipy, which will be used for implemented each of the proposed methods.

  from scipy.special import expi, erfi
  from scipy import exp, pi


  def goodwin_staton_1(x):
      mx2 = -x*x
      ans = 0.5 * exp(mx2) * (pi * erfi(x) - expi(-mx2))
      return ans

When $x$ is small the previous equation can be safely evaluated, but for $x \gtrsim 20$ the computation will overflow. This problem is avoidable if we consider some different approaches. A sensible approach consists of substituting the first term using a well-known relation between the imaginary error function and Dawson integral $F(x)$, $$\begin{equation}\frac{\pi \mathrm{erfi}(x)}{2e^{x^2}} = \sqrt{\pi}F(x) \end{equation}$$ On the other hand, the second term involving the exponential integral $E_1(-x^2)$ can be replaced by its asymptotic expansion as follows, $$\begin{equation} \frac{e^{-x^2} E_1(-x^2)}{2} \sim -\frac{1}{2x^2} \sum_{k=0}^{\infty} \frac{k!}{x^{2k}},\quad x \to \infty \end{equation}$$ Therefore, for large argument the Goodwin-Staton integral can be computed as $$\begin{equation} G(x) \sim \sqrt{\pi} F(x) -\frac{1}{2x^2} \sum_{k=0}^{\infty} \frac{k!}{x^{2k}},\quad x \to \infty \end{equation}$$

  from scipy.special import dawsn
  from scipy import pi, exp, sqrt


  def goodwin_staton_2(x):
      # compute Dawson's integral
      term1 = sqrt(pi) * dawsn(x)

      # compute series
      maxiter = 1000
      eps = 2.2e-16
      acc = False
      s = 1.0
      sp = s
      n = 1
      d = 1
      x2 = x*x

      for k in range(1, maxiter):
          n *= k
          d *= x2
          s += n / d
          if abs((s - sp) / sp) < eps:
              if acc:
                  break
              else:
                  acc = True
          else:
              sp = s

      ans = term1 - 0.5 * s / x2
      return ans

Another alternative is representing the second term as definite integral, $$\begin{equation} \frac{e^{-x^2} E_1(-x^2)}{2} = - \frac{1}{2}\int_0^1 \frac{dt}{x^2 + \log(t)} \end{equation}$$ Therefore, the Goodwin-Station is represented as follows, $$\begin{equation} G(x) = \sqrt{\pi}F(x) - \frac{1}{2}\int_0^1 \frac{dt} {x^2 + \log(t)} \end{equation}$$

from scipy.integrate import quad
from scipy.special import dawsn
from scipy import log, inf, sqrt, pi


def integrand(t, x):
    return 1 / (x * x + log(t))


def term2_integral(x):
    eps = 2.2e-16
    return quad(integrand, 0, 1, args=x, epsabs=eps, epsrel=eps, limit=100)


def goodwin_staton_3(x):
    # compute Dawson's integral
    term1 = sqrt(pi) * dawsn(x)

    # compute integral
    term2, err = term2_integral(x)

    return term1 - 0.5 * term2, err

Finally, we can evaluate the integral using an adaptive quadrature algorithm implemented in scipy.integrate.quad, based on the Fortran library QUADPACK.

from scipy.integrate import quad
from scipy import exp, inf


def integrand(t, x):
    return exp(-t**2) / (t + x)


def goodwin_staton_4(x):
    eps = 2.2e-16
    return quad(integrand, 0, inf, args=x, epsabs=eps, epsrel=eps, limit=100)

In order to perform a check of correctness, we evaluate the integral definition using the tanh-sinh or double exponential integration algorithm implemented in mpmath with 20 digits of precision. For example,

>>> from mpmath import *
>>> mp.dps = 20
>>> z = 10
>>> quad(lambda t: exp(-t*t)/(t + z), [0, inf])
mpf('0.084021593706602168771586')

Interior Point Method: Primal Affine scaling algorithm in Python

2015-07-01T12:44:00.003+02:00

On this post we implement The Primal Affine scaling algorithm, one of the Interior Point Methods. Let us start by the Problem Formulation:

$\begin{equation*}\label{a} \begin{aligned} & {\underset{x} min} & & c^T x\\ & \text{subject to} & & Ax= b \quad \quad \quad (P)\\ &&& x \ge 0\\ \end{aligned} \end{equation*}$
where $x\in \mathbb{R}^n$ is the vector of variables, $c \in \mathbb{R}^n$ is the vector of linear costs, $A \in \mathbb{R}^{m\times n}$ is the matrix of constraint and $b \in \mathbb{R}^m$ is the right-hand-side vector. In general, two main hypothesis are considered: Matrix $A$ is full row rank and $P$ is nondegenerate. However, as shown later, some techniques are available to deal with problems where $A$ is not full rank.

This algorithm includes a Big-M formulation and therefore the Big-M value must be chosen in advance, several approaches are possible:
A first approach is start $M$ with a large value, $10^{30}$, for instance. However, if $M$ is too large it can causes numerical problems, since it is known that the presence of coefficients larger than the rest of variables by many orders of magnitude have potentially undesirable numerical aspects. Therefore, ideally we should search for a $\textit{"sufficiently large"}$ $M$. A standard decision is choose a $M$ value which is at least $n$ times larger than the largest value of any of the variables in the specific Linear programming problem.
For this implementation, the second approach with $n=100$ was chosen. Nevertheless, the first approach has been also tested, obtaining bad results in several problems, with numerical instability when approaching to the optimal value, hence,

$\begin{equation}\label{big-m} M = n \cdot \mathrm{max}_{1 \le i \le n}\;|c_i| \end{equation}$
where $c$ is the cost vector and $n$ the number of variables. In any case, the presented implementation provides the option of supplying a custom $M$ in case the user previously knows an upper or lower bound on some quantity (e.g. a capacity limit).

We focus now on the algorithm implementation. The primal affine scaling algorithm has been implemented in a class called PrimalAffine() with several arguments. Furthermore, a function called _pre_check checks the input of the class: If the matrix and vectors are of the class ndarray. A relevant function is _check: This function validates that the vector $b$ and $c$ are column vectors and checks if the matrix $A$ is a full matrix, otherwise throw an exception and do not continue solving the problem. This check can be removed if the Cholesky decomposition for semidefinite matrices is implemented (we implement it later).

The following code is the main routine in the algorithm, where two functions $({\tt \_dual\_gap}$ and ${\tt \_compute\_y})$ are included. Basically, these two separated functions aim to improve readability.

import numpy as np
import scipy.io
from scipy.sparse import csr_matrix
import scipy
import pandas as pd
import time
import matplotlib.pyplot as plt
import math
import numba


class PrimalAffine(object):
    def __init__(self, c, a, b, opt_gap=1.0e-6, eps=1.0e-12, rho=0.95, m=None):
        self.c = c
        self.a = a
        self.b = b
        self.opt_gap = opt_gap
        self.eps = eps
        self.rho = rho
        self._alphas = []
        self._f_values = []
        self._gaps = []
        self._iter = []
        self._pre_check()

        if m is None:
            self.m = 100 * np.max(abs(self.c))
        else:
            self.m = m
        print("Max value = {}".format(np.max(abs(self.c))))
        print("Big-M = {}".format(self.m))

        self._optimize_primal()

    def _pre_check(self):
        if not (isinstance(self.a, np.ndarray) or isinstance(self.b, np.ndarray) or isinstance(self.c, np.ndarray)):
            raise Exception('Inputs must be a numpy arrays')

    def _check(self, ra, ca, rb, rc):
        if np.linalg.matrix_rank(self.a) != ra:
            raise Exception('Matrix A is not full rank')

        if rb != ra:
            if np.size(self.b) == ra:
                self.b = np.transpose(self.b)
            else:
                raise AttributeError('dimension of vector b is not correct')

        if rc != ca:
            if np.size(self.c) == ca:
                self.c = np.transpose(self.c)
            else:
                raise AttributeError('dimension of vector b is not correct')

    def _optimize_primal(self):

        r_a, c_a = np.shape(self.a)
        r_b = np.shape(self.b)[0]
        r_c = np.shape(self.c)[0]

        self._check(r_a, c_a, r_b, r_c)

        x = np.ones((c_a, 1))
        a = np.c_[self.a, self.b - np.dot(self.a, x)]
        x = np.ones((c_a + 1, 1))
        c = np.vstack([self.c, self.m])
        d = np.eye(c_a + 1)
        y = self._compute_y(a, d, c)

        i = 0
        while self._dual_gap(c, x, self.b, y) > self.opt_gap:
            z = c - np.dot(np.transpose(a), y)
            dx = np.dot(-d, z)
            if np.all(dx >= 0):
                print('Unbounded problem')
                break
            alpha = self.rho * np.min(-x[dx < self.eps] / dx[dx < self.eps].ravel())
            x += alpha*dx
            d = np.diag(x.ravel())**2
            y = self._compute_y(a, d, c)
            i += 1

            f_obj = np.dot(np.transpose(c), x)
            self._f_values.append(f_obj.ravel()[0])
            self._alphas.append(alpha.ravel()[0])
            self._gaps.append(self._dual_gap(c, x, self.b, y).ravel()[0])

        if x[c_a] > 1e-5:
            print('Infeasible problem. x[n+1] = {}'.format(x[c_a]))
        else:
            pass
            print('\n', x[0:c_a])

    def _compute_y(self, aa, dd, cc):
        term1 = np.dot(np.dot(aa, dd), np.transpose(aa))
        term2 = np.dot(np.dot(aa, dd), cc)
        L = self._chol_semidefinite(term1)
        sol = self._chol_solve(L, term2)
        return sol

    @staticmethod
    def _dual_gap(cc, xx, bb, yy):
        term1 = np.dot(np.transpose(cc), xx)
        term2 = term1 - np.dot(np.transpose(bb), yy)
        term3 = 1.0 + abs(term1)
        return abs(term2)/term3

    @staticmethod
    @numba.autojit
    def _chol_semidefinite(A):
        n = len(A)

        # create zero matrix for L
        L = np.zeros((n, n))

        # Perform the Cholesky decomposition
        for i in range(n):
            for k in range(i+1):
                tmp_sum = 0
                for j in range(k):
                    tmp_sum += L[i, j] * L[k, j]

                if i == k:  # diagonal
                    L[i, k] = math.sqrt(A[i, i] - tmp_sum)
                    if L[i, k] <= 0:
                        L[i, k] = 10**10

                else:
                    L[i, k] = (1.0 / L[k, k] * (A[i, k] - tmp_sum))
        return L

    @staticmethod
    def _chol_solve(L, b):
        L_inv = np.linalg.inv(L)
        z = np.dot(L_inv, b)
        y = np.dot(np.transpose(L_inv), z)
        return y


    def results(self):
        d = {'gap': self._gaps, 'fobj': self._f_values, 'alpha': self._alphas}
        df = pd.DataFrame(d)
        print(df)

A few comments about the above code:

An important point is to check whether the problem is infeasible or not. Therefore, if $x_{n+1}^* > \epsilon_{inf}$ then problem $(P)$ is infeasible, due to $M \gg 0$ in the optimal solution and hence, $x_{n+1}^* \approx 0$. We set $\epsilon_{inf}$ to $1{e-}5$.
We included the function to compute the Cholesky decomposition in a positive semidefinite matrix. Unfortunately, this implementation is not very efficient, even using good libraries like ${\tt Numpy}$, so in order to speed up the computation we use the library ${\tt Numba}$ so that the code can be just-in-time compiled by using the LLVM compiler. The ${\tt @numba.autojit}$ decorator tells ${\tt Numba}$ to compile the function. In average, by using a JIT compiler, it is possible to achieve $\times 3$ speed up. This function is based on the Cholesky-Banachiewicz algorithm applying a modification, so that if a value of the diagonal is 0 (zero pivot), it is replaced by a large value $L_{ik} = 10^{64}$ that is in a sense surrogates for $\infty$.

Furthermore, we can test our algorithm solving problems of the Netlib collection. The Netlib problems in format ${\tt .mat}$ have been downloaded from NETLIB problems. As the algorithm was implemented in Python a function has been included to read .mat files.

def mat_parse(file):
    mat_content = scipy.io.loadmat(file)
    mat_struct = mat_content['Problem']
    val = mat_struct[0, 0]
    A = csr_matrix(val['A']).todense()
    b = val['b']
    c = val['aux'][0][0][0]
    return A, b, c


if __name__ == '__main__':
    # http://www.cise.ufl.edu/research/sparse/mat/LPnetlib/
    A, b, c = mat_parse('lp_sctap3.mat')

    print(np.shape(A))
    print(np.shape(b))
    print(np.shape(c))

    t1 = time.perf_counter()
    primal = PrimalAffine(c, A, b, opt_gap=1e-6)
    elapsed = (time.perf_counter() - t1)
    print(elapsed)
    primal.results()

Finally, this is a small example that shows the behaviour of this algorithm: $\begin{equation} A = \begin{pmatrix} 2 & 1 & -1 & 0\\ 3& 4& 0 & 1\\ \end{pmatrix} \quad b = \begin{pmatrix} 2 \\ 12 \end{pmatrix} \quad c=\begin{pmatrix} 3 \\ 1 \\ 0 \\ 0 \end{pmatrix} \end{equation}$ And the following table and plot summarize the steps towards the optimal solution of the problem:

Iter	alpha	fobj	gap	x1	x2	x3	x4
0	0.0086	19.52871	0.76302	0.85813	1.95431	1.67057	1.40837
1	0.0636	5.15115	0.35592	0.78414	2.04874	1.61701	1.44264
2	0.7896	2.96509	0.08685	0.03921	2.54061	0.61903	1.71583
3	3.1027	2.27926	0.06842	0.03034	2.17289	0.23358	3.2172
4	4.4464	2.05176	0.01697	0.02578	1.96012	0.01168	4.082
5	36.9018	2.01468	0.00486	0.00129	2.00406	0.00664	3.97979
6	140.6546	2.00183	0.00061	0.00106	1.99833	0.00044	4.00352
7	900.1966	2.00055	0.00018	0.00005	2.00016	0.00027	3.9992
8	3579.9455	2.00009	0.00003	0.00004	1.99993	0.00001	4.00016
9	22196.8995	2.00002	0.00001	0	2.00001	0.00001	3.99997
10	101489.1179	2	0	0	2	0	4.00001
11	425757.3225	2	0	0	2	0	4

Figure 1: Steps towards the optimal solution. Optimal solution = 2. x1 = 0 and x2 = 2. Total number of iterations = 12.

AMPL reporting solver statistics

2015-04-24T03:25:00.000+02:00

This a simple script that I have been working on, which provides some Solver statistics to AMPL. At this moment, the script can handle the MINOS solver output and generate cool charts to summarize the process till the solver reaches a solution.
As usual, I have written a python script that is called by the shell command in AMPL.

option solver minos;
option minos_options "Superbasics_limit=100 summary_file=1 1=summary.txt summary_level = 1 summary_frequency=1";
  
printf "MINOS_report: generating AMPL results...\n";
shell 'C:\Python27\python.exe "MINOS_report.py" summary.txt';

As said, this script was written to work with MINOS, so after choosing MINOS as solver, we must provide to MINOS the output file which will contain the solver's statistics. The python script just parse the output file, clean it and generate the chart shown below. This code is just an example and does not aim to be very fast or optimum in terms of computational efficiency since is just a first test.

import pandas as pd
import matplotlib.pyplot as plt
import sys
import os


class Generate(object):
    def __init__(self, i_file):
        self.f_input = i_file
        self.f_output = "results_clean.csv"
        self.data = []
        self.output = []

    def clean_file(self):
        with open(self.f_input, "r") as f:
            it = 0
            for line in f:
                if "EXIT" in line:
                    break
                elif line[:7].strip().isdigit() and int(line[:7].strip()) > it \
                        and not "T" in line:
                    if len(line) > 10:
                        #print line
                        it += 1
                        self.output.append(line)
        f.close()
        w = open(self.f_output, "w")
        w.writelines(self.output)
        w.close()
        print "MINOS_report: new file > {0}".format(self.f_output)

    def get_data(self):
        cols = [(1, 7), (8, 16), (17, 22), (23, 32), (33, 48), (49, 54), (55, 60), (61, 65)]
        self.data = pd.read_fwf(self.f_output, colspecs=cols, header=None)

    def plots(self):

        itn = len(self.data)
        rgradient = self.data[1]
        ninf = self.data[2]
        sinf = self.data[3]
        objective = self.data[4]
        nobj = self.data[5]
        ncon = self.data[6]
        nsb = self.data[7]

        it = list(range(1, itn + 1))

        _min_obj = objective[itn - 1]
        _max_nsb = nsb[itn - 1]
        _max_nobj = nobj[itn - 1]
        _sinf = sinf[itn - 1]
        _ninf = ninf[itn - 1]
        _ncon = ncon[itn - 1]

        plt.figure("AMPL results")

        # plot 1
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ax1 = plt.subplot(211)
        ax1.set_title('Iterations={0}. Objective (optimal solution)={1}'
                      .format(itn, _min_obj), size=10)

        ax1.set_xlabel('iterations')
        ax1.set_ylabel('objective function')
        l_obj = ax1.plot(it, objective, '.-', color="#0E3EEB", label="obj")
        ax2 = ax1.twinx()
        ax2.set_ylabel('reduced gradient')
        l_rg = ax2.plot(it, rgradient, '.-', color="#EB710E", label="rg")

        lns = l_obj + l_rg
        labs = [l.get_label() for l in lns]

        ax1.grid()
        ax1.legend(lns, labs, loc=0)

        # plot 2
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ax3 = plt.subplot(223)
        ax3.set_title('number of superbasics={0}\n # objective and gradient calculations={1}'
                      .format(_max_nsb, _max_nobj), size=10)

        l_sb = ax3.plot(it, nsb, '.-', color="#0E3EEB", label="nsb")
        ax3.set_xlabel('iterations')
        ax3.set_ylabel('no.superbasics')
        ax4 = ax3.twinx()
        l_nobj = ax4.plot(it, nobj, '.-', color="#EB710E", label="nobj")
        ax4.set_ylabel('calculations obj. & gradient')

        lns = l_sb + l_nobj
        labs = [l.get_label() for l in lns]

        ax3.grid()
        ax3.legend(lns, labs, loc=4)

        # plot 3
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        plt.subplot(224)
        plt.title('no. of infeasibilities={0}\n sum of infeasibilities={1}\n '
        '# nonlinear constraints evaluations={2}'.format(_ninf, _sinf, _ncon), size=10)
        plt.plot(it, ncon, 'g-', label="ncon")
        plt.plot(it, ninf, 'y-', label="ninf")
        plt.plot(it, sinf, 'm-', label="sinf")
        plt.xlabel('iterations')
        plt.grid()
        plt.legend()

        plt.tight_layout()
        plt.show()


if __name__ == "__main__":
    summary_file = sys.argv[1]
    #summary_file = "C:/Users/Guillermo Navas/Desktop/MSc Statistics and Operations research/Continuous optimization/" \
    #                   "Modelling and solving optimization problems/Lab assignment 2 _ Network flow problems/solver_summary.txt"
    report = Generate(summary_file)
    report.clean_file()
    report.get_data()
    report.plots()

creating shared library AMPL

2014-10-31T11:45:00.001+01:00

On this post I will try to explain in detail how create a shared library for AMPL. At this moment the Extended Function Library extends AMPL with over 300 functions, all of them available on the GNU Scientific library. Hence, before creating your own library you should check if your function is already included in that library.

I started my custom library because I needed the sigmoid function for one project and that function was not available in the Extended Function Library.

This process has the next steps:

Download headers: funcadd.h and stdio.1.h

Visual Studio: create a new Win32 project -> DLL.

Build DLL and testing.

The headers can be downloaded from the pages below:

funcadd.h: download file

stdio1.h: download file

Once this is done, create a new project -> Win32 project in Visual Studio and select application type DLL. Add both headers to Header Files and create a new file for your functions (e.g mathFuncs.c). Now your solution explorer should look like this:

Figure 1. Solution Explorer - Visual Studio

As said, we will create the sigmoid function and an approximation, the C code below shows one possible implementation, surely can be improved to gain efficiency. But before getting into the code, let's explain some general concepts about solvers and sigmoid functions. In general, the solvers need the first and second order derivatives of the function, depending on the method implemented by the solver (quasi-newton, newton's method, etc.), they may require the first derivative or both, as example MINOS requires just the first order derivative since it approximates the second one. In this .dll we will provide both derivatives of the functions.

Sigmoid function: $\frac{1}{(1 +\exp(-x))}$

First order derivative: $\frac{1}{(1 +\exp(-x))} (1 - \frac{1}{(1 +\exp(-x))}) = sigmoid(x)(1-sigmoid(x))$

Second order derivative: $\frac{-(\exp(x)(\exp(x) - 1))}{(\exp(x) + 1)^3}$

A possible sigmoid function approximation ,computationally faster, is the "fast sigmoid" function. The evaluation of the exponential function has a relevant cost and in many situations the exact sigmoid function may not be needed. This post is an interesting resource about the topic with a good benchmarking between methods link.

Fast sigmoid function: $\frac{x}{(1 + |x|)}$

First order derivative: $\frac{(-\frac{x^2}{|x|} + |x| + 1)}{|x|(1 + |x|)^2}$

Second order derivative: $\frac{3x^2|x| - 3|x|^3 - 3|x|^2 + x^2}{|x|^3(|x| + 1)^3}$

#include "stdlib.h"
#include "math.h"
#include "funcadd.h"    /* includes "stdio1.h" */

#define MAX_VAL 1.0e+300

#ifdef __cplusplus
/* The #ifdef __cplusplus lines allow this to be compiled
 * either by an ANSI C compiler or by a C++ compiler.
 */
extern "C" {
#endif

 static real
sigmoid(arglist *al)
 {
     /*  Sigmoid function: sigmoid(x) = 1/(1 + exp(-x)) */
     
     real x = al->ra[0];
     AmplExports *ae = al->AE;

     if (x > MAX_VAL || x < -MAX_VAL) 
     {
         fprintf(Stderr, "x = %d exceeded the maximum allowed value.", x);
         return -1;
     }
     else
     {
        if (al->derivs)
        {
            //first derivative: sigmoid gradient
            //  domain: R (all real numbers)
            //  range: {y c R: 0 <= y <= 0.25}

            *al->derivs = (1.0f/(1.0f + exp(-x)))*(1.0f -1.0f/(1.0f + exp(-x)));
            if (al -> hes)
            {
                //second derivative (hessian)
                //  domain: R (all real numbers)
                //  range: {y c R: -1/(6*sqrt(3)) <= 1/(6*sqrt(3))}

                *al->hes = -(exp(x) * (exp(x) - 1.0f)/ (pow(exp(x) + 1.0f, 3)));
            }
        }
        return 1.0f/(1.0f + exp(-x));
     }   
 }
 
 static real
fast_sigmoid(arglist * al)
 {
    /* fast sigmoid: sigmoid function approximation.
       fast_sigmoid(x) = x/(1 + abs(x))

       Domain {x c R: x != 0} -> derivative discontinuity at x = 0
     */
     
     real x = al->ra[0];
     if (al->derivs)
     {
         // first derivative (gradient)
         *al->derivs = (-pow(x,2)/fabs(x) + fabs(x) + 1) /(fabs(x) * pow(1 + fabs(x), 2));
         
         if (al -> hes)
         {
             // second derivative (hessian)
             real d = 3 * pow(x,2) * fabs(x) - 3 * pow(fabs(x), 3);
             real d1 = 3 * pow(fabs(x),2) + pow(x,2);
             *al->hes = (d  - d1) / (pow(fabs(x), 3) * pow(fabs(x) + 1, 3));
         }
     }
     return x / (1 + fabs(x));  
 }

 void
funcadd(AmplExports *ae){
/* Insert calls on addfunc here... */

/* Arg 3, called type, must satisfy 0 <= type <= 6:
 * type&1 == 0: 0,2,4,6 ==> force all arguments to be numeric.
 * type&1 == 1: 1,3,5   ==> pass both symbolic and numeric arguments.
 * type&6 == 0: 0,1 ==> the function is real valued.
 * type&6 == 2: 2,3 ==> the function is char * valued; static storage
                suffices: AMPL copies the return value.
 * type&6 == 4: 4,5 ==> the function is random (real valued).
 * type&6 == 6: 6   ==> random, real valued, pass nargs real args,
 *              0 <= nargs <= 2.
 *
 * Arg 4, called nargs, is interpretted as follows:
 *  >=  0 ==> the function has exactly nargs arguments
 *  <= -1 ==> the function has >= -(nargs+1) arguments.
 */

    addfunc("sigmoid", (rfunc) sigmoid, 0, 1, 0);
    addfunc("f_sigmoid", (rfunc) fast_sigmoid, 0, 1, 0);
}

#ifdef __cplusplus
    }
#endif

When creating a custom functions there are two important points, write the first order derivative (*al ->derivs) and if there is second derivative (if the solver's method requires the second order derivative for computing the Hessian ; al -> hes) then (*al->hes).

A last comment regarding funcadd(); It is important to define appropiate names to call your function (first argument in addfunc()), in this case I chose "sigmoid" and "f_sigmoid", keep in mind this, because otherwise AMPL will not recognize your functions.

Now it is time to build the .dll and do some testing with AMPL. Paste your .dll into your AMPL directory:

Figure 2. AMPL directory (student-edition)

Open ampl.exe or amplide.exe and type the instructions below to test your functions; if you followed these steps your output should be identical.

Figure 3. AMPL console

That's all from my side, I hope you find it useful.

Sparse matrix: from categorical matrix to binary matrix

2014-08-14T00:01:00.000+02:00

On this post I will show a mini project I have been working on for the last few days. The goal is generate boolean or binary data from categorical data. For doing so, I will primarly use Pandas and Numpy libraries. This class was created for testing purposes, so performance has not been taken into account and therefore is not relevant at this moment, however I am convinced that the code can be much more optimized.

Firstly, a short explanation about what this code is doing:

The next table represents the input data, as shown, fields p1, p2, p3 and p4 store categorical data, in this case a product code. Each module consists of four products (p) at most. Note: field "modules" is the index (row names).

modules	p1	p2	p3	p4	price
m_1	A1	A2			9.90
m_2	A1		A3		10.50
m_3	B1	B2	A4	C1	20.30
m_4	B1		A1	C1	22.10

Now, let's imagine we want to solve a linear equation system (Ax = b) composed by field p1, p2, p3 and p4 as matrix "A" , unknows (A1, A2, A3, A4, B1, B2, C1) and price as vector "b":

$A=\begin{bmatrix}a_{11} & a_{12} & \cdots & a_{1n} \\a_{21} & a_{22} & \cdots & a_{2n} \\\vdots & \vdots & \ddots & \vdots \\a_{m1} & a_{m2} & \cdots & a_{mn}\end{bmatrix}, \quad x=\begin{bmatrix}x_1 \\x_2 \\\vdots \\x_n\end{bmatrix}, \quad b=\begin{bmatrix}b_1 \\b_2 \\\vdots \\b_m\end{bmatrix}$

The matrix equations shows that size(A) = m x n, size(x) = n and size(b) = m. On the above example columns(A) is not equal to size(x). Hence, we need a new matrix A with the correct shape, see below.

$A=\begin{bmatrix}1 & 1 & 0 & 0 & 0 & 0 & 0 \\1 & 0 & 1 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 1 & 1 & 1 & 1 \\1 & 0 & 0 & 0 & 1 & 0 & 1 \end{bmatrix}, \quad x=\begin{bmatrix}A1 \\A2 \\A3 \\A4 \\B1 \\B2 \\C1 \end{bmatrix}, \quad b=\begin{bmatrix}9.90 \\10.50 \\20.30 \\22.10 \end{bmatrix}$

If we get a matrix "A" like the previous one, then we can solve the system easier.
Using vectorial notation, m_1 = A(1, :).x (e.g. Matlab indexing)

Before we continue, I would like to remark that the above system can be easily solvead by isolating x from the equation or by using LU decomposition, especially for square matrices. When dealing with non square matrices ($m \neq n$) the next methods are normaly used:

Normal equation: $x = (A^{T} A)^{-1} A^{T} b $

QR factorization

SVD decomposition

Code:

    import numpy as np
    import pandas as pd
    import warnings
    import matplotlib.pyplot as plt


    class CatBin(object):
        def __init__(self, i_file, o_file=False, cat_cols="all", other_cols=None, op_nan=True):
            """
            Class to convert a matrix with categorical data into a matrix with binary data.

            Example: Convert the below matrix into a binary matrix

                initial_data =
                modules   p1    p2      p3      p4     price
                m_1       A1    A2      nan     nan    9.90
                m_2       A1    nan     A3      nan    10.50
                m_3       B1    B2      A4      C1     20.30
                m_4       B1    nan     A1      C1     22.10

                binary_data =
                modules   A1  A2  A3  A4  B1  B2  C1    price
                m_1       1   1   0   0   0   0   0     9.90
                m_2       1   0   1   0   0   0   0     10.50
                m_3       0   0   0   1   1   1   1     20.30
                m_4       1   0   0   0   1   0   1     22.10

            Uses:
            Solve a linear equation system: Ax = b
                A = matrix[0,1]
                x = [A1, A2, A3, A4, B1, B2, C1]
                b = price

                Note:
                if non-square matrix, we can solve the system by using some 
                of the well-known methods:

                    Normal equation
                    QR factorization
                    SVD decomposition

            Graphs - adjacency matrix - connectivity
            ...

            :param i_file: input_file in csv format
            :param o_file: if True a csv is generated. The path must be provided.
            :param cat_cols: slicing - e.g [1, 10] -> from column 1 to 10(not included - pandas)
            :param other_cols:

            if the new data frame must preserve others columns besides those
            with categorical data, position and name should be provided in a dictionary.
                e.g
                others = {'i' : 'modules', 'e', 'prices'}
                if 'i', that column will be inserted on the first position in the new data frame
                    data.insert(0, 'name', column_data)
                if 'e', that column will be inserted after the rest of columns
                    data['name'] = column_data

            :param op_nan: if true all NaN values will be removed.
            """

            self.o_file = o_file
            self.cat_cols = cat_cols
            self.other_cols = other_cols
            self.nan = op_nan
            self.list = []

            self.rows = None
            self.columns = None
            self.p_sparsity = -1

            self.input_data = pd.read_csv(i_file)
            self.boolean_m = None
            self.boolean_a = None

            if self.cat_cols is "all":
                if self.nan is True:
                    # if all data is categorical data to be converted, automatically 
                    # others_cols must be false to avoid inconsistencies.
                    
                    self.data = self.input_data.fillna(value='nan')
                else:
                    self.data = self.input_data
                if other_cols is not None:
                    other_cols = None
                    warnings.warn("other_cols has been set to False to avoid inconsistencies")
            else:
                # cat_cols is an array [start, end]
                if isinstance(self.cat_cols, list):
                    if len(self.cat_cols) == 2:
                        _start = self.cat_cols[0]
                        _end = self.cat_cols[1]
                        if self.nan is True:
                            self.data = self.input_data.fillna(value='nan').iloc[:, _start: _end]
                        else:
                            self.data = self.input_data.iloc[:, _start, _end]
                    else:
                        raise ValueError("'start'and 'end' position must be provided")
                else:
                    raise TypeError("cat_cols is not a list")

            if not self.data.empty:
                # extract unique values:
                self.list = self.unique_values()
            else:
                raise ImportError("errors or empty csv file. Check input file")

        def unique_values(self):
            """ return sorted list of unique values """

            u_values = pd.Series(self.data.values.ravel()).unique()
            if self.nan is True:
                return np.sort(u_values[u_values != 'nan'])
            else:
                return np.sort(u_values)

        def boolean_matrix(self):
            matrix = []
            for i in range(len(self.data)):
                v = self.data.iloc[i, :].values
                # avoid worrying about the argument op_nan
                intersect = set(v).intersection(self.list)

                ea = np.zeros(len(self.list))
                for n in intersect:
                    for index, v in enumerate(self.list):
                        if n == v:
                            ea[index] = 1
                matrix.append(ea)

            self.rows = len(matrix)
            self.columns = len(matrix[0])
            self.boolean_m = pd.DataFrame(matrix, columns=self.list)

            return self.boolean_m

        def boolean_augmented_matrix(self):
            """
            The goal of this function is create the augmented matrix:
            Ax = b --> augmented_matrix = [A|b]

            However, this function allows the option of inserting at least 
            one first column and many others columns after the boolean data, 
            in case that was required to suit your necessities.
            """
            if self.cat_cols is not 'all':
                if self.boolean_m is None:
                    self.boolean_a = self.boolean_matrix()
                else:
                    self.boolean_a = self.boolean_m.copy()

                if self.other_cols is not None and isinstance(self.other_cols, dict):
                    for key, value in self.other_cols.items():
                        if key[0] is 'i':
                            self.boolean_a.insert(0, value, self.input_data[value])
                        elif key[0] is 'e':
                            self.boolean_a[value] = self.input_data[value]
                        else:
                            raise ValueError("Incorrect key format")
                else:
                    warnings.warn("insufficient arguments")
                return self.boolean_a
            else:
                raise ValueError("cat_cols = 'all' the aug. matrix is not possible")

        def sparsity(self):
            """
            The sparsity or density is the fraction of non-zero elements in a matrix.
            :return: ratio non-zero/total
            """
            if self.rows and self.columns and self.boolean_m is not None:
                total = self.rows * self.columns
                non_zeros = sum(self.boolean_m[self.boolean_m == 1].count())
                self.p_sparsity = non_zeros / total
                return "sparsity = {0:.4}%".format(self.p_sparsity * 100)
            else:
                raise ValueError("incorrect data shape. A boolean transformation is required")

        def generate_csv(self, matrix_type='normal'):
            """
            generate a csv file with boolean matrix or augmented matrix
            note: removes first column (index)
            """
            if matrix_type in ['normal', 'augmented']:
                if matrix_type is 'normal':
                    if self.boolean_m is not None:
                        self.boolean_m.to_csv(self.o_file, index=False)
                else:
                    if self.boolean_a is not None:
                        self.boolean_a.to_csv(self.o_file, index=False)
            else:
                raise ValueError("select: normal or augmented")

        @property
        def clean_data(self):
            return self.data

        @property
        def raw_data(self):
            return self.input_data

        @property
        def initial_size(self):
            """ Only includes those fields with categorical data"""
            rows = self.data.shape[0]
            columns = self.data.shape[1]
            return "rows={0}, columns={1}, total={2} elements".\
                format(rows, columns, rows * columns)

        @property
        def final_size(self):
            rows = self.boolean_m.shape[0]
            columns = self.boolean_m.shape[1]
            return "rows={0}, columns={1}, total={2} elements".\
                format(rows, columns, rows * columns)

The code tries to be self-explanatory, however I consider that three functions need more details:

boolean_matrix: This function contains the basic code to convert the input data. Basically, it compares each row from the initial data "A" with its unique values and inserts one on those positions where matching. As example, this task can be done in Matlab much easier, Matlab can compare two arrays with different length by using the function "ismember()":
```
        
          p = {'A1', 'A2', 'A3', 'A4', 'B1', 'B2', 'C1'}; % character arrays or cell arrays
          p_test = {'A1', 'B1','C1'};
          ismember(p, p_test)

          ans =

             1     0     0     0     1     0     1
        
```
boolean_augmented_matrix: The augmented matrix is that obtained by appending the solution vector "b" to matrix "A", usually denoted as [A b] or (A|b). Having a solvable linear equation system, then the augmented matrix [A b] has the same rank as A. Besides of creating the augmented matrix, this function allows the possibility to append other columns if needed.

Sparsity: The sparsity or density is the fraction of non-zero elements in matrix. A typical example of sparse matrix is the adjacency matrix of an undirected graph.

Below a short example using the class CatBin:

    from others.catbin import CatBin
    import time

    file = "test data/i_simple.csv"
    output_file = "test data/clean_matrix.csv"

    others = {'i': 'modules', 'e': 'price'}

    t1 = time.time()
    ct = CatBin(i_file=file, o_file=output_file, cat_cols=[1, 5], other_cols=others, op_nan=True)
    
    unique_v = ct.unique_values()
    boolean_m = ct.boolean_matrix()
    boolean_a = ct.boolean_augmented_matrix()
    sparsity = ct.sparsity()
    t = time.time() - t1
    
    print(t)
    print(ct.initial_size)
    print(ct.final_size)

    ct.generate_csv()

The file "i_simple.csv" is exactly the first table on this post. This code uses all features of the class, producing both boolean matrix and boolean augmented matrix and calculates the sparsity ratio. The most simple example needs 5 lines (without verbose).

    from others.catbin import CatBin
    
    file = "test data/i_simple.csv"
    output_file = "test data/clean_matrix.csv"

    ct = CatBin(i_file=file, o_file=output_file, cat_cols=[1, 5], 
                other_cols={'i': 'modules', 'e': 'price'}, op_nan=True)

    boolean_a = ct.boolean_augmented_matrix()

Results:

    
      modules  A1  A2  A3  A4  B1  B2  C1  price
    0     m_1   1   1   0   0   0   0   0    9.9
    1     m_2   1   0   1   0   0   0   0   10.5
    2     m_3   0   0   0   1   1   1   1   20.3
    3     m_4   1   0   0   0   1   0   1   22.1

    0.010005950927734375
    rows=4, columns=4, total=16 elements
    rows=4, columns=7, total=28 elements

In my laptop(4GB, i5) the process took 0.01 seconds. The next step will be to test something bigger. The dimensiones of the new data is 2194 x 11:

    1.1577680110931396
    rows=2194, columns=9, total=19746 elements
    rows=2194, columns=294, total=645036 elements

Due to the number of unique values is much higher, the number of elements is almost 650.000 and the process now took 1.15 seconds, surely this can be improved, but for now it works for my purposes. As shown, the new dimensiones will be multiplied by a factor X(unique values / initial columns).

Hence, having large amount of binary data in a database can also have an impact on performance, so I will not recommend to store the new dataset in a DB better create a system file. There are other options of storing a sparse matrix much more efficient that save the complete data in a DB (CSR or CSC, for instance).

One more example: This dataset includes extra columns before and after the fields with categorical data (2194 x 13). If the order of the extra columns matters is convenient to use ordered dictionaries.

    from others.catbin import CatBin
    from collections import OrderedDict
    import time

    file = "test data/i_data_2.csv"
    output_file = "test data/clean_matrix.csv"

    others = {'i1': 'modules', 'i2': 'modules2', 'e1': 'price', 'e2': 'price2'}
    ordered_others = OrderedDict(sorted(others.items(), key=lambda t: t[0]))

    t1 = time.time()
    
    ct = CatBin(i_file=file, o_file=output_file, cat_cols=[2, 11],
                other_cols=others, op_nan=True)
    
    boolean_m = ct.boolean_matrix()
    boolean_a = ct.boolean_augmented_matrix()
    
    t = time.time() - t1
    
    sparsity = ct.sparsity()
    print(t)
    print(ct.initial_size)
    print(ct.final_size)
    print(boolean_a)

Finally, a function for plotting:

    def plot_sparse(self, ticks=False, x_ticks=False, y_ticks=False, y_label=None,
                    sparsity=False):
        """
        Notes: 
            if sparsity was previously calculated - this figure can be included.
            if data is extraordinarily large might be convenient to hide ticks.
        """
        if sparsity and self.p_sparsity >= 0:
            plt.title("sparsity ratio = {0:.4}%".format(self.p_sparsity * 100), fontsize=10)
            plt.suptitle("boolean matrix", fontsize=14, fontweight='bold')
        else:
            plt.title("boolean matrix", fontsize=14, fontweight='bold')
        plt.imshow(self.boolean_m, interpolation='nearest', cmap='gist_yarg')

        if ticks is True:
            plt.yticks(range(self.boolean_m.shape[0]))
            if x_ticks is True:
                plt.xticks(np.arange(len(self.list)), self.list)
            if y_label is not None and y_ticks is True:
                plt.yticks(np.arange(len(self.input_data[y_label])), self.input_data[y_label])
        else:
            plt.axis('off')

        plt.show()

Let's plot both matrices used on the previous examples:

    ct.plot_sparse(ticks=True, x_ticks=True, sparsity=True)

fig.1: small matrix

    ct.plot_sparse(ticks=False, sparsity=True)

fig.2: medium matrix

    ct.plot_sparse(ticks=True, x_ticks=True, sparsity=True)

fig.3: medium matrix - detail

That's all for this post, I hope you find it useful.

TSP / ATSP algorithm.

2014-03-29T20:17:00.000+01:00

This post will be dedicated to the Travelling salesman problem (TSP), one of the most known combinatorial optimization problem. The goal of this problem is simple: A salesman has to visit a number of places (nodes). Which route would visits them all in the shortest total distance (time, cost, etc...)
The simplest way to solve the problem is list and try all possible routes and them select the shortest one. The problem here is the computational time which will grow faster than exponential (Stirling's formula: $n! \approx \sqrt{2\pi n}(\frac{n}{e})^n$)

There are other solutions which can reduce this computational time, for example dynamic programming which solves the problem in exponential time E. The first implemented method using that more efficient method was the Held-Karp algorithm which finds the shortest route in $2^n n^2$
Other methods are: k-opt, Greedy algorithm, Greedy k-opt, Genetic algorithm, Simulating annealing, neuronal network approach among others.

On this post I will show an implementation for the simplest and naive method, an implementation for dealing with asymmetric input matrices and finally we will compare these methods with the genetic algorithm(this implementation has too many lines of code for one post)

import itertools
import numpy as np

INF = 1000000


class TSPExact:
    def __init__(self, matrix, redundant=False):
        """
        input: Adjacency matrix N x N
        len(nodes) = shape(input)
        @param matrix: N x N matrix
        @param redundant: false by default (n-1)!
        """
        # 1: Check if the adjacency matrix is symmetric
        self.cost_matrix = np.array(matrix)
        self.symmetry = None
        if self.symmetric() is False:
            self.symmetry = False
            self.as_symmetric()

        self.size = len(self.cost_matrix)
        self.nodes = range(self.size)
        self.all_tours = itertools.permutations

        if redundant is False:
            self.tours = self.non_redundant()
        else:
            self.tours = list(self.all_tours(self.nodes))

    def as_symmetric(self):
        """
        Reformulate an asymmetric TSP as a symmetric TSP: 
        "Jonker and Volgenant 1983"
        This is possible by doubling the number of nodes. For each city a dummy 
        node is added: (a, b, c) => (a, a', b, b', c, c')

        distance = "value"
        distance (for each pair of dummy nodes and pair of nodes is INF)
        distance (for each pair node and its dummy node is -INF)
        ------------------------------------------------------------------------
          |A'   |B'   |C'   |A    |B    |C    |
        A'|0    |INF  |INF  |-INF |dBA  |dCA  |
        B'|INF  |0    |INF  |dAB  |-INF |dCB  |
        C'|INF  |INF  |0    |dAC  |dBC  |-INF |
        A |-INF |dAB  |dAC  |0    |INF  |INF  |
        B |dBA  |-INF |dBC  |INF  |0    |INF  |
        C |dCA  |dCB  |-INF |INF  |INF  |0    |

        For large matrix an exact solution is infeasible
        if N > 5 (N = N*2 > 10) then use other techniques: 
        Heuristics and relaxation methods
        @return: new symmetric matrix

        [INF][A.T]
        [A  ][INF]
        """

        shape = len(self.cost_matrix)
        mat = np.identity(shape) * - INF + self.cost_matrix

        new_shape = shape * 2
        new_matrix = np.ones((new_shape, new_shape)) * INF
        np.fill_diagonal(new_matrix, 0)

        # insert new matrices
        new_matrix[shape:new_shape, :shape] = mat
        new_matrix[:shape, shape:new_shape] = mat.T
        # new cost matrix after transformation
        self.cost_matrix = new_matrix

    def total_tour(self, tour):
        total = sum(self.cost_matrix[tour[node], tour[(node + 1) % self.size]] for node in self.nodes)
        return total, tour

    def shortest(self):
        min_tour = min(self.total_tour(tour) for tour in self.tours)
        if self.symmetry is False:
            min_tour_real = min_tour[0] + self.size / 2 * INF
            min_cycle_real = min_tour[1][0::2]
            return min_tour_real, min_cycle_real
        else:
            return min_tour

    def non_redundant(self):
        start_node = self.nodes[0]
        new_set = list(self.nodes)[1:]
        return [list(np.append(start_node, tour)) for tour in self.all_tours(new_set)]

    def symmetric(self):
        #@return: boolean symmetric | asymmetric
        mat = [tuple(row) for row in self.cost_matrix]
        return mat == zip(*mat)

    @property
    def tour_iterations(self):
        return list(self.tours)

    @property
    def tour_size(self):
        return self.size

The code above tried to be self explanatory, but it might be some points unclear. The trickiest part is the transformation from asymmetric matrix to symmetric matrix. As shown in the code, we are using the implementation suggested by Jonker R, Volgenant T (1983) (“Transforming asymmetric into symmetric traveling salesman problems.” Operations Research Letters, 2, 161–163). The process behind is simple, doubling the number of nodes and transposing the original matrix it is possible get a symmetric matrix which can be solvable applying the current known methods. Finally this algorithm requires the removal of all the dummy nodes and a subtle operation in order to get the real shortest tour value.

$shortestRoute (real) = shortestRoute + \frac{Size(InputMatrix)}{2} \times INF.Value $
We can use slicing in Numpy for removing dummy nodes.

min_tour_real = min_tour[0] + self.size / 2 * INF
min_cycle_real = min_tour[1][0::2]

Unit testing

Now it's time to compare some algorithms. I have hard coded three different matrices for testing purpose. We must keep in mind that the matrix size will be multiplied by two if input matrix is asymmetric. In addition, one must point out that the basic algorithm becomes infeasible when the number of nodes is higher than 12 (6 - asymmetric matrix)

from Algorithms.scripts.TSP_Algorithm import *
import unittest
import time

MATRIX = [[0, 750, 43, 402], [750, 0, 736, 671], [43, 736, 0, 117], [402, 671, 117, 0]]
MATRIX5 = [[0, 750, 43, 402, 133], [750, 0, 736, 671, 85], [43, 736, 0, 117, 934], 
     [402, 671, 117, 0, 982], [133, 85, 934, 982, 0]]

MATRIX8 = [[0, 750, 431, 402, 133, 489, 310, 566], [750, 0, 736, 671, 85, 174, 870, 927],
           [431, 736, 0, 117, 934, 65, 939, 305], [402, 671, 117, 0, 982, 727, 911, 870],
           [133, 85, 934, 982, 0, 956, 834, 118], [489, 174, 65, 727, 956, 0, 322, 397],
           [310, 870, 939, 911, 834, 322, 0, 127], [566, 927, 305, 870, 118, 397, 127, 0]]


class TestAlgorithms(unittest.TestCase):
    @staticmethod
    def test_exact():
        start_time = time.time()
        tsp = TSPExact(MATRIX5)
        total_value, tour = tsp.shortest()
        print "__exact__"
        print total_value
        print tour
        print "Time (s)", time.time() - start_time
  
    @staticmethod
    def test_genetic():
        start_time = time.time()
        tsp = TSPGenetic(MATRIX5)
        total, route = tsp.calculate()
        print "__genetic__"
        print total
        print route
        print "Time (s)", time.time() - start_time

Some results

MATRIX(Asymmetric)
__exact__
1580.0
[0, 2, 3, 1]
Time (s) 0.130000114441
__genetic__
1580
[0, 2, 3, 1]
Time (s) 0.00100016593933

MATRIX5(Asymmetric)
__exact__
1049.0
[0, 2, 3, 1, 4]
Time (s) 10.1160001755
__genetic__
1049
[0, 2, 3, 1, 4]
Time (s) 0.000999927520752

MATRIX8(Asymmetric)
__exact__
16 nodes - time...
__genetic__
1398
[0, 3, 2, 5, 1, 4, 7, 6]
Time (s) 0.00399994850159

Comments

As shown on the previous tests, the problem becomes intractable for the Exact algorithm as the number of nodes increases. Using the Genetic algorithm we got the same results, however this algorithm might not obtain the exact solution in large graphs.

Python and Excel/VBA

2014-02-23T19:39:00.001+01:00

Today we are going to talk about how easy it could be run Python's scripts from Excel, which might be considered pretty interesting. Basically we will create a VBA code which launches a console running the selected script. This script creates an output which must be printed on the same excel file from where it was launched.

Continuing with the previous post, we will use the scraping script to print the extracted data into excel. First we will change a bit the previous code and we will create a class called "collect". UPDATE: The website has changed and the previous code doesn't work any more!

import numpy as np
import urllib2
import urlparse
import re
from bs4 import BeautifulSoup


class collect:

    def __init__(self, mainSites):
        self.sites = mainSites

    def connect(self, site):
        hdr = {'User-Agent': 'Mozilla/5.0 (X11; Linux x86_64)',
                'Accept': 'text/html,application/xhtml+xml,application/xml'
                ';q=0.9,*/*;q=0.8',
               'Accept-Charset': 'ISO-8859-1,utf-8;q=0.7,*;q=0.3',
               'Accept-Encoding': 'none',
               'Accept-Language': 'en-US,en;q=0.8',
               'Connection': 'keep-alive'}
        req = urllib2.Request(site, headers=hdr)
        page = urllib2.urlopen(req).read()
        soup = BeautifulSoup(page)
        return soup

    def extractUrls(self):
        FxUrls = []
        C1 = []
        C2 = []
        for mainSite in self.sites:
            soup = self.connect(mainSite)
            for link in soup.find_all('td', {'class': 'bold left noWrap'}):
                path = 'http://www.investing.com/currencies/'
                url = link.find(href=re.compile(path)).get('href')
                FxUrls.append(url)
                curr = urlparse.urlparse(url).path.split("/")[2].split("-")
                C1.append(curr[0])
                C2.append(curr[1])
        return np.vstack((FxUrls, C1, C2))

    def extractData(self):
        Value = []
        Data = self.extractUrls()
        for link in Data[0, :]:
            soup = self.connect(link)
            for span in soup.find_all('span', {'id': 'last_last'}):
                value = (span.text).replace(',', '')
                Value.append(value)
        return np.vstack((Data, Value))

Once we have the script is time to create the interaction with Excel. In order to do this operation we will use an excellent module: pyvot. I have previously prepared some special ranges for printing the data (I recommend Name ranges as a good Excel/VBA practice). The script below will create a connection with Excel, just selecting the desired ranges and sending the values. This script will also need some arguments, the workbook name from where the script is launched and our links.

import sys
import scraping
import xl

#arguments:
excelWb = sys.argv[1]
print "Excel workbook: ", excelWb

mainSites = []
for i in range(len(sys.argv) - 2):
    mainSites.append(str(sys.argv[i + 2]))
print "urls: ", mainSites

scr = scraping.collect(mainSites)
data = scr.extractData()

rngCurr1 = data[1, :]
rngCurr2 = data[2, :]
rngValue = data[3, :]

"""Excel connection"""
wb = xl.Workbook(excelWb)
wb.range("rngCurr1").set(rngCurr1)
wb.range("rngCurr2").set(rngCurr2)
wb.range("rngRate").set(rngValue)

Excel template

Ranges

Finally, we need a VBA code which can run a script through the Windows console. First we will be calling a function to clear ranges contents and the main function afterwards. The main function, called runScript, creates an object Shell and runs the script with its arguments, which are: workbook name and array of urls.

Option Explicit

Sub clean()

Dim rngForex As Range
Set rngForex = Sheet1.Range("rngForex")
rngForex.ClearContents

End Sub

Sub runScript()

Dim oShell As Object
Dim sScript As String, sScriptName As String, sWbName As String
Dim sArguments As String
Dim iRegions As Integer, i As Integer
Dim aRegion() As String, sLink As String

Call clean

Set oShell = CreateObject("WScript.Shell")

sWbName = ThisWorkbook.Name
sScriptName = "scrapingtest.py"
sLink = "http://www.investing.com/currencies/"

iRegions = Application.WorksheetFunction.CountA(Sheet1.Range("rngRegions"))
ReDim aRegion(iRegions - 1)
For i = 1 To iRegions
    aRegion(i - 1) = sLink & CStr(Sheet1.Range("rngRegions").Cells(i, 1).Value)
    sArguments = sArguments & aRegion(i - 1) & " "
Next i

sScript = sScriptName & " " & sWbName & " " & sArguments
oShell.Run ("cmd.exe /S /C " & sScript)

End Sub

Now, if you click "RUN" a console will appear showing a message similar to this one:

Conclusion:
My personal opinion is that at this moment there are not real alternatives to Excel. Is true that there are open source tools like LibreOffice or Apache OpenOffice but these are not even comparable to Excel. It is quite common create scripts in C# to automate processes which are not so fast in VBA. VS languages (VB /C#) are chosen because they are much more integrated with Excel and generally speaking with most of the Microsoft applications.

The aim of this post was show that we can use Python to accelerate the code and simply use Excel as a template where print our results, creating dashboards, reports and so on. Thus, for me Python and Excel/VBA is a good association and I will use it going forward.

Networks 1: Scraping + Data visualization + Graph stats

2014-01-25T00:54:00.000+01:00

These last weeks I have been reading about networks and optimization algorithms, I think is an interesting field with many applications, so my idea was write a new article (or series of articles) showing roughly how use some interesting python libraries like Networkx, for instance. In this article I also like to make a brief introduction to Scraping (a very good way to get data for your graphs).

And only to whet your appetite:

Figure 1. Undirected graph.

Web Scraping:

The first step to construct graphs like the previous one is get the data. As example, I wanted to visualize the whole foreign exchange market and the connections between each currency (Foreign exchange cross rates), but there are another good examples like user connections in social networks (Twitter, Facebook etc).
For this case study(Forex), there are several available websites where you can find financial data for free, nonetheless this data will be restricted to certain time periods, so do not expect 10 or 30 seconds candles for free.

List of Forex websites:

Once we have selected our source, we can think about how extract the info. Some websites like Yahoo Finance allow downloading end-of-day equities or forex data, normally in csv format. But let's image that our intention is to obtain the most possible current data. Here Scraping is the appropriate technique. Web scraping is the process of automatically collecting information from websites, usually simulating human interaction through a full-fledged web browser. Note: web scraping may be against the terms of use of several websites.

For this purpose, I used the next python modules:

urllib2

urlparse

BeautifulSoup

The goal then is extract all the foreign exchange cross rates. This data will be collected in three numpy arrays and finally we will create a dataframe with pandas to export a csv file. The tricky part of all this is parse the html file, however BeautifulSoup allows using regular expressions, which combined with urlparse creates a powerful and simple tool. So let's take a look at the script.

Script:

import numpy as np
import pandas as pd
import urllib2
import urlparse
import re
from bs4 import BeautifulSoup

#------------------------------------------------------------------------------
#Functions:
#------------------------------------------------------------------------------


def add_exchange(Data):
    Value = []
    for link in Data[0, :]:
        soup = connect(link)
        for span in soup.find_all('span', {'id': 'last_last'}):
            #remove thousands (",")
            value = (span.text).replace(',', '')
            Value.append(value)
    return np.vstack((Data, Value))


def connect(site):
    hdr = {'User-Agent': 'Mozilla/5.0 (X11; Linux x86_64)',
       'Accept': 'text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8',
       'Accept-Charset': 'ISO-8859-1,utf-8;q=0.7,*;q=0.3',
       'Accept-Encoding': 'none',
       'Accept-Language': 'en-US,en;q=0.8',
       'Connection': 'keep-alive'}
    req = urllib2.Request(site, headers=hdr)
    page = urllib2.urlopen(req).read()
    soup = BeautifulSoup(page)
    return soup

#------------------------------------------------------------------------------
#Extract list of currencies + Matrix
#------------------------------------------------------------------------------
mainSites = {'http://www.investing.com/currencies/middle-east',
    'http://www.investing.com/currencies/north-america',
    'http://www.investing.com/currencies/central-america',
    'http://www.investing.com/currencies/south-america',
    'http://www.investing.com/currencies/asia',
    'http://www.investing.com/currencies/pacific',
    'http://www.investing.com/currencies/europe',
    'http://www.investing.com/currencies/middle-east',
    'http://www.investing.com/currencies/africa',
    'http://www.investing.com/currencies/caribbean'}

FxUrls = []
C1 = []
C2 = []
for mainSite in mainSites:
    soup = connect(mainSite)
    for link in soup.find_all(href=re.compile("www.investing.com/currencies/")):
        #extract currencies:
        url = link.get('href')
        FxUrls.append(url)
        currencies = urlparse.urlparse(url).path.split("/")[2].split("-")
        C1.append(currencies[0])
        C2.append(currencies[1])

#Create Matrix
Forex = np.vstack((FxUrls, C1, C2))
#------------------------------------------------------------------------------

#Extract all foreign exchange cross rates:
Import = add_exchange(Forex)

#Pandas dataframe
data = {'currency_1': Import[1, :],
    'currency_2': Import[2, :],
    'value': Import[3, :]}

df = pd.DataFrame(data)
df.to_csv('Forex.csv', index=False, cols=['currency_1', 'currency_2', 'value'])

Summarizing what we have done in this script:

I wrote two functions, one in order to connect with website and capture the html page, and the other one to parse this html and select the required information. There are different ways to extract the data, select the one that best suits the structure of your html file. Run the script, just take into account that this script may take a while depending on your internet connection.
With this tool, you might consider a dynamic system to obtain frequently data. Well, that would be complicated and unlikely, your IP address would be blocked if the number of requests is considerable.

Networkx - Graphs:

Now we have collected the data, it's time to show what we got. For the creation and study of the graphs we will use the python package called Networkx. This is a pretty helpful package which provides many functions and algorithm that facilitate the study of complex networks. However, this may be not the best package to draw your graph, although we can use it as an interface to use Graphviz.

The first thing I did was write a new module, called "graphs", with the Networkx's functionalities I needed. This includes two classes, "create_graph"and "graphStats". The first one reads the produced csv file and creates the two elements that form a graph, nodes and edges. In this case we also have weights, the exchange rates (third column), so we need a list in order to create the labels, which looks like this:

example: {('D', 'E'):9.0, ('F', 'N'):6.0}...

The class "graphStats" includes two functions in order to calculate the degree centrality and betweenness centrality, we will talk about this in more detail later on.

Graph class:

import networkx as nx
from operator import itemgetters
import pandas as pd
import numpy as np


class create_graph:

    def __init__(self):
        data = 'csv/Forex.csv'
        headers = ['start', 'end', 'weights']
        self.Network = pd.read_csv(data, skiprows=1, names=headers)
        self.graph = np.asarray(self.Network[['start', 'end']])
        self.nodes = np.unique(self.graph)
        self.weights = list(map(float, self.Network['weights']))

    def networkList(self):
        G = nx.DiGraph()
        for node in self.nodes:
            G.add_node(node)

        self.graph_l = []
        for edge in self.graph:
            G.add_edge(edge[0], edge[1])
            self.graph_l.append((edge[0], edge[1]))

        labels = dict(list(zip(self.graph_l, self.weights)))
        return G, labels


class graphStats:

    def calculate_degree_centrality(self, graph):
        dgc_key = []
        dgc_value = []
        g = graph
        dc = nx.degree_centrality(g)
        nx.set_node_attributes(g, 'degree_cent', dc)
        degcent_sorted = sorted(dc.items(), key=itemgetter(1), reverse=True)
        for key, value in degcent_sorted:
            dgc_key.append(str(key))
            dgc_value.append(value)
        return dgc_key, dgc_value

    def calculate_betweenness_centrality(self, graph):
        btc_key = []
        btc_value = []
        g = graph
        bc = nx.betweenness_centrality(g)
        betcent_sorted = sorted(bc.items(), key=itemgetter(1), reverse=True)
        for key, value in betcent_sorted:
            btc_key.append(str(key))
            btc_value.append(value)
        return btc_key, btc_value

Now we have prepared the module, we will use matplotlib to draw the graph. Actually, this is not a huge network, it has only 133 nodes and 1918 edges, but it is difficult get a good resolution if your output file is .jpg or png, so I recommend a pdf file and some adjustments in the figure size.

Graph - plot:

import networkx as nx
import matplotlib.pyplot as plt
import graphs

#------------------------------------------------------------------------------
# Read data + create graph
#------------------------------------------------------------------------------
gr = graphs.create_graph()
graph = gr.graph
G, labels = gr.networkList()

#------------------------------------------------------------------------------
# Graph plot:
#------------------------------------------------------------------------------
fig = plt.figure()
#fig = plt.figure(figsize=(8,8), dpi=1000) - pdf file
pos = nx.graphviz_layout(G, prog='neato')
node_color = [float(G.degree(v)) for v in G]
nx.draw_networkx_nodes(G, pos, node_size=[float(G.degree(v)) * 2.5 for v in G],
             alpha=0.6, node_color=node_color)
nx.draw_networkx_edges(G, pos, alpha=0.25, edge_color='blue', width=0.5,
arrows=False)
nx.draw_networkx_labels(G, pos, font_size=2.5)

plt.axis('off')
plt.show()
#plt.savefig("directedGraph.pdf")

Figure 1.1. Undirected graph (detail).

Figure 2. Directed graph.

The second graph is which shows the real Forex network. The Forex network is a directed graph, a graph where its edges have a direction associated with them. This digraph is strongly connected due to it contains several pairs of vertices with both paths, (euro -> dolar and dolar -> euro, for instance).

Graph - stats:

Networkx is a great python package to study graphs, its main purpose. In this point we are going to look into some statistics which will help us to extract real information from our graphs with the aim of interpret the results properly.

1) Betweenness centrality:

Betweenness centrality quantifies the number of times a node acts as bridge along the shortest path between two other nodes. It is basically a measure of connectedness between components of the graph.

2) Degree centrality:

The centrality for a node measures its relative importance within the graph. The degree centrality for a node (n) is the number of nodes it is connected. The degree centrality values are normalized by the maximum degree in the graph. For a graph with self loops (a digraph) the maximum degree might be higher than 1 and the degree centrality as well.

In this article I will not expand the theory behind these measures.

import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
import graphs


def topTable(field1, field2, n_top):
    topM = max(field2) * 0.9
    right = len(field1) * 0.75
    plt.text(right, topM * 1.08, 'Top %s' % n_top, fontsize=12)
    for i in range(n_top):
        curr = field1[i]
        val = field2[i]
        plt.text(right, topM - i * topM / 20, '{}) {} = {}'.format(i + 1,
        curr.upper(), round(val, 3)), fontsize=8)

#------------------------------------------------------------------------------
# Read data + create graph
#------------------------------------------------------------------------------
gr = graphs.create_graph()
G, labels = gr.networkList()

#------------------------------------------------------------------------------
# Stats + plots:
#------------------------------------------------------------------------------
stats = graphs.graphStats()
betc_key, betc_value = stats.calculate_betweenness_centrality(G)
degc_key, degc_value = stats.calculate_degree_centrality(G)

#Average degree
N = G.order()
K = G.size()
avg_d = float(N) / K
avg_degree = 'Average degree: %.4f ' % (avg_d)


# Plot: Degree_centrality
plt.figure()

ax1 = plt.subplot(211)
plt.title('Degree centrality for nodes', fontsize=12)
a_lenght = np.arange(len(degc_value))
plt.bar(a_lenght, degc_value, color=cm.jet(degc_value), align='center')
plt.xticks(a_lenght, degc_key, size='small', rotation='vertical')
plt.tick_params(axis='x', labelsize=4)
plt.tick_params(axis='y', labelsize=8)
plt.autoscale(enable=True, axis='both', tight=None)

#Top degree centrality:
topTable(degc_key, degc_value, 10)
plt.text(len(degc_value) * 0.75, max(degc_value) * 0.4, avg_degree,
bbox={'facecolor': 'blue', 'alpha': 0.25, 'pad': 10}, fontsize=7)

# Plot: Betweenness_centrality
plt.subplot(212)
plt.title('Betweenness centrality for nodes', fontsize=12)
a_lenght = np.arange(len(betc_value))
plt.bar(a_lenght, betc_value, color=cm.jet(betc_value), align='center')
plt.xticks(a_lenght, betc_key, size='small', rotation='vertical')
plt.tick_params(axis='x', labelsize=4)
plt.tick_params(axis='y', labelsize=8)
plt.autoscale(enable=True, axis='both', tight=None)
plt.ylim(0, max(betc_value) * 1.1)
plt.plot(betc_value, '--b')

#Top degree betweenness:
topTable(betc_key, betc_value, 10)

plt.show()

The above plot shows the currencies which have more relevance into the network, and as we suspected these are: USD dolar, Euro and the pound sterling. We can make at least two significant considerations; The average degree centrality is quite low, what shows that there are currencies with a notable grade of popularity(Degree centrality : in-degree) and these same currencies play a broker role in the network due to cumulate most of trajectories for the shortest paths. For this particular example,this might mean that we can create a more simplistic network, which would represent pretty faithfully the Forex market, excluding those currencies with fewer number of connections.

Adjacency Matrix:

The adjacency matrix is a N x N matrix where the cells constitute edges. In this example the fields are coloured according the degree centrality and as the graph is directed, the matrix is not symmetric.

import numpy as np
import matplotlib.pyplot as plt
import graphs


gr = graphs.create_graph()
G = gr.networkList()[0]

n = len(gr.nodes)
oMapp = dict(list(zip(gr.nodes, np.arange(n))))

Mtx_w = (np.ones((n, n))) * np.inf
for i in range(len(gr.graph)):
    start = oMapp.get(gr.Network['start'][i])
    end = oMapp.get(gr.Network['end'][i])
    w = gr.Network['weights'][i]
    Mtx_w[start, end] = w

Mtx_nodes = []
Mtx_nodes = np.array([np.concatenate([Mtx_nodes, gr.nodes]) for _ in range(n)])
#-----------------------------------------------------------------------------

# Show matrix by colour:
Mtx_Image = Mtx_w.copy()
for i in range(n):
    count = ((Mtx_Image[i, :] != np.inf) & (Mtx_Image[i, :] != np.nan)).sum()
    Mtx_Image[i, :][Mtx_Image[i, :] != np.inf] = count

# Gradient based on degrees
plt.figure()
plt.title('Foreign Exchange cross rates', fontsize=16)
plt.imshow(Mtx_Image, interpolation='nearest')
plt.xticks(np.arange(len(gr.nodes)), gr.nodes, rotation=90)
plt.yticks(np.arange(len(gr.nodes)), gr.nodes)
plt.tick_params(axis='both', labelsize=4)

plt.tight_layout()
plt.show()
#plt.savefig("Forex_Matrix.pdf")

This first article related to networks has been just a brief introduction on this field. We will continue exploring graphs and graph algorithms in next posts. Finally, I would like to make some notes regarding the network images, these images have been modified reversing the RGB output (I preferred a black background), I used GIMP for that task.

Numerical methods: ODE and Finite difference method

2014-01-05T02:02:00.000+01:00

The goal of this post is show how solve an ordinary differential equation, numerically and using the finite difference method and compare the result with the analytic solution. The problem statement is as follow:

Numerically solve the ODE: $m\frac{\partial^{2}y}{\partial t^{2}} = -b\frac{\partial y}{\partial t} $ that represents an object whose velocity at $t=0$ is equal to $v_0$ and it is slowing due to a retardant force, which is proportional to its velocity.

The integration interval is $0 \leq t \leq 5$ seconds. For simplicity, we take $m = b = 1$.
Boundary conditions:

$y(0) = 0 m$

$y(5) = 1 m$

Analytic solution:
$Y(t) = \frac{m v_0}{b} (1 - e^{-\frac{b}{m}t})$ where $v_0 = 1$

Ordinary differential equation:
We will start solving the ordinary differential equation of second grade . I will explain the resolution step by step:

Step 1: Characteristic equation

Step 2: Finding constants of integration

$m\frac{\partial^{2}y}{\partial t^{2}} = -b\frac{\partial y}{\partial t} \Rightarrow mD^{2} = -bD \Rightarrow D(mD + b) = 0 $
Thus, we have two possible solutions: $s1 = 0; s2 = -\frac{b}{m}$. As these two roots are real and distinct ($s1 \ne s2$) we can use the general solution for this specific case:

$Y(x) = C1 e^{s_1 x} + C2 e^{s_2 x} \Rightarrow Y(t) = C1 e^{0 t} + C2 e^{-\frac{b}{m} t} = C1+C2 e^{-\frac{b}{m} t}$

Using the boundary conditions we can obtain the constants of integration C1 and C2:
$Y(t=0) \Rightarrow C1 = -C2$
$Y(t=5) \Rightarrow -C2 + C2e^{-\frac{b}{m} 5 } = 1 \Rightarrow (-1 + e^{-\frac{b}{m} 5 } ) C2 = 1 \Rightarrow C2 = \frac{1}{(-1 + e^{-\frac{b}{m} 5 } ) } $
Solution:
$Y(t) = -\frac{1}{(-1 + e^{-\frac{b}{m} 5 } ) } + \frac{1}{(-1 + e^{-\frac{b}{m} 5 } ) } e^{-\frac{b}{m} t}$

$m = b = 1 \Rightarrow Y(t) = 1.00678 - 1.00678 e^{-t}$

Finite difference method:
We use the centered difference at time t to obtain the central difference approximation in both equation sides.
\[y" = \frac{\partial^{2}y}{\partial t^{2}} = \frac{U_{i+1} - 2U_i + U_ {i-1}}{h^{2}} \enspace y' = \frac{\partial y}{\partial t} = \frac{U_{i+1}-U_{i-1}}{2h} \]
Now we can sustitude these above approximations in the initial expression:
\[\frac{m}{h^{2}}(U_{i+1} - 2U_i + U_ {i-1}) = -\frac{b}{2h}(U_{i+1}-U_{i-1})\] Isolate $U_i$:
$U_{i+1} - 2U_i + U_ {i-1} = -\frac{b}{2h}(U_{i+1}-U_{i-1})\frac{h^{2}}{m} \Rightarrow [-\frac{bh}{2m}(U_{i+1}-U_{i-1})-(U_{i+1}+U_{i-1})] \frac{1}{-2} = U_i$

Implementation:
For this simple example, we create a Python class called "ODE", which includes the analytic solution and the basic methodology based on the iterative scheme.

Iterations:

def solve(self, n, steps, m, b, h):
    v = np.zeros(n)
    #Boundary conditions
    #y(t = 0) = 0
    v[0] = 0
    #y(t = 5) = 1
    v[n - 1] = 1

    #Iterations
    for i in range(steps):
        for i in range(1, n - 1):
            a1 = (v[i - 1] - v[i + 1]) * (b * h / (2 * m))
            a2 = v[i + 1] + v[i - 1]
            v[i] = (a1 - a2) / (-2)
    return v

As shown, the first value in the array $= 0$, first boundary condition (we start the loop in the position 1) and the last value = 1 (position n - 1). Finally, the function returns the complete array.

The attached chart shows how the result is approximated by finite difference as the number of iterations increases.

Code:

import numpy as np
import matplotlib.pyplot as plt


class ODE(object):

    def __init__(self, m, b, v0, t_min, t_max, n):
        self.t = np.linspace(t_min, t_max, n)
        #Analytic solution
        self.s = m * v0 / b * (1 - np.exp(- b / m * self.t))
        #ODE
        C2 = 1 / (- 1 + np.exp(- b / m * t_max))
        C1 = -C2
        self.f = C1 + C2 * np.exp(- b / m * self.t)

    def solve(self, n, steps, m, b, h):
        v = np.zeros(n)
        #Boundary conditions
        #y(t = 0) = 0
        v[0] = 0
        #y(t = 5) = 1
        v[n - 1] = 1

        #Iterations
        for i in range(steps):
            for i in range(1, n - 1):
                a1 = (v[i - 1] - v[i + 1]) * (b * h / (2 * m))
                a2 = v[i + 1] + v[i - 1]
                v[i] = (a1 - a2) / (-2)
        return v

#inputs and solution:
m = 1
b = 1
v0 = b / m
t_min = 0.0
t_max = 5.0
n = 100
h = t_max / n

vel = ODE(m, b, v0, t_min, t_max, n)
plt.plot(vel.t, vel.f, 'b.-', label='ODE')
plt.plot(vel.t, vel.s, 'g.-', label='Analytic solution')

for i in range(250, 3000, 500):
    sol_i = vel.solve(n, i, m, b, h)
    plt.plot(vel.t, sol_i, '-.', label=r'Iteration = %i' % i)

plt.ylabel('Distance (m)')
plt.xlabel('Time (s)')

plt.legend(loc=2, prop={'size': 9})
plt.show()

Weibull analysis

2013-11-30T00:46:00.000+01:00

This first post in ComputSimu will explain the basic methodology to calculate Weibull distribution given a data set, in this case a data set provided by a three point bending test for specimens with circular cross section. The Weibull distribution is one of the most used distributions to model life data, due to its extreme flexibility to fit a wide range of data (Normal distribution or logarithmic distribution, for instance) and its applicability modelling different problems (weather forecasting, failure analysis, delivery times, etc...)

For creating the plot and calculate the Weibull distribution parameters we will use a python script (www.python.org). Python is widely extended as scripting language in computational science with an impressive support, which guarantees a constant number of new modules and improvements. Schema:

Three point bending test

Weibull distribution

Characteristic life

Weibull distribution characteristics

Extra: Gamma function

Code

1) Three point bending test
This test consists in an applied force (F) acting at the center of a beam to assess the fracture resistance. Required formulas:

Beam deflection (double pinned contact):
$Moment(max) =\dfrac{F}{2}\dfrac{L}{2} = \dfrac{FL}{4}$ where $L = span$ (mm) and $F = load$ (N)
$Stress(max) = \dfrac{M_t}{\dfrac{I_p}{r}} = \dfrac{M_t}{W_p}$ where $r = radius$ (mm); $W_p = W_y = W_z = \dfrac{I_y}{r} = \dfrac{\dfrac{\pi r^4}{4}}{r} = \dfrac{\pi r^3}{4}$
$Stress(max) = \dfrac{Moment(max)}{Resisting Bending Moment} = \dfrac{\dfrac{FL}{4}}{\dfrac{\pi r^3}{4}} = \dfrac{FL}{\pi r^3}$ (MPa)
We will use Numpy to import the data (there are many options, but this works pretty well). The completed code is on the bottom of the post, but I will explain quickly the most important parts. The first step is import the data into an array. These values are ordered ascendant, which is a must in order to calculate the failure probability, but you can also order your array once the values are imported. Example:

F(N)
333
440
440
...

 path = '../Data_Test.csv'
 aSamples = np.genfromtxt(path, delimiter=',', skip_header=1, dtype='f')

Now, we need two functions to create two new arrays with failure probabilities and the failure stresses. To calculate the probability we can choose one of the two methods below depending on the number of samples.

def FailureProb(n, N):
    if(N < 100):
        #Bernard's approximation
        return (n - 0.3) / (N + 0.4)
    else:
        #Mean Ranks
        return n / (N + 1)


def MaxStress(F, L, r):
    # Maximum Stress (cylinder 3 points)
    return (F * L) / (math.pi * r ** 3)

aSamples = np.genfromtxt(path, delimiter=',', skip_header=1, dtype='f')

aFailureProb = aSamples.copy()
aFailureStress = aSamples.copy()

for i in range(len(aFailureProb)):
    aFailureProb[i] = FailureProb((i + 1), float(len(aSamples) + 1))
    aFailureStress[i] = MaxStress(aSamples[i], L, r)

2) Weibull distribution
The probability density function of a random variable is: \[f(x; \alpha, \beta) = \begin{cases} \frac{\alpha}{\beta} (\frac{x}{\alpha})^{(\alpha - 1)} e^{-(\frac{x}{\beta})^{\alpha}} & \text{if } x \geq 0 \\ 0 & \text{if } x < 0 \end{cases}\]
and the cumulative distribution function is: \[F(x; \alpha, \beta) = \begin{cases} 1-e^{-(\frac{x}{\beta})^{\alpha}} & \text{if } x \geq 0 \\ 0 & \text{if } x < 0 \end{cases}\]
Shape parameter ($\alpha$):
$\alpha$ shows how the failure rate develops in time. We can describe three different situations for the parameter $\alpha$:

$\alpha$ < 1: This first situation is characterized by a decreasing failure rate or hazard function.

$\alpha$ = 1: The second situation is characterized by a constant failure rate (independent of time), this is called lack of memory. In this situation Weibull distribution is identical to the exponential distribution.

$\alpha$ > 1: The third situation is characterized by an increasing failure rate. This represents that the more older is the "element" the more likely it is to fail.

Scale parameter or characteristic life ($\beta$):
$\beta$ is the time/cycle/sample at which 63,2% are expected to fail. This is a constant for all Weibull distribution, regardless the shape parameter.

Plot:

Script:

"""
Weibull:
    1) Probability density function
    2) Cumulative distribution function
"""

import numpy as np
import matplotlib.pyplot as plt
from weibull import weibProbDist, weibCumDist

x = np.linspace(0, 5, 100)
a_alpha = np.array([0.5, 1, 2, 5])
label = r'${\alpha} = %s, {\beta} = %s$'

fig = plt.figure('Weibull')
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
for i in a_alpha:
    ax1.plot(x, weibProbDist(x, i, 1), label=label % (i, 1))
    ax2.plot(x, weibCumDist(x, i, 1), label=label % (i, 1))
ax1.set_title('Probability density function', size=12)
ax2.set_title('Cumulative distribution function', size=12)
ax1.legend(loc=1, prop={'size': 10})
ax2.legend(loc=7, prop={'size': 10})

fig.subplots_adjust(hspace=0.5)
plt.show()

3) Characteristic life
To obtain the Weibull parameters the first step is linearise the cumulative distribution functions: $F(x) = 1-e^{-(\frac{x}{\beta})^{\alpha}} \Rightarrow P_f = 1-e^{-(\frac{\sigma}{\sigma_0})^{m}}$ where m = slope and $\sigma$ = scale (characteristic stress)

$P_s = 1-P_f$
$ln(P_s^{-1}) = ln(e^{(\frac{\sigma}{\sigma_0})^{m}}) \Rightarrow ln[ln(P_s^{-1})] = m \enspace ln(\frac{\sigma}{\sigma_0}) \Rightarrow ln[ln(P_s^{-1})] = m \enspace ln(\sigma) - m \enspace ln(\sigma_0)$

Linear expression:
$\underbrace{ln[ln(P_s^{-1})]}_\text{f(x)} = \underbrace{ m \enspace ln(\sigma)}_\text{mx} - \underbrace{ m \enspace ln(\sigma_0)}_\text{b}$
Maximum likelihood estimation (MLE) - Linear regression:
$m (slope)= \dfrac{Cov(ln(\sigma), ln[ln(P_s^{-1})])}{Var(ln(\sigma))} = \dfrac{\dfrac{1}{n} \sum_{i=1}^{n}(ln(\sigma) - \overline{ln(\sigma)}) (ln[ln(P_s^{-1})] - \overline{ln[ln(P_s^{-1})])}}{\dfrac{1}{n} \sum_{i=1}^{n}(ln(\sigma) - \overline{ln(\sigma)})^2}$
$b (intercept)= \overline{ln[ln(P_s^{-1})])} - m \enspace \overline{ln(\sigma)}$In order to solve this simple equation we use Scipy and its module called stats which apply the described methodology above.

 #Survival probability = 1 - Failure Probability
 aSurvivalProb = 1 - aFailureProb

 #Linear regression: parameters (m, b)
 aLnPs = np.log(- np.log(aSurvivalProb))
 aLnFS = np.log(aFailureStress / 1000)

 slope, intercept, r_square, s, std_error = stats.linregress(aLnFS, aLnPs)
 regressLine = slope * aLnPs + intercept

 #Characteristic life : Stress
 Stress = np.exp(- (intercept / slope)) * 1000
 Pr_Stress = weibCumDist(Stress, slope, Stress)

Finally, we can calculate the parameter $\beta$ with the next expression:
$\beta, \sigma = e^{-(\frac{intercept(b)}{slope(m)})}$

Plotting analysis results:
As described above, the final plot shows how a stress = 125,08 MPa is the characteristic stress whereby the failure probability is equal to 63,2%. The plots below try to be self-explanatory.

4) Weibull distribution characteristics
You will find all characteristic parameters in the Weibull module. These parameters are exported into a .csv file.
$Mean = \beta \enspace \Gamma(1+\dfrac{1}{\alpha})$
$Median = \beta \enspace \ln(2)^{(\frac{1}{\alpha})}$
$Mode = \begin{cases} \beta \enspace (\frac{\alpha - 1}{\alpha})^{\frac{1}{\alpha}} & \text{if } \alpha > 1 \\ 0 & \text{if } \alpha = 1 \end{cases}$
$Variance = \beta^2 \enspace \Gamma(1+\dfrac{2}{\alpha}) - \mu^2$
$Skewness = \beta^2[\Gamma(1+\dfrac{3}{\alpha})+2(\Gamma(1+\dfrac{1}{\alpha}))^3 - 3[\Gamma(1+\dfrac{2}{\beta})][\Gamma(1+\dfrac{1}{\beta})]]$
Skewness

$\text{if } \alpha < 3.6: Skewness > 0 $
$\text{if }\alpha = 3.6: Skewness = 0 $
$\text{if }\alpha > 3.6: Skewness < 0 $

5) Extra: Gamma function

The Weibull module also include one of the possible approximations for the Gamma function, Lanczos approximation. You can find more information about different implementations on the links below:
C++ library + pseudocode
"Calculators and the Gamma Function", Viktor T. Toth (Recommended)

6) Code

import numpy as np
import matplotlib.pyplot as plt
import math
from scipy import stats
from weibull import *
import csv

#inputs:
L, r = float(40), float(4)
path = '../Data_Test.csv'
pathReport = '../WeibullReport.csv'


def FailureProb(n, N):
    if(N < 100):
        #Bernard's approximation
        return (n - 0.3) / (N + 0.4)
    else:
        #Mean Ranks
        return n / (N + 1)


def MaxStress(F, L, r):
    # Maximum Stress (cylinder 3 points)
    return (F * L) / (math.pi * r ** 3)

aSamples = np.genfromtxt(path, delimiter=',', skip_header=1, dtype='f')

aFailureProb = aSamples.copy()
aFailureStress = aSamples.copy()

for i in range(len(aFailureProb)):
    aFailureProb[i] = FailureProb((i + 1), float(len(aSamples) + 1))
    aFailureStress[i] = MaxStress(aSamples[i], L, r)

MeanForces = np.mean(aSamples)
MeanStress = np.mean(aFailureStress)

#Survival probability = 1 - Failure Probability
aSurvivalProb = 1 - aFailureProb

#Linear regression: parameters (m, b)
aLnPs = np.log(- np.log(aSurvivalProb))
aLnFS = np.log(aFailureStress / 1000)

slope, intercept, r_square, s, std_error = stats.linregress(aLnFS, aLnPs)
regressLine = slope * aLnPs + intercept

m = slope
b = intercept

#Characteristic life : Stress
Stress = np.exp(- (intercept / slope)) * 1000
Pr_Stress = weibCumDist(Stress, slope, Stress)

# export results to .csv
Results = [
            ['Parameters', 'Value', 'Units'],
            ['shape (alpha)', slope],
            ['scale (beta)', Stress, 'MPa'],
            ['Mean', weibMean(slope, Stress), 'MPa'],
            ['Median', weibMean(slope, Stress), 'MPa'],
            ['Mode', weibMode(slope, Stress), 'MPa'],
            ['Variance', weibVar(slope, Stress), 'MPa^2'],
            ['Standard Deviation', np.sqrt(weibVar(slope, Stress)), 'MPa'],
            ['Skewness', weibSkew(slope, Stress)]
          ]

with open(pathReport, 'wb') as f:
    writer = csv.writer(f)
    writer.writerows(Results)

#Plots:
fig = plt.figure('Weibull')
ax1 = fig.add_subplot(211)
scatt = ax1.scatter(aLnFS, aLnPs, label='samples', c='DarkTurquoise', s=15)
regres = ax1.plot(aLnPs, regressLine, 'b-', label='Linear regression (MLE)')
ax1.set_ylabel('$Ln[-Ln(Ps)]$')
ax1.set_xlabel('$Ln(\sigma) : Ln(GPa)$')
ax1.set_title('Weibull analysis')
ax1.set_ylim(min(aLnPs) * 1.1, max(aLnPs) * 1.25)
ax1.set_xlim(min(aLnFS) * 1.1, max(aLnFS) * 0.9)
ax1.legend(loc=2, prop={'size': 10})
ax1.text(min(aLnFS) * 1.05, -1, r'$y=%.4f m + %.4f$' % (m, b))
ax1.text(min(aLnFS) * 1.05, -1.5, r'$R^2=%.4f$' % r_square)

x = np.linspace(0, round(Stress * 2, 0), 100)
ax2 = fig.add_subplot(212)
ax2.set_title('Weibull distribution')
ax2.plot(x, weibProbDist(x, slope, Stress), label='Prob. Density function')
ax2.fill_between(x, weibProbDist(x, slope, Stress), 0, color='0.90')
ax2.set_ylabel(r'$f(x, {\alpha}, {\beta})$')
ax2.set_xlabel(r'${\sigma} (MPa)$')
ax3 = ax2.twinx()
ax3.plot(x, weibCumDist(x, slope, Stress), 'r', label='Cum. distr. function')
ax3.fill_between(x, weibCumDist(x, slope, Stress), alpha=0.15, color='r')
ax3.set_ylim(0, 1)
ax3.set_ylabel(r'$F(x) = 1-exp(-(x/{\beta})^{\alpha})$')
ax3.legend(loc=7, prop={'size': 10})
ax2.legend(loc=2, prop={'size': 10})

#Text:
ax3.text(5, 0.7, r'${\alpha} = %.4f, {\beta} = %.4f$' % (m, Stress))
ax3.text(5, 0.6, r'$F(%s) = %s$' % (round(Stress, 3), round(Pr_Stress, 3)),
color='m')

#Arrows:
ax3.arrow(Stress, 0, 0, Pr_Stress - 0.1, head_width=2, head_length=0.1,
ec='m', fc='m')
ax3.arrow(Stress, Pr_Stress, round(Stress * 2, 0) - Stress - 10, 0,
head_width=0.02, head_length=10, ec='m', fc='m')

fig.subplots_adjust(hspace=0.5)
plt.show()

Weibull module:

"""
Weibull distribution:
    1) Probability density function
    2) Cumulative distribution function
    3) Mean
    4) Median
    5) Mode
    6) Variance - Standard deviation
    7) Skewness
    8) Gamma function (Lanczos approximation)
"""

import numpy as np
import cmath as cm


def weibProbDist(x, a, b):
    return (a / b) * (x / b) ** (a - 1) * np.exp(-(x / b) ** a)


def weibCumDist(x, a, b):
    return 1 - np.exp(-(x / b) ** a)


def weibMean(a, b):
    return b * gamma(1 + 1 / a)


def weibMedian(a, b):
    return b * np.log(2) ** (1 / a)


def weibMode(a, b):
    if a > 1:
        return b * ((a - 1) / a) ** (1 / a)
    else:
        return 0


def weibVar(a, b):
    return (b ** 2) * gamma(1 + 2 / a) - weibMean(a, b) ** 2


def weibSkew(a, b):
    ax_1 = gamma(1 + 3 / a) + 2 * (gamma(1 + 1 / a) ** 3)
    ax_2 = 3 * (gamma(1 + 2 / a) * gamma(1 + 1 / a))
    return (b ** 3) * (ax_1 - ax_2)

"""
Gamma function:

Lanczos approximation implementation.
Numerical Recipes in C (2nd ed. Cambridge University Press 1992)

"""
g = 5

p0 = float(1.000000000190015)
p1 = float(76.18009172947146)
p2 = float(-86.50532032941677)
p3 = float(24.01409824083091)
p4 = float(-1.231739572450155)
p5 = float(1.208650973866179e-3)
p6 = float(-5.395239384953e-6)

Pn_Coef = np.array([p0, p1, p2, p3, p4, p5, p6])


def gamma(z):
    z = complex(z)
    if z.real < 0.5:
        #Recursion method:
        return cm.pi / (cm.sin(cm.pi * z) * gamma(1 - z))
    else:
        ax_1 = cm.sqrt(cm.pi * 2) / z
        ax2_sum = Pn_Coef[0]
        for i in range(1, g + 2):
            ax2_sum += Pn_Coef[i] / (z + i)
        t = z + g + 0.5
        ax_3 = (t ** (z + 0.5)) * cm.exp(-t)
        return ((ax_1 * ax2_sum) * ax_3).real