cola.h

/*
 * vim: ts=4 sw=4 et tw=0 wm=0
 *
 * libcola - A library providing force-directed network layout using the 
 *           stress-majorization method subject to separation constraints.
 *
 * Copyright (C) 2006-2015  Monash University
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * See the file LICENSE.LGPL distributed with the library.
 *
 * This library 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.
 *
 * Author(s):  Tim Dwyer
 *             Michael Wybrow
*/

#ifndef COLA_H
#define COLA_H

#include <utility>
#include <iterator>
#include <vector>
#include <valarray>
#include <algorithm>
#include <cmath>
#include <iostream>

#include "libcola/gradient_projection.h"
#include "libcola/cluster.h"
#include "libcola/straightener.h"
#include "libcola/exceptions.h"
#include "libcola/pseudorandom.h"

namespace vpsc { class Rectangle; }
namespace topology { 
    class ColaTopologyAddon;
}
namespace dialect {
    class Graph;
}


namespace cola {

class NonOverlapConstraints;
class NonOverlapConstraintExemptions;

typedef std::vector<unsigned> NodeIndexes;

typedef std::vector<NodeIndexes> ListOfNodeIndexes;

typedef std::pair<unsigned, unsigned> Edge;

typedef std::vector<double> EdgeLengths;
#define StandardEdgeLengths EdgeLengths()

class Lock {
public:
    Lock() {}
    Lock(unsigned id, double X, double Y) : id(id), x(X), y(Y) {
    }
    unsigned getID() const {
        return id;
    }
    double pos(vpsc::Dim dim) const {
        return dim==vpsc::HORIZONTAL?x:y;
    }
private:
    unsigned id;
    double x;
    double y;
};
typedef std::vector<cola::Lock> Locks;

class Resize {
public:
    Resize() {}
    Resize(unsigned id, double x, double y, double w, double h)
        : id(id), target(x,x+w,y,y+h) {}
    unsigned getID() const {
        return id;
    }
    const vpsc::Rectangle* getTarget() const {
        return &target;
    }
private:
    unsigned id;
    vpsc::Rectangle target;
};
typedef std::vector<cola::Resize> Resizes;

/*
 * Setting a desired position for a node adds a term to the goal function
 * drawing the node towards that desired position
 */
struct DesiredPosition {
    unsigned id;
    double x;
    double y;
    double weight;
};
typedef std::vector<cola::DesiredPosition> DesiredPositions;

/*
 * desired positions which should override those computed by applying forces
 * are passed in for a set of nodes.  The first entry is the Node->id, the
 * second is the desired position.
 */
typedef std::pair<unsigned,double> DesiredPositionInDim;
typedef std::vector<DesiredPositionInDim> DesiredPositionsInDim;

class PreIteration {
public:
    PreIteration(
            Locks& locks=__locksNotUsed,
            Resizes& resizes=__resizesNotUsed) 
        : locks(locks)
        , resizes(resizes) 
        , changed(true) {}
    PreIteration(Resizes& resizes) 
        : locks(__locksNotUsed)
        , resizes(resizes)
        , changed(true) {}

// To prevent C++ objects from being destroyed in garbage collected languages
// when the libraries are called from SWIG, we hide the declarations of the
// destructors and prevent generation of default destructors.
#ifndef SWIG
    virtual ~PreIteration() {}
#endif
    virtual bool operator()() { return true; }
    Locks& locks;
    Resizes& resizes;
    bool changed;
private:
    static Locks __locksNotUsed;
    static Resizes __resizesNotUsed;
};

class TestConvergence {
public:
    double old_stress;
    TestConvergence(const double tol = 1e-4, const unsigned maxiterations = 100)
        : tolerance(tol),
          maxiterations(maxiterations)
    { reset(); }
    virtual ~TestConvergence() {}

public:
    virtual bool operator()(const double new_stress, std::valarray<double> & X, std::valarray<double> & Y)
    {
        COLA_UNUSED(X);
        COLA_UNUSED(Y);

        iterations++;
        //std::cout<<"iteration="<<iterations<<", old_stress="<<old_stress<<", new_stress="<<new_stress<<std::endl;
        if (old_stress == DBL_MAX) {
            old_stress = new_stress;
            return iterations >= maxiterations;
        }
        // converged if relative decrease in stress falls below threshold
        // or if stress increases (shouldn't happen for straight majorization)
        bool converged = 
            (old_stress - new_stress) / (new_stress + 1e-10) < tolerance
             || iterations > maxiterations;
        old_stress = new_stress;
        return converged;
    }
    void reset() {
        old_stress = DBL_MAX;
        iterations = 0;
    }
    const double tolerance;
    const unsigned maxiterations;
    unsigned iterations;
};

class ConstrainedMajorizationLayout {
public:
    ConstrainedMajorizationLayout(
        vpsc::Rectangles& rs,
        std::vector<Edge> const & es,
        RootCluster* clusterHierarchy,
        const double idealLength,
        EdgeLengths eLengths = StandardEdgeLengths,
        TestConvergence *doneTest = nullptr,
        PreIteration* preIteration=nullptr,
        bool useNeighbourStress = false);
    void setConstraints(cola::CompoundConstraints* ccs) {
        constrainedLayout = true;
        this->ccs=ccs;
    }

