#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <unistd.h>
#include <math.h>
#include <algorithm>
using namespace std;
typedef vector<double> vd;
typedef vector<int> vi;
typedef vector< pair<int,int> > vpairs;

class Point2d
{
public:
    double x,y;
    Point2d (double X, double Y) { x=X; y=Y; }
};

class Point3d
{
public:
    double x,y,z;
    Point3d () {}
    Point3d (double X, double Y, double Z) { x=X; y=Y; z=Z; }
};

class Crossing
{
public:
    size_t arc1, arc2;
    Point2d point = Point2d(0,0), height = Point2d(0,0);
    bool over; // true if arc1 over arc2
    int index;
    Crossing (size_t a1, size_t a2, Point2d p, Point2d h, bool o, int i)
    {
        arc1=a1; arc2=a2; point=p; height=h; over=o; index=i;
    }
};

// declarations of functions
void Gauss_to_presentation(vi const& GaussCode, vector<vpairs>& Wirtinger);
bool sequence_to_Gauss (vector<Point3d>& Atoms, bool r1move, vi& GaussCode);
void shorten_sequence (vector<Point3d>& Atoms);
double volume (Point3d const& a, Point3d const& b, Point3d const& c, Point3d const& d);
double Det (Point3d const& A, Point3d const& B, Point3d const& C);
bool compare_dist (pair<size_t,double> a, pair<size_t,double> b);
void add_crossing (Point2d const& v, Crossing const& c, vector< vector<Crossing> >& vGaussCode);
void IntersectLines (Point2d const& v1, Point2d const& v2, Point2d const& v3, Point2d  const& v4, bool& flag, Point2d& v);
void read_pdb (string const& file_name, vector<Point3d>& Atoms);
void read_txt (string const& file_name, vector<Point3d>& Atoms);
void read_random_points (size_t const& nAtoms, vector<Point3d>& Atoms);
bool BoxesOverlap (vector<Point2d>const& FixPoly, vector<Point2d>const& VarPoly);
pair<Point2d, Point2d> BoundingBox (vector<Point2d>const& Poly);
void Det (Point2d const& u, Point2d const& v, double& det);
double min_coord (vd const& v);
double max_coord (vd const& v);
void Det (Point2d const& u, Point2d const& v, double& det);
double norm(Point2d v);
double sqdist(Point3d const& a,Point3d const& b);
void close_sequence(vd& xAtoms, vd& yAtoms, vd& zAtoms);
void PrintAtoms(vector<Point3d>const& Atoms);
void PrintWirtinger (vector<vpairs>const& Wirtinger);
void PrintWirtinger (vector<vpairs>const& Wirtinger, string& output);
void Distances ( vector<Point3d>const& Atoms, vector< pair<size_t,double> >& dAtoms );

int main()
{
    string folder = "/Users/kurlin/C++/KnottedProteins/pdbfiles/";
    string file_name, output;
    ofstream file;
    size_t nAtoms, nCrossings;
    vector<Point3d> Atoms;
    vector< pair<size_t,double> > dAtoms;
    string typeinput;
    bool r1move;
    vi GaussCode;
    vector<vpairs> Wirtinger; // set of relations with generators 1,...,m (full arcs)
    cout<<"Thank you for trying this KGG program computing the Knotted Graph Group.\n";
    cout<<"In line 1 of main(), please change the path to the folder on your computer.\n";
    cout<<"Input: ordered list of 3D coordinates of vertices of a polygonal chain K.\n";
    cout<<"Output: the group pi_1(R^3-K'), where K' is a closed simplification of K.\n";
    cout<<"Details are at http://kurlin.org/projects/invariants-knotted-graphs.html.\n";
    cout<<"If you have any questions, please e-mail Vitaliy, vitaliy.kurlin@gmail.com.\n";
    string welcome = "\nPlease type one of the following options for different types of input.\n";
    welcome += " demo: open polygonal trefoil\n";
    welcome += " pdb: protein backbone from PDB\n";
    welcome += " file: chain from your text file\n";
    welcome += " rand: random polygonal chain\n";
    welcome += " exit: finish the program \nInput: ";
    while (true)
    {
        cout<<welcome;
        cin>>typeinput;
        if (typeinput == "exit") break;
        if (typeinput == "demo")
        {
            Atoms.clear();
            Atoms.push_back( Point3d(-2,-2, 1) );
            Atoms.push_back( Point3d( 2, 2,-1) );
            Atoms.push_back( Point3d( 2,-1, 0) );
            Atoms.push_back( Point3d(-2,-1, 0) );
            Atoms.push_back( Point3d(-2, 2, 2) );
            Atoms.push_back( Point3d( 2,-2,-1) );
            cout<<"Vertices of the polygonal trefoil chain are\n";
            PrintAtoms(Atoms);
        }
        if (typeinput == "rand")
        {
            cout<<"\nA polygonal chain whose vertices are uniformly random in the unit cube.";
            cout<<"\nSuch a chain on more than 100 vertices is rarely simplifiable.";
            cout<<"\nPlease input the number of vertices.\nInput: ";
            cin>>nAtoms;
            srand( (int) time(0) );
            file_name = "random"+to_string(nAtoms);
            read_random_points (nAtoms, Atoms);
        }
        else if (typeinput == "pdb")
        {
            cout<<"\nPlease download a pdb file from http://www.rcsb.org/pdb to the same folder as main.cpp.";
            cout<<"\nExamples: 1V2X, 3OIL, 3OYS, 2RH3, 3NOU, 3NOT, 3NOP, 3ZQ5.";
            // 2LEN, 4JKJ, 1XD3, 1NS5, 1JS1, 1QMG, 1YRL, 4LRV, 3BJX, 4N2X, 3WJ8, 2EFV, 4UWE
            cout<<"\nThen please type the PDB id: with capital letters, but without .pdb, e.g. 1V2X.\nFile name: ";
            cin>>file_name;
            read_pdb( folder+file_name, Atoms);
        }
        else if (typeinput == "file")
        {
            cout<<"\nPlease put your input file with the extension .txt in the same folder as main.cpp.";
            cout<<"\nWrite down 3D coordinates of ordered vertices, in each line: x y z with spaces.";
            cout<<"\nThen please type the name of your file without extension, e.g. trefoil.\nFile name: ";
            cin>>file_name;
            read_txt( folder+file_name, Atoms);
        }
        else { cout<<"\nSorry, we couldn't understand your input\n"; continue; }
        // Stage 0: find initial crossings
        nAtoms = Atoms.size();
        if (nAtoms<3) { cout<<"Too few points: "<<nAtoms; continue; }
        r1move = false; // no simplification yet
        sequence_to_Gauss( Atoms, r1move, GaussCode );
        nCrossings = (size_t)( GaussCode.size() / 2 );
        file.open( folder + file_name + "_result.txt" );
        output = "The given chain has "+to_string( nAtoms ) + " vertices and " + to_string( nCrossings ) + " crossings.\n";
        cout<<output; file<<output;
        // Stage 1: shorten the given chain
        shorten_sequence( Atoms );
        // Stage 2: produce a Gauss code
        r1move = true; // try Reidemeister moves R1
        sequence_to_Gauss( Atoms, r1move, GaussCode );
        nAtoms = Atoms.size();
        nCrossings = (size_t)( GaussCode.size() / 2 );
        output="The simplified chain has "+to_string(nAtoms)+" vertices and "+to_string(nCrossings)+" crossings.\n";
        cout<<output; file<<output;
        output="The Gauss code has "+to_string(nCrossings)+" crossings: ";
        cout<<output; file<<output; output="";
        for ( size_t k = 0; k < GaussCode.size(); k++ ) output += to_string( GaussCode[k] ) + " ";
        file<<output; if (nCrossings>9) cout<<"too many to display"; else cout<<output;
        // Stage 3: write a presentation
        Gauss_to_presentation( GaussCode, Wirtinger );
        output="\nThe Wirtinger presentation has "+to_string( nCrossings ) + " generators and "+to_string(nCrossings-1)+" relations\n";
        cout<<output; file<<output;
        PrintWirtinger( Wirtinger, output );
        file<<output; if (nCrossings<=10) cout<<output;
        output = "The final vertices have the following coordinates.\n";
        for ( size_t k = 0; k < Atoms.size(); k++ )
            output += to_string( Atoms[k].x ) + " "+ to_string( Atoms[k].y ) + " "+ to_string( Atoms[k].z ) + "\n";
        file<<output; if (nCrossings<=10) cout<<output;
        cout<<"The output has been saved to the file "<<file_name<<"_result.txt"<<endl;
        file.close();
    }
    return 1;
}

