SM_StochasticBoltzmann.h

#ifndef SM_StochasticBoltzmann_H
#define SM_StochasticBoltzmann_H
 
//inherited class header
#include "SM_StochasticOutput.h"
 
//beam class header
#include "SM_Beam.h"
 
//suystem class header
#include "SM_System.h"
 
//runner header
#include "CORE_Run.h"
 
//timer include
#include "CORE_Chrono.h"
 
//system headers
#include "SM_MonteCarloSystem.h"
#include "SM_LandauLifschitzSystem.h"
 
class SM_StochasticBoltzmann  : public SM_StochasticOutput<SM_StochasticBoltzmann>  {
    
    // ATTRIBUTES
    
public:
    
    
    
private:
    //class types definition
    typedef SM_StochasticBoltzmann Self;
    typedef SM_StochasticOutput<Self> SuperClass;
    
   
    //easy axis
    std::array<tReal,SM_Constants::DIM> mU;
 
    //last step index registered 
    tIndex mStepIndex;
    
    //number of subdivisions of the support of the probability
    tInteger mN;
    
    //amplitude of one subdivision
    tReal mDTheta;
 
    //Effective temperature
    tReal mTeff;
    
    //boltzmann files
    tString mBoltzmannFileName;
    tString mSBoltzmannFileName;
    
public:
    // METHODS
  
    // CONSTRUCTORS 
  
    SM_StochasticBoltzmann(void) {
        mN=100;
        mDTheta=M_PI/mN;
        mU={0,0,1};
        setOutputDescription({
                "<theta_i>",
                "<number of angles of S with easy axes in [theta_i,theta_i+1[>",
                "<total number of computed angles>",
                "<normalized propability>",
                "<theorical value>"
            });
        //no output values for cycle
        std::valarray<tReal>& outputs=getOutputValues();
        outputs.resize(0);
    }
 
  
  
    // DESTRUCTORS 
public:
    
    virtual ~SM_StochasticBoltzmann(void) {
    }
 
 
 
public:
    //MEMORY
    
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+getContentsMemorySize();
    }
    
    virtual tMemSize getContentsMemorySize() const {
        tMemSize mem=SuperClass::getContentsMemorySize();
        mem+=mU.size()*sizeof(tIndex);
        mem+=sizeof(tChar)*mBoltzmannFileName.length();
        mem+=sizeof(tChar)*mSBoltzmannFileName.length();
        return mem;
    }
 
    inline static CORE_UniquePointer<Self> New() {
        return CORE_UniquePointer<Self>(new Self(),
                                        SM_Object::Delete());
    }
    inline static CORE_UniquePointer<Self> New(const Self& c) {
        CORE_UniquePointer<Self> p=New();
        p->copy(c);
        return p;
    }
 
    //SET & GET Methods
    //=================
    
    virtual void adimensionize(const SM_Material& material) override {
      
    }
    virtual void copy(const SM_StochasticOutputComponent& c) override {
        SuperClass::copy(c);
        const Self *cc=dynamic_cast<const Self*>(&c);
        if (cc!=null) {
            setSegmentsNumber(cc->getSegmentsNumber());
            setAxis(cc->getAxis());
        }
    }
 
 
    inline void setSegmentsNumber(const tInteger& n) {
        mN=n;
    }
    inline const tInteger&  getSegmentsNumber() const {
        return mN;
    }
     
    inline void setAxis(const std::array<tReal,SM_Constants::DIM>& axis) {
        mU=axis;
    } 
    inline void setAxis(const std::vector<tReal>& axis) {
        if (axis.size()<SM_Constants::DIM) {
            //mU=0;
            memset(mU.data(),0,sizeof(tReal)*SM_Constants::DIM);
            //mU=axis
            tReal *iU=mU.data();
            for(const auto& Uk:axis) {
                (*iU)=Uk;
                iU++;
            }
        } else {
            std::vector<tReal>::const_iterator iAxis=axis.begin();
            for(auto& Uk:mU) {
                Uk=(*iAxis);
                iAxis++;
            }
        }
    }
    inline const std::array<tReal,SM_Constants::DIM>&  getAxis() const {
        return mU;
    }
 
    
 
    
    //computed methods
    //================
    
    inline void printHeader(std::ofstream& file,const SM_Beam& beam) {
    
        //header is the data
        
        file<<"# Boltzmann normalized probability for 1 spin  with "<<getPackedSimulationsIndex()<<"/"<<getPackedSimulationsNumber()<<" packed simulations of size "<<beam.getBeamSize()<<" simulations generated by soft:"<<CORE_Run::GetSoftName()<<" Version:"<<CORE_Run::GetVersion()<<"\n";
        
        //print beam output
        tString str=beam.toString();
        functions_string::replaceAll("\n","\n#",str);
        file<<"#"<<str<<"\n";
        
        //print strochastic output
        str=this->toString();
        functions_string::replaceAll("\n","\n#",str);
        file<<"#"<<str<<"\n";
 
        //Teff=Ku/(T.k_B)
        mTeff=beam.getSystem().getMaterial().getAnisotropyEnergyFactor()/(beam.getSystem().getNoiseTemperature()*SM_Constants::K_B);
        file<<"#Teff:"<<mTeff<<"\n";
        file<<"#";
        for(const auto& d:getOutputDescription()) file<<d<<"\t";
        file<<"\n";
        
    }
 
  
    
    
    //implemented templated methods
    //=============================
    
    inline tBoolean open(const SM_Beam& beam) {
        
       
        //number of raws values number is null
        tIndex &nRaws=getPackedRawValuesNumber();
        nRaws=0;
 
        //set the size of the distribution
        std::valarray<tReal>& rawValues=getPackedRawValues();
        rawValues.resize(mN);
        
        //initialize to 0
        memset(&rawValues[0],0,rawValues.size()*sizeof(tReal));
        
        //set the amplitude of one segment
        mDTheta=M_PI/mN;
        
        //open the Boltzmann Propability file for writting
        mBoltzmannFileName=CORE_IO::GetAbsolutePath(this->getOutputPath()+"/"+this->getPrefix());
        if (getPackedSimulationsNumber()>1) {
            mBoltzmannFileName+="-"+std::to_string(getPackedSimulationsIndex())+"p";
        }
        //sysmmetric boltzamnn file
        mSBoltzmannFileName=mBoltzmannFileName;
        
        mBoltzmannFileName+=".bd";
        mSBoltzmannFileName+=".sbd";
        if (getIndex()==0) {//only one core print the output file 
            if (hasLogPerPackedSimulations() || (getPackedSimulationsNumber()==1)) {
                std::ofstream file(mBoltzmannFileName.c_str(),
                                   std::ios::out);
                if (file) {
                    printHeader(file,beam);
                    file.close();
                } else {
                    std::cout<<"Fatal Error : "<<mBoltzmannFileName<<" can not be opened \n";
                    return false;
                }
                
                std::ofstream sfile(mSBoltzmannFileName.c_str(),
                                    std::ios::out);
                if (sfile) {
                    printHeader(sfile,beam);
                    sfile.close();
                } else {
                    std::cout<<"Fatal Error : "<<mSBoltzmannFileName<<" can not be opened \n";
                    return false;
                }
            }
        }
 
        //set step index to an unreachable step value (infinity)
        mStepIndex=functions_numeric::getMax<tIndex>();
        
      
        
        return true;
    }
 
  
    
    
    inline tBoolean open(const tIndex& s,const SM_System& system) {
        tBoolean ok=true;
        mStepIndex=functions_numeric::getMax<tIndex>();
        return ok;
    }
    
    inline tBoolean  store(const SM_System& system) {
 
        if ((getRootIndex()==-1) || (getIndex()==0)) {
            if (system.getStepIndex()!=mStepIndex) {
                mStepIndex=system.getStepIndex();
                
                //get magnetic moment St
                const SM_RealField& S=system.getMagneticMomentDirections();
                
                //begin iterator on S
                const tReal* iS=&S[0];
                //end iterator on S
                const tReal* eS=iS;eS+=S.getSize();
                
                //sSU=\sum_p <Sp,U> for  particle
                tReal sSU=0;
                const tReal *iU=mU.data();
                const tReal *eU=iU;eU+=mU.size();
                while (iS!=eS) {
                    iU=mU.data();
                    while (iU!=eU) {
                        sSU+=(*iS)*(*iU);   
                        iS++;
                        iU++;
                    }
                }
                //arcos(<S,U>) in [0,PI[
                sSU=acos(sSU);
                
                //register the angle in the distribution array
                tIndex i=sSU/mDTheta;
                getPackedRawValues()[i]++;
                
                //number of raw data
                getPackedRawValuesNumber()++;
            }
        }
        
        return true;
    }
    
    inline tBoolean close(const tIndex& s,const SM_System& system,const tBoolean& hasSucceeded) {
        tBoolean ok=true;
        if (hasSucceeded) store(system);
        return ok;
    }
 
    inline tBoolean close(const SM_Beam& beam) {
        if (getRootIndex()==-1) {
            tIndex& rawValuesNumber=getPackedRawValuesNumber();
            std::valarray<tReal>& rawValues=getPackedRawValues();
            beam.getSystem().sumOverNetworks(0,rawValuesNumber,rawValues.size(),&rawValues[0]);
        }
        
        if (getIndex()==0) {
            if (hasLogPerPackedSimulations() || (getPackedSimulationsNumber()==1)) {
                std::ofstream file(mBoltzmannFileName.c_str(),std::ios::app);
                if (file) {
                    //std::cout<<"print Disribution "<<getIndex()<<"\n";
                    printDistribution(file,
                                      getPackedRawValuesNumber(),
                                      getPackedRawValues());
                    file.close();
                }
                
                std::ofstream sFile(mSBoltzmannFileName.c_str(),std::ios::app);
                if (sFile) {
                    printSymmetricDistribution(sFile,
                                               getPackedRawValuesNumber(),
                                               getPackedRawValues());
                    sFile.close();
                }
            }
        }
        return true;
            
    }
    
    virtual void closePackedSimulations(const SM_Beam& beam,
                                        const tIndex& rawValuesNumber,
                                        const std::valarray<tReal>& rawValues) override {
        
        if (getIndex()==0) {
            std::ofstream file(CORE_IO::GetAbsolutePath(this->getOutputPath()+"/"+this->getPrefix()+".bd").c_str(),
                               std::ios::out);
            if (file) {
                file<<"#Closed Packed Simulations\n";
                printHeader(file,beam);
                printDistribution(file,
                                  rawValuesNumber,
                                  rawValues);
                file.close();
            }
       
            
            std::ofstream sFile(CORE_IO::GetAbsolutePath(this->getOutputPath()+"/"+this->getPrefix()+".sbd").c_str(),
                                std::ios::out);
            if (sFile) {
                sFile<<"#Closed Symmetric Packed Simulations\n";
                printHeader(sFile,beam);
                printSymmetricDistribution(sFile,
                                           rawValuesNumber,
                                           rawValues);
                sFile.close();
            }
        }
    }
 
    
 
    inline static tBoolean LoadDistributionFile(const tString& fn,
                                                std::vector<tReal>& intervals,
                                                std::vector<tIndex>& intervalSamplesNumber,
                                                tIndex& samplesNumber,
                                                std::vector<tReal>& probability,
                                                std::vector<tReal>& boltzmannProbability)  {
        std::ifstream f(fn.c_str(),std::ios::in);
        if (f) {
            std::string line;
            
            intervals.clear();
            intervalSamplesNumber.clear();
            probability.clear();
            boltzmannProbability.clear();
            
            samplesNumber=0;
            tReal x;
            tIndex n,samplesNumber;
            tReal Pk;
            tReal thPk;
            while ((!f.eof()) && (std::getline(f,line))) {
                
                //ignore empty lines & comment lines
                if ((line.length()==0) || (line[0]=='#')) continue;
                
                //parser of line
                std::stringstream sstr;
                
                //read the line
                sstr<<line;
                
               
                sstr>>x;//theta_i
                sstr>>n;//number of samples n [theta_i,theta_{i+1})
                sstr>>samplesNumber;//total number of samples
                sstr>>Pk;//normalized probability
                sstr>>thPk;//boltzmann probability
                
                //get the line values
                //std::cout<<"line:"<<line<<" of "<< fn<<" -->"<<intervals.size()<<" theta="<<x<<" nSamples="<<n<<" nTotalSamplesNumber="<<samplesNumber<<" Pk:"<<Pk<<" \n";
                //register the values
                intervals.push_back(x);
                intervalSamplesNumber.push_back(n);
                probability.push_back(Pk);
                boltzmannProbability.push_back(thPk);
                
            }
            f.close();
            
            return true;
        }
        return false;
    }
    virtual tString toString() const override {
        std::stringstream ret;
        ret<<SuperClass::toString()<<"\n";
        ret<<"\t Number of subdivisions :"<<mN<<"\n";
        ret<<"\t Easy axis :"<<functions_array::toString(mU)<<"\n";
        return ret.str();
    }
 