    void setConstraintsVector(cola::CompoundConstraints& ccs) {
        constrainedLayout = true;
        cola::CompoundConstraints *ccsp = new cola::CompoundConstraints;
        for (size_t i = 0; i < ccs.size(); ++i) {
            ccsp->push_back(ccs.at(i));
        }
        this->ccs=ccsp;
    }

    void setUnsatisfiableConstraintInfo(
            UnsatisfiableConstraintInfos *unsatisfiableX,
            UnsatisfiableConstraintInfos *unsatisfiableY) {
        this->unsatisfiableX = unsatisfiableX;
        this->unsatisfiableY = unsatisfiableY;
    }
    void setStickyNodes(const double stickyWeight, 
            std::valarray<double> const & startX,
            std::valarray<double> const & startY);

    void setScaling(bool scaling) {
        this->scaling=scaling;
    }
    void setExternalSolver(bool externalSolver) {
        this->externalSolver=externalSolver;
    }
    void setAvoidOverlaps(bool horizontal = false) {
        constrainedLayout = true;
        this->avoidOverlaps = horizontal ? Horizontal : Both;
    }
    void setNonOverlappingClusters() {
        constrainedLayout = true;
        nonOverlappingClusters = true;
    }
    void setStraightenEdges(std::vector<straightener::Edge*>* straightenEdges, 
            double bendWeight = 0.01, double potBendWeight = 0.1,
            bool xSkipping = true) {
        for(std::vector<straightener::Edge*>::const_iterator e=straightenEdges->begin();
                e!=straightenEdges->end();e++) {
            (*e)->rerouteAround(boundingBoxes);
        }
        constrainedLayout = true;
        this->xSkipping = xSkipping;
        this->straightenEdges = straightenEdges;
        this->bendWeight = bendWeight;
        this->potBendWeight = potBendWeight;
    }
    void moveBoundingBoxes() {
        for(unsigned i=0;i<n;i++) {
            boundingBoxes[i]->moveCentre(X[i],Y[i]);
        }
    }

