#
##############################################################################
##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`');
##############################################################################