/* -*- mode: C -*-  */
/* vim:set ts=2 sw=2 sts=2 et: */
/* 
   IGraph R package.
   Copyright (C) 2014  Gabor Csardi <csardi.gabor@gmail.com>
   334 Harvard street, Cambridge, MA 02139 USA
   
   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 2 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, write to the Free Software
   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 
   02110-1301 USA

*/

#include "igraph_layout.h"
#include "igraph_interface.h"
#include "igraph_random.h"
#include "igraph_math.h"

/**
 * \ingroup layout
 * \function igraph_layout_gem
 *
 * The GEM layout algorithm, as described in Arne Frick, Andreas Ludwig,
 * Heiko Mehldau: A Fast Adaptive Layout Algorithm for Undirected Graphs,
 * Proc. Graph Drawing 1994, LNCS 894, pp. 388-403, 1995.
 * \param graph The input graph. Edge directions are ignored in
 *        directed graphs.
 * \param res The result is stored here. If the \p use_seed argument
 *        is true (non-zero), then this matrix is also used as the
 *        starting point of the algorithm.
 * \param use_seed Boolean, whether to use the supplied coordinates in
 *        \p res as the starting point. If false (zero), then a
 *        uniform random starting point is used.
 * \param maxiter The maximum number of iterations to
 *        perform. Updating a single vertex counts as an iteration.
 *        A reasonable default is 40 * n * n, where n is the number of 
 *        vertices. The original paper suggests 4 * n * n, but this 
 *        usually only works if the other parameters are set up carefully.
 * \param temp_max The maximum allowed local temperature. A reasonable
 *        default is the number of vertices.
 * \param temp_min The global temperature at which the algorithm
 *        terminates (even before reaching \p maxiter iterations). A
 *        reasonable default is 1/10.
 * \param temp_init Initial local temperature of all vertices. A
 *        reasonable default is the square root of the number of
 *        vertices.
 * \return Error code.
 * 
 * Time complexity: O(t * n * (n+e)), where n is the number of vertices,
 * e is the number of edges and t is the number of time steps
 * performed.
 */

