/* This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */

package etomica.virial;

import etomica.atom.IAtomList;
import etomica.box.Box;
import etomica.integrator.mcmove.MCMoveAtom;
import etomica.space.Space;
import etomica.space.Vector;
import etomica.util.random.IRandom;

/**
 * Class that samples molecule positions based on a tree of hard spheres of
 * diameter sigma.  The tree structure is chosen randomly as a Prüfer sequence.
 * <p>
 * https://en.wikipedia.org/wiki/Pr%C3%BCfer_sequence#Algorithm_to_convert_a_Pr.C3.BCfer_sequence_into_a_tree
 *
 * @author Andrew
 */
public class MCMoveClusterAtomHSTree extends MCMoveAtom {

    public MCMoveClusterAtomHSTree(IRandom random, Space _space, double sigma) {
        super(random, null, _space);
        this.sigma = sigma;
    }

    public void setBox(Box box) {
        super.setBox(box);
        int n = box.getLeafList().size();
        degree = new int[n];
        a = new int[n-2];
        inserted = new int[n];
        bonds = new int[n*(n-1)/2][2];
    }

    public boolean doTrial() {

        IAtomList leafAtoms = box.getLeafList();
        int n = leafAtoms.size();
        for (int i=0; i<n; i++) {
            degree[i] = 1;
        }
        for (int i=0; i<n-2; i++) {
            a[i] = random.nextInt(n);
            degree[a[i]]++;
        }

        int numBonds = 0;
        for (int i=0; i<n-2; i++) {
            int ii = a[i];
            for (int j=0; j<n; j++) {
                if (degree[j] == 1) {
                    bonds[numBonds][0] = ii;
                    bonds[numBonds][1] = j;
                    numBonds++;
                    degree[ii]--;
                    degree[j]--;
                    break;
                }
            }
        }
        int u = -1, v = -1;
        for (int i=0; i<n; i++) {
            if (degree[i] == 1) {
                if (u==-1) {
                    u = i;
                } else {
                    v = i;
                }
            }
        }
        bonds[numBonds][0] = u;
        bonds[numBonds][1] = v;
        numBonds++;

        leafAtoms.get(0).getPosition().E(0);
        inserted[0] = 0;
        int numInserted = 1;
        // inserted is a list of points that have inserted, but not coordinated
        // coordinated is a list of points that have been inserted and coordinated
        int coordinatedMask = 0;
        while (numInserted > 0) {
            int nbr = inserted[numInserted-1];
            numInserted--;
            coordinatedMask |= 1<<nbr;
            for (int iBond = 0; iBond<numBonds; iBond++) {
                int[] b = bonds[iBond];
                int nbr2 = -1;
                if (b[0] == nbr) {
                    nbr2 = b[1];
                } else if (b[1] == nbr) {
                    nbr2 = b[0];
                } else {
                    continue;
                }
                if ((coordinatedMask & (1<<nbr2)) != 0) {
                    // already inserted nbr2, move along
                    continue;
                }
                // insert nbr2 around nbr
                Vector pos = leafAtoms.get(nbr2).getPosition();

                pos.setRandomInSphere(random);
                double sig = getSigma(nbr, nbr2);
                if (sig < 0) {
                    // we want to force the position to be in the well (between 1 and sigma)
                    sig = -sig;
                    while (pos.squared() < 1 / (sig * sig)) {
                        pos.setRandomInSphere(random);
                    }
                }
                pos.TE(sig);
                pos.PE(leafAtoms.get(nbr).getPosition());
                inserted[numInserted] = nbr2;
                numInserted++;
            }
        }

        ((BoxCluster)box).trialNotify();
        return true;
    }

    // override this to do a mixture
    protected double getSigma(int i, int j) {
        return sigma;
    }

    public double getChi(double temperature) {
        return 1;
    }

    public void rejectNotify() {
        throw new RuntimeException("nope");
    }

    public void acceptNotify() {
        ((BoxCluster)box).acceptNotify();
    }

    protected final double sigma;
    protected int[][] bonds;
    protected int[] degree, a;
    protected int[] inserted;
}

@@ -3048,18 +3048,20 @@
 ask & (1
+ 
 %3C%3C
+ 
 nbr2)) !
@@ -3278,17 +3278,99 @@
 tion();%0A
-%0A
+                boolean insideBox = false;%0A                while(!insideBox) %7B%0A    
         
@@ -3424,16 +3424,20 @@
         
+   
+ 
 double s
@@ -3474,24 +3474,28 @@
             
+    
 if (sig %3C 0)
@@ -3493,24 +3493,28 @@
 (sig %3C 0) %7B%0A
+    
             
@@ -3614,16 +3614,20 @@
         
+    
 sig = -s
@@ -3630,16 +3630,20 @@
 = -sig;%0A
+    
         
@@ -3712,32 +3712,36 @@
                 
+    
 pos.setRandomInS
@@ -3767,32 +3767,36 @@
                 
+    
 %7D%0A              
@@ -3789,34 +3789,42 @@
                 
+ 
-%7D%0A
+   %7D%0A    
                 
@@ -3836,16 +3836,384 @@
 E(sig);%0A
+                    insideBox = true;%0A                    for (int i = 0; i %3C pos.getD(); i++) %7B%0A                        if (Math.abs(pos.getX(i)) %3E box.getBoundary().getBoxSize().getX(i)/2) %7B%0A                            insideBox = false;%0A                            break;%0A                        %7D%0A                    %7D%0A                    if (insideBox) %7B%0A        
         
@@ -4270,32 +4270,40 @@
                 
+      
+  
 inserted%5BnumInse
@@ -4324,32 +4324,40 @@
                 
+        
 numInserted++;%0A 
@@ -4351,24 +4351,64 @@
 Inserted++;%0A
+                    %7D%0A                %7D%0A
             

Handle case where part of tree is out of box