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

class CDF {
  int numData;
  float datamin, datamax;
  float sortData[];
  CachedCDF cached = null;
  
  CDF(float[] data, int dim) {
    DataIterator idata;

    idata = new Iterator1d(data, dim);
    numData = dim;
    sortData = idata.getSortedData();
    //for (int i=0; i < numData; ++i) {
    //  println(sortData[i]);
    //}
    datamin = sortData[0];
    datamax = sortData[sortData.length-1];
  }
  
  CDF(float[] data) {
    this(data, data.length);
  }
  
  CDF(float[][][] data, int dim1, int dim2, int dim3) {
    DataIterator idata;

    idata = new Iterator3d(data, dim1, dim2, dim3);
    numData = dim1*dim2*dim3;
    sortData = idata.getSortedData();
    datamin = sortData[0];
    datamax = sortData[sortData.length-1];
  }

  CDF(ScalarField scalars, int arr) {
    DataIterator idata;

    idata = new IteratorScalarField(scalars, arr);
    numData = idata.getNumData();
    sortData = idata.getSortedData();
    datamin = sortData[0];
    datamax = sortData[sortData.length-1];
  }
  
  CDF(ScalarField scalars) {
    this(scalars, scalars.getActiveScalars());
  }
  
  float getProbability(float x) {
    int count = 0, below = 0, above = numData;
    
    if (x <= datamin) return 0.0;
    if (x >= datamax) return 1.0;
    
    if (cached == null) {
      cached = new CachedCDF();
    }
    else {
      if (x>=cached.value) {
        count = cached.index;
        below = cached.index;
      }
      else {
        count = 0;
        below = 0;
        above = min(cached.index+1, numData);;
      }
    }
    
    for (int i=below; i < above; ++i) {
      if (x >= sortData[i]) {
        ++count;
      }
      else {
        break;
      }
    }
    
    cached.value = x;
    cached.index = count;

    return ((float) count / (float) numData);
  }
}

class CachedCDF {
  int index;
  float value;
  
  CachedCDF() {
  }
}

abstract class DataIterator {

  abstract int getNumData();
  abstract void reset();
  abstract float getNextData();
  
  float[] getSortedData() {
    int numData;
    int numSwapped;
    
    numData = getNumData();
    float[] sorted = new float[numData];

    reset();
    for (int i=0; i < numData; ++i) {
      sorted[i] = getNextData();
    }

    while(true) {
      numSwapped = 0;
      for (int i=0; i < numData-1; ++i) {
        float tmp;
        if (sorted[i] > sorted[i+1]) {
          tmp = sorted[i+1];
          sorted[i+1] = sorted[i];
          sorted[i] = tmp;
          ++numSwapped;
        }
      }
      if (numSwapped == 0) break;
    }
    return sorted;
  }
}

class Iterator1d extends DataIterator {
  float[] data;
  int dim;
  int counter;
  
  Iterator1d(float[] data, int dim) {
    this.data = data;
    this.dim = dim;
  }
  
  void reset() {
    counter = 0;
  }
  
  float getNextData() {
    return data[counter++];
  }
  
  int getNumData() {
    return dim;
  }
}

class Iterator3d extends DataIterator {
  float[][][] data;
  int dim1, dim2, dim3;
  int counter;

  Iterator3d(float[][][] data, int dim1, int dim2, int dim3) {
    this.data = data;
    this.dim1 = dim1;
    this.dim2 = dim2;
    this.dim3 = dim3;
  }

  void reset() {
    counter = 0;
  }
  
  float getNextData() {
    return getData(counter++);
  }
  
  float getData(int n) {
    int i, j, k, tmp;
    i = n / (dim2*dim3);
    tmp = n - (i*dim2*dim3);
    j = tmp / dim3;
    k = tmp - j * dim3;
    
    return data[i][j][k];
  }
  
  int getNumData() {
    return dim1*dim2*dim3;
  }
}

class IteratorScalarField extends DataIterator {
  ScalarField sc;
  int arr;
  int last;
  int numData;
  int counter;
  
  IteratorScalarField(ScalarField sc, int arr) {
    this.sc = sc;
    this.arr = arr;
    this.numData = computeNumData();
  }

  void reset() {
    counter = 0;
  }
  
  int computeNumData() {
    int nData = 0;
    if (sc.hasInvalid) {
      for (int i=0; i < sc.xRes; ++i) {
        for (int j=0; j < sc.yRes; ++j) {
          for (int k=0; k < sc.zRes; ++k) {
            if (!sc.isInvalid(i, j, k)) {
              ++nData;
            }
          }
        }
      }
      return nData;
    }
    return sc.xRes*sc.yRes*sc.zRes;
  }

  int getNumData() {
    return numData;
  }
  
  float getNextData() {
    int i, j, k, tmp, maxData;
    float value;
    maxData = sc.xRes*sc.yRes*sc.zRes;
    
    while(true) {
      i = counter / (sc.yRes*sc.zRes);
      tmp = counter - (i*sc.yRes*sc.zRes);
      j = tmp / sc.zRes;
      k = tmp - j * sc.zRes;
      ++counter;
      if (!sc.isInvalid(i, j, k)) {
        return sc.getData(i, j, k);
      }
      if (counter >= maxData) break;
    }
    return 0.0;
  }
}

/* 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)
*/