// Copyright (C) 1996 DIMACS Center, Rutgers, The State University of New Jersey
// Author(s): Michael Dineen (Los Alamos National Laboratory)

// This software is copyrighted by the DIMACS Center at Rutgers, The State
// University of New Jersey.  IT IS PROVIDED AS IS, AND THE AUTHORS, DIMACS, AND
// RUTGERS, THE STATE UNIVERSITY OF NEW JERSEY  DISCLAIM
// ALL LIABILITY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL
// DAMAGES ARISING OUT OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY
// DERIVATIVES THEREOF, EVEN IF THE AUTHORS HAVE BEEN ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.

// THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES,
// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE, AND NON-INFRINGEMENT.  THIS SOFTWARE
// IS PROVIDED ON AN "AS IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE
// NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
// MODIFICATIONS.

// The authors hereby grant permission to use, copy, modify, distribute,
// and license this software and its documentation for any purpose, provided
// that existing copyright notices are retained in all copies and that this
// notice is included verbatim in any distributions. No written agreement,
// license, or royalty fee is required for any of the authorized uses.
// Modifications to this software may be copyrighted by their authors
// and need not follow the licensing terms described here, provided that
// the new terms are clearly indicated on the first page of each file where
// they apply.

// Last File Update: 31-Jul-1996
// 

///////////////////////////////////////////////////////////////////////////
// LINK: Generic Graph Tool and Class Library
//                               
//      Function name: SpringLayout(Graph*)
//
//		       Standard 'distance matrix' spring layout of a graph
//                                                                           
//      Creation:  02-Dec-94, Michael J. Dinneen, mjd@lanl.gov
//                                                                            
//      History:

/*************************************************************************
 *
 *  You might check this prog which is a straight implementation
 *  of an iterative solution of the spring model: compute the shortest
 *  path for all pairs of vertices and place springs of appropriate lengths
 *  between them. Then find the equilibrium state of this system of springs
 *  and points. - This model is standard, I don't have any reference handy.
 *
 *  The prog is in Turbo Pascal, but should be easily modifyable.
 *  Hope this helps. Johannes - joe@inuo30.mathematik.uni-jena.de}
 *
 *  Lazy conversion to C/C++ by 'mjd@csr.uvic.ca' (Feb 93)
 * ************************************************************************/

//   **** See the followup algorithm in layout/IncArrangement.cc ****
                                                                            
//////////////////////////////////////////////////////////////////////////////

#include <LINK/basic/Array.h>  

#include <stdlib.h>
#include <iostream.h>
#include <math.h>
#include <LINK/basic/Matrix.h>
#include <LINK/graph/Attribute.h>
#include <LINK/graph/Graph.h>
#include <LINK/graph/Vertex.h>
#include <LINK/algorithm/Algorithms.h>
#include <LINK/layout/Layouts.h>

#define START_RANDOM 0

inline double sqr(double x)     { return x*x; }
const double GS_eps = 0.0000005;
const double Viewport_Buffer = 0.15;
static int n;    		// order of the graph


static void prepare(Graph *G, AsymMatrix<double> &D, double **f,
		    double *x, double *y, int start_random)
{
	int i, j;

          FloydWarshall(G, D);

          // { compute spring constants }
          //
          for (i=0; i<n; i++) {
              for ( j = 0; j<n; j++ ) {
                  if (i != j) {
		     //cout << "about to compute" << endl;
                     f[i][j] = 1.0 / sqr (D(i,j));
		     //cout << "computed" << endl;
                  }
              }
          }
          for (i=0; i<n; i++)
          {
		if (start_random) {
            		x[i] = Link_randomDouble0_1();
            		y[i] = Link_randomDouble0_1();
		} else {
              		x[i] = cos (i * 2 * M_PI / n);
              		y[i] = sin (i * 2 * M_PI / n);
		}
          }
}

static void forces(const AsymMatrix<double> &D, double *x, double *y,
		   double *fx, double *fy, double *fs, double **f) 
{
int i, j;
double sx, sy, dx, dy, tmp;

       for (i=0; i<n; i++) {
              sx = 0; sy = 0;
              for ( j = 0; j<n; j++ ) {
                  if (i != j) {
                     dx = x[i] - x[j]; dy = y[i] - y[j];
                     tmp = f[i][j] * (1 - D.at(i,j)
                                            / sqrt (sqr (dx) + sqr (dy)));
                     sx = sx + tmp * dx; sy = sy + tmp * dy;
                  }
              }
              fx[i] = sx; fy[i] = sy;
              fs[i] = sqrt (sqr (sx) + sqr (sy));
          }
}

