//
// GriddedSet.java
//

/*
VisAD system for interactive analysis and visualization of numerical
data.  Copyright (C) 1996 - 2017 Bill Hibbard, Curtis Rueden, Tom
Rink, Dave Glowacki, Steve Emmerson, Tom Whittaker, Don Murray, and
Tommy Jasmin.

This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public
License as published by the Free Software Foundation; either
version 2 of the License, or (at your option) any later version.

This library 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
Library General Public License for more details.

You should have received a copy of the GNU Library General Public
License along with this library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA
*/

package visad;

/**
   GriddedSet is implemented by those Set sub-classes whose samples
   lie on a rectangular grid topology (but note the geometry need
   not be rectangular).  It is a M-dimensional array of points in
   R^N where ManifoldDimension = M <= N = DomainDimension.


   The order of the samples is the rasterization of the orders of
   the 1D components, with the first component increasing fastest.
   For more detail, see the example in Linear2DSet.java.


   Grid coordinates are zero-based and in the range -0.5 to length-0.5
   i.e., there are "length" intervals of length 1.0,
   the first centered on 0.0 and the last centered on length-1.0
   points outside this range are indicated by the grid coordinate
   Double.NaN.

*/
public class GriddedSet extends SampledSet implements GriddedSetIface {

  int[] Lengths;
  // error tolerance for Newton's method solvers
  // not static because may become data dependent
  float EPS = (float) 1.0E-15;
  // Pos records the grid's orientation
  // (i.e., the sign of the cross-products of the grid edges)
  boolean Pos;

  /** construct a GriddedSet with samples */
  public GriddedSet(MathType type, float[][] samples, int[] lengths)
         throws VisADException {
    this(type, samples, lengths, null, null, null, true);
  }

  /** construct a GriddedSet with samples and non-default CoordinateSystem */
  public GriddedSet(MathType type, float[][] samples, int[] lengths,
                    CoordinateSystem coord_sys, Unit[] units,
                    ErrorEstimate[] errors) throws VisADException {
    this(type, samples, lengths, coord_sys, units, errors, true);
  }

  public GriddedSet(MathType type, float[][] samples, int[] lengths,
             CoordinateSystem coord_sys, Unit[] units,
             ErrorEstimate[] errors, boolean copy)
             throws VisADException {
    super(type, lengths.length, coord_sys, units, errors);
    init_lengths(lengths);
    if (samples == null ) {
	setMySamples(null);
    }
    else {
      init_samples(samples, copy);
    }
  }