void PrintWirtinger (vector<vpairs>const& Wirtinger)
{
    for (size_t i=1; i<Wirtinger.size(); i++)
    {
        for (size_t j=0; j<Wirtinger[i].size(); j++)
        {
            cout << " x_" << Wirtinger[i][j].first;
            if (Wirtinger[i][j].second == -1) cout << "^{-1}";
        }
        cout<<" = 1\n";
    }
}

void PrintWirtinger (vector<vpairs>const& Wirtinger, string& output)
{
    output = "";
    for (size_t i=1; i<Wirtinger.size(); i++)
    {
        for (size_t j=0; j<Wirtinger[i].size(); j++)
        {
            output += " x_" + to_string(Wirtinger[i][j].first);
            if (Wirtinger[i][j].second == -1) output += "^{-1}";
        }
        output += " = 1\n";
    }
}

void Gauss_to_presentation(vi const& GaussCode, vector<vpairs>& Wirtinger)
{
    size_t m = size_t ( GaussCode.size() / 2 );
    //cout << "\nRelations: " << m << endl;
    // build generators for overcrossing arcs between negative crossings
    vector<size_t> generators(2*m);
    size_t nGen=1;
    for (size_t k=0; k<2*m; k++)
    {
        if (GaussCode[k]<0) nGen++;
        generators[k] = nGen;
    }
    // find generator for each overcrossing arc
    vector<size_t> overcrossgen(m+1);
    for (size_t k=0; k<2*m; k++)
        if (GaussCode[k]>0)
            overcrossgen[GaussCode[k]] = generators[k];
    // write relation for each negative crossing
    Wirtinger.clear();
    for (size_t k=0; k<2*m; k++)
        if (GaussCode[k]<0)
        {
            size_t after = generators[k], before = after-1, over = overcrossgen[-GaussCode[k]];
            Wirtinger.push_back( { make_pair(over,1), make_pair(before,1), make_pair(over,-1), make_pair(after,-1)} );
            //cout << "x_" << over << " x_" << before << " x_" << over << "^{-1}" << " x_" << after << "^{-1} = 1 \n";
        }
    for (size_t i=0; i<m; i++)
        for (size_t j=0; j<Wirtinger[i].size(); j++)
            if (Wirtinger[i][j].first == m+1 ) Wirtinger[i][j].first = 1;
}