    ~ConstrainedMajorizationLayout() {
        if (using_default_done)
        {
            delete done;
        }

        if(constrainedLayout) {
            delete gpX;
            delete gpY;
        }
    }
    void run(bool x=true, bool y=true);
    void runOnce(bool x=true, bool y=true);
    void straighten(std::vector<straightener::Edge*>&, vpsc::Dim);
    void setConstrainedLayout(bool c) {
        constrainedLayout=c;
    }
    double computeStress();
private:
    double euclidean_distance(unsigned i, unsigned j) {
        return sqrt(
            (X[i] - X[j]) * (X[i] - X[j]) +
            (Y[i] - Y[j]) * (Y[i] - Y[j]));
    }
    double compute_stress(std::valarray<double> const & Dij);
    void majorize(std::valarray<double> const & Dij,GradientProjection* gp, std::valarray<double>& coords, std::valarray<double> const & startCoords);
    void newton(std::valarray<double> const & Dij,GradientProjection* gp, std::valarray<double>& coords, std::valarray<double> const & startCoords);
    unsigned n; //< number of nodes
    //std::valarray<double> degrees;
    std::valarray<double> lap2; //< graph laplacian
    std::valarray<double> Q; //< quadratic terms matrix used in computations
    std::valarray<double> Dij; //< all pairs shortest path distances
    double tol; //< convergence tolerance
    TestConvergence *done; //< functor used to determine if layout is finished
    bool using_default_done; // Whether we allocated a default TestConvergence object.
    PreIteration* preIteration; //< client can use this to create locks on nodes
    vpsc::Rectangles boundingBoxes; //< node bounding boxes
    /*
     * stickyNodes controls whether nodes are attracted to their starting
     * positions (at time of ConstrainedMajorizationLayout instantiation)
     * stored in startX, startY
     */
    std::valarray<double> X, Y;
    bool stickyNodes;
    double stickyWeight;
    std::valarray<double> startX;
    std::valarray<double> startY;
    double edge_length;
    bool constrainedLayout;
    bool nonOverlappingClusters;
    /*
     * A cluster is a set of nodes that are somehow semantically grouped
     * and should therefore be kept together a bit more tightly than, and
     * preferably without overlapping, the rest of the graph.
     *
     * We achieve this by augmenting the L matrix with stronger attractive
     * forces between all members of a cluster (other than the root)
     * and by maintaining a (preferably convex) hull around those 
     * constituents which, using constraints and dummy variables, is 
     * prevented from overlapping other parts of the graph.
     *
     * Clusters are defined over the graph in a hierarchy starting with
     * a single root cluster.
     *
     * Need to:
     *  - augment Lap matrix with intra cluster forces
     *  - compute convex hull of each cluster
     *  - from convex hull generate "StraightenEdges"
     */
    RootCluster *clusterHierarchy;
    GradientProjection *gpX, *gpY;
    cola::CompoundConstraints *ccs;
    UnsatisfiableConstraintInfos *unsatisfiableX, *unsatisfiableY;
    NonOverlapConstraintsMode avoidOverlaps;
    std::vector<straightener::Edge*>* straightenEdges;
    
    double bendWeight, potBendWeight;
    /* 
     * determines whether we should leave some overlaps to be resolved
     * vertically when generating straightening constraints in the x-dim
     */
    bool xSkipping;
    /* 
     * when using the gradient projection optimisation method, the following
     * controls whether the problem should be preconditioned by affine scaling
     */
    bool scaling;
    /* 
     * if the Mosek quadratic programming environment is available it may be
     * used to solve each iteration of stress majorization... slow but useful
     * for testing 
     */
    bool externalSolver;
    bool majorization;
};

vpsc::Rectangle bounds(vpsc::Rectangles& rs);

class ConstrainedFDLayout;

class TopologyAddonInterface 
{
    public:
        TopologyAddonInterface()
        {
        }

        virtual ~TopologyAddonInterface()
        {
        }

        virtual TopologyAddonInterface *clone(void) const
        {
            return new TopologyAddonInterface(*this);
        }
        
        virtual void freeAssociatedObjects(void)
        {
        }

