SM_StochasticJs.h

#ifndef SM_StochasticJs_H
#define SM_StochasticJs_H
 
//inherited class header
#include "SM_StochasticOutput.h"
 
//beam class header
#include "SM_LandauLifschitzBeam.h"
 
//suystem class header
#include "SM_System.h"
 
//runner header
#include "CORE_Run.h"
 
//string utility header
#include "core_string.h"
#include "core_numeric.h"
 
 
 
class SM_StochasticJs  : public SM_StochasticOutput<SM_StochasticJs>  {
  
    // ATTRIBUTES
  
public:
   
  
   
private:
 
    //class types definition
    typedef SM_StochasticJs Self;
    typedef SM_StochasticOutput<Self> SuperClass;
 
    tBoolean mIsLLBeam;
 
    tString mCTPackedFileName;
    tString mCTFileName;
   
    tIndex mStepIndex;
    tReal mE;//esperance
 
 
    //S at equilibrium
    SM_RealField mSeq;
    std::array<tReal,SM_Constants::DIM> mBarS;//mean of S over the particles
    tReal mNormS;
 
    tUCInt mStochasticIndex;
    tUCInt mMeanSIndex;
    tUCInt mMeanJIndex;
    tUCInt mNSamplesIndex;
    tUCInt mNSimsIndex;
 
    tUCInt mESRawIndex;
    tUCInt mEJRawIndex;
    tUCInt mNRawIndex;
    
 
    tString mJFunctionName;
    std::function<void(const tReal&,const tReal& ,const tReal&,tReal& )> mJFunction;
    
    std::ofstream *mJsFile;
    
public:
    // METHODS
  
    // CONSTRUCTORS 
  
    SM_StochasticJs(void) {
      
        setOutputDescription({
                "<stochastic noise>",
                "<Mean of |S|>",
                "<mean of J>",
                "<samples number>",
                "<number of simulations>",
                });
        
        mStochasticIndex=0;
        mMeanSIndex=1;
        mMeanJIndex=2;
        mNSamplesIndex=3;
        mNSimsIndex=4;
 
        mESRawIndex=0;
        mEJRawIndex=1;
        mNRawIndex=2;
        std::valarray<tReal>& outputs=getOutputValues();
        outputs.resize(5);
        
        memset(&outputs[0],0,outputs.size()*sizeof(tReal));
 
        setJFunction("increase");
        
        mJsFile=null;
    }
 
  
  
    // DESTRUCTORS 
public:
    
    virtual ~SM_StochasticJs(void) {
    }
 
 
 
public:
    //MEMORY
    
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+getContentsMemorySize();
    }
    
    virtual tMemSize getContentsMemorySize() const {
        tMemSize mem=SuperClass::getContentsMemorySize();
        mem+=mCTFileName.length()*sizeof(tChar);
        mem+=mCTPackedFileName.length()*sizeof(tChar);
        mem+=mBarS.size()*sizeof(tReal);
        mem+=mSeq.getContentsMemorySize();
        mem+=mJFunctionName.length()*sizeof(tChar);
        mem+=mCTPackedFileName.length()*sizeof(tChar);
        mem+=mCTFileName.length()*sizeof(tChar);
        mem+=(mJsFile==null)?0:sizeof(*mJsFile);
        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 copy(const SM_StochasticOutputComponent& c) override {
        SuperClass::copy(c);
    }
 
 
    inline void setJFunction(const tString& FName) {
        mJFunctionName=FName;
        switch(FName[0]) {
       
        case 'd'://decrease
            mJFunction=([](const tReal& J,const tReal& eps,const tReal& normal,tReal& v) {
                            //v=J/(1+|eps.N|)
                            v=normal;
                            v*=eps;
                            v=fabs(v);
                            v++;
                            v=1./v;
                            v*=J;
                        });
            break;
        case 'c'://cubic-decrease
            mJFunction=([](const tReal& J,const tReal& eps,const tReal& normal,tReal& v) {
                            //v=J/(1+|eps.N|)^3
                            v=normal;
                            v*=eps;
                            v=fabs(v);
                            v++;
                            v=v*v*v;
                            v=1./v;
                            v*=J;
                        });
            break;
        case 'i'://increase
        default:
            mJFunction=([](const tReal& J,const tReal& eps,const tReal& normal,tReal& v) {
                            //v=std::max(0,J*(1+eps.N))
                            v=normal;
                            v*=eps;
                            v++;
                            v*=J;
                            v=std::max((tReal)0,v);
                        });
            break;
        }
        
 
    }
    inline const std::function<void(const tReal&,const tReal& ,const tReal&,tReal& )>& getJFunction() const {
        return mJFunction;
    }
    
    inline void setEquilibriumState(const SM_RealField& Seq) {
        mSeq=Seq;
    } 
    inline const SM_RealField&  getEquilibriumState() const {
        return mSeq;
    }
    inline SM_RealField&  getEquilibriumState()  {
        return mSeq;
    }
    
    virtual void adimensionize(const SM_Material& material) override {
    }
    
    inline void printHeader(std::ofstream& file,const SM_Beam& beam) {
        //header is the data
        //get the system
        file<<"# Stochastic J Curie Temperature   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:"+beam.toString();
        core_string::replaceAll("\n","\n#",str);
        file<<"#"<<str<<"\n";
 
        //print strochastic output
        str="StochasticOutput:"+toString();
        core_string::replaceAll("\n","\n#",str);
        file<<"#"<<str<<"\n";
        
        //Tc=J*\epsilon * z / 3.kB
        tReal Tc=beam.getSystem().getMaterial().getExchangeEnergyFactor();
        Tc*=beam.getSystem().getMaterial().getSpinWaveMFCorrection();
        Tc*=beam.getSystem().getMaterial().getCrystalStructure().getCoordinationNumber();
        Tc/=(3.0*SM_Constants::K_B);
        
        file<<"# T: temperature N : simulations number\n";
        file<<"# for all simulation s in [0,N[:\n";
        //file<<"# \t for all step t in [t_{eq},T[, bar S_k(eps,s,t)=frac{1}{P} sum_{p=0}^{p=P-1} <S^s_k(eps,t,p),S^0_k(0,t_{eq},p)>\n";
        file<<"# \t for all step t in [t_{eq},T[, bar S_k(s,t)=frac{1}{P} sum_{p=0}^{p=P-1} S^s_k(t,p) \n"; 
        file<<"# \t E_s=frac{1}{T} sum_{t=0}^{t=T-1} |bar S_k(eps,s,t)|\n";
        file<<"# E=frac{1}{N} sum_{s=0}^{s=N-1} E_s\n";
        file<<"\n";
        file<<"# Curie Temperature: "<<Tc<<"\n";
        file<<"#";
        for(const auto& d:getOutputDescription()) file<<d<<"\t";
        file<<"\n";
    }
     
  
    
    //implemented templated methods
    //=============================
    
    inline tBoolean open(const SM_Beam& beam) {
        tString routineName="SM_StochasticJs::open()";
        CORE_Profiler::StartCallRoutine(routineName);
  
        //set step index to an unreachable step value (infinity)
        mStepIndex=core_numeric::getInfinity<tIndex>();
 
            
        if (getIndex()==0) {
            std::valarray<tReal>& outputs=getOutputValues();
            
            //temperature
            outputs[mStochasticIndex]=beam.getSystem().getStochasticNoise();
            
            //open the core & global CT Files for writting
            //--------------------------------------------
            tString prefix=this->getPrefix();
            
            if ((prefix.length()>0)&&(getIndex()==0)) {
                //save the stochastic output for simulations in file
                mCTFileName=CORE_IO::GetAbsolutePath(this->getOutputPath()+"/"+prefix);
                mCTPackedFileName=mCTFileName;
                if (getPackedSimulationsNumber()>1) {
                    mCTPackedFileName+="-"+std::to_string(getPackedSimulationsIndex())+"p";
                }
                mCTPackedFileName+=".ct";
                mCTFileName+=".ct";
            } else {
                //prefix is null : no saving
                mCTFileName="";
                mCTPackedFileName="";
            }
            //std::cout<<"JCT file name:"<< mCTFileName<<" index:"<<getIndex()<<" prefix:"<<prefix<<"\n";
            if (mCTPackedFileName.length()>0) {
                if (hasLogPerPackedSimulations() || (getPackedSimulationsNumber()==1) ) {
                    
                    std::ofstream file(mCTPackedFileName.c_str(),std::ios::out);
                    if (file) {
                        printHeader(file,beam);
                        
                    } else {
                        std::cout<<"Fatal Error : "<<mCTPackedFileName<<" can not be opened \n";
                        return false;
                    }
                }
            }
            
            mJsFile=new std::ofstream(CORE_IO::GetAbsolutePath(this->getOutputPath()+"/"+prefix+".j0s").c_str(),
                                      std::ios::out);
            (*mJsFile)<<"#simulation epsilon time <J for the particle 0>\n";
            
        
            
            
            //iinit the stochastic computings
            //---------------------------------
            //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(3);//E(|S|),E(J),nSimulations per pack
            
            //initialize to 0
            memset(&rawValues[0],0,rawValues.size()*sizeof(tReal));
            
            //init the class state
            //----------------------
          
        }
        
        CORE_Profiler::EndCallRoutine(routineName);
        return true;
    }
    
    
    
 
    inline tBoolean open(const tIndex& s,const SM_System& system) {
        
        tString routineName="SM_StochasticJs::open(s,system)";
        CORE_Profiler::StartCallRoutine(routineName);
 
        tBoolean ok=true;
        
        //initialize the step index
        mStepIndex=core_numeric::getInfinity<tIndex>();
 
   
        
        CORE_Profiler::EndCallRoutine(routineName);
        
        return ok;
    }
    
    inline tBoolean  store(const SM_System& system) {
        tString routineName="SM_StochasticJs::store()";
        CORE_Profiler::StartCallRoutine(routineName);
 
      
        //compute the esperance at time t
        if ((getIndex()==0) || (getRootIndex()==-1)) {
            computeEsperanceContribution(system.getStepIndex(),system);
        }
        if (getIndex()==0) {
            //print J for the connections of the particles 0
            tIndex i,n=system.getNetwork().getNeighboringParticlesNumber(0);
            const SM_LandauLifschitzSystem* lls=dynamic_cast<const SM_LandauLifschitzSystem*>(&system);
            
            if (mJsFile!=null) {
                (*mJsFile)<<lls->getStochasticNoise()<<"\t"<<lls->getTime()<<"\t";
                for (i=0;i<n;i++) {
                    (*mJsFile)<<lls->getNetwork().getConnectionValues()[i]<<"\t";
                }
                (*mJsFile)<<"\n";
            }
        }
        
        CORE_Profiler::EndCallRoutine(routineName);
        return true;
    }
    
 
    
    inline tBoolean close(const tIndex& s,const SM_System& system,const tBoolean& hasSucceeded) {
        tBoolean ok=true;
        tString routineName="SM_StochasticJs::close(s,system,hasSucceeded)";
        CORE_Profiler::StartCallRoutine(routineName);
        if (hasSucceeded) {
            if ((getIndex()==0) || (getRootIndex()==-1)) { 
                //compute the esperance at last time 
                computeEsperanceContribution(system.getStepIndex(),system);
            }
        }
    
        CORE_Profiler::EndCallRoutine(routineName);
        return ok;
    }
 
    inline tBoolean close(const SM_Beam& beam) {
        tString routineName="SM_StochasticJs::close(beam)";
        CORE_Profiler::StartCallRoutine(routineName);
 
        if (getIndex()==0) {
            //number of simulation of the beam
            getPackedRawValues()[mNRawIndex]=beam.getBeamSize();
            
            if ( hasLogPerPackedSimulations() || (getPackedSimulationsNumber()==1) ) {
                registerStochasticOutput(getPackedRawValuesNumber(),getPackedRawValues());
                if (mCTPackedFileName.length()>0) {
                    //\f$ E(\varepsilon)= \displaystyle \sqrt { \sum_{k=0}^{k=d-1} E_k(\bar mu_k(\varepsilon,T))^2 
                    std::ofstream file(mCTPackedFileName.c_str(),std::ios::app);
                    //file<<"# raw values number:"<<getPackedRawValuesNumber()<<"\n";
                    //file<<"# raw values:";
                    //for(const auto& v:getPackedRawValues()) file<<"\t"<<v;
                    //file<<"\n";
                    flushStochasticOutput(file);
                    file.close();
                }
            }
            if (mJsFile!=null) {
                mJsFile->close();
                mJsFile=null;
            }
        }
        CORE_Profiler::EndCallRoutine(routineName);
        return true;
    
    }
 
 
    virtual void closePackedSimulations(const SM_Beam& beam,
                                        const tIndex& rawValuesNumber,
                                        const std::valarray<tReal>& rawValues) override {
        tString routineName="SM_StochasticJs::closePackedSimulations(beam,nRaws,raws)";
        CORE_Profiler::StartCallRoutine(routineName);
 
        if (getIndex()==0) {
            registerStochasticOutput(getPackedRawValuesNumber(),getPackedRawValues());
            if (mCTFileName.length()>0) {
                std::ofstream file(mCTFileName,std::ios::out);
                if (file) {
                    //file<<"#Closed Packed Simulations\n";
                    //file<<"# raw values number:" <<rawValuesNumber<<"\n";
                    //file<<"# raw values:";
                    //for(const auto& v:rawValues) file<<"\t"<<v;
                    //file<<"\n";
                    printHeader(file,beam);
                    flushStochasticOutput(file);
                    file.close();
                }
            }
        }
        
        CORE_Profiler::EndCallRoutine(routineName);
    }
    
    virtual tString toString() const override {
        std::stringstream ret;
        ret<<SuperClass::toString();
        ret<<"\t JFunctionName:"<<mJFunctionName<<"\n";
        return ret.str();
    }
 
    
private:
 
    inline void computeEsperanceContribution(const tIndex& stepIndex,const SM_System& system) {
 
        
        if (mStepIndex!=stepIndex) {
            mStepIndex=stepIndex;
            tReal sumJij=0;
            tIndex nConnections=0;
 
            if ((getIndex()==0) || (getRootIndex()==-1)) {
                
                //compute barS=1/P \sum S(p)
                //for OMP/MPIC version: barS=1/P sum(P) for p of the whole domain
                //for MPI version :  barS=1/P sum(P) for p of the whole domain
                system.computeMagneticMomentDirectionsMeanOverNetworks(0,mBarS);
 
                //compute  mean J
                const std::valarray<tReal>& Js=system.getNetwork().getConnectionValues();
                sumJij=0;
                for(const auto& Jij:Js) sumJij+=Jij;
                nConnections=Js.size();
                system.sumOverNetworks(0,nConnections,1,&sumJij);
            }
            
            if (getIndex()==0) {
                //compute norm of the mean of S
                mNormS=0;
                for(const auto& Sk:mBarS) {
                    mNormS+=Sk*Sk;
                }
                mNormS=sqrt(mNormS);
                //std::cout<<"|bar S|="<<mNormS<<"\n";
                
                //esperance contribution
                std::valarray<tReal> &Es=getPackedRawValues();
                //Es[0]+=mE;
                Es[mESRawIndex]+=mNormS;
                
                
                //mean of J
                sumJij/=nConnections;
                Es[mEJRawIndex]+=sumJij;
                
                //the s-th simulation has succeeded : mEs[k]+=\bar \mu_s[k]
                getPackedRawValuesNumber()++;
                //std::cout<<"stepIndex:"<<mStepIndex<<" nRawValues:"<<getPackedRawValuesNumber()<<"\n";
            }
        }
    }
 
    inline void registerStochasticOutput(const tIndex& rawValuesNumber,
                                         const std::valarray<tReal>& rawValues) {
        std::valarray<tReal>& outputs=getOutputValues();
        outputs[mMeanSIndex]=rawValues[mESRawIndex]/rawValuesNumber;//mean of |S|
        outputs[mMeanJIndex]=rawValues[mEJRawIndex]/rawValuesNumber;//mean of |S|
        outputs[mNSamplesIndex]=rawValuesNumber;
        outputs[mNSimsIndex]=rawValues[mNRawIndex];//sum of simulations of packed simulations
    }
 
    
    inline void flushStochasticOutput(std::ofstream& file) {
        std::valarray<tReal>& outputs=getOutputValues();
        for(const auto v:outputs) {
            file<<v<<"\t";
        }
        file<<"\n";
    }
  
 
 
  
};
 
#endif