  private void init_lengths(int[] lengths) throws VisADException {
    Lengths = new int[ManifoldDimension];
    Length = 1;
    for (int j=0; j
   *                             lengths[i} contains the number of points
   *                             in the manifold for dimension i.
   * @throws VisADException      problem creating the set
   */
  public static GriddedSet create(MathType type, float[][] samples,
                                  int[] lengths) throws VisADException {
    return create(type, samples, lengths, null, null, null);
  }

  /**
   * General Factory method for creating the proper gridded set
   * (Gridded1DSet, Gridded2DSet, etc.).
   *
   * @param type                 MathType for the returned set
   * @param samples              Set samples
   * @param lengths              The dimensionality of the manifold.  
   *                             lengths[i} contains the number of points
   *                             in the manifold for dimension i.
   * @param coord_sys            CoordinateSystem for the GriddedSet
   * @param units                Unit-s of the values in samples
   * @param errors               ErrorEstimate-s for the values
   * @throws VisADException      problem creating the set
   */
  public static GriddedSet create(MathType type, float[][] samples,
                                  int[] lengths, CoordinateSystem coord_sys,
                                  Unit[] units, ErrorEstimate[] errors)
         throws VisADException {
    return create(type, samples, lengths, coord_sys, units, errors,
                  true, true);
  }

  /**
   * General Factory method for creating the proper gridded set
   * (Gridded1DSet, Gridded2DSet, etc.).
   *
   * @param type                 MathType for the returned set
   * @param samples              Set samples
   * @param lengths              The dimensionality of the manifold.  
   *                             lengths[i} contains the number of points
   *                             in the manifold for dimension i.
   * @param coord_sys            CoordinateSystem for the GriddedSet
   * @param units                Unit-s of the values in samples
   * @param errors               ErrorEstimate-s for the values
   * @param copy                 make a copy of the samples
   * @throws VisADException      problem creating the set
   */
  public static GriddedSet create(MathType type, float[][] samples,
                                  int[] lengths, CoordinateSystem coord_sys,
                                  Unit[] units, ErrorEstimate[] errors, 
                                  boolean copy)
         throws VisADException {
    return create(type, samples, lengths, coord_sys, units, errors,
                  copy, true);
  }

  /**
   * General Factory method for creating the proper gridded set
   * (Gridded1DSet, Gridded2DSet, etc.).
   *
   * @param type                 MathType for the returned set
   * @param samples              Set samples
   * @param lengths              The dimensionality of the manifold.  
   *                             lengths[i} contains the number of points
   *                             in the manifold for dimension i.
   * @param coord_sys            CoordinateSystem for the GriddedSet
   * @param units                Unit-s of the values in samples
   * @param errors               ErrorEstimate-s for the values
   * @param copy                 make a copy of the samples
   * @param test                 test to make sure samples are valid.  Used
   *                             for creating Gridded*DSets where the 
   *                             manifold dimension is equal to the domain
   *                             dimension
   * @throws VisADException      problem creating the set
   */
    public static GriddedSet create(MathType type, float[][] samples,
                                   int[] lengths, CoordinateSystem coord_sys,
                           Unit[] units, ErrorEstimate[] errors,
                           boolean copy, boolean test)
          throws VisADException {
    int domain_dimension = samples.length;
    int manifold_dimension = lengths.length;
    if (manifold_dimension > domain_dimension) {
      throw new SetException("GriddedSet.create: manifold_dimension " +
                             manifold_dimension + " is greater than" +
                             " domain_dimension " + domain_dimension);
    }

    switch (domain_dimension) {
      case 1:
        return new Gridded1DSet(type, samples,
                                lengths[0],
                                coord_sys, units, errors, copy);
      case 2:
        if (manifold_dimension == 1) {
          return new Gridded2DSet(type, samples,
                                  lengths[0],
                                  coord_sys, units, errors, copy);
        }
        else {
          return new Gridded2DSet(type, samples,
                                  lengths[0], lengths[1],
                                  coord_sys, units, errors, copy, test);
        }
      case 3:
        if (manifold_dimension == 1) {
          return new Gridded3DSet(type, samples,
                                  lengths[0],
                                  coord_sys, units, errors, copy);
        }
        else if (manifold_dimension == 2) {
          return new Gridded3DSet(type, samples,
                                  lengths[0], lengths[1],
                                  coord_sys, units, errors, copy);
        }
        else {
          return new Gridded3DSet(type, samples,
                                  lengths[0], lengths[1], lengths[2],
                                  coord_sys, units, errors, copy, test);
        }
      default:
        return new GriddedSet(type, samples,
                              lengths,
                              coord_sys, units, errors, copy);
    }
  }

  public Set makeSpatial(SetType type, float[][] samples) throws VisADException {
    return create(type, samples, Lengths, null, null, null, false, false);
  }

  public int getLength(int i) {
    return Lengths[i];
  }

  public int[] getLengths() {
    int[] lens = new int[Lengths.length];
    for (int i=0; i=0; i--) {
            wedge[k] = wedge[i] + base;
            k++;
          }
        }
        else {
          for (i=0; i=0; j--) {
        if (Double.isNaN(grid[j][i])) {
          k = -1;
          break;
        }
        k = ((int) (grid[j][i] + 0.5)) + Lengths[j] * k;
      }
      index[i] = k;
    }
    return index;
  }

  /** transform an array of non-integer grid coordinates to an array
      of values in R^DomainDimension */
  public float[][] gridToValue(float[][] grid) throws VisADException {
    if (Length > 1) {
      for (int j=0; j 1");
        }
      }
    }
    throw new UnimplementedException("GriddedSet.gridToValue");
  }

