ForwardFlux.h

#pragma once 

#include "Method.h"
#include <fstream>
#include <random>
#include <deque>
#include "sys/stat.h"

namespace SSAGES
{

    class ForwardFlux : public Method
    {
    protected:
        class FFSConfigID
        {
           public:
            unsigned int l; 
            unsigned int n;      
            unsigned int a;      
            unsigned int lprev;      
            unsigned int nprev;      
            unsigned int aprev;      

            FFSConfigID(const unsigned int l, 
                        const unsigned int n, 
                        const unsigned int a, 
                        const unsigned int lprev, 
                        const unsigned int nprev, 
                        const unsigned int aprev): 
             l(l),n(n),a(a),lprev(lprev),nprev(nprev),aprev(aprev)
            {}

            FFSConfigID():
             l(0),n(0),a(0),lprev(0),nprev(0),aprev(0)
            {}
        };

        double _ninterfaces;

        std::vector<double> _interfaces;

        bool _interfaces_increase;

        double _cvvalue_previous;

        double _cvvalue;

        double _rate;

        std::vector<FFSConfigID> Lambda0ConfigLibrary;

        double _N0TotalSimTime;
        
        unsigned int _N0Target;

        double _fluxA0;

        std::vector<unsigned int> _M;

        std::vector<unsigned int> _A;

        //bool _computefluxA0;

        std::vector<double> _P;

        std::vector<unsigned int> _S;

        std::vector<unsigned int> _N ;
        
        bool _pop_tried_but_empty_queue;

        bool _initialFluxFlag;

        bool initializeQueueFlag;

        FFSConfigID myFFSConfigID;

        bool _saveTrajectories;

        unsigned int _nfailure_total;

        std::vector<std::vector<double>> _pB;

        unsigned int _current_interface;

        std::default_random_engine _generator;

        std::deque<FFSConfigID> FFSConfigIDQueue; 

        std::string _output_directory;

        std::ofstream _trajectory_file;

        uint iteration_;

        //-----------------------------------------------------------------
        // Private Functions
        //-----------------------------------------------------------------
        
        void CheckInitialStructure(const CVList&);

        void WriteInitialFlux();

        void AddNewIDsToQueue();

        bool HasReturnedToA(double);


        int HasCrossedInterface(double, double, unsigned int interface);

        void WriteFFSConfiguration(Snapshot *snapshot,FFSConfigID& ffsconfig, bool wassuccess);

        void ReadFFSConfiguration(Snapshot *,FFSConfigID&,bool);
       
        void ComputeInitialFlux(Snapshot*, const CVList&);

        virtual void CheckForInterfaceCrossings(Snapshot*, const class CVManager&) = 0;

        virtual void InitializeQueue(Snapshot*, const CVList&) =0;


        void ComputeTransitionProbabilities();

        void PrintQueue();

        void PopQueueMPI(Snapshot*, const CVList&, unsigned);
        

        void FluxBruteForce(Snapshot*, const CVList&);
        
        void ReconstructTrajectories(Snapshot *);

        void ComputeCommittorProbability(Snapshot *);
        
        void AppendTrajectoryFile(Snapshot*, std::ofstream&);

        void OpenTrajectoryFile(std::ofstream&);

    public:

