# <PLAINTEXT>
##############################################################################
##PROCEDURE(notest)     SpringIntegrate
##AUTHOR                Allan Wittkopf
##DATE                  December 3, 2005
##
##CALLINGSEQ
##- \`dsolve/numeric/SpringIntegrate\`(
##      reledge::list,
##      ini::{integer,hfarray},
##      dim::{identical(2),identical(3)},
##      nframes::nonnegint,
##      maxiter::{posint,name}
##  )
##
##PARAMETERS
##- 'reledge' : `list`; this parameter is an ordered list with each entry
##               containing a list of the vertices that the entry corresponding
##               to the current vertex is connected to.
##               For example, for a complete graph of 3 vertices this is
##               [[2,3],[1,3],[2,3]], meaning vertex 1 is connected to 2,3,
##               vertex 2 is connected to 1,3, etc.
##- 'ini' : This parameter has 3 main values; 0 means to have the code
##          generate structured (grid) initial positions for the vertices;
##          1 means to have the code generate random initial positions for
##          the vertices; and the third input form is an array of initial
##          values for the vertices.
##          For the third form, the data must be a #vertex x #dim C_order
##          float[8] rectangular storage rtable (Matrix or Array).
##          Note that the code is designed so that the final values for an
##          animation approx. fill the range [-1..1$dim], so initial positions
##          should be specified to fill the same range, as otherwise the
##          graph will either grow or shrink in an undesirable way during the
##          animation process.
##- 'dim' : This is the number of dimensions in which the graph is to be
##          built, and must be either 2 or 3.
##- 'nframes' : This should be 0 if only the final values are desired
##              (and these are returned in a #vertex x #dim array). For
##              values > 0, the computation proceeds as an animation
##              (3 times as expensive) and returns the data in a
##              #frames x #vertex x #dim array.
##- 'maxiter' : This is an iteration maximum - it serves to limit the amount
##              of work performed for an animation. A value of 20000 is
##              suitable for any graphs up to the 'extremely difficult'
##              level. This can be specified in 2 ways:
##              1) Just pass a value.
##                 This is then a limit.
##              2) Pass a value through an assigned quoted variable.
##                 This method also sets the number of iterations performed
##                 at output time.
##                 **If** the graph did not complete, then this value will
##                 be 1 greater than the input max iter count on output.
##              The max allowed value is 2^31-1.
##
##RETURNS
##- the procedure returns depend upon the parameters. There are the following
##  cases:
##- nframes=0: returns a #vertex x #dim C_order, rectangular, float[8] rtable
##      that represents the final value computed for the graph.
##- nframes=1: returns a #frames x #vertex x #dim C_order, rectangular,
##      float[8] rtable that represents a smooth transition from the initial
##      value to the final value computed for the graph.
##- For both cases if maxiter is a name, and is assigned to a value no greater
##  than the input maxiter value, then the graph has stabilized, and the answer
##  is final. If the value is greater, then the graph did not stabilize in
##  the specified number of steps.
##      
##
##DESCRIPTION
##- The \`dsolve/numeric/SpringIntegrate\` procedure uses numerical methods,
##  combined with a specialized model system, to integrate a graph from
##  either specified, structured, or random initial positions to a final
##  state where the model system is in a steady state (stable).
##- The model system has been designed so that the steady state is visually
##  appealing based on the following 2 criteria: 1) No overlapping vertices;
##  2) When possible, connected vertices are close to each other.
##
##- Note that if the graph is not connected, a steady state can *never* be
##  reached, as vertices repel one another, and edges counter-balance the
##  effect.
##- Note also that the code is not designed to account-for or ignore
##  connections from a vertex to itself, so the code will simply fail in
##  this case.
##
##SUBSECTION The Model
##- The model is basically a second order ODE system force balance model
##  in 'dim' dimensions. The following are the rules:
##- Each vertex exerts an attractive force to any other vertex to which
##  it is attached by an edge. This force is proportional to the distance
##  between the vertices.
##- Each vertex exerts a constant repulsive force to all other vertices,
##  regardless of whether they are connected by an edge.
##- All vertex motion proceeds in the presence of friction that serves to
##  damp oscillations in the model.
##- Finally, there is one other small force used to 'center' the model
##  (a fictitious force that attempts to center the average positions in
##  the model at the origin.
##
##- The differential equation for each vertex can be described as follows:
##  x[i]\'\'=-a\*sum((x[i]-x[j])/(n[i]+n[j]),j=edges[i])-b\*x[i]\'
##         +c\*sum((x[i]-x[j])/|x[i]-x[j]|,j=1..i-1,i+1..n)-0.001\*xav,
##  where the x[i] is a vector for the position of vertex i, n[i]
##  counts the number of edges connected to vertex i, edges[i] gives
##  the other vertex for each edge connected to vertex i, and xav is the
##  average positions over all nodes.
##- The values of the constants 'a','b','c' are internally managed by the
##  code, but have been optimized (via a very lengthy derivation) to provide
##  stabilization of the model in as little time as possible, and to keep
##  the scale uniform when possible.
##- Note that one may question the use of a constant force repulsion, but
##  this is critical for several reasons. One of the two major reasons is that
##  a constant force model allows for cross-overs in 2-d problems (i.e. one
##  vertex can pass right through another). This makes the dependence of the
##  model steady state much less dependent on the initial positions chosen for
##  the model. Clearly an advantage when the initial values are random.
##  Use of an inverse or inverse square repulsion law does not allow this (as
##  it prevents nodes from getting close enough to allow a crossing).
##  The other major reason is efficiency. It was found that for large
##  graphs, use of an inverse law took many times longer to stabilize.
##  The constant force model is closer to a linear model, and as a result
##  has better behavior with respect to integration to a steady state.
##ENDSUBSECTION
##
##SUBSECTION Final Fix-up
##- The model in use has one fatal flaw.
##  Specifically since the repulsive forces between two vertices are
##  constant, it is possible to have overlapping vertices.
##  This is clearly undesirable.
##- To fix this, there is a final phase of the integration that, after
##  the primary steady state is achieved, checks the state for vertices that
##  overlap, or are too close. The criteria for this is a little complex, but
##  is basically concerned with the average of the minimum vertex distances
##  above a fixed bar, and enforcing that the minimum node distance is at
##  least 0.6 times that average.
##- If these are detected, then the constant repulsive force between any pair
##  of close vertices is increased, and a few more integration steps are taken,
##  until a steady state is achieved with no close/overlapping vertices.
##- It should again be noted that an attempt was made to change the model from
##  the current model to an inverse law model *after* the steady state was
##  achieved, but it was found that the desirable steady state was an unstable
##  stationary state for some models of interest under an inverse law, so this
##  was abandoned in favor of the current criteria.
##ENDSUBSECTION
##
##SUBSECTION Numerical Solution of the Model
##- The model is numerically solved using an order 4 Adams solver, utilizing
##  an rk2 method as a startup for the multi-step integration.
##  This was chosen to minimize the number of function evaluations performed
##  by the algorithm (multi-step req. only 1 function eval per step).
##- The step size used depends upon the stability of the model, but given the
##  choices made for optimization of a,b,c, it turns out that the optimal
##  step size is nearly always in the 0.4-0.6 range.
##- Basically the model proceeds by integrating the system forward in time
##  until the distance traveled by the vertices in a single step drops below
##  a fixed minimum value of 1e-8. Of course the first few steps are excluded
##  from the test, as the initial velocities in the model are always zero,
##  and the system takes a little time to get 'up to speed' as it were.
##- Also the second phase of the numerical integration proceeds in much the
##  same manner.
##ENDSUBSECTION
##
##SUBSECTION Animation Details
##- Animations are computed a little bit differently than just the steady
##  state. The algorithm needs to perform up to 3 passes before producing the
##  final animation:
##- Step 1: Integrate system forward in time to reach the steady state, and
##          note the boundaries of the region for the final state.
##          If these boundaries are within 2% of -1..1, -1..1, -1..1 then
##          proceed to step 3, otherwise adjust 'a','b','c' and proceed to
##          step 2.
##- Step 2: Integrate for the adjusted constants.
##- Step 3: Use the computed total distance traveled in the course of the
##          computation to make the animation uniform in the sense:
##          \"the sum of the absolute values of the distance changed between
##           any two consecutive frames is equal\".
##         Note that this is quite relevant, as the rate of change in the
##         model varies massively, and an equi-spaced time sampling does
##         not produce an esthetically pleasing animation.
##- There are disadvantages to the 'total distance traveled' approach when
##  dealing with small graphs containing very few nodes. The animation may
##  not appear very smooth, as initially the model is stabilizing with
##  respect to relative location of the nodes, and the small force that is
##  used to center the graph is much slower acting.
##  This means that the graph may appear to stabilize first by relative
##  position, then slowly shift to being centered.
##  This is something of a judgment call, as for the equal spaced time
##  approach the first few frames of animation will contain nearly all
##  the movement, and the last 20 will be small shifts - this is also ugly
##  but it is ugly for nearly all graphs.
##ENDSUBSECTION
##

`dsolve/numeric/SpringIntegrate` := proc()
local lib;
global `dsolve/numeric/SpringIntegrate`;
	try
		lib := ExternalCalling:-ExternalLibraryName("ode2",'HWFloat');
		`dsolve/numeric/SpringIntegrate` :=
			ExternalCalling:-DefineExternal('hw_SpringIntegrate',lib);
   catch:
      error("external library %1 could not be found/used",lib);
   end try;
   `dsolve/numeric/SpringIntegrate`(args[1..nargs]);
end proc:

#savelib('`dsolve/numeric/SpringIntegrate`');
##############################################################################