/*
 *                    BioJava development code
 *
 * This code may be freely distributed and modified under the
 * terms of the GNU Lesser General Public Licence.  This should
 * be distributed with the code.  If you do not have a copy,
 * see:
 *
 *      http://www.gnu.org/copyleft/lesser.html
 *
 * Copyright for this code is held jointly by the individual
 * authors.  These should be listed in @author doc comments.
 *
 * For more information on the BioJava project and its aims,
 * or to join the biojava-l mailing list, visit the home page
 * at:
 *
 *      http://www.biojava.org/
 *
 */
package org.biojava.nbio.structure.contact;


import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collection;
import java.util.List;

import javax.vecmath.Point3d;

import org.biojava.nbio.structure.Atom;
import org.biojava.nbio.structure.Calc;


/**
 * A grid to be used for calculating atom contacts through a spatial hashing algorithm.
 * <p>
 * The grid is composed of cells of size of the cutoff so that the distances that need to be calculated
 * are reduced to those within each cell and to the neighbouring cells.
 * <p>
 * Usage, for generic 3D points:
 * <pre>
 *  Point3d[] points = ...;
 *  Grid grid = new Grid(8.0);
 *  grid.addCoords(points);
 *  List&lt;Contact&gt; contacts = getIndicesContacts();
 * </pre>
 * Usage, for atoms:
 * <pre>
 *  Atom[] atoms = ...;
 *  Grid grid = new Grid(8.0);
 *  grid.addCoords(atoms);
 *  AtomContactSet contacts = getAtomContacts();
 * </pre> 
 *
 * @author Jose Duarte
 *
 */
public class Grid {

        /**
         * The scale: we use units of hundredths of Angstroms (thus cutoffs can be specified with a maximum precision of 0.01A)
         */
        private static final int SCALE=100;

        private GridCell[][][] cells;

        private double cutoff;
        private int cellSize;

        private Point3d[] iAtoms;
        private Point3d[] jAtoms;
        
        private Atom[] iAtomObjects;
        private Atom[] jAtomObjects;

        // the bounds in int grid coordinates
        private int[] bounds;

        // the i and j bounding boxes in original double coordinates
        private BoundingBox ibounds;
        private BoundingBox jbounds;

        private boolean noOverlap; // if the 2 sets of atoms are found not to overlap then this is set to true

        /**
         * Creates a <code>Grid</code>, the cutoff is in the same units as the coordinates
         * (Angstroms if they are atom coordinates) and can
         * be specified to a precision of 0.01.
         * @param cutoff
         */
        public Grid(double cutoff) {
                this.cutoff = cutoff;
                this.cellSize = (int) Math.floor(cutoff*SCALE);
                this.noOverlap = false;
        }

        private int getFloor(double number) {
                return (cellSize*((int)Math.floor(number*SCALE/cellSize)));
        }

        private int xintgrid2xgridindex(int xgridDim) {
                return (xgridDim-bounds[0])/cellSize;
        }

        private int yintgrid2ygridindex(int ygridDim) {
                return (ygridDim-bounds[1])/cellSize;
        }

        private int zintgrid2zgridindex(int zgridDim) {
                return (zgridDim-bounds[2])/cellSize;
        }

        /**
         * Adds the i and j atoms and fills the grid. Their bounds will be computed.
         * Subsequent call to {@link #getIndicesContacts()} or {@link #getAtomContacts()} will produce the interatomic contacts.
         * @param iAtoms
         * @param jAtoms
         */
        public void addAtoms(Atom[] iAtoms, Atom[] jAtoms) {
                addAtoms(iAtoms, null, jAtoms, null);
        }

        /**
         * Adds the i and j atoms and fills the grid, passing their bounds (array of size 6 with x,y,z minima and x,y,z maxima)
         * This way the bounds don't need to be recomputed.
         * Subsequent call to {@link #getIndicesContacts()} or {@link #getAtomContacts()} will produce the interatomic contacts.
         * @param iAtoms
         * @param icoordbounds
         * @param jAtoms
         * @param jcoordbounds
         */
        public void addAtoms(Atom[] iAtoms, BoundingBox icoordbounds, Atom[] jAtoms, BoundingBox jcoordbounds) {
                this.iAtoms = Calc.atomsToPoints(iAtoms);
                this.iAtomObjects = iAtoms;

                if (icoordbounds!=null) {
                        this.ibounds = icoordbounds;
                } else {
                        this.ibounds = new BoundingBox(this.iAtoms);
                }

                this.jAtoms = Calc.atomsToPoints(jAtoms);
                this.jAtomObjects = jAtoms;

                if (jAtoms==iAtoms) {
                        this.jbounds=ibounds;
                } else {
                        if (jcoordbounds!=null) {
                                this.jbounds = jcoordbounds;
                        } else {
                                this.jbounds = new BoundingBox(this.jAtoms);

                        }
                }

                fillGrid();
        }

