//  $Id: Octree.java,v 1.1 2009/05/08 20:40:58 agramada Exp $
//
//  Copyright 2000-2004 The Regents of the University of California.
//  All Rights Reserved.
//
//  Permission to use, copy, modify and distribute any part of this
//  Molecular Biology Toolkit (MBT)
//  for educational, research and non-profit purposes, without fee, and without
//  a written agreement is hereby granted, provided that the above copyright
//  notice, this paragraph and the following three paragraphs appear in all
//  copies.
//
//  Those desiring to incorporate this MBT into commercial products
//  or use for commercial purposes should contact the Technology Transfer &
//  Intellectual Property Services, University of California, San Diego, 9500
//  Gilman Drive, Mail Code 0910, La Jolla, CA 92093-0910, Ph: (858) 534-5815,
//  FAX: (858) 534-7345, E-MAIL:invent@ucsd.edu.
//
//  IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY FOR
//  DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING
//  LOST PROFITS, ARISING OUT OF THE USE OF THIS MBT, EVEN IF THE
//  UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
//  THE MBT PROVIDED HEREIN IS ON AN "AS IS" BASIS, AND THE
//  UNIVERSITY OF CALIFORNIA HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT,
//  UPDATES, ENHANCEMENTS, OR MODIFICATIONS. THE UNIVERSITY OF CALIFORNIA MAKES
//  NO REPRESENTATIONS AND EXTENDS NO WARRANTIES OF ANY KIND, EITHER IMPLIED OR
//  EXPRESS, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
//  MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, OR THAT THE USE OF THE
//  MBT WILL NOT INFRINGE ANY PATENT, TRADEMARK OR OTHER RIGHTS.
//
//  For further information, please see:  http://mbt.sdsc.edu
//
//  History:
//  $Log: Octree.java,v $
//  Revision 1.12 2024/11/27 cdedman-rollet
//  This is a modified version of the original Octree class which removed some 
//  unused method for our project purposes.
//
//  Revision 1.1  2009/05/08 20:40:58  agramada
//  This is a modified version of the MBT toolkit that I need for charge profile comparison.
//
//  Revision 1.10  2004/04/09 00:15:21  moreland
//  Updated copyright to new UCSD wording.
//
//  Revision 1.9  2004/01/29 17:14:52  agramada
//  Removed General Atomics from copyright
//
//  Revision 1.8  2004/01/14 19:25:51  agramada
//  Remove useless code, commented public methods. Also, the build methods now
//  throw an ExcessiveDivisionException when the size of the box becomes too
//  "small".
//
//  Revision 1.7  2003/12/15 21:34:42  moreland
//  Corrected open comment.
//
//  Revision 1.6  2003/12/15 21:33:27  moreland
//  Commented out debug print statments.
//
//  Revision 1.5  2003/10/17 18:19:32  moreland
//  Fixed a javadoc comment.
//
//  Revision 1.4  2003/10/01 21:11:51  agramada
//  Added methods needed by DerivedInformation class for derivation of
//  secondary structures.
//
//  Revision 1.3  2003/07/17 16:40:32  agramada
//  Removed some comments.
//
//  Revision 1.2  2003/07/11 20:23:05  agramada
//  Made get*Bonds methods return Bond objects instead of BondInfo. Removed
//  some coments.
//
//  Revision 1.1  2003/07/11 18:17:53  moreland
//  Modifed Apostol's Octree classes to genate Bonds from the BondFactory
//  and in turn the StructureMap class.
//
//  Revision 1.3  2003/07/07 23:02:47  agramada
//  Added a method to return a vector of covalent bonds rather than an array.
//
//  Revision 1.2  2003/06/24 22:19:46  agramada
//  Reorganized the geometry package. Old classes removed, new classes added.
//
//  Revision 1.1  2003/04/24 18:46:36  agramada
//  First version of the Octree class to be used for an efficient search of
//  the set of bonds between a set of atoms.
//
//  Revision 1.2  2003/02/20 17:08:22  agramada
//  First working version.
//
//  Revision 1.0  2003/02/20 23:38:39  agramada
//
package src;

import java.util.*;

/**
 * Octree class constructed with recursion in mind.
 * Each child is itself a tree.
 * <P>
 * @author      Apostol Gramada
 */