bool sequence_to_Gauss (vector<Point3d>& Atoms, bool r1move, vi& GaussCode)
{
    bool flag;
    size_t n=0;
    size_t nR1moves=0;
    size_t nAtoms = Atoms.size();
    vector< vector<Crossing> > vGaussCode(nAtoms-1);
    for (size_t i=0; i<nAtoms-1; i++)
        for (size_t j=i+1; j<nAtoms-1; j++)
        {
            Point2d v(0,0),v1(Atoms[i].x,Atoms[i].y), v2(Atoms[i+1].x,Atoms[i+1].y), v3(Atoms[j].x,Atoms[j].y), v4(Atoms[j+1].x,Atoms[j+1].y);
            vector<Point2d> v12 = { v1, v2 }, v34 = { v3, v4 };
            if ( !BoxesOverlap(v12,v34) ) continue; //
            IntersectLines(v1,v2,v3,v4,flag,v);
            if (!flag) continue;
            n++; // new intersection
            double t12 = (v.x-v1.x)/(v2.x-v1.x);
            double t34 = (v.x-v3.x)/(v4.x-v3.x);
            double z12 = (1-t12)*Atoms[i].z+t12*Atoms[i+1].z, z34 = (1-t34)*Atoms[j].z+t34*Atoms[j+1].z;
            bool over = true;
            if (z12 < z34) over = false;
            add_crossing ( v1, Crossing(i,j,v,Point2d(z12,z34),over,0), vGaussCode );
            add_crossing ( v3, Crossing(j,i,v,Point2d(z34,z12),!over,0), vGaussCode );
        }
    GaussCode.clear();
    int ind = 0; // abs index of crossing in Gauss code
    for (size_t i=0; i<Atoms.size()-1;i++)
        for (size_t k=0; k<vGaussCode[i].size();k++)
        {
            Crossing c = vGaussCode[i][k];
            if ( c.index != 0 ) // c was passed
            {
                GaussCode.push_back(c.index);
                int s = (int)GaussCode.size();
                if ( r1move && ( GaussCode[s-1] + GaussCode[s-2] == 0 ) ) // last 2 crossings can be removed
                {
                    nR1moves++;
                    size_t i=c.arc1, j=c.arc2, mn = min(i,j), mx=max(i,j);
                    Atoms[mn+1].x = c.point.x;
                    Atoms[mn+1].y = c.point.y;
                    Atoms[mx].x = c.point.x;
                    Atoms[mx].y = c.point.y;
                    if (mn==i)
                    {
                        Atoms[mn+1].z = c.height.x;
                        Atoms[mx].z = c.height.y;
                    }
                    else
                    {
                        Atoms[mn+1].z = c.height.y;
                        Atoms[mx].z = c.height.x;
                    }
                    //cout<<"r1move: i="<<i<<" j="<<j<<" mn="<<mn<<" Zmn="<<Atoms[mn+1].z<<" mx="<<mx<<" Zmx="<<Atoms[mx].z<<endl;
                    for (size_t l=mn+2; l<mx; l++) Atoms.erase(Atoms.begin()+l);
                    GaussCode.erase (GaussCode.begin()+s-1);
                    GaussCode.erase (GaussCode.begin()+s-2);
                }
            }
            else
            {
                ind++; // c is a new crossing
                int coef = 1;
                if (!c.over) coef = -1;
                c.index = ind * coef; // signed index
                size_t j = c.arc2, l=0;
                for (; l<vGaussCode[j].size();l++)
                    if (vGaussCode[j][l].arc2==i) break;
                if (l==vGaussCode[j].size())
                    cout << "symmetric crossing not found! \n";
                vGaussCode[j][l].index = -c.index;
                GaussCode.push_back(c.index);
            }
        }
    /* applying R1 to terminal codes
    int s = (int)GaussCode.size();
    if ( r1move && ( GaussCode[0] + GaussCode[s-1] == 0 ) )
    {
        Crossing c = vGaussCode[0][0];
        Atoms[0].x = c.point.x;
        Atoms[0].y = c.point.y;
        Atoms[s-1].x = c.point.x;
        Atoms[s-1].y = c.point.y;
        if (c.over)
        {
            Atoms[0].z = c.height.x;
            Atoms[s-1].z = c.height.y;
        }
        else
        {
            Atoms[0].z = c.height.y;
            Atoms[s-1].z = c.height.x;
        }
        GaussCode.erase (GaussCode.begin()+s-1);
        GaussCode.erase (GaussCode.begin());
        nR1moves++;
    }*/
    //cout<<"Number of Reidemeister moves R1 applied: "<<nR1moves;
    return r1move;
}

void Distances ( vector<Point3d>const& Atoms, vector< pair<size_t,double> >& dAtoms )
{
    dAtoms.resize( Atoms.size()-2 );
    for (size_t i=1; i<Atoms.size()-1; i++)
        dAtoms[i-1] = make_pair( i, sqdist( Atoms[i-1],Atoms[i+1] ) );
    sort (dAtoms.begin(), dAtoms.end(), compare_dist);
}