        /**
         * Adds a set of atoms, subsequent call to {@link #getIndicesContacts()} or {@link #getAtomContacts()} will produce the interatomic contacts.
         * The bounding box of the atoms will be computed based on input array
         * @param atoms
         */
        public void addAtoms(Atom[] atoms) {
                addAtoms(atoms, (BoundingBox) null);
        }

        /**
         * Adds a set of atoms, subsequent call to {@link #getIndicesContacts()} or {@link #getAtomContacts()} will produce the interatomic contacts.
         * The bounds calculated elsewhere can be passed, or if null they are computed.
         * @param atoms
         * @param bounds
         */
        public void addAtoms(Atom[] atoms, BoundingBox bounds) {
                this.iAtoms = Calc.atomsToPoints(atoms);
                this.iAtomObjects = atoms;

                if (bounds!=null) {
                        this.ibounds = bounds;
                } else {
                        this.ibounds = new BoundingBox(iAtoms);
                }

                this.jAtoms = null;
                this.jAtomObjects = null;
                this.jbounds = null;

                fillGrid();
        }
        
        /**
         * Adds the i and j coordinates and fills the grid. Their bounds will be computed.
         * Subsequent call to {@link #getIndicesContacts()} will produce the 
         * contacts, i.e. the set of points within distance cutoff.
         * 
         * Subsequent calls to method {@link #getAtomContacts()} will produce a NullPointerException 
         * since this only adds coordinates and no atom information.
         * @param iAtoms
         * @param jAtoms
         */
        public void addCoords(Point3d[] iAtoms, Point3d[] jAtoms) {
                addCoords(iAtoms, null, jAtoms, null);
        }

        /**
         * Adds the i and j coordinates and fills the grid, passing their bounds (array of size 6 with x,y,z minima and x,y,z maxima)
         * This way the bounds don't need to be recomputed.
         * Subsequent call to {@link #getIndicesContacts()} will produce the 
         * contacts, i.e. the set of points within distance cutoff.
         * 
         * Subsequent calls to method {@link #getAtomContacts()} will produce a NullPointerException 
         * since this only adds coordinates and no atom information.
         * @param iAtoms
         * @param icoordbounds
         * @param jAtoms
         * @param jcoordbounds
         */
        public void addCoords(Point3d[] iAtoms, BoundingBox icoordbounds, Point3d[] jAtoms, BoundingBox jcoordbounds) {
                this.iAtoms = iAtoms;
                this.iAtomObjects = null;

                if (icoordbounds!=null) {
                        this.ibounds = icoordbounds;
                } else {
                        this.ibounds = new BoundingBox(this.iAtoms);
                }

                this.jAtoms = jAtoms;
                this.jAtomObjects = null;

                if (jAtoms==iAtoms) {
                        this.jbounds=ibounds;
                } else {
                        if (jcoordbounds!=null) {
                                this.jbounds = jcoordbounds;
                        } else {
                                this.jbounds = new BoundingBox(this.jAtoms);

                        }
                }

                fillGrid();
        }

        /**
         * Adds a set of coordinates, subsequent call to {@link #getIndicesContacts()} will produce the 
         * contacts, i.e. the set of points within distance cutoff.
         * The bounding box of the atoms will be computed based on input array.
         * Subsequent calls to method {@link #getAtomContacts()} will produce a NullPointerException 
         * since this only adds coordinates and no atom information.
         * @param atoms
         */
        public void addCoords(Point3d[] atoms) {
                addCoords(atoms, (BoundingBox) null);
        }

        /**
         * Adds a set of coordinates, subsequent call to {@link #getIndicesContacts()} will produce the 
         * contacts, i.e. the set of points within distance cutoff.
         * The bounds calculated elsewhere can be passed, or if null they are computed.
         * Subsequent calls to method {@link #getAtomContacts()} will produce a NullPointerException 
         * since this only adds coordinates and no atom information.
         * @param atoms
         * @param bounds
         */
        public void addCoords(Point3d[] atoms, BoundingBox bounds) {
                this.iAtoms = atoms;
                this.iAtomObjects = null;

                if (bounds!=null) {
                        this.ibounds = bounds;
                } else {
                        this.ibounds = new BoundingBox(iAtoms);
                }

                this.jAtoms = null;
                this.jAtomObjects = null;
                this.jbounds = null;

                fillGrid();
        }

