/* 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: SpringEmbedderBase.tcc 7579 2007-01-22 11:34:46Z gawrilow $"

#include <vector>
#include <cstdlib>

namespace polymake { namespace graph {

template <typename Top, typename Traits>
void SpringEmbedderBase<Top, Traits>::init_params(const argv_option *options)
{
   if (options[opt_eps])
      epsilon=std::atof(options[opt_eps]);
   else
      epsilon=1e-4;
   epsilon_2=epsilon*epsilon;

   if (options[opt_viscosity])
      viscosity=std::atof(options[opt_viscosity]);
   else
      viscosity=1;

   if (options[opt_inertion])
      inertion=std::atof(options[opt_inertion]);
   else
      inertion=1;
}

static inline
void calc_internal_constants(double& a, double& b, double& c, double& d, double viscosity, double inertion)
{
   a=std::exp(-viscosity/inertion);
   b=(1-a)/viscosity;
   c=inertion*b;
   d=(1-c)/viscosity;
}

template <typename Top, typename Traits>
bool SpringEmbedderBase<Top, Traits>::calculate(int max_iterations)
{
   const int n_nodes=G.nodes();
   std::vector<vector3> forces(n_nodes);

   const double viscosity_increasing_factor=1.2;

   int incr_viscosity=0, decr_viscosity=0;
   double a,b,c,d;
   calc_internal_constants(a,b,c,d,viscosity,inertion);

   for (int iter=0; iter<max_iterations; ++iter) {
      this->top().calculate_forces(forces);
#ifdef POLY_DEBUG
      cout << "iteration " << iter << endl;
#endif

      int oscillated=0, moved=0;
      typename std::vector<vector3>::iterator f=forces.begin();
      for (typename Entire< Nodes<typename Traits::graph_type> >::iterator this_node=entire(nodes(G));
	   !this_node.at_end();  ++this_node, ++f) {
#ifdef POLY_DEBUG
	 cout << '[' << this_node.index() << "]: x=(" << this_node->x << "); f=(" << *f << "); v=(" << this_node->v << ")\n";
#endif
	 vector3 v1=a*this_node->v + b*(*f),
	         dx=c*this_node->v + d*(*f);
	 if (sqr(dx) >= epsilon_2) {
	    if (this_node->v * v1 < 0)
	       ++oscillated;
	    else
	       ++moved;
	 }
	 this_node->v=v1;
	 this_node->x+=dx;
	 if (gravity) this_node->x-=barycenter;
      }

      if (oscillated*2 >= n_nodes) {
	 if (++incr_viscosity>=2) {
	    viscosity *= viscosity_increasing_factor;
	    calc_internal_constants(a,b,c,d,viscosity,inertion);
	    incr_viscosity=0;
#ifdef POLY_DEBUG
	    cout << "++ viscosity ++\n";
#endif
	 }
	 decr_viscosity=0;
      } else {
	 if (!moved) return true;
	 if (!oscillated) {
	    if (++decr_viscosity>=10 && viscosity>1 ||
		decr_viscosity>=200 && viscosity>0.75) {
	       viscosity /= viscosity_increasing_factor;
	       calc_internal_constants(a,b,c,d,viscosity,inertion);
	       decr_viscosity=0;
#ifdef POLY_DEBUG
	    cout << "-- viscosity --\n";
#endif
	    }
	 } else {
	    decr_viscosity=0;
	 }
	 incr_viscosity=0;
      }	    
#ifdef POLY_DEBUG
      cout << "-------------------------------------------\n"
	      "moved=" << moved << " oscillated=" << oscillated << " viscosity=" << viscosity << '\n';
#endif
   }
   return false;
}

} }

// Local Variables:
// mode:C++
// c-basic-offset:3
// End: