SM_StochasticCurieTemperature.h

#ifndef SM_StochasticCurieTemperature_H
#define SM_StochasticCurieTemperature_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_StochasticCurieTemperature  : public SM_StochasticOutput<SM_StochasticCurieTemperature>  {
  
    // ATTRIBUTES
  
public:
   
  
   
private:
 
    //class types definition
    typedef SM_StochasticCurieTemperature Self;
    typedef SM_StochasticOutput<Self> SuperClass;
 
    tBoolean mIsLLBeam;
 
    tString mCTPackedFileName;
    tString mCTFileName;
   
   
    tIndex mStepIndex;
    std::array<tReal,SM_Constants::DIM> mBarS;//mean of S over the particles
    tReal mNormS;
 
    tUCInt mTemperatureIndex;
    tUCInt mMeanSIndex;
    tUCInt mNSamplesIndex;
    tUCInt mNSimsIndex;
public:
    // METHODS
  
    // CONSTRUCTORS 
  
    SM_StochasticCurieTemperature(void) {
      
        setOutputDescription({
                "<Temperature (K)>",
                "<Mean of |S|>",
                "<samples number>",
                "<number of simulations>",
                });
        
        mTemperatureIndex=0;
        mMeanSIndex=1;
        mNSamplesIndex=2;
        mNSimsIndex=3;
        
        std::valarray<tReal>& outputs=getOutputValues();
        outputs.resize(4);
        
        memset(&outputs[0],0,outputs.size()*sizeof(tReal));
 
    }
 
  
  
    // DESTRUCTORS 
public:
    
    virtual ~SM_StochasticCurieTemperature(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);
        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);
    }
   
    
    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<<"# 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 [1,T[, bar S_k(s,t)=frac{1}{P} sum_{p=0}^{p=P-1} S^s_k(t,p) \n"; file<<"# \t for all step t in [1,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(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_StochasticCurieTemperature::open()";
        CORE_Profiler::StartCallRoutine(routineName);
    
        
        //init the class state
        //----------------------
        
        //set step index to an unreachable step value (infinity)
        mStepIndex=core_numeric::getInfinity<tIndex>();
        
        if (getIndex()==0) {
            std::valarray<tReal>& outputs=getOutputValues();
            
            //temperature
            outputs[mTemperatureIndex]=beam.getSystem().getNoiseTemperature();
            
            //open the core & global CT Files for writting
            //--------------------------------------------
            mCTFileName=CORE_IO::GetAbsolutePath(this->getOutputPath()+"/"+this->getPrefix());
            mCTPackedFileName=mCTFileName;
            if (getPackedSimulationsNumber()>1) {
                mCTPackedFileName+="-"+std::to_string(getPackedSimulationsIndex())+"p";
            }
            mCTPackedFileName+=".ct";
            mCTFileName+=".ct";
            
            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;
                }
            }
            
        
            //init 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(2);//E(|S|),nSimulations per pack
            
            //initialize to 0
            memset(&rawValues[0],0,rawValues.size()*sizeof(tReal));
            
         
        }
        CORE_Profiler::EndCallRoutine(routineName);
        return true;
    }
    
  
  
 
    inline tBoolean open(const tIndex& s,const SM_System& system) {
        
        tString routineName="SM_StochasticCurieTemperature::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_StochasticCurieTemperature::store()";
        CORE_Profiler::StartCallRoutine(routineName);
        //compute the esperance at time t
        //std::cout<<"SM_StochasticCurieTemperature::store() Index:"<<getIndex()<<" rootIndex:"<<getRootIndex()<<"\n";
        if ((getIndex()==0)||(getRootIndex()==-1)) computeEsperanceContribution(system.getStepIndex(),system);
        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_StochasticCurieTemperature::close(s,system,hasSucceeded)";
        CORE_Profiler::StartCallRoutine(routineName);
        
        if ((getIndex()==0)||(getRootIndex()==-1)){
            if (hasSucceeded) {
                //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_StochasticCurieTemperature::close(beam)";
        CORE_Profiler::StartCallRoutine(routineName);
        
        
        if (getIndex()==0) {
            //number of simulation of the beam
            getPackedRawValues()[1]=beam.getBeamSize();
            
            if ( hasLogPerPackedSimulations() || (getPackedSimulationsNumber()==1) ) {
                //\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,getPackedRawValuesNumber(),getPackedRawValues());
                file.close();
            }
        }
        
        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_StochasticCurieTemperature::closePackedSimulations(beam,nRaws,raws)";
        CORE_Profiler::StartCallRoutine(routineName);
 
        if (getIndex()==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,rawValuesNumber,rawValues);
                file.close();
            }
        }
        CORE_Profiler::EndCallRoutine(routineName);
    }
    
    virtual tString toString() const override {
        std::stringstream ret;
        ret<<SuperClass::toString();
        return ret.str();
    }
 
    
private:
 
    inline void computeEsperanceContribution(const tIndex& stepIndex,const SM_System& system) {
 
        
        if (mStepIndex!=stepIndex) {
            
            mStepIndex=stepIndex;
 
            //compute mean of S over the networks store in index 0
            system.computeMagneticMomentDirectionsMeanOverNetworks(0,mBarS);
        
            if (getIndex()==0) {
                //compute norm of the mean of S
                mNormS=0;
                for(const auto& Sk:mBarS) {
                    mNormS+=Sk*Sk;
                }
                mNormS=sqrt(mNormS);
                
                //esperance contribution
                std::valarray<tReal> &Es=getPackedRawValues();
                Es[0]+=mNormS;
            
                
                //the s-th simulation has succeeded : mEs[k]+=\bar \mu_s[k]
                getPackedRawValuesNumber()++;
                //std::cout<<"stepIndex:"<<mStepIndex<<" nRawValues:"<<getPackedRawValuesNumber()<<"\n";
            }
        }
        
    }
 
    inline void flushStochasticOutput(std::ofstream& file,
                                      const tIndex& rawValuesNumber,
                                      const std::valarray<tReal>& rawValues) {
        std::valarray<tReal>& outputs=getOutputValues();
        outputs[mMeanSIndex]=rawValues[0]/rawValuesNumber;//mean of |S|
        outputs[mNSamplesIndex]=rawValuesNumber;
        outputs[mNSimsIndex]=rawValues[1];//sum of simulations of packed simulations
        for(const auto v:outputs) {
            file<<v<<"\t";
        }
        file<<"\n";
    }
  
 
 
  
};
 
#endif