        virtual void handleResizes(const Resizes& resizeList, unsigned n,
                std::valarray<double>& X, std::valarray<double>& Y, 
                cola::CompoundConstraints& ccs, 
                vpsc::Rectangles& boundingBoxes,
                cola::RootCluster* clusterHierarchy)
        {
            COLA_UNUSED(resizeList);
            COLA_UNUSED(n);
            COLA_UNUSED(X);
            COLA_UNUSED(Y);
            COLA_UNUSED(ccs);
            COLA_UNUSED(boundingBoxes);
            COLA_UNUSED(clusterHierarchy);
        }
        virtual void computePathLengths(unsigned short** G)
        {
            COLA_UNUSED(G);
        }
        virtual double computeStress(void) const
        {
            return 0;
        }
        virtual bool useTopologySolver(void) const
        {
            return false;
        }
        virtual void makeFeasible(bool generateNonOverlapConstraints, 
                vpsc::Rectangles& boundingBoxes, 
                cola::RootCluster* clusterHierarchy)
        {
            COLA_UNUSED(generateNonOverlapConstraints);
            COLA_UNUSED(boundingBoxes);
            COLA_UNUSED(clusterHierarchy);
        }
        virtual void moveTo(const vpsc::Dim dim, vpsc::Variables& vs, 
                vpsc::Constraints& cs, std::valarray<double> &coords,
                cola::RootCluster* clusterHierarchy)
        {
            COLA_UNUSED(dim);
            COLA_UNUSED(vs);
            COLA_UNUSED(cs);
            COLA_UNUSED(coords);
            COLA_UNUSED(clusterHierarchy);
        }
        virtual double applyForcesAndConstraints(ConstrainedFDLayout *layout, 
                const vpsc::Dim dim, std::valarray<double>& g, 
                vpsc::Variables& vs, vpsc::Constraints& cs, 
                std::valarray<double> &coords,
                DesiredPositionsInDim& des, double oldStress)
        {
            COLA_UNUSED(layout);
            COLA_UNUSED(dim);
            COLA_UNUSED(g);
            COLA_UNUSED(vs);
            COLA_UNUSED(cs);
            COLA_UNUSED(coords);
            COLA_UNUSED(des);
            COLA_UNUSED(oldStress);
            return 0;
        }
};


class ConstrainedFDLayout {
public:
    ConstrainedFDLayout(
        const vpsc::Rectangles& rs,
        const std::vector<cola::Edge>& es,
        const double idealLength,
        const EdgeLengths& eLengths = StandardEdgeLengths, 
        TestConvergence* doneTest = nullptr,
        PreIteration* preIteration = nullptr);
    ~ConstrainedFDLayout();
  
    void run(bool x=true, bool y=true);

    void runOnce(bool x=true, bool y=true);

    void setConstraints(const cola::CompoundConstraints& ccs);

    void setAvoidNodeOverlaps(bool avoidOverlaps,
            ListOfNodeIndexes listOfNodeGroups =
            ListOfNodeIndexes());

    void setTopology(TopologyAddonInterface *topology); 
    TopologyAddonInterface *getTopology(void);
    
    void setDesiredPositions(DesiredPositions *desiredPositions);

    void setClusterHierarchy(RootCluster *hierarchy)
    {
        clusterHierarchy = hierarchy;
    }
    void setUnsatisfiableConstraintInfo(
            UnsatisfiableConstraintInfos *unsatisfiableX,
            UnsatisfiableConstraintInfos *unsatisfiableY)
    {
        unsatisfiable.resize(2);
        unsatisfiable[0]=unsatisfiableX;
        unsatisfiable[1]=unsatisfiableY;
    }
    void makeFeasible(double xBorder=1, double yBorder=1);

    void freeAssociatedObjects(void);
    
    void outputInstanceToSVG(std::string filename = std::string());

    void setUseNeighbourStress(bool useNeighbourStress);

    std::vector<double> readLinearD(void);

    std::vector<unsigned> readLinearG(void);