        /**
         * Creates the grid based on the boundaries defined by all atoms given (iAtoms and jAtoms)
         * and places the atoms in their corresponding grid cells.
         * Checks also if the i and j grid overlap, i.e. the enclosing bounds of
         * the 2 grids (i and j) are no more than one cell size apart. If they don't
         * overlap then they are too far apart so there's nothing to calculate, we set
         * the noOverlap flag and then {@link #getIndicesContacts()} will do no calculation at all.
         */
        private void fillGrid() {

                if (jbounds!=null && !ibounds.overlaps(jbounds, cutoff)) {
                        //System.out.print("-");
                        noOverlap = true;
                        return;
                }

                findFullGridIntBounds();

                cells = new GridCell[1+(bounds[3]-bounds[0])/cellSize]
                                    [1+(bounds[4]-bounds[1])/cellSize]
                                    [1+(bounds[5]-bounds[2])/cellSize];

                int i = 0;
                for (Point3d atom:iAtoms) {

                        int xind = xintgrid2xgridindex(getFloor(atom.x));
                        int yind = yintgrid2ygridindex(getFloor(atom.y));
                        int zind = zintgrid2zgridindex(getFloor(atom.z));
                        if (cells[xind][yind][zind]==null) {
                                cells[xind][yind][zind] = new GridCell(this);
                        }
                        cells[xind][yind][zind].addIindex(i);
                        i++;
                }

                if (jAtoms==null) return;

                int j = 0;
                for (Point3d atom:jAtoms) {

                        int xind = xintgrid2xgridindex(getFloor(atom.x));
                        int yind = yintgrid2ygridindex(getFloor(atom.y));
                        int zind = zintgrid2zgridindex(getFloor(atom.z));
                        if (cells[xind][yind][zind]==null) {
                                cells[xind][yind][zind] = new GridCell(this);
                        }
                        cells[xind][yind][zind].addJindex(j);
                        j++;
                }

        }

        /**
         * Calculates an int array of size 6 into member variable bounds:
         * - elements 0,1,2: minimum x,y,z of the iAtoms and jAtoms
         * - elements 3,4,5: maximum x,y,z of the iAtoms and jAtoms
         */
        private void findFullGridIntBounds() {
                int[] iIntBounds = getIntBounds(ibounds);

                bounds = new int[6];
                if (jbounds==null) {
                        bounds = iIntBounds;
                } else {
                        int[] jIntBounds = getIntBounds(jbounds);
                        bounds[0] = Math.min(iIntBounds[0],jIntBounds[0]);
                        bounds[1] = Math.min(iIntBounds[1],jIntBounds[1]);
                        bounds[2] = Math.min(iIntBounds[2],jIntBounds[2]);
                        bounds[3] = Math.max(iIntBounds[3],jIntBounds[3]);
                        bounds[4] = Math.max(iIntBounds[4],jIntBounds[4]);
                        bounds[5] = Math.max(iIntBounds[5],jIntBounds[5]);
                }
        }

        /**
         * Returns an int array of size 6 :
         * - elements 0,1,2: minimum x,y,z (in grid int coordinates) of the given atoms
         * - elements 3,4,5: maximum x,y,z (in grid int coordinates) of the given atoms
         * @return
         */
        private int[] getIntBounds(BoundingBox coordbounds) {
                int[] bs = new int[6];
                bs[0] = getFloor(coordbounds.xmin);
                bs[1] = getFloor(coordbounds.ymin);
                bs[2] = getFloor(coordbounds.zmin);
                bs[3] = getFloor(coordbounds.xmax);
                bs[4] = getFloor(coordbounds.ymax);
                bs[5] = getFloor(coordbounds.zmax);
                return bs;
        }

        /**
         * Returns all contacts, i.e. all atoms that are within the cutoff distance.
         * If both iAtoms and jAtoms are defined then contacts are between iAtoms and jAtoms,
         * if jAtoms is null, then contacts are within the iAtoms.
         * @return
         */
        public AtomContactSet getAtomContacts() {

                AtomContactSet contacts = new AtomContactSet(cutoff);

                List<Contact> list = getIndicesContacts();

                if (jAtomObjects == null) {
                        for (Contact cont : list) {
                                contacts.add(new AtomContact(new Pair<Atom>(iAtomObjects[cont.getI()],iAtomObjects[cont.getJ()]),cont.getDistance()));
                        }

                } else {
                        for (Contact cont : list) {
                                contacts.add(new AtomContact(new Pair<Atom>(iAtomObjects[cont.getI()],jAtomObjects[cont.getJ()]),cont.getDistance()));
                        }
                }

                return contacts;
        }
        