        ForwardFlux(const MPI_Comm& world,
                    const MPI_Comm& comm, 
                    double ninterfaces, std::vector<double> interfaces,
                    unsigned int N0Target, std::vector<unsigned int> M,
                    bool initialFluxFlag, bool saveTrajectories,
                    unsigned int currentInterface, std::string output_directory, unsigned int frequency) : 
         Method(frequency, world, comm), _ninterfaces(ninterfaces), _interfaces(interfaces), _N0Target(N0Target), 
         _M(M), _initialFluxFlag(initialFluxFlag), _saveTrajectories(saveTrajectories), _current_interface(currentInterface),
          _output_directory(output_directory), _generator(1), iteration_(0) 
        {
            //_output_directory = "FFSoutput";
            mkdir(_output_directory.c_str(),S_IRWXU); //how to make directory?

            //_current_interface = 0;
            _N0TotalSimTime = 0;

            if (!_initialFluxFlag)
                initializeQueueFlag = true;

            _A.resize(_ninterfaces);
            _P.resize(_ninterfaces);
            _S.resize(_ninterfaces);
            _N.resize(_ninterfaces);

            //_N[_current_interface] = _N0Target;
            /*_M.resize(_ninterfaces);

            //code for setting up simple simulation and debugging
            _ninterfaces = 5;
            int i;
            _interfaces.resize(_ninterfaces);
            _interfaces[0]=-1.0;
            _interfaces[1]=-0.95;
            _interfaces[2]=-0.8;
            _interfaces[3]= 0;
            _interfaces[4]= 1.0;

            _saveTrajectories = true;

            _initialFluxFlag = true;

            for(i=0;i<_ninterfaces;i++) _M[i] = 50;

            _N0Target = 100;*/
            //_N0Target = 100;
            _nfailure_total = 0;


            //check if interfaces monotonically increase or decrease
            bool errflag = false;
            if (_interfaces[0]< _interfaces[1]) _interfaces_increase = true;
            else if (_interfaces[0] > _interfaces[1]) _interfaces_increase = false;
            else errflag = true;
            for (unsigned int i=0;i<_ninterfaces-1;i++)
            {
                if((_interfaces_increase) && (_interfaces[i] >= _interfaces[i+1]))
                    errflag = true;
                else if ((!_interfaces_increase) && (_interfaces[i] <= _interfaces[i+1]))
                    errflag = true;
            }
            if (errflag)
            {
                std::cerr << "Error! The interfaces are poorly defined. They must be monotonically increasing or decreasing and cannot equal one another! Please fix this.\n";
                for (auto interface : _interfaces)
                    std::cerr << interface << " ";
                std::cerr << "\n";
                MPI_Abort(world_, EXIT_FAILURE);
            }

            // This is to generate an artificial Lambda0ConfigLibrary, Hadi's code does this for real
            // THIS SHOULD BE SOMEWHERE ELSE!!! 
            Lambda0ConfigLibrary.resize(_N0Target);
            std::normal_distribution<double> distribution(0,1);
            for (unsigned int i = 0; i < _N0Target ; i++)
            {
                Lambda0ConfigLibrary[i].l = 0;
                Lambda0ConfigLibrary[i].n = i;
                Lambda0ConfigLibrary[i].a = 0;
                Lambda0ConfigLibrary[i].lprev = 0;
                Lambda0ConfigLibrary[i].nprev = i;
                Lambda0ConfigLibrary[i].aprev = 0;
                //FFSConfigID ffsconfig = Lambda0ConfigLibrary[i];
            }
            /*// Write the dump file out
            std::ofstream file;
            std::string filename = _output_directory + "/l" + std::to_string(ffsconfig.l) + "-n" + std::to_string(ffsconfig.n) + ".dat";
            file.open(filename.c_str());

            //first line gives ID of where it came from
            file << ffsconfig.lprev << " " << ffsconfig.nprev << " " << ffsconfig.aprev << "\n";
            //write position and velocity
            file << "1 -1 0 0 " << distribution(_generator) << " "  << distribution(_generator) << " 0\n";

            }
            */
            _pop_tried_but_empty_queue = false;
        }


        void PreSimulation(Snapshot* snapshot, const class CVManager& cvmanager) override;


        virtual void PostIntegration(Snapshot* snapshot, const class CVManager& cvmanager) = 0;


        void PostSimulation(Snapshot* snapshot, const class CVManager& cvmanager) override;

        static ForwardFlux* Build(const Json::Value& json, 
                                      const MPI_Comm& world,
                                      const MPI_Comm& comm,
                                      const std::string& path);
    };
}

/*
File Formats:
_indexfile
interface(some integer) dump_file_name(a string that contains interface and trial number)
example: 1 dump_1_10.xyz

dumpfile
atomid posx posy posz vx vy vz


*/