  /** Transform an array of values in R^DomainDimension to an array
      of non-integer grid coordinates. 
      A guess for for the first value may be supplied: useful when making
      repeated calls with one value on this set if the resulting grid coords are
      in proximity with one another. The last determined grid location becomes the guess.
  */   
  public float[][] valueToGrid(float[][] value, int[] guess) throws VisADException {
    if (Length > 1) {
      for (int j=0; j 1");
        }
      }
    }
    throw new UnimplementedException("GriddedSet.valueToGrid with supplied guess");
  }
  
  /** transform an array of values in R^DomainDimension to an array
      of non-integer grid coordinates */
  public float[][] valueToGrid(float[][] value) throws VisADException {
    if (Length > 1) {
      for (int j=0; j 1");
        }
      }
    }
    throw new UnimplementedException("GriddedSet.valueToGrid");      
  }
  
  /** for each of an array of values in R^DomainDimension, compute an array
      of 1-D indices and an array of weights, to be used for interpolation;
      indices[i] and weights[i] are null if i-th value is outside grid
      (i.e., if no interpolation is possible) 
      A guess for the first value may be supplied */
  
  public void valueToInterp(float[][] value, int[][] indices, float[][] weights)
              throws VisADException {
    valueToInterp(value, indices, weights, null);
  }
  
  public void valueToInterp(float[][] value, int[][] indices, float[][] weights, int[] guess)
              throws VisADException {
    if (value.length != DomainDimension) {
      throw new SetException("GriddedSet.valueToInterp: value dimension " +
                             value.length + " not equal to Domain dimension " +
                             DomainDimension);
    }
    int length = value[0].length; // number of values
    if (indices.length != length) {
      throw new SetException("GriddedSet.valueToInterp: indices length " +
                             indices.length +
                             " doesn't match value[0] length " +
                             value[0].length);
    }
    if (weights.length != length) {
      throw new SetException("GriddedSet.valueToInterp: weights length " +
                             weights.length +
                             " doesn't match value[0] length " +
                             value[0].length);
    }
    if (guess != null && guess.length != DomainDimension) {
      throw new SetException("GriddedSet.valueToInterp: guess length " +
                             guess.length +
                             " not equal to Domain dimension " +
                             DomainDimension);
    }
    
    float[][] grid;
    // convert value array to grid coord array
    if (guess != null) {
      grid = valueToGrid(value, guess);
    }
    else {
      grid = valueToGrid(value);
    }
    
    int i, j, k; // loop indices
    int lis; // temporary length of is & cs
    int length_is; // final length of is & cs, varies by i
    int isoff; // offset along one grid dimension
    float a, b; // weights along one grid dimension; a + b = 1.0
    int[] is; // array of indices, becomes part of indices
    float[] cs; // array of coefficients, become part of weights

    int base; // base index, as would be returned by valueToIndex
    int[] l = new int[ManifoldDimension]; // integer 'factors' of base
    float[] c = new float[ManifoldDimension]; // fractions with l; -0.5 <= c <= 0.5

    // array of index offsets by grid dimension
    int[] off = new int[ManifoldDimension];
    off[0] = 1;
    for (j=1; j= 0.0))) {
          // only interp along ManifoldDimension-1 if between two valid grid coords
          length_is *= 2;
        }
        base = l[ManifoldDimension-1];
        if (base >= Lengths[ManifoldDimension-1]) base = -1;
      }
      for (j=ManifoldDimension-2; j>=0 && base>=0; j--) {
        if (Double.isNaN(grid[j][i])) {
          base = -1;
        }
        else {
          l[j] = (int) (grid[j][i] + 0.5);
          if (l[j] == Lengths[j]) l[j]--; // WLH 23 Dec 99
          c[j] = grid[j][i] - ((float) l[j]);
          if (!((l[j] == 0 && c[j] <= 0.0) ||
                (l[j] == Lengths[j] - 1 && c[j] >= 0.0))) {
            // only interp along dimension j if between two valid grid coords
            length_is *= 2;
          }
          base = l[j] + Lengths[j] * base;
          if (l[j] < 0 || l[j] >= Lengths[j]) base = -1;
        }
      }

      if (base < 0) {
        // value is out of grid so return null
        is = null;
        cs = null;
      }
      else {
        // create is & cs of proper length, and init first element
        is = new int[length_is];
        cs = new float[length_is];
        is[0] = base;
        cs[0] = 1.0f;
        lis = 1;

        for (j=0; j= 0.0))) {
            // only interp along dimension j if between two valid grid coords
            if (c[j] >= 0.0) {
              // grid coord above base
              isoff = off[j];
              a = 1.0f - c[j];
              b = c[j];
            }
            else {
              // grid coord below base
              isoff = -off[j];
              a = 1.0f + c[j];
              b = -c[j];
           }
            // float is & cs; adjust new offsets; split weights
            for (k=0; kneighbors.length should be at least getLength().
   * On return, neighbors[i][j] will be the index of the j
   * -th neighboring point of point i.  This method
   * will allocate and set the array neighbors[i] for all
   * i.  For points in the interior of the set, the number of
   * neighboring points is equal to 2*getManifoldDimension().  The
   * number of neighboring points decreases, however, if the point in question
   * is an exterior point on a hyperface, hyperedge, or hypercorner.
   *
   * @param neighbors                The array to contain the indexes of the
   *                                 neighboring points.
   * @throws NullPointerException    if the array is null.
   * @throws ArrayIndexOutOfBoundsException
   *                                 if neighbors.length < getLength()
   *                                 .
   * @throws VisADException          if a VisAD failure occurs.
   */
  public void getNeighbors( int[][] neighbors )
              throws VisADException
  {
    int ii, ix, iy, iz, ii_R, ii_L, ii_U, ii_D, ii_F, ii_B;
    int a, b, c, d;
    int LengthX;
    int LengthY;
    int LengthZ;
    int LengthXY;
    int LengthXYZ;

    switch( ManifoldDimension )
    {
       case 1:
         neighbors[0] = new int[1];
         neighbors[Length - 1] = new int[1];
         neighbors[0][0] = 1;
         neighbors[Length - 1][0] = Length - 2;

         for ( ii = 1; ii < (Length - 1); ii++ ) {
           neighbors[ii] = new int[2];
           neighbors[ii][0] = ii - 1;
           neighbors[ii][1] = ii + 1;
         }
         break;
       case 2:

         LengthX = Lengths[0];
         LengthY = Lengths[1];

         //- corners:
         ii = 0;
         neighbors[ii] = new int[2];
         neighbors[ii][0] = ii + LengthX;
         neighbors[ii][1] = ii + 1;

         ii = Length - 1;
         neighbors[ii] = new int[2];
         neighbors[ii][0] = ii - LengthX;
         neighbors[ii][1] = ii - 1;

         ii = Length - LengthX;
         neighbors[ii] = new int[2];
         neighbors[ii][0] = ii - LengthX;
         neighbors[ii][1] = ii + 1;

         ii = LengthX - 1;
         neighbors[ii] = new int[2];
         neighbors[ii][0] = ii + LengthX;
         neighbors[ii][1] = ii - 1;

         //- edge points, not corners:
         for ( iy = 1; iy < (LengthY - 1); iy++ )
         {
           ii = iy*LengthX;
           neighbors[ii] = new int[3];
           neighbors[ii][0] = ii + LengthX;
           neighbors[ii][1] = ii - LengthX;
           neighbors[ii][2] = ii + 1;

           ii = iy*LengthX + LengthX - 1;
           neighbors[ii] = new int[3];
           neighbors[ii][0] = ii + LengthX;
           neighbors[ii][1] = ii - LengthX;
           neighbors[ii][2] = ii - 1;
         }

         for ( ix = 1; ix < (LengthX - 1); ix++ )
         {
           ii = ix;
           neighbors[ii] = new int[3];
           neighbors[ii][0] = ii - 1;
           neighbors[ii][1] = ii + 1;
           neighbors[ii][2] = ii + LengthX;

           ii = (LengthY-1)*LengthX + ix;
           neighbors[ii] = new int[3];
           neighbors[ii][0] = ii - 1;
           neighbors[ii][1] = ii + 1;
           neighbors[ii][2] = ii - LengthX;
         }
         //- interior points:
         for ( iy = 1; iy < (LengthY - 1); iy++) {
           for ( ix = 1; ix < (LengthX - 1); ix++) {

             ii   = iy*LengthX + ix;
             ii_R = ii + 1;
             ii_L = ii - 1;
             ii_U = ii + LengthX;
             ii_D = ii - LengthX;
             neighbors[ii] = new int[4];
             neighbors[ii][0] = ii_R;
             neighbors[ii][1] = ii_L;
             neighbors[ii][2] = ii_U;
             neighbors[ii][3] = ii_D;
           }
         }
         break;

       case 3:

         LengthX = Lengths[0];
         LengthY = Lengths[1];
         LengthZ = Lengths[2];
         LengthXY = LengthX*LengthY;
         LengthXYZ = LengthX*LengthY*LengthZ;

         //- corners:
         ii = 0;
         neighbors[ii] = new int[3];
         neighbors[ii][0] = ii + LengthX;
         neighbors[ii][1] = ii + 1;
         neighbors[ii][2] = ii + LengthXY;

         ii = LengthX*LengthY*( LengthZ - 1);
         neighbors[ii] = new int[3];
         neighbors[ii][0] = ii + 1;
         neighbors[ii][1] = ii + LengthX;
         neighbors[ii][2] = ii - LengthXY;

         ii = LengthX*LengthY - 1;
         neighbors[ii] = new int[3];
         neighbors[ii][0] = ii - LengthX;
         neighbors[ii][1] = ii - 1;
         neighbors[ii][2] = ii + LengthXY;

         ii = LengthX*LengthY*LengthZ - 1;
         neighbors[ii] = new int[3];
         neighbors[ii][0] = ii - 1;
         neighbors[ii][1] = ii - LengthX;
         neighbors[ii][2] = ii - LengthXY;

         ii = LengthX*(LengthY - 1);
         neighbors[ii] = new int[3];
         neighbors[ii][0] = ii - LengthX;
         neighbors[ii][1] = ii + 1;
         neighbors[ii][2] = ii + LengthXY;

         ii = LengthXY*(LengthZ - 1) + LengthX*(LengthY - 1);
         neighbors[ii] = new int[3];
         neighbors[ii][0] = ii + 1;
         neighbors[ii][1] = ii - LengthX;
         neighbors[ii][2] = ii - LengthXY;

         ii = LengthX - 1;
         neighbors[ii] = new int[3];
         neighbors[ii][0] = ii + LengthX;
         neighbors[ii][1] = ii - 1;
         neighbors[ii][2] = ii + LengthXY;

         ii = LengthXY*(LengthZ - 1) + (LengthX - 1);
         neighbors[ii] = new int[3];
         neighbors[ii][0] = ii - 1;
         neighbors[ii][1] = ii + LengthX;
         neighbors[ii][2] = ii - LengthXY;

         //- edges:
         for ( iy = 1; iy < (LengthY - 1); iy++ )
         {
           a = iy*LengthX;
           b = a + LengthX-1;

           ii = a;
           neighbors[ii] = new int[4];
           neighbors[ii][0] = ii + LengthX;
           neighbors[ii][1] = ii - LengthX;
           neighbors[ii][2] = ii + 1;
           neighbors[ii][3] = ii + LengthXY;

           ii = a + LengthX-1;
           neighbors[ii] = new int[4];
           neighbors[ii][0] = ii + LengthX;
           neighbors[ii][1] = ii - LengthX;
           neighbors[ii][2] = ii - 1;
           neighbors[ii][3] = ii + LengthXY;

           ii = a + (LengthZ-1)*LengthXY;
           neighbors[ii] = new int[4];
           neighbors[ii][0] = ii + LengthX;
           neighbors[ii][1] = ii - LengthX;
           neighbors[ii][2] = ii - 1;
           neighbors[ii][3] = ii - LengthXY;

           ii = b + (LengthZ-1)*LengthXY;
           neighbors[ii] = new int[4];
           neighbors[ii][0] = ii + LengthX;
           neighbors[ii][1] = ii - LengthX;
           neighbors[ii][2] = ii - 1;
           neighbors[ii][3] = ii - LengthXY;
         }
         for ( ix = 1; ix < (LengthX - 1); ix++ )
         {
           ii = ix;
           neighbors[ii] = new int[4];
           neighbors[ii][0] = ii - 1;
           neighbors[ii][1] = ii + 1;
           neighbors[ii][2] = ii + LengthX;
           neighbors[ii][3] = ii + LengthXY;

           ii = (LengthY-1)*LengthX + ix;
           neighbors[ii] = new int[4];
           neighbors[ii][0] = ii - 1;
           neighbors[ii][1] = ii + 1;
           neighbors[ii][2] = ii - LengthX;
           neighbors[ii][3] = ii + LengthXY;

           ii = (LengthZ-1)*LengthXY + ix;
           neighbors[ii] = new int[4];
           neighbors[ii][0] = ii - 1;
           neighbors[ii][1] = ii + 1;
           neighbors[ii][2] = ii + LengthX;
           neighbors[ii][3] = ii - LengthXY;

           ii = (LengthZ-1)*LengthXY + (LengthY-1)*LengthX + ix;
           neighbors[ii] = new int[4];
           neighbors[ii][0] = ii - 1;
           neighbors[ii][1] = ii + 1;
           neighbors[ii][2] = ii - LengthX;
           neighbors[ii][3] = ii - LengthXY;
         }

         for ( iz = 1; iz < (LengthZ - 1); iz++ )
         {
           a = iz*LengthXY;
           b = a + (LengthX-1);

           ii = a;
           neighbors[ii] = new int[4];
           neighbors[ii][0] = ii + 1;
           neighbors[ii][1] = ii + LengthX;
           neighbors[ii][2] = ii + LengthXY;
           neighbors[ii][3] = ii - LengthXY;

           ii = a + (LengthX-1);
           neighbors[ii] = new int[4];
           neighbors[ii][0] = ii - 1;
           neighbors[ii][1] = ii + LengthX;
           neighbors[ii][2] = ii + LengthXY;
           neighbors[ii][3] = ii + LengthXY;

           ii = a + LengthX*(LengthY-1);
           neighbors[ii] = new int[4];
           neighbors[ii][0] = ii + 1;
           neighbors[ii][1] = ii - LengthX;
           neighbors[ii][2] = ii + LengthXY;
           neighbors[ii][3] = ii - LengthXY;

           ii = a + (LengthXY-1);
           neighbors[ii] = new int[4];
           neighbors[ii][0] = ii - 1;
           neighbors[ii][1] = ii - LengthX;
           neighbors[ii][2] = ii + LengthXY;
           neighbors[ii][3] = ii - LengthXY;
         }

         //- sides:

         for ( iy = 1; iy < (LengthY - 1); iy++ )
         {
           for ( ix = 1; ix < ( LengthX - 1); ix++ )
           {
             ii = iy*LengthX + ix;
             neighbors[ii] = new int[5];
             neighbors[ii][0] = ii + 1;
             neighbors[ii][1] = ii - 1;
             neighbors[ii][2] = ii + LengthX;
             neighbors[ii][3] = ii - LengthX;
             neighbors[ii][4] = ii + LengthXY;

             ii = (LengthZ - 1)*LengthXY + iy*LengthX + ix;
             neighbors[ii] = new int[5];
             neighbors[ii][0] = ii + 1;
             neighbors[ii][1] = ii - 1;
             neighbors[ii][2] = ii + LengthX;
             neighbors[ii][3] = ii - LengthX;
             neighbors[ii][4] = ii - LengthXY;
           }
         }

         for ( iz = 1; iz < (LengthZ - 1); iz++ )
         {
           for ( ix = 1; ix < ( LengthX - 1); ix++ )
           {
             ii = iz*LengthXY + ix;
             neighbors[ii] = new int[5];
             neighbors[ii][0] = ii + 1;
             neighbors[ii][1] = ii - 1;
             neighbors[ii][2] = ii + LengthX;
             neighbors[ii][3] = ii + LengthXY;
             neighbors[ii][4] = ii - LengthXY;

             ii = iz*LengthXY + LengthX*(LengthY - 1) + ix;
             neighbors[ii] = new int[5];
             neighbors[ii][0] = ii + 1;
             neighbors[ii][1] = ii - 1;
             neighbors[ii][2] = ii - LengthX;
             neighbors[ii][3] = ii + LengthXY;
             neighbors[ii][4] = ii - LengthXY;
           }
         }

         for ( iz = 1; iz < (LengthZ - 1); iz++ )
         {
           for ( iy = 1; iy < ( LengthY - 1); iy++ )
           {
             ii = iz*LengthXY + iy*LengthX;
             neighbors[ii] = new int[5];
             neighbors[ii][0] = ii + 1;
             neighbors[ii][1] = ii + LengthX;
             neighbors[ii][2] = ii - LengthX;
             neighbors[ii][3] = ii + LengthXY;
             neighbors[ii][4] = ii - LengthXY;

             ii = iz*LengthXY + iy*LengthX + ( LengthX - 1);
             neighbors[ii] = new int[5];
             neighbors[ii][0] = ii - 1;
             neighbors[ii][1] = ii + LengthX;
             neighbors[ii][2] = ii - LengthX;
             neighbors[ii][3] = ii + LengthXY;
             neighbors[ii][4] = ii - LengthXY;
           }
         }

         //- interior points:
         for ( iz = 1; iz < (LengthZ - 1); iz++) {
           for ( iy = 1; iy < (LengthY - 1); iy++) {
             for ( ix = 1; ix < (LengthX - 1); ix++) {

               ii   = iz*LengthXY + iy*LengthX + ix;
               ii_R = ii + 1;
               ii_L = ii - 1;
               ii_F = ii + LengthX;
               ii_B = ii - LengthX;
               ii_U = ii + LengthXY;
               ii_D = ii - LengthXY;
               neighbors[ii] = new int[6];
               neighbors[ii][0] = ii_R;
               neighbors[ii][1] = ii_L;
               neighbors[ii][2] = ii_F;
               neighbors[ii][3] = ii_B;
               neighbors[ii][4] = ii_U;
               neighbors[ii][5] = ii_D;
             }
           }
         }
         break;

       default:
         throw new UnimplementedException("getNeighbors(): ManifoldDimension >"+
                                 ManifoldDimension+" not currently implemented" );

    }
  }

  /**
   * Returns the indexes of the neighboring points along a manifold
   * dimesion for every point in the set. Elements [i][0]
   * and [i][1] of the returned array are the indexes of the
   * neighboring sample points in the direction of decreasing and increasing
   * manifold index, respectively, for sample point i.  If a sample
   * point doesn't have a neighboring point (because it is an exterior point,
   * for example) then the value of the corresponding index will be -1.
   *
   * @param manifoldIndex          The index of the manifold dimension along
   *                               which to return the neighboring points.
   * @throws ArrayIndexOutOfBoundsException
   *                               if manifoldIndex < 0 || 
   *                               manifoldIndex >= getManifoldDimension()
   *                               .
   * @see #getManifoldDimension()
   */
  public int[][] getNeighbors( int manifoldIndex )
  {
    int[][] neighbors = new int[ Length ][2];
    int[] m_coords = new int[ ManifoldDimension ];
    int[][] indeces = new int[2][ ManifoldDimension ];
    int ii_tmp, idx_u, idx_d, mm, tt, ii, jj, kk, k;

    for ( ii = 0; ii < Length; ii++ )
    {
      //- get the manifold coordinates for each index  -*
      ii_tmp = ii;
      for ( jj = 0; jj < (ManifoldDimension-1); jj++ ) {
        m_coords[jj] = ii_tmp % Lengths[jj];
        ii_tmp /= Lengths[jj];
      }
      m_coords[ManifoldDimension-1] = ii_tmp;

      //- get coordinates of two neighbors on each side  -*
      for ( kk = 0; kk < 2; kk++) {
        for ( jj = 0; jj < ManifoldDimension; jj++) {
          indeces[kk][jj] = m_coords[jj];
        }
      }

      idx_u = m_coords[ manifoldIndex ] + 1;
      idx_d = m_coords[ manifoldIndex ] - 1;
      indeces[1][manifoldIndex] = idx_u;
      indeces[0][manifoldIndex] = idx_d;
      if ( idx_u >= Lengths[manifoldIndex] ) {
        indeces[1][manifoldIndex] = -1;
      }
      else if ( idx_d < 0 ) {
        indeces[0][manifoldIndex] = -1;
      }

      //- compute index for each new coordinates  -*
      for ( kk = 0; kk < 2; kk++ )
      {
        if ( indeces[kk][manifoldIndex] != -1 ) {
          ii_tmp = 0;
          for ( mm = (ManifoldDimension-1); mm>=0; mm--) {
            k = indeces[kk][mm];
            for ( tt = 0; tt < mm; tt++ ) {
              k = k*Lengths[tt];
            }
            ii_tmp += k;
          }
          neighbors[ii][kk] = ii_tmp;
        }
        else {
          neighbors[ii][kk] = -1;
        }
      }
    }
    return neighbors;
  }

  public boolean equals(Object set) {
    if (!(set instanceof GriddedSet) || set == null ||
        set instanceof LinearSet ||
        set instanceof Gridded1DDoubleSet) return false;
    if (this == set) return true;
    if (testNotEqualsCache((Set) set)) return false;
    if (testEqualsCache((Set) set)) return true;
    if (!equalUnitAndCS((Set) set)) return false;
    try {
      int i, j;
      if (DomainDimension != ((GriddedSet) set).getDimension() ||
          ManifoldDimension != ((GriddedSet) set).getManifoldDimension() ||
          Length != ((GriddedSet) set).getLength()) return false;
      for (j=0; j