[Vegan-commits] r1565 - in pkg/vegan: . R inst man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Apr 1 15:15:12 CEST 2011
Author: jarioksa
Date: 2011-04-01 15:15:12 +0200 (Fri, 01 Apr 2011)
New Revision: 1565
Added:
pkg/vegan/R/monoMDS.R
pkg/vegan/man/monoMDS.Rd
pkg/vegan/src/monoMDS.f
Modified:
pkg/vegan/DESCRIPTION
pkg/vegan/inst/ChangeLog
pkg/vegan/man/vegan-package.Rd
Log:
add monoMDS: Peter Minchin's NMDS code
Modified: pkg/vegan/DESCRIPTION
===================================================================
--- pkg/vegan/DESCRIPTION 2011-03-31 18:34:49 UTC (rev 1564)
+++ pkg/vegan/DESCRIPTION 2011-04-01 13:15:12 UTC (rev 1565)
@@ -1,10 +1,10 @@
Package: vegan
Title: Community Ecology Package
-Version: 1.18-27
-Date: March 31, 2011
+Version: 1.18-28
+Date: April 1, 2011
Author: Jari Oksanen, F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre,
- R. B. O'Hara, Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens,
- Helene Wagner
+ Peter R. Minchin, R. B. O'Hara, Gavin L. Simpson, Peter Solymos,
+ M. Henry H. Stevens, Helene Wagner
Maintainer: Jari Oksanen <jari.oksanen at oulu.fi>
Suggests: MASS, mgcv, lattice, cluster, scatterplot3d, rgl, tcltk
Description: Ordination methods, diversity analysis and other
Added: pkg/vegan/R/monoMDS.R
===================================================================
--- pkg/vegan/R/monoMDS.R (rev 0)
+++ pkg/vegan/R/monoMDS.R 2011-04-01 13:15:12 UTC (rev 1565)
@@ -0,0 +1,73 @@
+monoMDS <-
+ function(dist, y, k = 2,
+ model = c("global", "local", "hybrid"), threshold = 0.8,
+ maxit = 200, tol = 0.0001, ...)
+{
+ model <- match.arg(model)
+ if (model == "global") {
+ ## global NMDS: lower triangle
+ mat <- as.matrix(dist)
+ dist <- mat[lower.tri(mat)]
+ iidx <- row(mat)[lower.tri(mat)]
+ jidx <- col(mat)[lower.tri(mat)]
+ iregn <- 1
+ ngrp <- 1
+ nobj <- nrow(mat)
+ istart <- 1
+ } else if (model == "local") {
+ ## local NMDS: whole matrix without the diagonal, and rows in
+ ## a row (hence transpose)
+ mat <- t(as.matrix(dist))
+ dist <- mat[col(mat) != row(mat)]
+ iidx <- col(mat)[col(mat) != row(mat)] # transpose!
+ jidx <- row(mat)[col(mat) != row(mat)]
+ iregn <- 1
+ nobj <- nrow(mat)
+ ngrp <- nobj
+ istart <- seq(1, length(dist), by = (nobj-1))
+ } else if (model == "hybrid") {
+ ## Hybrid NMDS: two lower triangles, first a complete one,
+ ## then those with dissimilarities below the threshold
+ mat <- as.matrix(dist)
+ dist <- mat[lower.tri(mat)]
+ iidx <- row(mat)[lower.tri(mat)]
+ jidx <- col(mat)[lower.tri(mat)]
+ ## second group: dissimilarities below threshold
+ ngrp <- 2
+ istart <- c(1, length(dist) + 1)
+ take <- dist < threshold
+ dist <- c(dist, dist[take])
+ iidx <- c(iidx, iidx[take])
+ jidx <- c(jidx, jidx[take])
+ iregn <- 3
+ nobj <- nrow(mat)
+ }
+ ## ndis: number dissimilarities
+ ndis <- length(dist)
+ ## starting configuration
+ if (missing(y)) {
+ y <- matrix(runif(nobj*k, -1, 1), nobj, k)
+ ## centre
+ y <- sweep(y, 2, colMeans(y), "-")
+ }
+ ## y to vector
+ y <- as.vector(as.matrix(y))
+ ## Fortran call
+ sol <- .Fortran("monoMDS", nobj = as.integer(nobj), nfix=as.integer(0),
+ ndim = as.integer(k), ndis = as.integer(ndis),
+ ngrp = as.integer(ngrp), diss = as.double(dist),
+ iidx = as.integer(iidx), jidx = as.integer(jidx),
+ xinit = as.double(y), istart = as.integer(istart),
+ isform = as.integer(1), ities = as.integer(1),
+ iregn = as.integer(iregn), iscal = as.integer(1),
+ maxits = as.integer(maxit),
+ sratmx = as.double(0.99999), strmin = as.double(tol),
+ sfgrmn = as.double(1e-7), dist = double(ndis),
+ dhat = double(ndis), points = double(k*nobj),
+ stress = double(1), iters = integer(1),
+ icause = integer(1), PACKAGE = "vegan")
+ sol$points <- matrix(sol$points, nobj, k)
+ class(sol) <- "monoMDS"
+ sol
+}
+
Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog 2011-03-31 18:34:49 UTC (rev 1564)
+++ pkg/vegan/inst/ChangeLog 2011-04-01 13:15:12 UTC (rev 1565)
@@ -2,8 +2,22 @@
VEGAN DEVEL VERSIONS at http://r-forge.r-project.org/
-Version 1.18-27 (opened March 31, 2011)
+Version 1.18-28 (opened April 1, 2011)
+ * Peter Minchin joined the vegan team.
+
+ * monoMDS: a new function with Peter Minchin's Fortran90 code for
+ NMDS. The full Fortran90 code has global, local, linear and hybrid
+ MDS, configurable and valid tie treatment, handles missing values,
+ and allows adding new points to existing ordinations. The Fortran
+ code is highly tuned, and much faster than other alternatives in
+ R. The current R interface allows only access to a subset of the
+ capabilities of the Fortran90 code, but the missing pieces will be
+ added gradually. The function will eventually replaces isoMDS() of
+ the MASS package as the main NMDS engine in metaMDS().
+
+Version 1.18-27 (closed April 1, 2011)
+
* orderingKM: Kurt Hornik found a problem when inspecting the
Fortran subroutines with gcc 4.6 tools when checking vegan release
1.17-9.
Added: pkg/vegan/man/monoMDS.Rd
===================================================================
--- pkg/vegan/man/monoMDS.Rd (rev 0)
+++ pkg/vegan/man/monoMDS.Rd 2011-04-01 13:15:12 UTC (rev 1565)
@@ -0,0 +1,65 @@
+\name{monoMDS}
+\alias{monoMDS}
+
+\title{
+ Kruskal's Non-metric Multidimensional Scaling
+}
+\description{
+ Function implements Kruskal's (1964a,b) non-metric multidimensional scaling
+ using monotone regression and primary (\dQuote{weak}) treatment of ties.
+}
+\usage{
+monoMDS(dist, y, k = 2, model = c("global", "local", "hybrid"),
+ threshold = 0.8, maxit = 200, tol = 1e-04, ...)
+}
+
+\arguments{
+ \item{dist}{Input dissimilarities.}
+ \item{y}{Starting configuration. A random configuration will be
+ generated if this is missing.}
+ \item{k}{Number of dimensions.}
+ \item{model}{MDS model: \code{"global"} is normal non-metric MDS
+ with a monotone regression, \code{"local"} is non-metric MDS with
+ separate regressions for each point, and \code{"hybrid"} uses
+ linear MDS for dissimilarities below the threshold, and global
+ non-metric MDS for dissimilarities above the threshold. }
+ \item{threshold}{Dissimilarity above which monotone regression is
+ used instead of linear regression. }
+ \item{maxit}{Maximum number of iterations.}
+ \item{tol}{Convergence tolerance.}
+ \item{\dots}{Other parameters to the functions (ignored).}
+}
+
+\details{
+ This is the first version of the function, and still under construction.
+}
+
+\value{
+ A complicated object (under construction).
+}
+
+\references{
+ Kruskal, J.B. 1964a. Multidimensional scaling by optimizing
+ goodness-of-fit to a nonmetric hypothesis. \emph{Psychometrika}
+ 29, 1--28.
+
+ Kruskal, J.B. 1964b. Nonmetric multidimensional scaling: a numerical
+ method. \emph{Psychometrika} 29, 115--129.
+}
+
+\author{
+Peter R. Michin (Fortran core) and Jari Oksanen (R interface).
+}
+
+\note{
+ Still under development.
+}
+
+\seealso{ \code{\link[vegan]{metaMDS}}, \code{\link[MASS]{isoMDS}},
+ \code{\link[smacof]{smacofSym}}, \code{\link[ecodist]{nmds}}. }
+
+\examples{
+##
+}
+\keyword{ multivariate }
+
Modified: pkg/vegan/man/vegan-package.Rd
===================================================================
--- pkg/vegan/man/vegan-package.Rd 2011-03-31 18:34:49 UTC (rev 1564)
+++ pkg/vegan/man/vegan-package.Rd 2011-04-01 13:15:12 UTC (rev 1565)
@@ -45,10 +45,10 @@
}
\author{ The \pkg{vegan} development team is Jari Oksanen,
-F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre, R. B. O'Hara,
-Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens, Helene
-Wagner. Many other people have contributed to individual functions: see
-credits in function help pages.
+F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre, Peter
+R. Minchin, R. B. O'Hara, Gavin L. Simpson, Peter Solymos, M. Henry
+H. Stevens, Helene Wagner. Many other people have contributed to
+individual functions: see credits in function help pages.
The maintainers at the R-Forge are Jari Oksanen <jari.oksanen at oulu.fi>
and Gavin Simpson <gavin.simpson at ucl.ac.uk>.
Added: pkg/vegan/src/monoMDS.f
===================================================================
--- pkg/vegan/src/monoMDS.f (rev 0)
+++ pkg/vegan/src/monoMDS.f 2011-04-01 13:15:12 UTC (rev 1565)
@@ -0,0 +1,1117 @@
+ SUBROUTINE monoMDS (NOBJ, NFIX, NDIM, NDIS, NGRP,
+ . DISS, IIDX, JIDX, XINIT, ISTART,
+ . ISFORM, ITIES, IREGN, ISCAL, MAXITS, SRATMX,
+ . STRMIN, SFGRMN,
+ . DIST, DHAT, X, STRESS, ITERS, ICAUSE)
+C
+C Subroutine for multidimensional scaling.
+C
+C Written by Dr. Peter R. Minchin
+C Department of Biological Sciences
+C Southern Illinois University Edwardsville
+C PO Box 1651
+C Edwardsville, IL 62026-1651, U.S.A.
+C Phone: +1-618-650-2975 FAX: +1-618-650-3174
+C Email: pminchi at siue.edu
+C
+C Starting from a supplied initial configuarion, uses steepest descent
+C to minimize Kruskal's stress, a measure of badness-of-fit of one
+C or more regressions of distances onto the supplied dissimilarities.
+C
+C This routine is very flexible in the form of the dissimilarity
+C matrices it can handle because dissimilarities are input as a
+C vector (DISS) with the associated i- and j-indices in parallel
+C vectors (IIDX, JIDX). There is no need for missing values, as
+C these index vectors specify which pair of points each dissimilarity
+C belongs to and missing dissimilarities are simply omitted. Another
+C element of flexibility is that there can be more than one group of
+C dissimilarities and a separate regression of distance on
+C dissimilarity is done for each group. Local or "row conditional"
+C scaling can be performed by inputting each row of the dissimilarity
+C matrix (minus the diagonals) as a separate group. Hybrid scaling
+C of Faith et. al (1987) can be done by inputting the complete
+C dissimilarity matrix as group 1 and only those dissimilarities
+C below a specified threshold as group 2 (and also setting input
+C parameter IREGN to 3). New points can be inserted into an
+C existing ordination by inputting dissimilarities from the
+C rectangular matrix in which the existing points are rows and the
+C new points are columns. The parameter NFIX is then set to the
+C number of existing points (rows) and these will be the first
+C NFIX points in XINIT, which supplies their coordinates. Normally,
+C when NFIX>0, the scaling parameter, ISCAL, should be 0, so that
+C the coordinates of the fixed points remain unchanged in the
+C output ordination scores, X.
+C
+C========INPUT ARGUMENTS:
+C
+C NOBJ = number of objects being ordinated
+C
+C NFIX = number of fixed objects (< NOBJ - if NFIX>0 the input
+C dissim matrix is assumed to be rectangular, with the
+C first NFIX objects, rows, remaining fixed and the rest
+C able to move)
+C
+C NDIM = number of dimensions for the MDS ordination
+C
+C NDIS = total number of dissimilarities among objects
+C
+C NGRP = number of groups of dissimilarities (a separate regression
+C of distances on dissimilarities is performed for each
+C group)
+C
+C DISS(NDIS) = vector of dissimilarities among objects
+C
+C IIDX(NDIS) = i-indices of objects in the dissimilarity vector
+C
+C JIDX(NDIS) = j-indices of objects in the dissimilarity vector
+C
+C XINIT(NOBJ,NDIM) = initial coordinates for NOBJ objects on NDIM
+C dimensions
+C
+C ISTART(NGRP) = index of starting position of each group of
+C dissimilarities within vectors DISS, IIDX and JIDX
+C
+C ISFORM = Kruskal stress formula to be used: 1 or 2
+C
+C ITIES = Kruskal tie treatment: 1 = weak/primary, 2 = strong/secondary
+C
+C IREGN = Regression type: 1 = Kruskal monotone, 2 = linear, 3 = hybrid
+C (Kruskal monotone for the first NGRP/2 groups, linear for
+C the remainder of the groups - when this is selected it is
+C assumed that NGRP>1 and NGRP is even)
+C
+C ISCAL = scaling of final ordination: 0 = preserved (output raw scores
+C with no standardization), 1 = normalize (center and adjust
+C RMS distance of points to centroid to 1.0)
+C
+C MAXITS = maximum number of iterations - at least 200 recommended
+C
+C SRATMX = maximum average stress ratio (iterations stop when average
+C stress ratio exceeds this value but is still < 1) - a value
+C of 0.9999 or higher is recommended
+C
+C STRMIN = minimum stress to stop (iterations stop if stress falls
+C below this value) - a value of 0.001 or lower is
+C recommended
+C
+C SFGRMN = minimum scale factor of the gradient (iterations stop if
+C the scale factor of the gradient drops below this)
+C
+C========OUTPUT ARGUMENTS:
+C
+C DIST(NDIS) = distances among objects in the final ordination
+C
+C DHAT(NDIS) = fitted distances in final regression of distance
+C on dissimilarity
+C
+C X(NOBJ,NDIM) = final ordination coordinates of NOBJ objects in NDIM
+C dimensions
+C
+C STRESS = final value of stress
+C
+C ITERS = number of iterations performed
+C
+C ICAUSE = reason for termination of iterations: 1 = max iterations
+C used, 2 = stress fell below STRMIN, 3 = stress ratio
+C exceeded SRATMX, 4 = scale factor of gradient fell below
+C SFGRMN
+C
+C---INPUT ARGUMENTS
+ INTEGER, INTENT(IN) :: NOBJ, NFIX, NDIM, NDIS, NGRP,
+ . ISFORM, ITIES, IREGN, ISCAL, MAXITS
+ INTEGER, INTENT(IN) :: IIDX(NDIS), JIDX(NDIS), ISTART(NGRP)
+ DOUBLE PRECISION, INTENT(IN) :: XINIT(NOBJ,NDIM), DISS(NDIS),
+ . SRATMX, STRMIN, SFGRMN
+C---OUTPUT ARGUMENTS
+ INTEGER, INTENT(OUT) :: ITERS, ICAUSE
+ DOUBLE PRECISION, INTENT(OUT) :: X(NOBJ,NDIM), DIST(NDIS),
+ . DHAT(NDIS), STRESS
+C---ALLOCATABLE TEMPORARY ARRAYS
+ INTEGER, ALLOCATABLE :: IWORK(:)
+ DOUBLE PRECISION, ALLOCATABLE :: GRAD(:,:), GRLAST(:,:), STRS(:)
+C
+ DOUBLE PRECISION :: STRLST, SQRTN, SRATF1, SRATF2, FNGRP,
+ . STRINC, COSAV, ACOSAV, SRATAV, STEP, FNDIM, SFGR, SRATIO,
+ . SFACT, TFACT, DMEAN, CAGRGL, SFGLST, STPNEW, SSFACB, SSFACT,
+ . PARAM(2)
+C
+C ALLOCATE THE TEMPORARY ARRAYS NEEDED
+C
+ ALLOCATE (IWORK(NDIS), GRAD(NOBJ,NDIM), GRLAST(NOBJ,NDIM),
+ . STRS(NGRP))
+C
+C INITIALIZE SOME PARAMETERS
+C
+ SQRTN=SQRT(REAL(NDIS))
+ SRATF1=0.5*(1.0+SRATMX)
+ SRATF2=1.0-SRATF1
+ FNGRP=REAL(NGRP)
+ STRINC=1.1
+ STRESS=1.0D0
+ COSAV=0.0
+ ACOSAV=0.0
+ SRATAV=0.8
+ STEP=0.0
+ FNDIM=REAL(NDIM)
+ SFGR=SQRTN
+C
+C PREPARE DISSIMILARITIES FOR USE IN REGRESSION(S)
+C
+ DO IGRP=1,NGRP
+ I1=ISTART(IGRP)
+ IF (IGRP.LT.NGRP) THEN
+ N=ISTART(IGRP+1)-I1
+ ELSE
+ N=NDIS+1-I1
+ ENDIF
+ IF (N.GT.0) THEN
+C
+C PRE-SORT IF NECESSARY. NOT REQUIRED FOR LINEAR REGRESSION (IREGN=2)
+C NOR FOR DISSIM GROUPS > NGRP/2 IN HYBRID REGRESSION (IREGN=3)
+C
+ IF (IREGN.EQ.1.OR.(IREGN.EQ.3.AND.IGRP.LE.NGRP/2)) THEN
+ CALL ASORT4 (DISS(I1),N,IIDX(I1),JIDX(I1))
+ ENDIF
+ ENDIF
+ ENDDO
+C
+C COPY INITIAL CONFIGURATION TO CURRENT CONFIGURATION
+C
+ CALL MACOPY (XINIT,NOBJ,NOBJ,NDIM,X,NOBJ)
+C
+C INITALIZE GRADIENT (WILL BE USED AS THE FIRST "LAST GRADIENT")
+C
+ CALL MAINIT (GRAD,NOBJ,NDIM,NOBJ,SQRT(1.0/FNDIM))
+C=======================================================================
+C
+C START OF ITERATION LOOP
+C
+ DO IT=0,MAXITS
+C---Reset backup counter
+ NBACK=0
+C
+C SAVE LAST STRESS VALUE AND SET STRESS BACK TO ZERO
+C
+C---Jump back to here when a backup if required
+ 10 STRLST=STRESS
+ STRESS=0.0D0
+C
+C Normalize the current configuration (unless ISCAL=0)
+C
+ IF (ISCAL.GT.0) THEN
+ CALL NRMCON (X,NOBJ,NDIM,NOBJ,SSFACT)
+ ELSE
+ SSFACT=1.0D0
+ ENDIF
+C
+C MOVE CURRENT GRADIENT TO LAST GRADIENT (WITH APPROPRIATE SCALING)
+C AND SAVE LAST GRADIENT SCALING FACTOR
+C
+ CALL MACOPY (GRAD,NOBJ,NOBJ,NDIM,GRLAST,NOBJ)
+ CALL MAMAS (GRLAST,NOBJ,NOBJ,NDIM,SSFACT)
+ SFGLST=SFGR*SSFACT
+C
+C CLEAR CURRENT GRADIENT
+C
+ CALL MAINIT (GRAD,NOBJ,NDIM,NOBJ,0.0D0)
+C
+C COMPUTE DISTANCES
+C
+ CALL CLCDIS (X,NOBJ,NDIM,DIST,IIDX,JIDX,NDIS)
+C-----------------------------------------------------------------------
+C
+C LOOP OVER THE NGRP GROUPS OF DISSIMILARITIES
+C
+ DO IGRP=1,NGRP
+ I1=ISTART(IGRP)
+ IF (IGRP.LT.NGRP) THEN
+ N=ISTART(IGRP+1)-I1
+ ELSE
+ N=NDIS+1-I1
+ ENDIF
+ IF (N.GT.0) THEN
+C
+C PERFORM REGRESSION OF DISTANCES ON DISSIMILARITIES
+C
+ IF (IREGN.EQ.1.OR.(IREGN.EQ.3.AND.IGRP.LE.NGRP/2)) THEN
+ CALL MONREG (DISS(I1),DIST(I1),DHAT(I1),IIDX(I1),
+ . JIDX(I1),IWORK(I1),N,ITIES)
+ ELSEIF (IREGN.EQ.2.OR.(IREGN.EQ.3.AND.IGRP.GT.NGRP/2)) THEN
+ CALL LINREG (DISS(I1),DIST(I1),DHAT(I1),N,PARAM)
+ ENDIF
+C
+C COMPUTE AND ACCUMULATE STRESS
+C
+ CALL CLCSTR (DIST(I1),DHAT(I1),N,SFACT,TFACT,STRS(IGRP),
+ . ISFORM,DMEAN)
+ STRESS=STRESS+SFACT/TFACT
+C
+C ACCUMULATE THE NEGATIVE GRADIENT
+C
+ CALL CLCGRD (X,GRAD,NOBJ,NDIM,DIST(I1),DHAT(I1),
+ . IIDX(I1),JIDX(I1),N,STRS(IGRP),SFACT,TFACT,ISFORM,DMEAN)
+ ENDIF
+ ENDDO
+C
+C END OF DISSIMILARITY GROUP LOOP
+C
+C-----------------------------------------------------------------------
+C
+C COMPUTE OVERALL STRESS AND ITS REDUCTION RATIO
+C
+ STRESS=SQRT(STRESS/FNGRP)
+ IF (IT.EQ.0) THEN
+ SRATIO=0.8D0
+ ELSE
+ SRATIO=STRESS/STRLST
+ ENDIF
+C
+C TO PREVENT FIXED POINTS BEING MOVED, ZERO THEIR GRADIENT ELEMENTS
+C
+ IF (NFIX.GT.0) THEN
+ DO I=1,NFIX
+ DO IDIM=1,NDIM
+ GRAD(I,IDIM)=0.0D0
+ ENDDO
+ ENDDO
+ ENDIF
+C
+C CALCULATE SCALE FACTOR OF NEW GRADIENT AND ITS DIRECTION COSINE
+C WITH THE LAST GRADIENT
+C
+ CALL CLCSFA (GRAD,GRLAST,NOBJ,NDIM,NOBJ,SFGR,CAGRGL,
+ . SFGLST)
+C
+C IF STRESS HAS INCREASED APPRECIABLY, OR THE ANGLE BETWEEN THE NEW
+C GRADIENT AND THE LAST GRADIENT IS APPROACHING 180 DEGREES, REDUCE
+C THE STEP SIZE AND BACK UP ALONG THE LAST GRADIENT
+C (BACKUP PROCEDURE IS PERFORMED A MAXIMUM OF 3 TIMES PER ITERATION)
+C
+ IF ((SRATIO.GT.STRINC.OR.CAGRGL.LT.(-0.95)).AND.NBACK.LT.3)
+ . THEN
+ STPNEW=STEP*0.1D0
+ CALL BACKUP (X,GRAD,GRLAST,NOBJ,NDIM,NOBJ,NBACK,SSFACT,
+ . SSFACB,STRESS,STRLST,SFGR,SFGLST,STEP,STPNEW)
+ GO TO 10
+ ENDIF
+C
+C UPDATE AVERAGES OF STRESS RATIO AND PREVIOUS DIRECTION COSINES
+C
+ SRATAV=SRATIO**0.33334D0 * SRATAV**0.66666D0
+ COSAV=CAGRGL*0.66666D0 + COSAV*0.33334D0
+ ACOSAV=ABS(CAGRGL)*0.66666D0 + ACOSAV*0.33334D0
+C
+C COMPUTE NEW STEP SIZE
+C
+ CALL CLCSTP (STEP,IT,SFGR,STRESS,COSAV,ACOSAV,SRATIO,
+ . SRATAV)
+C
+C CHECK WHETHER OR NOT TO KEEP ITERATING
+C
+C IF STRESS IN SMALL ENOUGH, TERMINATE
+C
+ ITERS=IT
+ IF (STRESS.LT.STRMIN) THEN
+ ICAUSE=2
+ EXIT
+C
+C IF STRESS IS DECREASING SLOWLY, TERMINATE.
+C
+ ELSEIF (ABS(SRATAV-SRATF1).LE.SRATF2.AND.
+ . ABS(SRATIO-SRATF1).LE.SRATF2) THEN
+ ICAUSE=3
+ EXIT
+C
+C IF SCALE FACTOR OF GRADIENT IS SMALL ENOUGH, TERMINATE
+C
+ ELSEIF (SFGR.LE.SFGRMN) THEN
+ ICAUSE=4
+ EXIT
+C
+C IF THIS IS THE FINAL ITERATION, TERMINATE
+C
+ ELSEIF (IT.EQ.MAXITS) THEN
+ ICAUSE=1
+ EXIT
+ ENDIF
+C
+C COMPUTE NEW CONFIGURATION
+C
+ CALL NEWCON (X,GRAD,NOBJ,NDIM,NOBJ,STEP,SFGR)
+ ENDDO
+C
+C END OF ITERATION LOOP
+C
+C=======================================================================
+C
+C DEALLOCATE THE TEMPORARY ARRAYS AND RETURN
+C
+ DEALLOCATE (IWORK, GRAD, GRLAST, STRS)
+ RETURN
+ END SUBROUTINE monoMDS
+
+ SUBROUTINE ASORT4 (X,N,I1,I2)
+C
+C SORTS THE REAL*8 VECTOR X ASCENDING, SIMULTANEOUSLY PERMUTING THE
+C INTEGER VECTORS I1 AND I2 IN PARALLEL.
+C
+C ADAPTED FROM SUPERCHARGED SHELL SORT ALGORITHM PUBLISHED ON PAGE 488
+C OF THE MAY 1983 EDITION OF BYTE MAGAZINE.
+C
+C Adapted by Dr. Peter R. Minchin
+C Department of Biological Sciences
+C Southern Illinois University Edwardsville
+C PO Box 1651
+C Edwardsville, IL 62026-1651, U.S.A.
+C Phone: +1-618-650-2975 FAX: +1-618-650-3174
+C Email: pminchi at siue.edu
+C
+ INTEGER I1(N), I2(N), I1TEMP, I2TEMP
+ DOUBLE PRECISION X(N), XTEMP
+C
+ IF (N.LT.2) RETURN
+ FN=REAL(N)
+ NLOOPS=MAX(NINT(LOG(FN)/LOG(2.)),1)
+ M=2**(NLOOPS-1)
+ DO II=1,NLOOPS
+ FM=M
+ DO I=1,MAX(1,N-M)
+ IF (X(I).GT.X(I+M)) THEN
+ XTEMP=X(I+M)
+ I1TEMP=I1(I+M)
+ I2TEMP=I2(I+M)
+ X(I+M)=X(I)
+ I1(I+M)=I1(I)
+ I2(I+M)=I2(I)
+ IF (I.LE.M) THEN
+ X(I)=XTEMP
+ I1(I)=I1TEMP
+ I2(I)=I2TEMP
+ ELSE
+ DO J=I-M,1,-M
+ IF (XTEMP.GE.X(J)) EXIT
+ X(J+M)=X(J)
+ I1(J+M)=I1(J)
+ I2(J+M)=I2(J)
+ ENDDO
+ X(J+M)=XTEMP
+ I1(J+M)=I1TEMP
+ I2(J+M)=I2TEMP
+ ENDIF
+ ENDIF
+ ENDDO
+ M=INT(FM/2.)
+ ENDDO
+ RETURN
+ END SUBROUTINE ASORT4
+
+ SUBROUTINE BACKUP (X,GRAD,GRLAST,NOBJ,NDIM,MAXOBJ,NBACK,SSFACT,
+ . SSFACB,STRESS,STRLST,SFGR,SFGLST,STEP,STPNEW)
+C
+C BACKS THE CONFIGURATION UP ALONG THE LAST GRADIENT
+C
+C Written by Dr. Peter R. Minchin
+C Department of Biological Sciences
+C Southern Illinois University Edwardsville
+C PO Box 1651
+C Edwardsville, IL 62026-1651, U.S.A.
+C Phone: +1-618-650-2975 FAX: +1-618-650-3174
+C Email: pminchi at siue.edu
+C
+ DOUBLE PRECISION X(MAXOBJ,NDIM), GRAD(MAXOBJ,NDIM),
+ . GRLAST(MAXOBJ,NDIM), SSFACT, SSFACB, STRESS, STRLST, SFGR,
+ . SFGLST,STEP,STPNEW,FACTOR
+C
+ NBACK=NBACK+1
+ IF (NBACK.EQ.1) THEN
+ SSFACB=1.0D0
+ ELSE
+ SSFACB=SSFACB*SSFACT
+ ENDIF
+ FACTOR=SSFACB*(STEP-STPNEW)/SFGLST
+ DO IDIM=1,NDIM
+ DO IOBJ=1,NOBJ
+ X(IOBJ,IDIM)=X(IOBJ,IDIM)-FACTOR*GRLAST(IOBJ,IDIM)
+ GRAD(IOBJ,IDIM)=GRLAST(IOBJ,IDIM)
+ ENDDO
+ ENDDO
+ STEP=STPNEW
+ SFGR=SFGLST
+ STRESS=STRLST
+ RETURN
+ END SUBROUTINE BACKUP
+
+ SUBROUTINE CLCDIS (X,MAXOBJ,NDIM,DIST,IIDX,JIDX,NDIS)
+C
+C COMPUTES EUCLIDEAN DISTANCES BETWEEN EACH PAIR OF THE NOBJ POINTS
+C WHOSE CO-ORDINATES IN NDIM DIMENSIONS ARE IN X(NOBJ,NDIM).
+C
+C THE DISTANCES ARE STORED IN THE VECTOR DIST AND THE INDEX VECTORS
+C IIDX AND JIDX HOLD THE I,J INDICES OF THE POINTS IN THE ORDER IN
+C WHICH THE COMPARISONS ARE TO BE MADE.
+C
+C Written by Dr. Peter R. Minchin
+C Department of Biological Sciences
+C Southern Illinois University Edwardsville
+C PO Box 1651
+C Edwardsville, IL 62026-1651, U.S.A.
+C Phone: +1-618-650-2975 FAX: +1-618-650-3174
+C Email: pminchi at siue.edu
+C
+ INTEGER IIDX(NDIS), JIDX(NDIS)
+ DOUBLE PRECISION X(MAXOBJ,NDIM), DIST(NDIS)
+C
+C INITIALIZE DISTANCES
+C
+ DO I=1,NDIS
+ DIST(I)=0.0
+ ENDDO
+C
+C ACCUMULATE SUMS OF SQUARED DIFFERENCES ON EACH DIMENSION
+C
+ DO IDIM=1,NDIM
+ DO K=1,NDIS
+ DIST(K)=DIST(K)+(X(IIDX(K),IDIM)-X(JIDX(K),IDIM))**2
+ ENDDO
+ ENDDO
+C
+C TAKE SQUARE ROOTS OF TOTALS
+C
+ DO I=1,NDIS
+ DIST(I)=SQRT(DIST(I))
+ ENDDO
+ RETURN
+ END SUBROUTINE CLCDIS
+
+ SUBROUTINE CLCGRD (X,GRAD,MAXOBJ,NDIM,DIST,DHAT,IIDX,JIDX,
+ . NDIS,STRESS,SFACT,TFACT,ISFORM,DMEAN)
+C
+C ACCUMULATES THE NEGATIVE GRADIENT IN GRAD(NOBJ,NDIM).
+C
+C Written by Dr. Peter R. Minchin
+C Department of Biological Sciences
+C Southern Illinois University Edwardsville
+C PO Box 1651
+C Edwardsville, IL 62026-1651, U.S.A.
+C Phone: +1-618-650-2975 FAX: +1-618-650-3174
+C Email: pminchi at siue.edu
+C
+ INTEGER IIDX(NDIS), JIDX(NDIS)
+ DOUBLE PRECISION X(MAXOBJ,NDIM), GRAD(MAXOBJ,NDIM), DIST(NDIS),
+ . DHAT(NDIS), STRESS, SFACT, TFACT, DMEAN, SOTSQ, RECIPT,
+ . DELTA
+C
+ IF (STRESS.LE.0.0D0) RETURN
+ SOTSQ=SFACT/(TFACT**2)
+ RECIPT=1.0/TFACT
+ DO IDIM=1,NDIM
+ IF (ISFORM.LE.1) THEN
+C---Kruskal's stress formula 1
+ DO K=1,NDIS
+ IF (DIST(K).GT.0.0) THEN
+ DELTA=(SOTSQ-RECIPT*(DIST(K)-DHAT(K))/DIST(K))*
+ . (X(IIDX(K),IDIM)-X(JIDX(K),IDIM))
+ GRAD(IIDX(K),IDIM)=GRAD(IIDX(K),IDIM)+DELTA
+ GRAD(JIDX(K),IDIM)=GRAD(JIDX(K),IDIM)-DELTA
+ ENDIF
+ ENDDO
+ ELSE
+C---Kruskal's stress formula 2
+ DO K=1,NDIS
+ IF (DIST(K).GT.0.0) THEN
+ DELTA=(SOTSQ*(DIST(K)-DMEAN)/DIST(K)-
+ . RECIPT*(DIST(K)-DHAT(K))/DIST(K))*
+ . (X(IIDX(K),IDIM)-X(JIDX(K),IDIM))
+ GRAD(IIDX(K),IDIM)=GRAD(IIDX(K),IDIM)+DELTA
+ GRAD(JIDX(K),IDIM)=GRAD(JIDX(K),IDIM)-DELTA
+ ENDIF
+ ENDDO
+ ENDIF
+ ENDDO
+ RETURN
+ END SUBROUTINE CLCGRD
+
+ SUBROUTINE CLCSFA (GRAD,GRLAST,NOBJ,NDIM,MAXOBJ,SFGR,CAGRGL,
+ . SFGLST)
+C
+C COMPUTES SCALE FACTOR OF THE NEGATIVE GRADIENT IN GRAD(NOBJ,NDIM)
+C AND ALSO FINDS ITS DIRECTION COSINE WITH THE PREVIOUS GRADIENT,
+C STORED IN GRLAST(NOBJ,NDIM).
+C
+C Written by Dr. Peter R. Minchin
+C Department of Biological Sciences
+C Southern Illinois University Edwardsville
+C PO Box 1651
+C Edwardsville, IL 62026-1651, U.S.A.
+C Phone: +1-618-650-2975 FAX: +1-618-650-3174
+C Email: pminchi at siue.edu
+C
+ DOUBLE PRECISION GRAD(MAXOBJ,NDIM), GRLAST(MAXOBJ,NDIM),
+ . SFGR, CAGRGL, SFGLST, FN, FACTOR
+C
+ FN=DBLE(NOBJ)
+ SFGR=0.0D0
+ CAGRGL=0.0D0
+ DO IDIM=1,NDIM
+ DO IOBJ=1,NOBJ
+ SFGR=SFGR+GRAD(IOBJ,IDIM)**2
+ CAGRGL=CAGRGL+GRAD(IOBJ,IDIM)*GRLAST(IOBJ,IDIM)
+ ENDDO
+ ENDDO
+ SFGR=SQRT(SFGR/FN)
+ FACTOR=SFGR*SFGLST*FN
+ IF (FACTOR.NE.0.0) CAGRGL=CAGRGL/FACTOR
+ RETURN
+ END SUBROUTINE CLCSFA
+
+ SUBROUTINE CLCSTP (STEP,IT,SFGR,STRESS,COSAV,ACOSAV,SRATIO,
+ . SRATAV)
+C
+C COMPUTES NEW STEP SIZE. BASED ON J.B. KRUSKAL'S STEP SIZE ALGORITHM
+C IN THE BELL LABORATORIES PROGRAM "KYST".
+C
+C Written by Dr. Peter R. Minchin
+C Department of Biological Sciences
+C Southern Illinois University Edwardsville
+C PO Box 1651
+C Edwardsville, IL 62026-1651, U.S.A.
+C Phone: +1-618-650-2975 FAX: +1-618-650-3174
+C Email: pminchi at siue.edu
+C
+ DOUBLE PRECISION STEP, SFGR, STRESS, COSAV, ACOSAV, SRATIO,
+ . SRATAV, FACTR1, FACTR2, FACTR3
+C
+C MAKE A WILD GUESS FOR INITIAL STEP SIZE (IT=0).
+C
+ IF (IT.EQ.0) THEN
+ STEP=25.0*STRESS*SFGR
+C
+C OTHERWISE, UPDATE STEP SIZE.
+C
+ ELSE
+ FACTR1=4.0**COSAV
+ FACTR2=1.6/( (1.0+(MIN(1.0,SRATAV))**5) *
+ . (1.0+ACOSAV-ABS(COSAV)) )
+ FACTR3=SQRT(MIN(1.0,SRATIO))
+ STEP=STEP*FACTR1*FACTR2*FACTR3
+ ENDIF
+ RETURN
+ END SUBROUTINE CLCSTP
+
+ SUBROUTINE CLCSTR (DIST,DHAT,N,SNUMER,SDENOM,STRESS,ISFORM,DMEAN)
+C
+C COMPUTES KRUSKAL'S STRESS (FORMULA 1 or 2)
+C
+C Written by Dr. Peter R. Minchin
+C Department of Biological Sciences
+C Southern Illinois University Edwardsville
+C PO Box 1651
+C Edwardsville, IL 62026-1651, U.S.A.
+C Phone: +1-618-650-2975 FAX: +1-618-650-3174
+C Email: pminchi at siue.edu
+C
+ DOUBLE PRECISION DIST(N), DHAT(N), SNUMER, SDENOM, STRESS,
+ . DMEAN
+C
+ SNUMER=0.0D0
+ SDENOM=0.0D0
+ DMEAN=0.0D0
+ IF (ISFORM.GE.2) THEN
+ DO I=1,N
+ DMEAN=DMEAN+DIST(I)
+ ENDDO
+ DMEAN=DMEAN/DBLE(N)
+ DO I=1,N
+ SNUMER=SNUMER+(DIST(I)-DHAT(I))**2
+ SDENOM=SDENOM+(DIST(I)-DMEAN)**2
+ ENDDO
+ ELSE
+ DO I=1,N
+ SNUMER=SNUMER+(DIST(I)-DHAT(I))**2
+ SDENOM=SDENOM+DIST(I)**2
+ ENDDO
+ ENDIF
+ STRESS=SQRT(SNUMER/SDENOM)
+ RETURN
+ END SUBROUTINE CLCSTR
+
+ SUBROUTINE LINREG (DISS,DIST,DHAT,N,C)
+C
+C PERFORMS LINEAR REGRESSION OF DIST ON DISS, PLACING THE FITTED VALUES
+C IN DHAT. THE REGRESSION COEFFICIENTS ARE RETURNED IN C:
+C C(1) = INTERCEPT, C(2) = SLOPE.
+C
+C Written by Dr. Peter R. Minchin
+C Department of Biological Sciences
+C Southern Illinois University Edwardsville
+C PO Box 1651
+C Edwardsville, IL 62026-1651, U.S.A.
+C Phone: +1-618-650-2975 FAX: +1-618-650-3174
+C Email: pminchi at siue.edu
+C
+ DOUBLE PRECISION DISS(N), DIST(N), DHAT(N),
+ . AVDISS, SSDISS, AVDIST, SPROD, C(2), FN
+C
+C COMPUTE MEANS, SUMS OF SQUARES AND SUM OF PRODUCTS
+C
+ FN=DBLE(N)
+ AVDIST=0.0D0
+ AVDISS=0.0D0
+ DO K=1,N
+ AVDIST=AVDIST+DIST(K)
+ AVDISS=AVDISS+DISS(K)
+ ENDDO
+ AVDIST=AVDIST/FN
+ AVDISS=AVDISS/FN
+ SSDISS=0.0D0
+ SPROD=0.0D0
+ DO K=1,N
+ SSDISS=SSDISS+(DISS(K)-AVDISS)**2
+ SPROD=SPROD+(DIST(K)-AVDIST)*(DISS(K)-AVDISS)
+ ENDDO
+C
+C COMPUTE SLOPE AND INTERCEPT
+C
+ C(2)=SPROD/SSDISS
+ C(1)=AVDIST-C(2)*AVDISS
+C
+C CALCULATE FITTED VALUES
+C
+ DO K=1,N
+ DHAT(K)=C(1)+C(2)*DISS(K)
+ ENDDO
+ RETURN
+ END SUBROUTINE LINREG
+
+ SUBROUTINE MACOPY (A,MAXN1,N,M,B,MAXN2)
+C
+C COPY MATRIX A(N,M) TO B(N,M).
+C
+C CAN BE USED TO PACK AN ARRAY BY PASSING THE SAME ARRAY AS A AND B,
+C WITH MAXN2 SET TO N.
+C
+C Written by Dr. Peter R. Minchin
+C Department of Biological Sciences
+C Southern Illinois University Edwardsville
+C PO Box 1651
+C Edwardsville, IL 62026-1651, U.S.A.
+C Phone: +1-618-650-2975 FAX: +1-618-650-3174
+C Email: pminchi at siue.edu
+C
+ DOUBLE PRECISION A(MAXN1,M), B(MAXN2,M)
+C
+ DO J=1,M
+ DO I=1,N
+ B(I,J)=A(I,J)
+ ENDDO
+ ENDDO
+ RETURN
+ END SUBROUTINE MACOPY
+
+ SUBROUTINE MAINIT (X,M,N,MAXM,CONST)
+C
+C INITIALIZES MATRIX X(M,N) WITH THE VALUE CONST.
+C
+C Written by Dr. Peter R. Minchin
+C Department of Biological Sciences
+C Southern Illinois University Edwardsville
+C PO Box 1651
+C Edwardsville, IL 62026-1651, U.S.A.
+C Phone: +1-618-650-2975 FAX: +1-618-650-3174
+C Email: pminchi at siue.edu
+C
+ DOUBLE PRECISION X(MAXM,N), CONST
+C
+ DO J=1,N
+ DO I=1,M
+ X(I,J)=CONST
+ ENDDO
+ ENDDO
+ RETURN
+ END SUBROUTINE MAINIT
+
+ SUBROUTINE MAMAB (A,B,L,M,N,MAXL1,MAXM,C,MAXL2)
+C
+C MATRIX MULTIPLICATION C=A.B
+C
+C MATRIX A(L,M) is multiplied by B(M,N) to give C(L,N)
+C
+C Written by Dr. Peter R. Minchin
+C Department of Biological Sciences
+C Southern Illinois University Edwardsville
+C PO Box 1651
+C Edwardsville, IL 62026-1651, U.S.A.
+C Phone: +1-618-650-2975 FAX: +1-618-650-3174
+C Email: pminchi at siue.edu
+C
+ DOUBLE PRECISION A(MAXL1,M), B(MAXM,N), C(MAXL2,N)
+C
+ DO I=1,L
+ DO J=1,N
+ C(I,J)=0.0
+ DO K=1,M
+ C(I,J)=C(I,J)+A(I,K)*B(K,J)
+ ENDDO
+ ENDDO
+ ENDDO
+ RETURN
+ END SUBROUTINE MAMAB
+
+ SUBROUTINE MAMAS (A,MAXL,L,M,S)
+C
+C MULTIPLY MATRIX A(L,M) BY THE SCALAR S
+C
+C Written by Dr. Peter R. Minchin
+C Department of Biological Sciences
+C Southern Illinois University Edwardsville
+C PO Box 1651
+C Edwardsville, IL 62026-1651, U.S.A.
+C Phone: +1-618-650-2975 FAX: +1-618-650-3174
+C Email: pminchi at siue.edu
+C
+ DOUBLE PRECISION A(MAXL,M), S
+C
+ DO I=1,L
+ DO J=1,M
+ A(I,J)=A(I,J)*S
+ ENDDO
+ ENDDO
+ RETURN
+ END SUBROUTINE MAMAS
+
+ SUBROUTINE MAMATA (A,L,M,MAXL,B,MAXM)
+C
+C MATRIX MULTIPLICATION B=A'.A
+C
+C The transpose of matrix A(L,M) is multiplied A(L,M) to give B(M,M)
+C
+C Written by Dr. Peter R. Minchin
+C Department of Biological Sciences
+C Southern Illinois University Edwardsville
+C PO Box 1651
+C Edwardsville, IL 62026-1651, U.S.A.
+C Phone: +1-618-650-2975 FAX: +1-618-650-3174
+C Email: pminchi at siue.edu
+C
+ DOUBLE PRECISION A(MAXL,M), B(MAXM,M)
+C
+C COMPUTE B
+C
+ DO I=1,M
+ DO J=1,I
+ B(I,J)=0.0
+ DO K=1,L
+ B(I,J)=B(I,J)+A(K,I)*A(K,J)
+ ENDDO
+ B(J,I)=B(I,J)
+ ENDDO
+ ENDDO
+ RETURN
+ END SUBROUTINE MAMATA
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/vegan -r 1565
More information about the Vegan-commits
mailing list