public class Octree {
    // Fields
    private OctreeDataItem[] dataItems = null;
    private Octree[] children = null;
    private Octree parent = null;
    private Octree root = null;
    private int dimension; // Keep it general as long as practical. In practice, maximum 3.
    private int maxNumberOfChildren;
    private int numberOfChildren;
    private int childId; // Number from 0 to 2^dimension
    private int weight;
    private int leafCutOff = 1;
    private int majorAxis; // Axis with the maximum spatial extension
    private int minorAxis; // Axis with the minimal spatial extension
    private int type1Count, type2Count;
    private double cellRadius;
    private double[] firstCorner;
    private double[] secondCorner;
    private double[] geometricalCenter;
    private double[] margin;
    private double[] cellCenter;
    private double[] size;
    private double[] mid; // Center of the cell
    private String path; // Sequence of digits showing the path from root to this child
    private boolean leaf = false;
    private Hashtable<OctreeDataItem, Octree> leafHash = null;
    private static int countAppBondOp = 0;
    private static int rejectedPaths = 0;
    private static double searchCutOff;

    /**
     * A counter for the number of children instantiated in the tree
     */
    public static int countChildren = 0;

    // Constructor for child nodes
    public Octree(Octree root, Octree parent, int id, OctreeDataItem[] data,
            double[] firstCorner, double[] secondCorner) {
        this.root = root;
        this.childId = id;
        this.path = parent.path + id;
        this.dimension = root.getDimension(); // Assume dimension passed is the same as dimension of data
        // this.dataItems = data;
        this.firstCorner = firstCorner;
        this.secondCorner = secondCorner;
        this.parent = parent;
        this.maxNumberOfChildren = root.getMaxNumberOfChildren();

        // Set the size in each direction.
        double maxSize = 0.0;
        double minSize = 1.0E6;

        size = new double[dimension];
        mid = new double[dimension];

        cellRadius = 0.0;

        for (int j = 0; j < dimension; j++) {
            size[j] = Math.abs(secondCorner[j] - firstCorner[j]);
            mid[j] = (secondCorner[j] + firstCorner[j]) / 2.0;

            cellRadius += size[j] * size[j];

            if (size[j] > maxSize) {
                maxSize = size[j];
                majorAxis = j;
            }

            if (size[j] < minSize) {
                minSize = size[j];
                minorAxis = j;
            }
        }

        cellRadius = Math.sqrt(cellRadius);
        cellRadius /= 2.0;
    }

    /**
     * Constructor that is mostly used for instantiating the root tree.
     */
    public Octree(int spaceDimension, OctreeDataItem[] data, double[] offset) {
        // Derive the domain spaned by the coordinates of the elements.
        initialize(spaceDimension, data, offset);
    }

    // Constructor for root node
    public Octree(int spaceDimension, double[][] coords1, double[][] coords2, double[] offset) {
        type1Count = coords1.length;
        type2Count = coords2.length;

        int treeDataLength = coords1.length + coords2.length;

        OctreeCoordinateItem[] treeData = new OctreeCoordinateItem[treeDataLength];

        for (int i = 0; i < coords1.length; i++) {
            treeData[i] = new OctreeCoordinateItem(coords1[i], i, 1);
        }

        for (int i = coords1.length; i < treeDataLength; i++) {
            treeData[i] = new OctreeCoordinateItem(coords2[i - coords1.length], i - coords1.length, 2);
        }

        initialize(spaceDimension, treeData, offset);
    }

    // Initializes the root node
    public void initialize(int spaceDimension, OctreeDataItem[] data, double[] offset) {
        // Derive the domain spaned by the coordinates of the elements.
        dimension = spaceDimension; // Assume dimension passed is the same as dimension of data
        dataItems = data;

        weight = data.length;
        margin = new double[dimension];
        // leafHash = new Hashtable();

        // System.out.println( "Root tree data length: " + dataItems.length );
        firstCorner = new double[dimension];
        secondCorner = new double[dimension];

        for (int i = 0; i < dimension; i++) {
            margin[i] = offset[i];
            firstCorner[i] = 1.0E6;
            secondCorner[i] = -1.0E6;
        }

        maxNumberOfChildren = 1;
        for (int i = 0; i < dimension; i++) {
            maxNumberOfChildren <<= 1;
        }

        // Determine the two corners from data coordinates
        double x;
        for (int i = 0; i < dataItems.length; i++) {

            double coordinate[] = dataItems[i].getCoordinate();

            for (int j = 0; j < dimension; j++) {
                x = coordinate[j];

                if (x <= firstCorner[j]) {
                    firstCorner[j] = x;
                }

                if (x >= secondCorner[j]) {
                    secondCorner[j] = x;
                }
            }
        }

        // Add margins to the corners to make the box somewhat larger than
        // the container of spatial data and set the size in each direction.
        double maxSize = 0.0;
        double minSize = 1.0E6;

        size = new double[dimension];
        mid = new double[dimension];
        cellRadius = 0.0;

        for (int j = 0; j < dimension; j++) {
            firstCorner[j] -= margin[j];
            secondCorner[j] += margin[j];

            size[j] = Math.abs(secondCorner[j] - firstCorner[j]);
            mid[j] = (secondCorner[j] + firstCorner[j]) / 2.0;

            cellRadius += size[j] * size[j];

            if (size[j] > maxSize) {
                maxSize = size[j];
                majorAxis = j;
            }

            if (size[j] < minSize) {
                minSize = size[j];
                minorAxis = j;
            }
        }

        cellRadius = Math.sqrt(cellRadius);
        cellRadius /= 2.0;

        // We have a root of the octree now
        childId = 0;
        path = "" + childId;
        // levels = 0;
        setRoot(this);
        // if ( data.length > 0 ) build( );
    }

