Genetic algorithms

Got very little time on this weekend. I was again interested in techniques for solving optimisation problems. Many of the problems in real-life are optimisation problems, where unknows are more than the information (equations) at hand. Examples: fitting of surface patches / curves on points, closest path - path planning computations, Nesting (manufacturing) for optimising the material cut from the raw material sheets, inverse kinematics, fluid mechanics, medicine, many other areas...

Optimisation techniques is a huge research area, and it would require considerable amount of time to study every technique. I have used basic optimization methods before, but I am starting to learn about advanced techniques. Many problems can be formulated and solved by calculus methods, but they become exceedingly complicated with more number of unknowns. A simple minima, maxima problem in single variable can be solved  by calculus, where f`(x)=0... then solving by Newton Raphson method, other numerical methods, but when there are multiple variables we have to find partial derivatives which become complex with more variables. Some of the calculus approaches can only find out the local maximas for the functions, in that case optimization results heavily depend on the initial guesses. There are advanced numerical methods to solve non-linear optimization problems to get global roots.

GA uses a heuristic approach based on natural selection phenomenon. Got reading this book again about GA by David Goldberg, I am fascinated by the GA technique. Its a searching technique which narrows down to the answer in solution (search) space using different combination functions on the existing pool of intermediate solutions which they call as chromosomes. It can be used for problems that are NP-Hard as well. Most important part of this techniques is the problem formulation, we have to represent the problem in a particular form so that we can use GA. To get familiar with this technique I want to try and solve a problem with this, possibly Travelling salesman problem or best fit circle in points or both.

Useful links:
http://hyperphysics.phy-astr.gsu.edu/hbase/math/maxmin.html
http://en.wikipedia.org/wiki/Nonlinear_programming
http://www.ai-junkie.com/ga/intro/gat1.html

More later...

Quickhull algo in javascript - pt4

Finally got the hang of OOP in javascript.. the public functions for the classes are prototype methods.

The quickhull works! Don't use the code for some serious stuff.. its come out of casual coding and might have some errors. Some escape sequences and optimizations remain, but as for the javascript experiment it looks good enough.. I had to put quite a bit of code in for even this simple algo, CVector, CPoint, CEdge and CHull classes. Perhaps a reusable good js geometry util can be written for general usage :-)
More later..

Have a go at it below. Enjoy :-) [Only works on Google Chrome]

    /* This program demonstrates the creating of Convex Hull of points
    * using the Quickhull algorithm 
    * Copyright (C) 2013 Rishikesh Parkhe
    *
    * This program is free software: you can redistribute it and/or modify
    * it under the terms of the GNU General Public License as published by
    * the Free Software Foundation, either version 3 of the License, or
    * (at your option) any later version.
    *
    * This program is distributed in the hope that it will be useful,
    * but WITHOUT ANY WARRANTY; without even the implied warranty of
    * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    * GNU General Public License for more details.
    *
    * You should have received a copy of the GNU General Public License
    * along with this program. If not, see .
    */

    /*
    ******************************************
    * Point class
    *******************************************
    */
    function CPoint(x, y) {
        this.x = x;
        this.y = y;

        this.inside = false;
    }

    CPoint.prototype.GetX = function () {
        return this.x;
    }

    CPoint.prototype.GetY = function () {
        return this.y;
    }

    CPoint.prototype.SetInside = function () {
        this.inside = true;
    }

    CPoint.prototype.GetInside = function () {
        return this.inside;
    }

    /* 
    ******************************************
    * Convex hull class
    *******************************************
    */
    function CHull() {
        this.edges = new Array();
        this.numOfEdges = 0;

        this.ptsInside = new Array();
        this.numInside = 0;

        this.ptsOutside = new Array();
        this.numOutside = 0;

        //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        // Add edge function
        //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        this.AddEdge = function (edge) {
            edge.SetNext(this.numOfEdges + 1);
            this.edges[this.numOfEdges] = edge;
            this.numOfEdges++;
        }

        //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        // Split edge function
        //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        this.SplitEdge = function (index, point) {
            var newArray = new Array();
            var numEdges = 0;

            for (var i = 0; i < index; i++) {
                newArray[numEdges] = this.edges[i];
                numEdges++;
            }

            ///////////////////////////////////////////////////////////////////
            // Actual split of edge
            var splt1 = new CEdge(this.edges[index].GetA(), point);
            splt1.SetNext(index + 1);
            newArray[numEdges] = splt1;
            numEdges++;

            var splt2 = new CEdge(point, this.edges[index].GetB());
            splt2.SetNext(index + 2);
            newArray[numEdges] = splt2;
            numEdges++;

            ///////////////////////////////////////////////////////////////////

            if (index + 1 < this.numOfEdges) {
                for (var i = index + 1; i < this.numOfEdges; i++) {
                    newArray[numEdges] = this.edges[i];
                    numEdges++;
                }
            }

            this.edges = newArray;
            this.numOfEdges = this.numOfEdges + 1;

            if (this.numOfEdges == 2) {
                var closure = new CEdge(this.edges[1].GetB(), this.edges[0].GetA());
                closure.SetNext(0);
                newArray[numEdges] = closure;
                numEdges++;
                this.numOfEdges = this.numOfEdges + 1;
            }

            this.edges = newArray;
        }


        //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        // Render function
        //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        this.Render = function (ctx) {
            ctx.fillStyle = '#f00';
            ctx.strokeStyle = '#f00';
            ctx.lineWidth = 2;
            ctx.beginPath();

            for (var i = 0; i < this.numOfEdges; i++) {
                ctx.lineTo(this.edges[i].GetA().GetX(), this.edges[i].GetA().GetY());
                ctx.lineTo(this.edges[i].GetB().GetX(), this.edges[i].GetB().GetY());
            }
            ctx.stroke();
            ctx.closePath();

            ctx.fillStyle = "rgba(255, 255, 255, 0.8)";
            for (var i = 0; i < this.numOfEdges; i++) {
                ctx.fillRect(this.edges[i].GetB().GetX() - 1, this.edges[i].GetB().GetY() - 1, 2, 2);
            }
        }

        //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        // Compute hull function
        //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        this.ComputeHull = function (points, numPts, minIndex, maxIndex) {
            this.edges = new Array();
            this.numOfEdges = 0;

            if (numPts >= 2) {

                var edge = new CEdge(points[minIndex], points[maxIndex]);
                points[minIndex].SetInside();
                points[maxIndex].SetInside();

                this.AddEdge(edge);

                if (numPts > 2) {
                    var i = 0;
                    var count = 0;
                    while (i < this.numOfEdges || count > 1000) {
                        // compute farthest point to this edge
                        var index = this.edges[i].CalcFarthestPoint(points, numPts);
                        var liesAbove = this.edges[i].IsPointAboveEdge(points[index]);

                        if (liesAbove) {
                            // construct triangle exculde points inside triangle
                            // split edge
                            this.SplitEdge(i, points[index]);
                            this.MarkPointsInsideNewTri(this.edges[i].GetA(), this.edges[i].GetB(), points[index], points, numPts);

                            points[index].SetInside();
                            i = 0; // do again..
                        }
                        else i++;

                        count++;
                    }
                }
            }
        }

        this.MarkPointsInsideNewTri = function (A, B, C, points, numPts) {
            var tri = new CTriangle(A, B, C);
            for (var i = 0; i < numPts; i++) {
                if (points[i].GetInside() == false) {   // only for points outside the hull
                    var inside = tri.CheckPtIsInside(points[i]);
                    if (inside == true) points[i].SetInside();
                }
            }
        }
    }


    /*
    ******************************************
    * Edge class
    *******************************************
    */
    function CEdge(ptA, ptB) {
        this.A = ptA;
        this.B = ptB;
        this.index = -1;

        // Getter and setter for next index
        this.SetNext = function (index) {
            this.index = index;
        }
        this.GetNext = function () {
            return this.index;
        }

    }

    // Getter and setter for ptA
    CEdge.prototype.SetA = function (pt) {
        this.A = pt;
    }
    CEdge.prototype.GetA = function () {
        return this.A;
    }

    // Getter and setter for ptB
    CEdge.prototype.SetB = function (pt) {
        this.B = pt;
    }
    CEdge.prototype.GetB = function () {
        return this.B;
    }

    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Check if point is above the edge, using Right hand rule
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    CEdge.prototype.IsPointAboveEdge = function point_above_line(point) {
        var cross = ((this.B.GetX() - this.A.GetX()) * (point.GetY() - this.A.GetY())
                        - (this.B.GetY() - this.A.GetY()) * (point.GetX() - this.A.GetX())
                    );

        if (cross > 0) // right hand rule
            return true;

        return false;
    }

    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Compute distance of a point to the edge
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    CEdge.prototype.CalcDistToPoint = function (pt) {
        var numer =
            (this.A.GetY() - this.B.GetY()) * pt.GetX() +
            (this.B.GetX() - this.A.GetX()) * pt.GetY() +
            (this.A.GetX() * this.B.GetY() - this.B.GetX() * this.A.GetY());

        var edge = new CVector(0, 0);
        edge = edge.Create(this.A, this.B);
        var length = edge.Mag();

        var dist = -99999999999.0;
        if (length != 0)
            dist = numer / length;

        return dist;
    }

    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Compute farthest point to the edge
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    CEdge.prototype.CalcFarthestPoint = function (pts, numPts) {
        if (numPts > 0) {
            var max_index = -1;
            var max_dist = -99999999999;

            for (var i = 0; i < numPts; i++) {
                //if (pts[i].GetInside() == false) 
                {  // for all outside points
                    var dist = this.CalcDistToPoint(pts[i]);

                    if (dist > max_dist) {
                        max_index = i;
                        max_dist = dist;
                    }
                }
            }
        }
        return max_index;
    }

    /*
    ******************************************
    * Vector class
    *******************************************
    */
    function CVector(x, y) {
        this.x = x;
        this.y = y;

        this.Create = function (ptA, ptB) {
            this.x = ptB.x - ptA.x;
            this.y = ptB.y - ptA.y;

            return this;
        }
    }

    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Compute magnitude of this vector
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    CVector.prototype.Mag = function () {
        return Math.sqrt(this.x * this.x + this.y * this.y);
    }

    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Normalise this vector
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    CVector.prototype.Normalise = function () {
        var mag = this.Mag();
        if (mag != 0) {
            this.x = this.x / mag;
            this.y = this.y / mag;
        }
    }

    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Dot product of two vectors
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    CVector.prototype.Dot = function (other) {
        var product = this.x * other.x + this.y * other.y;
        return product;
    }

    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Compute angle between two vectors
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    CVector.prototype.ComputeAngle = function (other) {
        var nv1 = new CVector(this.x, this.y);
        nv1.Normalise();
        var nv2 = other.Normalise();

        var cosTheta = nv1.Dot(nv2);
        var angle = Math.acos(cosTheta);

        return angle;
    }


    /*
    ******************************************
    * Triangle class
    *******************************************
    */
    function CTriangle(ptA, ptB, ptC) {
        this.A = ptA;
        this.B = ptB;
        this.C = ptC;

        //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        // Check if point lies inside triangle
        //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        this.CheckPtIsInside = function(pt) {

            var e0 = new CVector();
            e0.Create(this.A, this.C);

            var e1 = new CVector();
            e1.Create(this.A, this.B);

            var ap = new CVector();
            ap.Create(this.A, pt);

            var dot00 = e0.Dot(e0) // mag of edge 0
            var dot01 = e0.Dot(e1) // projection of edge1 on edge 0
            var dot02 = e0.Dot(ap) // projection of ap on edge 0
            var dot11 = e1.Dot(e1) // mag of edge 1
            var dot12 = e1.Dot(ap) // projection of ap on edge1

            // Compute barycentric 
            invDenom = 1 / (dot00 * dot11 - dot01 * dot01)
            u = (dot11 * dot02 - dot01 * dot12) * invDenom
            v = (dot00 * dot12 - dot01 * dot02) * invDenom

            // Check if point lies inside triangle by checking the two projections
            return (u >= 0) && (v >= 0) && (u + v < 1)
        }
    }

    var ctx;
    var m_num_pts = 0;
    var m_num_hull_pts = 0;
    var m_pts = new Array();
    var m_hull = new CHull();

    var mouse_x = 0;
    var mouse_y = 0;

    var m_minIndex = 0;
    var m_maxIndex = 0;

    var m_AABB = {
        "min": {
            "x": 0,
            "y": 0
        },
        "max": {
            "x": 0,
            "y": 0
        }
    };

    function init() {
        var canvas = document.getElementById("canvas");
        ctx = canvas.getContext("2d");
        canvas.addEventListener('mousemove', mouse_move_handler, false);
        canvas.addEventListener('mouseup', mouse_up_handler, false);

        draw();
    }

    function draw() {
        ctx.fillStyle = "rgba(0, 10, 30, 1)";
        ctx.fillRect(0, 0, 400, 400);

        ctx.fillStyle = "rgba(200, 100, 10, 0.8)";
        ctx.fillRect(mouse_x - 10, mouse_y - 10, 10, 10);

        if (m_num_pts > 1) {
            ctx.moveTo(m_pts[0].GetX(), m_pts[0].GetY());
            ctx.fillStyle = '#00f';
            ctx.strokeStyle = "rgba(230, 220, 250, 0.2)";
            ctx.lineWidth = 1;
            ctx.beginPath();
            for (var i = 0; i < m_num_pts; i++) {
                ctx.lineTo(m_pts[i].GetX(), m_pts[i].GetY());
            }
            ctx.stroke();
            ctx.closePath();
        }

        ctx.fillStyle = "rgba(250, 20, 50, 0.5)";
        for (var i = 0; i < m_num_pts; i++) {
            ctx.fillRect(m_pts[i].GetX() - 1, m_pts[i].GetY() - 1, 2, 2);
        }

        if (m_num_pts > 0) {
            ctx.fillStyle = "rgba(255, 170, 0, 1)";

            draw_ABB();
            m_hull.Render(ctx);
        }

        ctx.fillStyle = "rgba(170, 170, 170, 0.4)";
        ctx.fillText("RISHIKESH PARKHE", 10, 370);
        ctx.fillText("Quick Hull", 10, 385);
    }

    function mouse_move_handler(evt) {
        // Get the mouse position relative to the canvas element.
        if (evt.layerX || evt.layerX == 0) { // for Firefox
            mouse_x = evt.layerX;
            mouse_y = evt.layerY;
        } else if (evt.offsetX || evt.offsetX == 0) { // for Opera
            mouse_x = evt.offsetX;
            mouse_y = evt.offsetY;
        }

        draw();
    }

    function mouse_up_handler(evt) {
        mouse_move_handler(evt);
        var pt = new CPoint(mouse_x, mouse_y);
        m_pts[m_num_pts] = pt;
        m_num_pts++;

        compute_quickhull();
        draw();
    }

    function compute_quickhull() {
        // compute the bounding box
        compute_AABB();

        // store points lying inside the polygon
        computeMinMax(m_pts, m_num_pts);
        m_hull.ComputeHull(m_pts, m_num_pts, m_minIndex, m_maxIndex);
    }

    function computeMinMax(pts, numPts) {
        for (var i = 0; i < m_num_pts; i++) {
            m_AABB.min.x = m_pts[0].GetX();
            m_AABB.max.x = m_pts[0].GetX();
            m_AABB.min.y = m_pts[0].GetY();
            m_AABB.max.y = m_pts[0].GetY();

            for (var i = 0; i < m_num_pts; i++) {
                if (m_pts[i].x < m_AABB.min.x) {
                    m_AABB.min.x = m_pts[i].GetX();
                    m_minIndex = i;
                }

                if (m_pts[i].x > m_AABB.max.x) {
                    m_AABB.max.x = m_pts[i].GetX();
                    m_maxIndex = i;
                }

                if (m_pts[i].y < m_AABB.min.y) {
                    m_AABB.min.y = m_pts[i].GetY();
                    m_minIndex = i;
                }

                if (m_pts[i].y > m_AABB.max.y) {
                    m_AABB.max.y = m_pts[i].GetY();
                    m_maxIndex = i;
                }
            }
        }
    }

        

    function compute_AABB() {
        if (m_num_pts > 0) {

            m_AABB.min.x = m_pts[0].GetX();
            m_AABB.max.x = m_pts[0].GetX();
            m_AABB.min.y = m_pts[0].GetY();
            m_AABB.max.y = m_pts[0].GetY();

            for (var i = 0; i < m_num_pts; i++) {
                if (m_pts[i].x < m_AABB.min.x) {
                    m_AABB.min.x = m_pts[i].GetX();
                }

                if (m_pts[i].x > m_AABB.max.x) {
                    m_AABB.max.x = m_pts[i].GetX();
                }

                if (m_pts[i].y < m_AABB.min.y) {
                    m_AABB.min.y = m_pts[i].GetY();
                }

                if (m_pts[i].y > m_AABB.max.y) {
                    m_AABB.max.y = m_pts[i].GetY();
                }
            }
        }
    }


    function draw_ABB() {

        ctx.moveTo(m_AABB.min.x, m_AABB.min.y);
        ctx.fillStyle = '#00f';
        ctx.strokeStyle = "rgba(0, 100, 22, 0.2)";;
        ctx.lineWidth = 1;
        ctx.beginPath();
        ctx.lineTo(m_AABB.min.x, m_AABB.min.y);
        ctx.lineTo(m_AABB.min.x, m_AABB.max.y);
        ctx.lineTo(m_AABB.max.x, m_AABB.max.y);
        ctx.lineTo(m_AABB.max.x, m_AABB.min.y);
        ctx.lineTo(m_AABB.min.x, m_AABB.min.y);

        ctx.stroke();
        ctx.closePath();
    }

-->

Quickhull into js and html5- p3

After this, I will do an implementation for all the popular convex hull finding techniques.

It would be so easy to write just a recursive function for the quickhull in c++/c#, but putting the algorithm into javascript is not so straight forward.. but its fun :-)
More vector logic and triangle point inside checks added, the logic now checks for extreme points for each edge, I just need to add them into the edge and recourse through the whole logic for new edge definitions.

Link to source: https://skydrive.live.com/redir?resid=700BE2A366278EA0!133&authkey=!AATEd_SmxpiGDW8

More later..

Quickhull experiment in js & html5 - part2

Cleaning up the code for Quickhull. Added more math routines in it. Created classes inside javascript. Its simple to create classes(objects only) in javascript. One has to use the Java Script Object Notation (JSON)

Example:


var bounding_box = {
    "min": {
        "x":0,
        "y":0
    },
    "max": {
        "x":0,
        "y":0
    }
};

The source code so far: QuickHull java script code. 

Needs to recursively call on the edge subdivision. That will be when I get time next :-)