/* Copyright (c) 1997-2007
   Ewgenij Gawrilow, Michael Joswig (Technische Universitaet Berlin, Germany)
   http://www.math.tu-berlin.de/polymake,  mailto:polymake@math.tu-berlin.de

   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, or (at your option) any
   later version: http://www.gnu.org/licenses/gpl.txt.

   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.
*/

#ident "$Project: polymake $$Id: SpringEmbedder.cc 7579 2007-01-22 11:34:46Z gawrilow $"

#include <cstdlib>
#include <cmath>
#include <SpringEmbedder.h>

namespace polymake { namespace graph {

void SpringEmbedder::init_params(const argv_option *options)
{
   if (options[opt_scale]) {
      scale=std::atof(options[opt_scale]);
      epsilon_2*=scale*scale;
   } else {
      scale=1;
   }

   if (options[opt_repulsion])
      rep=std::atof(options[opt_repulsion]);
   else
      rep=1;

   if (options[opt_z_factor])
      z_factor=std::atof(options[opt_z_factor]);
   else
      z_factor=1;
}

void SpringEmbedder::adjust_edge_lengths(bool different_edge_lengths)
{
   double min_edge_weight, avg_edge_weight;
   if (different_edge_lengths) {
      min_edge_weight=std::numeric_limits<double>::infinity();
      avg_edge_weight=0;
      for (Entire< Edges<graph_type> >::const_iterator e=entire(edges(G)); !e.at_end(); ++e) {
	 if (e->wanted_length <= 0)
	    throw std::runtime_error("non-positive edge length encountered");
	 if (e->wanted_length < min_edge_weight)
	    min_edge_weight=e->wanted_length;
	 avg_edge_weight += e->wanted_length;
      }
      avg_edge_weight/=min_edge_weight;
   } else {
      avg_edge_weight=scale;
      min_edge_weight=1/scale;
   }

   int n_edges=0;
   for (Entire< Edges<graph_type> >::iterator e=entire(edges(G)); !e.at_end(); ++e, ++n_edges)
      e->inv_wanted_length=min_edge_weight/e->wanted_length;
   if (different_edge_lengths)
      avg_edge_weight/=n_edges;

   const double n_nodes=G.nodes();
   scale=avg_edge_weight/4*std::sqrt(n_nodes);
   const double avg_deg=2*n_edges/n_nodes;
   if (avg_deg>=3) {
      const double x=2*M_PI/avg_deg;
      scale*=std::sqrt(std::sin(x)/x);
   }
#ifdef POLY_DEBUG
   cout << "initial adjustment:"
           "\n  min_edge_weight=" << min_edge_weight
	<< "\n  avg_edge_weight=" << avg_edge_weight
	<< "\n  avg_degree=" << avg_deg
	<< "\n  n_nodes=" << n_nodes
	<< "\n  scale=" << scale
	<< endl;
#endif
}

void SpringEmbedder::set_z_ordering(Poly& p, const argv_option *options)
{
   Vector<double> Obj=p.give(options[opt_z_ordering]);
   if (!strcmp(options[opt_z_ordering], "LINEAR_OBJECTIVE")) {
      const Matrix<double> V=p.give("VERTICES");
      Obj=V*Obj;

      const Set<int> far_vertices=p.give("FAR_FACE");
      if (!far_vertices.empty()) {
	 if (options[opt_skip_node_attr]) {
	    // a coarse hack: currently this option is used with BOUNDED_SUBGRAPH only
	    Obj=Obj.slice(~far_vertices);
	 } else {
	    const double max_obj=accumulate(Obj, operations::max()),
		         min_obj=accumulate(Obj, operations::min());
	    for (Entire< Set<int> >::const_iterator f=entire(far_vertices); !f.at_end(); ++f)
	       if (Obj[*f]>0)
		  Obj[*f]=2*max_obj-min_obj;
	       else
		  Obj[*f]=2*min_obj-max_obj;
	 }
      }
   } else if (options[opt_skip_node_attr]) {
      throw std::runtime_error("don't know how to select vertices");
   }

   Entire< Vector<double> >::const_iterator obj=entire(Obj);
   double z_order_min=*obj, z_order_max=z_order_min;
   while (!(++obj).at_end())
      pm::assign_min_max(z_order_min, z_order_max, *obj);

   const double z_mid=(z_order_max+z_order_min)/2;
   z_order_max-=z_order_min;
   if (z_order_max > 1e-3*scale) {
      z_ordering.resize(Obj.size());
      obj=entire(Obj);
      for (Entire< Vector<double> >::iterator z=entire(z_ordering); !z.at_end(); ++z, ++obj)
	 *z=(*obj-z_mid)/z_order_max;
   }
}

void SpringEmbedder::calculate_forces(std::vector<vector3>& forces)
{
   std::vector<vector3>::iterator f=forces.begin();
   Entire< Vector<double> >::const_iterator z=entire(z_ordering);
   double new_z_min=G.node(0).x[2], new_z_max=new_z_min;
   z_max-=z_min;
   if (gravity) barycenter.fill(0.);

   for (Entire< Nodes<graph_type> >::const_iterator this_node=entire(nodes(G));  
	!this_node.at_end(); ++this_node, ++f) {

      f->fill(0.);
      if (gravity) barycenter += this_node->x;

      if (!z.at_end()) {
	 const double this_z=this_node->x[2];
	 (*f)[2] = ((*z)*z_max - this_z) * z_factor;
	 pm::assign_min_max(new_z_min, new_z_max, this_z);
	 ++z;
      }

      int n=this_node.index();
      graph_type::out_edge_list::const_iterator edge=this_node.out_edges().begin();

      for (int n2=0; n2<n; ++n2) {
	 const vector3 delta=G.node(n2).x - this_node->x;
	 const double delta_sqr=std::max(sqr(delta),epsilon_2), delta_abs=std::sqrt(delta_sqr);
	    
	 if (!edge.at_end() && edge.index()==n2) {
	    // the nodes are neighbors
	    const vector3 attraction= (edge->inv_wanted_length - 1/delta_abs) * delta;
	    *f += attraction;
	    forces[n2] -= attraction;
	    ++edge;

	 } else {
	    // the nodes are not neighbors
	    const vector3 repulsion= (rep/delta_sqr/delta_abs) * delta;
	    *f -= repulsion;
	    forces[n2] += repulsion;
	 }
      }
   }

   z_min=new_z_min;  z_max=new_z_max;
   if (gravity) {
      barycenter /= G.nodes();
      if (!z_ordering.empty()) barycenter[2]=0;
   }

   for (Entire< Set<int> >::const_iterator fi = entire(fixed_vertices); !fi.at_end(); ++fi)
      forces[*fi].fill(0);
}

} }

// Local Variables:
// c-basic-offset:3
// End: