/* Copyright (c) 2010, Chris Want 
 Please see BSD style licence at the end of this file */

abstract class Search {
  boolean enabled;

  abstract KrigData search(KrigData in, float x, float y, float z);
  abstract Search duplicate();
  
  void initialize(KrigData kd) {
    // Do nothing, let subclasses do things if they like
  }
  
  boolean isEnabled() {
    return enabled;
  }
  void enable() {
    enabled=true;
  }
  void disable() {
    enabled=false;
  }
  
}

/*
class RadiusSearch extends Search {
  RotMatrix rot;
  float radius, r2;

  RadiusSearch(float radius) {
    setRadius(radius);
  }

  RadiusSearch(float range1,float range2, float range3,
               float angle1, float angle2, float angle3) {
    setRange(range1, range2, range3, angle1, angle2, angle3);
    enabled=false;
  }
  void setRadius(float radius) {
    this.radius = radius;
    r2 = radius*radius;
  }

  void setRange(float range1,float range2, float range3,
               float angle1, float angle2, float angle3) {
      this.rot = new RotMatrix(range1, range2, range3, angle1, angle2, angle3);
  }
  
  KrigData search(KrigData in, float x, float y, float z) {
    if (in==null) return null;

    KrigData out = new KrigData();
    for (int i=0; i<in.ndata; ++i) {
      if (d2(in.x[i], in.y[i], in.z[i], x, y, z) <= r2) {
        out.addPoint(in.x[i], in.y[i], in.z[i], in.value[i]);
      }
    }
    if (out.ndata == 0) {
      return null;
    }
    return out;
  }
}

class RadiusNumSearch extends Search {
  RadiusSearch radiusSearch;
  NumSearch numSearch;
  
  RadiusNumSearch(float radius, int numdata) {
    radiusSearch = new RadiusSearch(radius);
    radiusSearch.disable();
    numSearch = new NumSearch(numdata);
    numSearch.disable();
    enabled=false;
  }
  
  KrigData search(KrigData in, float x, float y, float z) {
    KrigData out = in;
    if (radiusSearch.isEnabled()) {
      out = radiusSearch.search(out, x, y, z);
    }
    if (numSearch.isEnabled()) {
      out = numSearch.search(out, x, y, z);
    }
    return out;
  }
}
*/

class EllipsoidSearch extends Search {
  RotMatrix rot;
  int minPoints, maxPoints;
  
  EllipsoidSearch(RotMatrix rot, int minPoints, int maxPoints) {
    setValues(rot, minPoints, maxPoints);
  }
  EllipsoidSearch(Variogram vg, int minPoints, int maxPoints) {
    if (vg != null) {
      if (vg.structures != null) {
        if (vg.structures.length > 0) {
          setValues(vg.structures[vg.structures.length-1].rot, minPoints, maxPoints);
        }
      }
    }
  }
  void setValues(RotMatrix rot, int minPoints, int maxPoints) {
    this.rot = rot.deepCopy();
    this.minPoints = minPoints;
    this.maxPoints = maxPoints;
  }
  
  KrigData search(KrigData in, float x, float y, float z) {
    KrigData out;
    float[] dist = new float[in.ndata];
    float[] u, v;
    int swapped = 0, inRange = 0;
    float tmp;
    v = new float[] {x ,y, z};
    out = in.deepCopy();
    for (int i=0; i < out.ndata; ++i) {
      u = new float[] { out.x[i], out.y[i], out.z[i] };
      
      dist[i] = rot.dist(u, v)/rot.anis[0];
    }

    while (true) {
      swapped = 0;
      for (int i=0; i < out.ndata-1; ++i) {
        if (dist[i] > dist[i+1]) {
          tmp = dist[i];
          dist[i] = dist[i+1];
          dist[i+1] = tmp;
          out.swapPoints(i, i+1);
          ++swapped;
        }
      }
      if (swapped == 0) break;
    }
    for (int i=0; i < out.ndata-1; ++i) {
      if (dist[i] <= 1.0) {
        ++inRange;
      }
      else {
        break;
      }
    }

    if (inRange >= maxPoints) {
      return out.deepCopy(maxPoints);
    }
    if (inRange < minPoints) {
      // don't estimate, discard
      return null;
    }
    return out.deepCopy(inRange);
  }
  
  Search duplicate() {
    EllipsoidSearch dup;
    dup = new EllipsoidSearch(rot, minPoints, maxPoints);
    if (isEnabled()) {
      dup.enable();
    }
    else {
      dup.disable();
    }
    return dup;
  }
}

float d2(float x0, float y0, float z0, float x1, float y1, float z1) {
  float dx,dy,dz;
  dx = x1 - x0;
  dy = y1 - y0;
  dz = z1 - z0;
  return(dx*dx + dy*dy + dz*dz);
}

/* Copyright (c) 2010, Chris Want, Research Support Group,
 AICT, University of Alberta. All rights reserved.
 
 Redistribution and use in source and binary forms, with or without 
 modification, are permitted provided that the following conditions are met:
 
 1) Redistributions of source code must retain the above copyright 
 notice, this list of conditions and the following disclaimer.
 2) Redistributions in binary form must reproduce the above copyright 
 notice, this list of conditions and the following disclaimer in the 
 documentation and/or other materials provided with the distribution.
 
 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF  
 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF 
 THE POSSIBILITY OF SUCH DAMAGE.
 
 Contributors: Chris Want (University of Alberta),
 Jeff Boisvert (University of Alberta)
 
 */