creating shared library AMPL

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.

TSP / ATSP algorithm.

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.