void shorten_sequence (vector<Point3d>& Atoms)
{
    if (Atoms.size()<=2) return;
    vector< pair<size_t,double> > dAtoms;
    Distances( Atoms, dAtoms );
    //for (int k=0; k<nAtoms-2; k++) cout << " k= " << k << " dist = " << dAtoms[k].second << " ind = " << dAtoms[k].first << endl;
    //cout << "\nErased: ";
    for ( int k = 0; k < dAtoms.size(); k++ )
    {
        size_t i = dAtoms[k].first; // index of vertex with k-th shortest distance
        bool flag = false; // no intersections yet
        for ( size_t j=0; j<Atoms.size()-1; j++ )
        {
            if ( (j+2>=i) && (j<=i+1) ) continue; // too close segments
            double vol_up = volume(Atoms[i-1],Atoms[i],Atoms[i+1],Atoms[j]);
            double vol_down = volume(Atoms[i-1],Atoms[i],Atoms[i+1],Atoms[j+1]);
            if ( vol_up * vol_down >=0 ) continue; // segment j on one side of the i-th plane
            double vol1 = volume(Atoms[i-1],Atoms[i],Atoms[j],Atoms[j+1]);
            double vol2 = volume(Atoms[i+1],Atoms[i],Atoms[j],Atoms[j+1]);
            double vol3 = volume(Atoms[i-1],Atoms[i+1],Atoms[j],Atoms[j+1]);
            if ( fabs( fabs(vol1)+fabs(vol2)+fabs(vol3) - fabs(vol_up)-fabs(vol_down) ) > 1e-10 ) continue; // intersection outside triangle
            flag = true; // intersection inside triangle
            break;
        }
        if (flag) continue; // we can't shorten triangle with intersection
        // now remove i-th vertex
        Atoms.erase( Atoms.begin() + i );
        if (Atoms.size()<=2) break;
        dAtoms.resize(Atoms.size()-2);
        for (size_t i=1; i<Atoms.size()-1; i++)
            dAtoms[i-1] = make_pair( i, sqdist(Atoms[i-1],Atoms[i+1]) );
        sort( dAtoms.begin(), dAtoms.end(), compare_dist );
        k=-1; // start from shortest triangle again
    }
}

double volume (Point3d const& a, Point3d const& b, Point3d const& c, Point3d const& d)
{
    Point3d A(b.x-a.x,b.y-a.y,b.z-a.z), B(c.x-a.x,c.y-a.y,c.z-a.z), C(d.x-a.x,d.y-a.y,d.z-a.z);
    return Det(A,B,C);
}

double Det (Point3d const& A, Point3d const& B, Point3d const& C)
{
    //cout << " A=(" << A.x << "," << A.y << "," << A.z <<") B=(" << B.x << "," << B.y << "," << B.z <<") C=(" << C.x << "," << C.y << "," << C.z;
    double vol = A.x*(B.y*C.z-B.z*C.y)-A.y*(B.x*C.z-B.z*C.x)+A.z*(B.x*C.y-B.y*C.x);
    //cout << " vol = "<< vol <<endl;
    return vol;
}

bool compare_dist (pair<size_t,double> a, pair<size_t,double> b) { return (a.second<b.second); }

void add_crossing(Point2d const& v, Crossing const& c, vector< vector<Crossing> >& vGaussCode)
{
    size_t i = c.arc1;
    Point2d p = c.point;
    size_t k = 0;
    for (; k<vGaussCode[i].size(); k++)
    {
        double d = (v.x-p.x)*(v.x-p.x) + (v.y-p.y)*(v.y-p.y);
        Point2d pk = vGaussCode[i][k].point;
        double dk = (v.x-pk.x)*(v.x-pk.x) + (v.y-pk.y)*(v.y-pk.y);
        if (d<dk) break;
    }
    vGaussCode[i].insert(vGaussCode[i].begin()+k,c);
}