    public void build() throws ExcessiveDivisionException {
        countChildren++;

        children = new Octree[maxNumberOfChildren];

        Vector[] dataSets = new Vector[maxNumberOfChildren]; // Collects data for each child

        // double[] mid = new double[ dimension ];
        double[] childSize = new double[dimension];
        double[] corner2;
        int setFirstBit = 1;

        // Determine the center of the parallelepiped.
        for (int i = 0; i < dimension; i++) {
            childSize[i] = size[i] / 2.0;
            // mid[i] = ( firstCorner[i] + secondCorner[i] ) / 2;
        }

        // Split the data according to the cell.
        for (int i = 0; i < maxNumberOfChildren; i++) {
            dataSets[i] = new Vector();
        }

        for (int i = 0; i < dataItems.length; i++) {
            int child = 0;
            double coordinate[] = dataItems[i].getCoordinate();

            for (int j = 0; j < dimension; j++) {
                child <<= 1;

                if (coordinate[j] <= mid[j]) {
                    child |= setFirstBit;
                }
            }

            dataSets[child].add(dataItems[i]);
        }

        OctreeDataItem[] tmpData = null;
        Iterator dataIterator = null;

        double[] secondCorner = null;
        int dataSize = 0;
        int l = 0;
        int testBit;

        numberOfChildren = 0;

        if (dataItems.length > leafCutOff) {
            for (int i = 0; i < maxNumberOfChildren; i++) {
                dataSize = dataSets[i].size();

                if (dataSize > 0) {
                    l = 0;
                    tmpData = new OctreeDataItem[dataSize];
                    dataIterator = dataSets[i].iterator();

                    while (dataIterator.hasNext()) {
                        tmpData[l] = (OctreeDataItem) dataIterator.next();
                        l++;
                    }

                    // Calculate the second corner of the child
                    corner2 = new double[dimension];
                    testBit = 1;

                    for (int j = dimension - 1; j >= 0; j--) {
                        if ((i & testBit) > 0) {
                            corner2[j] = mid[j] - childSize[j];
                        }

                        else {
                            corner2[j] = mid[j] + childSize[j];
                        }

                        testBit <<= 1;
                    }

                    children[i] = new Octree(this.root, this, i, null, mid, corner2);
                    children[i].build(tmpData);
                }

                numberOfChildren++;
            }
        }

        else {
            setLeafFlag(true);
            /*
             * for (int j = 0; j < dataItems.length; j++) {
             * leafHash.put((OctreeDataItem) dataItems[j], this);
             * }
             */
            // System.out.println( "Warning: Root is also leaf " );
            // TODO: Perhaps this should throw an exception?
        }
    }