    double computeStress() const;

private:
    unsigned n; // number of nodes
    std::valarray<double> X, Y;
    vpsc::Rectangles boundingBoxes;
    double applyForcesAndConstraints(const vpsc::Dim dim,const double oldStress);
    double computeStepSize(const SparseMatrix& H, const std::valarray<double>& g,
            const std::valarray<double>& d) const;
    void computeDescentVectorOnBothAxes(const bool xaxis, const bool yaxis,
            double stress, std::valarray<double>& x0, std::valarray<double>& x1);
    void moveTo(const vpsc::Dim dim, std::valarray<double>& target);
    double applyDescentVector(
            const std::valarray<double>& d,
            const std::valarray<double>& oldCoords,
            std::valarray<double> &coords, 
            const double oldStress, 
            double stepsize
            /*,topology::TopologyConstraints *s=nullptr*/);
    void computePathLengths(
            const std::vector<Edge>& es, std::valarray<double> eLengths);
    void generateNonOverlapAndClusterCompoundConstraints(
            vpsc::Variables (&vs)[2]);
    void handleResizes(const Resizes&);
    void setPosition(std::valarray<double>& pos);
    void moveBoundingBoxes();
    bool noForces(double, double, unsigned) const;
    void computeForces(const vpsc::Dim dim, SparseMap &H, 
            std::valarray<double> &g);
    void recGenerateClusterVariablesAndConstraints(
            vpsc::Variables (&vars)[2], unsigned int& priority, 
            cola::NonOverlapConstraints *noc, Cluster *cluster, 
            cola::CompoundConstraints& idleConstraints);
    std::vector<double> offsetDir(double minD);

    void computeNeighbours(std::vector<Edge> es);
    std::vector<std::vector<unsigned> > neighbours;
    std::vector<std::vector<double> > neighbourLengths;
    TestConvergence *done;
    bool using_default_done; // Whether we allocated a default TestConvergence object.
    PreIteration* preIteration;
    cola::CompoundConstraints ccs;
    double** D;
    unsigned short** G;
    double minD;
    PseudoRandom random;

    TopologyAddonInterface *topologyAddon;
    std::vector<UnsatisfiableConstraintInfos*> unsatisfiable;
    bool rungekutta;
    DesiredPositions *desiredPositions;
    cola::CompoundConstraints extraConstraints;
    
    RootCluster *clusterHierarchy;
    double rectClusterBuffer;
    double m_idealEdgeLength;
    bool m_generateNonOverlapConstraints;
    bool m_useNeighbourStress;
    const std::valarray<double> m_edge_lengths;

    NonOverlapConstraintExemptions *m_nonoverlap_exemptions;

    friend class topology::ColaTopologyAddon;
    friend class dialect::Graph;
};

struct ProjectionResult {
    int errorLevel;
    std::string unsatinfo;
};

ProjectionResult projectOntoCCs(vpsc::Dim dim, vpsc::Rectangles &rs, cola::CompoundConstraints ccs,
                                bool preventOverlaps, int accept=0, unsigned debugLevel=0);

ProjectionResult solve(vpsc::Variables &vs, vpsc::Constraints &cs, vpsc::Rectangles &rs,
                        unsigned debugLevel=0);


ConstrainedMajorizationLayout* simpleCMLFactory(
        vpsc::Rectangles& rs,
        std::vector<Edge> const & es,
        RootCluster* clusterHierarchy,
        const double idealLength,
        bool useNeighbourStress = false
    );

/*
 * find shortest path lengths from node s to all other nodes.
 * @param s starting node
 * @param n total number of nodes
 * @param d n vector of path lengths
 * @param es edge pairs
 * @param eLengths edge weights
 */
void dijkstra(const unsigned s, const unsigned n, double* d, 
              const std::vector<cola::Edge>& es,
              const std::valarray<double> & eLengths);

#if 0
void removeClusterOverlapFast(RootCluster& clusterHierarchy, vpsc::Rectangles& rs, Locks& locks);
#endif

void setupVarsAndConstraints(unsigned n, const CompoundConstraints& ccs,
        const vpsc::Dim dim, vpsc::Rectangles& boundingBoxes,
        RootCluster *clusterHierarchy,
        vpsc::Variables& vs, vpsc::Constraints& cs, 
        std::valarray<double> &coords);
void setVariableDesiredPositions(vpsc::Variables& vs, vpsc::Constraints& cs,
        const DesiredPositionsInDim& des, std::valarray<double>& coords);

}
#endif // COLA_H