private:
    
    tReal computeBoltzmannDistribution(const tReal& theta) const {
        tReal s=sin(theta);
        tReal ret=s;
        ret*=s;
        ret*=mTeff;
        ret*=-1.0;
        ret=exp(ret);
        ret*=s;
        return ret;
    }
    tReal computeBoltzmannDistributionMaximum() const {
        tReal Pmax=exp(-mTeff);
        if (mTeff>=0.5) {
            tReal Pmax2=sqrt(1./(2.*mTeff*exp(1.)));
            Pmax=(Pmax2>Pmax)?Pmax2:Pmax;
        }
        return Pmax;
    }
    
    inline void printDistribution(std::ofstream& file,
                                  const tInteger& rawDataSize,
                                  const std::valarray<tReal>& rawData) const {
        
        
        tReal PBmax=computeBoltzmannDistributionMaximum();
        //get the max value
        tReal maxPk=rawData[0];
        for(const auto& Pk:rawData) {
            maxPk=(maxPk<Pk)?Pk:maxPk;
        }
        
        if (maxPk>0) {
            tReal theta=0;
            for(const auto& Pk:rawData) {
                //file<<"#theta_i N_i N=sim N_i p=N_i/max(N_i)";   
                file<<theta<<"\t"<<Pk<<"\t"<<rawDataSize<<"\t"<< (Pk/maxPk)<<"\t"<<(computeBoltzmannDistribution(theta)/PBmax)<<"\n";
                theta+=mDTheta;
            }
        } else {
            tReal theta=0;
            for(const auto& Pk:rawData) {
                //file<<"#theta_i N_i N=sim N_i p=N_i/max(N_i)";   
                file<<theta<<"\t"<<Pk<<"\t"<<rawDataSize<<"\t"<< Pk<<"\t"<<(computeBoltzmannDistribution(theta)/PBmax)<<"\n";
                theta+=mDTheta;
            }
        }
        
    }
    inline void printSymmetricDistribution(std::ofstream& file,
                                           const tInteger& rawDataSize,
                                           const std::valarray<tReal>& rawData) const {
        
        tInteger i,n=rawData.size();
        
        //compute the max
        
        const tReal *iRawData=&rawData[0];
        const tReal *eRawData=iRawData;eRawData+=rawData.size();
        
        const tReal *iEnd=eRawData;iEnd--;
        
        n/=2;
        i=0;
        tReal nMaxSamples=0,nSamples;
        while (i<n) {
            
            nSamples=(*iRawData)+(*iEnd);
            nMaxSamples=(nMaxSamples<nSamples)?nSamples:nMaxSamples;
 
            //next i
            i++;
            iRawData++;
            iEnd--;
        }
        tReal PBmax=computeBoltzmannDistributionMaximum();
        
        if (nMaxSamples>0) {
            tReal theta=0;
           
            iRawData=&rawData[0];
            iEnd=eRawData;iEnd--;
            while (iRawData!=eRawData) {
                //file<<"#theta_i N_i N=sim N_i p=N_i/max(N_i)";
                if (iRawData!=iEnd) nSamples=(*iRawData)+(*iEnd);
                else nSamples=(*iRawData);
                
                file<<theta<<"\t"<<nSamples<<"\t"<<rawDataSize<<"\t"<< (nSamples/nMaxSamples)<<"\t"<<(computeBoltzmannDistribution(theta)/PBmax)<<"\n";
                
                theta+=mDTheta;
                iRawData++;
                iEnd--;
            }
        } else {
            tReal theta=0;
            for(const auto& Pk:rawData) {
                //file<<"#theta_i N_i N=sim N_i p=N_i/max(N_i)";   
                file<<theta<<"\t"<<Pk<<"\t"<<rawDataSize<<"\t"<< Pk<<"\t"<<(computeBoltzmannDistribution(theta)/PBmax)<<"\n";
                theta+=mDTheta;
            }
        }
        
    }
 
   
 
};
 
#endif