    // Method to build children recursively
    public void build(OctreeDataItem[] data) throws ExcessiveDivisionException {
        countChildren++;
        weight = data.length;

        if (data.length <= leafCutOff) {
            setData(data);
            /*
             * for (int j = 0; j < data.length; j++) {
             * if (leafHash.get(data[j]) == null) {
             * leafHash.put((OctreeDataItem) data[j], this);
             * }
             * }
             */
            leaf = true;
            return;
        }

        children = new Octree[maxNumberOfChildren];

        Vector[] dataSets = new Vector[maxNumberOfChildren]; // Collects data for each child

        // double[] mid = new double[ dimension ];
        double[] childSize = new double[dimension];
        double[] corner2;
        int setFirstBit = 1;

        // Determine the center of the parallelepiped.
        for (int i = 0; i < dimension; i++) {
            childSize[i] = size[i] / 2.0;
            // mid[i] = ( firstCorner[i] + secondCorner[i] ) / 2;
        }

        boolean throwException = true;
        for (int i = 0; i < dimension; i++) {
            if (childSize[i] > 0.001) {
                throwException = false;
                break;
            }
        }

        if (throwException) {

            for (int i = 0; i < data.length; i++) {
                double[] coord = data[i].getCoordinate();
                System.err.println("Point " + i + "  " + coord[0] + "  "
                        + coord[1] + "  " + coord[2] + " Index "
                        + data[i].getIndex() + "  Type " + data[i].getType());
            }

            setData(data);
            /*
             * for (int j = 0; j < data.length; j++) {
             * if (leafHash.get(data[j]) == null) {
             * leafHash.put((OctreeDataItem) data[j], this);
             * }
             * }
             */
            leaf = true;
            return;

            // throw new ExcessiveDivisionException(
            // "Excessive division of the data bounding box");
        }

        // Split the data according to the cell.
        for (int i = 0; i < maxNumberOfChildren; i++) {
            dataSets[i] = new Vector();
        }

        for (int i = 0; i < data.length; i++) {
            int child = 0;
            double coordinate[] = data[i].getCoordinate();

            for (int j = 0; j < dimension; j++) {
                child <<= 1;

                if (coordinate[j] <= mid[j]) {
                    child |= setFirstBit;
                }
            }

            dataSets[child].add(data[i]);
        }

        OctreeDataItem[] tmpData = null;
        Iterator dataIterator = null;

        double[] secondCorner = null;
        int dataSize = 0;
        int l = 0;
        int testBit;

        numberOfChildren = 0;

        /*
         * for( int i=0; i<maxNumberOfChildren; i++ )
         * {
         * System.out.println( "Data Sizes: " + i + "  " + dataSets[i].size() );
         * }
         */

        for (int i = 0; i < maxNumberOfChildren; i++) {
            dataSize = dataSets[i].size();

            if (dataSize > 0) {
                l = 0;
                tmpData = new OctreeDataItem[dataSize];
                dataIterator = dataSets[i].iterator();

                while (dataIterator.hasNext()) {
                    tmpData[l] = (OctreeDataItem) dataIterator.next();
                    l++;
                }

                // Calculate the second corner of the child.
                corner2 = new double[dimension];
                testBit = 1;

                for (int j = dimension - 1; j >= 0; j--) {

                    if ((i & testBit) > 0) {
                        corner2[j] = mid[j] - childSize[j];
                    }

                    else {
                        corner2[j] = mid[j] + childSize[j];
                    }

                    testBit <<= 1;
                }

                children[i] = new Octree(this.root, this, i, null, mid, corner2);
                children[i].build(tmpData);
                numberOfChildren++;
            }
        }
    }

    // Get partition based on nearest items
    public int[][] getPartition() {
        int[][] part = new int[type1Count][2];

        OctreeDataItem initialOutpt = null;

        for (int i = 0; i < type1Count; i++) {
            Octree.searchCutOff = 2 * cellRadius;
            part[i][0] = i;
            part[i][1] = ((OctreeCoordinateItem) nearest(dataItems[i].getCoordinate(), initialOutpt, 2)).getIndex();
        }

        return part;
    }

    // Find nearest item in the octree
    public Object nearest(double[] point, Object initialOutpt, int type) {
        Object outpt = initialOutpt;
        double dist = getDistance(point, mid) - cellRadius;
        double[] tmp = null;
        Object tmpOutpt = null;

        if (dist > Octree.searchCutOff) {
            return outpt;
        }

        else {
            // System.err.println("Number of dataitems " + dataItems.length);
            if (leaf) {
                for (int i = 0; i < dataItems.length; i++) {
                    tmp = dataItems[i].getCoordinate();
                    dist = getDistance(point, dataItems[i].getCoordinate());

                    if ((dist <= Octree.searchCutOff) && (type == dataItems[i].getType())) {
                        Octree.searchCutOff = dist;
                        outpt = dataItems[i];
                    }
                }

                return outpt;
            }

            else {
                for (int i = 0; i < maxNumberOfChildren; i++) {
                    if (children[i] != null) {
                        tmpOutpt = children[i].nearest(point, outpt, type);

                        if (tmpOutpt != null) {
                            outpt = tmpOutpt;
                        }

                    }
                }

                return outpt;
            }

            // System.out.println(
            // "This doesn't seem to be a root tree, No data found " );
            // TODO: Perhaps this should throw an exception?
        }
    }

    // Calculate distance between two points
    private final double getDistance(double[] coord1, double[] coord2) {
        double distance = 0.0;

        for (int i = 0; i < coord1.length; i++) {
            distance += (coord2[i] - coord1[i]) * (coord2[i] - coord1[i]);
        }

        distance = Math.sqrt(distance);

        return distance;
    }

    // Getters
    public int getDimension() {
        return dimension;
    }

    public int getMaxNumberOfChildren() {
        return maxNumberOfChildren;
    }

    public void setRoot(Octree rootTree) {
        this.root = rootTree;
    }

    public void setData(OctreeDataItem[] data) {
        this.dataItems = data;
    }

    private void setLeafFlag(boolean leafFlag) {
        this.leaf = leafFlag;
    }

}