        /**
         * Returns all contacts, i.e. all atoms that are within the cutoff distance.
         * If both iAtoms and jAtoms are defined then contacts are between iAtoms and jAtoms,
         * if jAtoms is null, then contacts are within the iAtoms.
         * @return
         * @deprecated use {@link #getAtomContacts()} instead
         */     
        @Deprecated
        public AtomContactSet getContacts() {
                return getAtomContacts();
        }
        
        /**
         * Returns all contacts, i.e. all atoms that are within the cutoff distance, as simple Contact objects containing the atom indices pairs and the distance.
         * If both iAtoms and jAtoms are defined then contacts are between iAtoms and jAtoms,
         * if jAtoms is null, then contacts are within the iAtoms.
         * @return
         */
        public List<Contact> getIndicesContacts() {

                List<Contact> list = new ArrayList<>();

                // if the 2 sets of atoms are not overlapping they are too far away and no need to calculate anything
                // this won't apply if there's only one set of atoms (iAtoms), where we would want all-to-all contacts
                if (noOverlap) return list;


                for (int xind=0;xind<cells.length;xind++) {
                        for (int yind=0;yind<cells[xind].length;yind++) {
                                for (int zind=0;zind<cells[xind][yind].length;zind++) {
                                        // distances of points within this cell
                                        GridCell thisCell = cells[xind][yind][zind];
                                        if (thisCell==null) continue;

                                        list.addAll(thisCell.getContactsWithinCell());

                                        // distances of points from this box to all neighbouring boxes: 26 iterations (26 neighbouring boxes)
                                        for (int x=xind-1;x<=xind+1;x++) {
                                                for (int y=yind-1;y<=yind+1;y++) {
                                                        for (int z=zind-1;z<=zind+1;z++) {
                                                                if (x==xind && y==yind && z==zind) continue;

                                                                if (x>=0 && x<cells.length && y>=0 && y<cells[x].length && z>=0 && z<cells[x][y].length) {
                                                                        if (cells[x][y][z] == null) continue;

                                                                        list.addAll(thisCell.getContactsToOtherCell(cells[x][y][z]));
                                                                }
                                                        }
                                                }
                                        }
                                }
                        }
                }

                return list;
        }

        /**
         * Fast determination of whether any atoms from a given set fall within
         * the cutoff of iAtoms. If {@link #addAtoms(Atom[], Atom[])} was called
         * with two sets of atoms, contacts to either set are considered.
         * @param atoms
         * @return
         */
        public boolean hasAnyContact(Point3d[] atoms) {
                return hasAnyContact(Arrays.asList(atoms));
        }
        public boolean hasAnyContact(Collection<Point3d> atoms) {
                for(Point3d atom : atoms) {
                        // Calculate Grid cell for the atom
                        int xind = xintgrid2xgridindex(getFloor(atom.x));
                        int yind = yintgrid2ygridindex(getFloor(atom.y));
                        int zind = zintgrid2zgridindex(getFloor(atom.z));

                        // Consider 3x3x3 grid of cells around point
                        for (int x=xind-1;x<=xind+1;x++) {
                                if( x<0 || cells.length<=x) continue;
                                for (int y=yind-1;y<=yind+1;y++) {
                                        if( y<0 || cells[x].length<=y ) continue;
                                        for (int z=zind-1;z<=zind+1;z++) {
                                                if( z<0 || cells[x][y].length<=z ) continue;
                                                
                                                GridCell cell = cells[x][y][z];
                                                // Check for contacts in this cell
                                                if(cell != null && cell.hasContactToAtom(iAtoms, jAtoms, atom, cutoff)) {
                                                        return true;
                                                }
                                        }
                                }
                        }
                }
                return false;
        }

        public double getCutoff() {
                return cutoff;
        }

        /**
         * Tells whether (after having added atoms to grid) the i and j grids are not overlapping.
         * Overlap is defined as enclosing bounds of the 2 grids being no more than one cell size apart.
         * @return true if the 2 grids don't overlap, false if they do
         */
        public boolean isNoOverlap() {
                return noOverlap;
        }
        
        protected Point3d[] getIAtoms() {
                return iAtoms;
        }
        
        protected Point3d[] getJAtoms() {
                return jAtoms;
        }

}