int igraph_layout_gem(const igraph_t *graph, igraph_matrix_t *res,
		      igraph_bool_t use_seed, igraph_integer_t maxiter,
		      igraph_real_t temp_max, igraph_real_t temp_min,
		      igraph_real_t temp_init) {

  igraph_integer_t no_nodes = igraph_vcount(graph);
  igraph_vector_int_t perm;
  igraph_vector_float_t impulse_x, impulse_y, temp, skew_gauge;
  igraph_integer_t i;
  float temp_global;
  igraph_integer_t perm_pointer = 0;
  float barycenter_x = 0.0, barycenter_y = 0.0;
  igraph_vector_t phi;
  igraph_vector_t neis;
  const float elen_des2 = 128 * 128;
  const float gamma = 1/16.0;
  const float alpha_o = M_PI;
  const float alpha_r = M_PI / 3.0;
  const float sigma_o = 1.0 / 3.0;
  const float sigma_r = 1.0 / 2.0 / no_nodes;
  
  if (maxiter < 0) {
    IGRAPH_ERROR("Number of iterations must be non-negative in GEM layout",
		 IGRAPH_EINVAL);
  }
  if (use_seed && (igraph_matrix_nrow(res) != no_nodes ||
		   igraph_matrix_ncol(res) != 2)) {
    IGRAPH_ERROR("Invalid start position matrix size in GEM layout",
		 IGRAPH_EINVAL);
  }
  if (temp_max <= 0) {
    IGRAPH_ERROR("Maximum temperature should be positive in GEM layout",
		 IGRAPH_EINVAL);
  }
  if (temp_min <= 0) {
    IGRAPH_ERROR("Minimum temperature should be positive in GEM layout",
		 IGRAPH_EINVAL);
  }
  if (temp_init <= 0) {
    IGRAPH_ERROR("Initial temperature should be positive in GEM layout",
		 IGRAPH_EINVAL);
  }
  if (temp_max < temp_init || temp_init < temp_min) {
    IGRAPH_ERROR("Minimum <= Initial <= Maximum temperature is required "
		 "in GEM layout", IGRAPH_EINVAL);
  }

  if (no_nodes == 0) { return 0; }

  IGRAPH_CHECK(igraph_vector_float_init(&impulse_x, no_nodes));
  IGRAPH_FINALLY(igraph_vector_float_destroy, &impulse_x);
  IGRAPH_CHECK(igraph_vector_float_init(&impulse_y, no_nodes));
  IGRAPH_FINALLY(igraph_vector_float_destroy, &impulse_y);
  IGRAPH_CHECK(igraph_vector_float_init(&temp, no_nodes));
  IGRAPH_FINALLY(igraph_vector_float_destroy, &temp);
  IGRAPH_CHECK(igraph_vector_float_init(&skew_gauge, no_nodes));
  IGRAPH_FINALLY(igraph_vector_float_destroy, &skew_gauge);
  IGRAPH_CHECK(igraph_vector_int_init_seq(&perm, 0, no_nodes-1));
  IGRAPH_FINALLY(igraph_vector_int_destroy, &perm);
  IGRAPH_VECTOR_INIT_FINALLY(&phi, no_nodes);
  IGRAPH_VECTOR_INIT_FINALLY(&neis, 10);

  RNG_BEGIN();

  /* Initialization */
  igraph_degree(graph, &phi, igraph_vss_all(), IGRAPH_ALL, IGRAPH_LOOPS);
  if (!use_seed) {
    const igraph_real_t width_half=no_nodes*100, height_half=width_half;
    IGRAPH_CHECK(igraph_matrix_resize(res, no_nodes, 2));
    for (i=0; i<no_nodes; i++) {
      MATRIX(*res, i, 0) = RNG_UNIF(-width_half, width_half);
      MATRIX(*res, i, 1) = RNG_UNIF(-height_half, height_half);
      barycenter_x += MATRIX(*res, i, 0);
      barycenter_y += MATRIX(*res, i, 1);
      VECTOR(phi)[i] *= (VECTOR(phi)[i] / 2.0 + 1.0);
    }
  } else {
    for (i=0; i<no_nodes; i++) {
      barycenter_x += MATRIX(*res, i, 0);
      barycenter_y += MATRIX(*res, i, 1);
      VECTOR(phi)[i] *= (VECTOR(phi)[i] / 2.0 + 1.0);
    }
  }
  igraph_vector_float_fill(&temp, temp_init);
  temp_global = temp_init * no_nodes;
  
  while (temp_global > temp_min * no_nodes && maxiter > 0) {
    
    /* choose a vertex v to update */
    igraph_integer_t u, v, nlen, j;
    float px, py, pvx, pvy;
    if (!perm_pointer) { 
      igraph_vector_int_shuffle(&perm); 
      perm_pointer=no_nodes-1;
    }
    v=VECTOR(perm)[perm_pointer--];
    
    /* compute v's impulse */
    px = (barycenter_x/no_nodes - MATRIX(*res, v, 0)) * gamma * VECTOR(phi)[v];
    py = (barycenter_y/no_nodes - MATRIX(*res, v, 1)) * gamma * VECTOR(phi)[v];
    px += RNG_UNIF(-32.0, 32.0);
    py += RNG_UNIF(-32.0, 32.0);

    for (u = 0; u < no_nodes; u++) {
      float dx, dy, dist2;
      if (u == v) { continue; }
      dx=MATRIX(*res, v, 0) - MATRIX(*res, u, 0);
      dy=MATRIX(*res, v, 1) - MATRIX(*res, u, 1);
      dist2=dx * dx + dy * dy;
      if (dist2 != 0) {
	px += dx * elen_des2 / dist2;
	py += dy * elen_des2 / dist2;
      }
    }

    IGRAPH_CHECK(igraph_neighbors(graph, &neis, v, IGRAPH_ALL));
    nlen=igraph_vector_size(&neis);
    for (j = 0; j < nlen; j++) {
      igraph_integer_t u=VECTOR(neis)[j];
      float dx=MATRIX(*res, v, 0) - MATRIX(*res, u, 0);
      float dy=MATRIX(*res, v, 1) - MATRIX(*res, u, 1);
      float dist2= dx * dx + dy * dy;
      px -= dx * dist2 / (elen_des2 * VECTOR(phi)[v]);
      py -= dy * dist2 / (elen_des2 * VECTOR(phi)[v]);
    }

    /* update v's position and temperature */
    if (px != 0 || py != 0) {
      float plen = sqrtf(px * px + py * py);
      px *= VECTOR(temp)[v] / plen;
      py *= VECTOR(temp)[v] / plen;
      MATRIX(*res, v, 0) += px;
      MATRIX(*res, v, 1) += py;
      barycenter_x += px;
      barycenter_y += py;
    }
    
    pvx=VECTOR(impulse_x)[v]; pvy=VECTOR(impulse_y)[v];
    if (pvx != 0 || pvy != 0) {
      float beta = atan2f(pvy - py, pvx - px);
      float sin_beta = sinf(beta);
      float sign_sin_beta = (sin_beta > 0) ? 1 : ((sin_beta < 0) ? -1 : 0);
      float cos_beta = cosf(beta);
      float abs_cos_beta = fabsf(cos_beta);
      float old_temp=VECTOR(temp)[v];
      if (sin(beta) >= sin(M_PI_2 + alpha_r / 2.0)) {
	VECTOR(skew_gauge)[v] += sigma_r * sign_sin_beta;
      }
      if (abs_cos_beta >= cosf(alpha_o / 2.0)) {
	VECTOR(temp)[v] *= sigma_o * cos_beta;
      }
      VECTOR(temp)[v] *= (1 - fabsf(VECTOR(skew_gauge)[v]));
      if (VECTOR(temp)[v] > temp_max) { VECTOR(temp)[v] = temp_max; }
      VECTOR(impulse_x)[v] = px;
      VECTOR(impulse_y)[v] = py;
      temp_global += VECTOR(temp)[v] - old_temp;
    }

    maxiter--;

  } /* while temp && iter */
  

  RNG_END();
    
  igraph_vector_destroy(&neis);
  igraph_vector_destroy(&phi);
  igraph_vector_int_destroy(&perm);
  igraph_vector_float_destroy(&skew_gauge);
  igraph_vector_float_destroy(&temp);
  igraph_vector_float_destroy(&impulse_y);
  igraph_vector_float_destroy(&impulse_x);
  IGRAPH_FINALLY_CLEAN(7);
  
  return 0;
}