// { compute offset for point m }
//
static void handle(int m, const AsymMatrix<double> &D, double *x, double *y,
		   double *fx, double *fy, double **f) 
{     
	double fxx=0.0, fxy=0.0, fyy=0.0, dx, dy, tmp;
	int i;

          for (i=0; i<n; i++) {
              if (i != m) {
                 dx  = x[m] - x[i]; dy = y[m] - y[i];
                 tmp = sqr (dx) + sqr (dy);
                 tmp = 1.0 / tmp / sqrt (tmp);
                 fxx = fxx + f[m][i] * (1 - D.at(m,i) * sqr (dy) * tmp);
                 fxy = fxy + f[m][i] * D.at(m,i) * dx * dy * tmp;
                 fyy = fyy + f[m][i] * (1 - D.at(m,i) * sqr (dx) * tmp);
              }
          }
          tmp   = sqr (fxy) - fxx * fyy;
          x[m] = x[m] + (fx[m] * fyy - fy[m] * fxy) / tmp;
          y[m] = y[m] + (fy[m] * fxx - fx[m] * fxy) / tmp;
}

static void loop(const AsymMatrix<double> &D, double *x, double *y,
		 double *fx, double *fy, double *fs, double **f) 
{     
	Bool leaving = FALSE;
	int maxind, i;
	double max;
	int maxTries = 1000;

          //repeat  forces;
          //
          while (leaving==FALSE && maxTries--)
           {
                  forces(D, x, y, fx, fy, fs, f);

                  maxind = 0; max = fs[0];
                  for (i = 1; i<n; i++) {
                      if (fs[i] > max) {
                         maxind = i; max = fs[i];
                      }
                  }

                  if (max < GS_eps) 
			leaving = TRUE;
                  else 	  	 
			handle(maxind, D, x, y, fx, fy, f); 
           } 
}

static void getSpringerLocations(Graph *G, Array<double>& xcoords,
					   Array<double>& ycoords)
{
    n = G->order();  //assert(n<=GS_index);

    double **f = new double*[n]; 	// spring constants
    f[0] = new double[n*n];
    for (int nloop =1; nloop<n; nloop++) f[nloop] = f[0] + nloop*n;
    double *x = new double[n]; 		// x coordinates (world)
    double *y = new double[n];		// y coordinates (world)
    double *fx = new double[n]; 	// x component of force at point
    double *fy = new double[n];		// y component of force at point
    double *fs = new double[n];		// = fx^2 + fy^2

    AsymMatrix<double> D(n,n);
    prepare(G, D, f, x, y, START_RANDOM);
    loop(D, x, y, fx, fy, fs, f); 

    int i;
    for (i=0; i<n ; i++) {
	xcoords[i] = x[i];
	ycoords[i] = y[i];
    }

    delete [] f[0]; delete [] f;
    delete [] x; delete [] y;
    delete [] fx; delete [] fy;
    delete [] fs;

    return;
}

void SpringLayout(Graph* graph)
{
    if (!graph->binaryQ()) {
	HyperSpringLayout(graph);
	return;
    }
    int order = graph->order();
    Array<double> wxcoords(order); 
    Array<double> wycoords(order); 

    getSpringerLocations(graph, wxcoords, wycoords);

    double xmin = wxcoords.min() - Viewport_Buffer;
    double xmax = wxcoords.max() + Viewport_Buffer;
    double ymin = wycoords.min() - Viewport_Buffer;
    double ymax = wycoords.max() + Viewport_Buffer;
    double delta_x = xmax-xmin;
    double delta_y = ymax-ymin;
    //double delta_max = (delta_x > delta_y) ? delta_x : delta_y;


    const SortedArray<Vertex*>& verts = graph->vertices();
    for (int i=0; i<order; i++)
    {
	setAttribute((GraphObject*) verts[i], "x", 
			(wxcoords[i] - xmin)/delta_x);
	setAttribute((GraphObject*) verts[i], "y", 
			(wycoords[i] - ymin)/delta_y);
    }
}