void IntersectLines(Point2d const& v1, Point2d const& v2, Point2d const& v3, Point2d  const& v4, bool& flag, Point2d& v)
{
    flag = false; // we haven't found intersection
    Point2d v21(v2.x-v1.x,v2.y-v1.y), v43(v4.x-v3.x,v4.y-v3.y), v31(v3.x-v1.x,v3.y-v1.y);
    Point2d v32(v3.x-v2.x,v3.y-v2.y), v41(v4.x-v1.x,v4.y-v1.y), v42(v4.x-v2.x,v4.y-v2.y);
    // check if 2 points are one side of the line through other points
    double det123, det124;
    Det(v21, v31, det123);
    Det(v21, v41, det124);
    if (((det123 > 0) && (det124 > 0)) || ((det123 < 0) && (det124 < 0)))
        return; // v3, v4 on one side of the line through v1,v2
    double det341, det342;
    Det(v43, v31, det341);
    Det(v43, v32, det342);
    if (((det341 > 0) && (det342 > 0)) || ((det341 < 0) && (det342 < 0)))
        return; // v1, v2 are on one side of the line through v3,v4
    // now segments [v1,v2] and [v3,v4] are in the same line or meet, possibly at an endpoint
    double l43 = norm(v43);
    // check if v1 is in the line through v3, v4
    if (det341 == 0) // v1 is in the line through v3,v4
    {
        if ( (norm(v31) > l43) || (norm(v41) > l43) ) // v1 is outside segment [v3,v4]
        {
            if (det342 == 0) // v2 is also in the line through v3,v4
            {
                if ((norm(v32) > l43) || (norm(v42) > l43)) return; // v2 is outside segment [v3,v4]
                else return; // v2 is inside the segment [v3,v4], it should be considered earlier
            }
            else return; // v2 is outside the line through v3, v4, so no intersection
        }
        else return; // v1 is inside the segment [v3,v4], it should be considered earlier
    }
    // check if v2 is in the line through v3, v4
    if (det342 == 0) // v2 is in the line through v3,v4
    {
        if ((norm(v32) > l43) || (norm(v42) > l43)) // v2 is outside segment [v3,v4]
        {
            if (det341 == 0)  return; // v1 is also in the line through v3,v4
            else return; // v1 is outside the line through v3, v4, so no intersection
        }
        else return; // v2 is inside the segment [v3,v4], it should be considered earlier
    }
    double l21 = norm(v21);
    // check if v3 is on the line through v1, v2
    if (det123 == 0)  // v3 is in the line through v1,v2
    {
        if ((norm(v31) > l21) || (norm(v32) > l21)) // v3 is outside segment [v1,v2]
        {
            if (det124 == 0) // v4 is also in the line through v1,v2
            {
                if ((norm(v41) > l21) || (norm(v42) > l21)) return; // v4 is outside segment [v1,v2]
                else return; // v4 is inside the segment [v1,v2], it should be considered earlier
            }
            else return; // v4 is outside the line through v1, v2, so no intersection
        }
        else return; // v3 is inside the segment [v1,v2], it should be considered earlier
    }
    // check if v4 is on the line through v1, v2
    if (det124 == 0) // v4 is in the line through v1,v2
    {
        if ((norm(v41) > l21) || (norm(v42) > l21)) // v4 is outside segment [v1,v2]
        {
            if (det123 == 0) return; // v3 is also in the line through v1,v2
            else return; // v3 is outside the line through v1, v2, so no intersection
        }
        else return; // v4 is inside the segment [v1,v2], it should be considered earlier
    }
    // now segments [v1,v2] and [v3,v4] should meet at an interior point
    flag = true;
    double det; // = v21.x*v43.y - v21.y*v43.x;
    Det(v21, v43, det); // depends on variables x1, y1, x2, y2
    if (det == 0) // parallel lines, this case should be considered earlier
    {
        cout << " det=0 " << endl;
        return; // mistake
    }
    double top; // = v31.x*v43.y - v31.y*v43.x;
    Det(v31, v43, top); // depends on variables x1, y1
    double t = top / det; // scale t when 2 lines meet, depends on variables x1, y1, x2, y2
    if ((t <= 0) || (t >= 1)) // no internal intersection, this case should be considered earlier
    {
        cout << " t<=0 or t=>1 " << endl;
        return; // mistake
    }
    // intersection point of line segments (v1,v2) and (v3,v4)
    v.x = v1.x + t*v21.x;
    v.y = v1.y + t*v21.y;
}

void read_pdb(string const& file_name, vector<Point3d>& Atoms)
{
    Atoms.clear();
    ifstream file;
    file.open(file_name+".pdb");
    if (file.is_open()) cout << "File ./" << file_name << " is open.\n";
    else { cout << "Error opening " << file_name << ".\n"; return; }
    string line, col1, col2, col3, col4, col5, col6;
    double x,y,z;
    while(std::getline(file, line))  // read one line from ifs
    {
        std::istringstream iss(line); // access line as a stream
        iss >> col1;
        if (col1!="ATOM") continue; // no ATOM yet
        iss >> col2 >> col3;
        if (col3!="C") continue; // non-carbon atom
        iss >> col4 >> col5 >> col6; // irrelevant columns
        iss >> x >> y >> z;
        Atoms.push_back(Point3d(x,y,z));
    }
}

void read_txt(string const& file_name, vector<Point3d>& Atoms)
{
    Atoms.clear();
    ifstream file;
    file.open( file_name+".txt" );
    if (file.is_open()) cout << "File ./" << file_name << " is open.\n";
    else { cout << "Error opening " << file_name << ".\n"; return; }
    string line;
    double x,y,z;
    while(std::getline(file, line))  // read one line from ifs
    {
        std::istringstream iss(line); // access line as a stream
        iss >> x >> y >> z;
        Atoms.push_back(Point3d(x,y,z));
    }
}

void read_random_points ( size_t const& nAtoms, vector<Point3d>& Atoms)
{
    Atoms.clear();
    Atoms.push_back ( Point3d (-1,-1,-1) );
    for (size_t i=2; i<nAtoms; i++) {
        double rx = (double) rand() / RAND_MAX;
        double ry = (double) rand() / RAND_MAX;
        double rz = (double) rand() / RAND_MAX;
        Atoms.push_back ( Point3d (rx,ry,rz) );
    }
    Atoms.push_back ( Point3d (2,2,2) );
}

bool BoxesOverlap(vector<Point2d>const& FixPoly, vector<Point2d>const& VarPoly)
{
    pair<Point2d, Point2d> boxF = BoundingBox(FixPoly);
    pair<Point2d, Point2d> boxV = BoundingBox(VarPoly);
    //printPoint(boxF.first); printPoint(boxF.second);
    if (boxF.first.x >= boxV.second.x) return false;
    if (boxF.second.x <= boxV.first.x) return false;
    if (boxF.first.y >= boxV.second.y) return false;
    if (boxF.second.y <= boxV.first.y) return false;
    //std::cout << "overlap true"<<endl;
    return true;
}

pair<Point2d, Point2d> BoundingBox(vector<Point2d>const& Poly)
{
    double xmin = Poly[0].x;
    double xmax = Poly[0].x;
    double ymin = Poly[0].y;
    double ymax = Poly[0].y;
    for (auto it = Poly.begin(); it != Poly.end(); ++it) {
        double x = it->x;
        if (xmin > x) xmin = x;
        if (xmax < x) xmax = x;
        double y = it->y;
        if (ymin > y) ymin = y;
        if (ymax < y) ymax = y;
    }
    //printPoint(Point2d(xmin, ymin));
    //printPoint(Point2d(xmax, ymax));
    return make_pair(Point2d(xmin, ymin), Point2d(xmax, ymax));
}

double min_coord (vd const& v)
{
    size_t n = v.size();
    double m = 1e+64;
    for (size_t i=0;i<n;i++)
        if (v[i]<m) m=v[i];
    return m;
}

double max_coord (vd const& v)
{
    size_t n = v.size();
    double m = -1e+64;
    for (size_t i=0;i<n;i++)
        if (v[i]>m) m=v[i];
    return m;
}

void Det(Point2d const& u, Point2d const& v, double& det) { det = u.x * v.y - u.y * v.x; }

double norm(Point2d v) { return sqrt(v.x*v.x+v.y*v.y); }

double sqdist(Point3d const& a, Point3d const& b) {
    return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z);
}

void close_sequence(vd& xAtoms, vd& yAtoms, vd& zAtoms)
{
    size_t nAtoms = xAtoms.size();
    if (nAtoms<=1) return;
    double Xmin = min_coord(xAtoms);
    double Xmax = max_coord(xAtoms);
    double Ymin = min_coord(yAtoms);
    double Ymax = max_coord(yAtoms);
    double Zmin = min_coord(zAtoms);
    double Zmax = max_coord(zAtoms);
    
    // move the first point away from the bounding box
    double vx = xAtoms[0]-xAtoms[1];
    double vy = yAtoms[0]-yAtoms[1];
    double vz = zAtoms[0]-zAtoms[1];
    double vl = sqrt(vx*vx+vy*vy+vz*vz);
    if (vl==0) { cout << "vl=0" << endl; return; }
    vx/=vl; vy/=vl; vz/=vl;
    xAtoms[0] += vx * (Xmax-Xmin);
    yAtoms[0] += vy * (Ymax-Ymin);
    zAtoms[0] += vz * (Zmax-Zmin);
    
    // move the start and end points away from the bounding box
}

void PrintAtoms (vector<Point3d>const& Atoms)
{
    for (size_t i=0; i<Atoms.size(); i++)
        cout<<"v["<<i<< "] = ( " << Atoms[i].x << " , " << Atoms[i].y << " , " << Atoms[i].z << " )\n";
}