SM_StochasticDecreasingNoise.h

#ifndef SM_StochasticDecreasingNoise_H
#define SM_StochasticDecreasingNoise_H
 
//inhereited class header
#include "SM_StochasticOutput.h"
 
//beam class header
#include "SM_Beam.h"
 
//LL system header
#include "SM_LandauLifschitzSystem.h"
 
//IO header
#include "CORE_IO.h"
 
class SM_StochasticDecreasingNoise  : public SM_StochasticOutput<SM_StochasticDecreasingNoise>  {
  
    // ATTRIBUTES
  
public:
   
  
   
private:
    
    //type class
    typedef SM_StochasticDecreasingNoise SelfClass;
    typedef SM_StochasticOutput<SelfClass> SuperClass;
 
    //mean of S overt the particles
    std::array<tReal,SM_Constants::DIM> mBarS;
 
    //normalized solution
    std::array<tReal,SM_Constants::DIM> mSinf;
    
    //normalized times  at t
    std::valarray<tReal> mTimes;
    tReal *miTime;//iterator on time
    
    
    
    tReal* miE;//iterator on mE
    
    tString mPackedFileName;
    tString mFileName;
  
    
    //number of registered times number
    tIndex mRegisteredTimesNumber;
   
    
public:
    // METHODS
  
    // CONSTRUCTORS 
  
    SM_StochasticDecreasingNoise(void) {
       
        setOutputDescription({
                "time",
                "E(|Sinf-S|^2/eps^2)"});
    }
  
  
    // DESTRUCTORS 
    virtual ~SM_StochasticDecreasingNoise(void) {
    }
 
 
 
public:
    //MEMORY
    
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+getContentsMemorySize();
    }
    
    virtual tMemSize getContentsMemorySize() const {
        tMemSize mem=SuperClass::getContentsMemorySize();
        mem+=sizeof(tReal)*mTimes.size();
        mem+=mSinf.size()*sizeof(tReal);
        mem+=mBarS.size()*sizeof(tReal);
        return mem;
    }
 
    inline static CORE_UniquePointer<SM_StochasticDecreasingNoise> New() {
        return CORE_UniquePointer<SM_StochasticDecreasingNoise>(new SM_StochasticDecreasingNoise(),
                                                              SM_StochasticDecreasingNoise::Delete());
    }
 
    //SET & GET Methods
    //=================
    
    virtual void adimensionize(const SM_Material& material) override {
    }
    inline void setSpinDirectionsAtRelaxationTime(const std::array<tReal,SM_Constants::DIM>& Hext) {
 
        const tReal *iHext=Hext.data();
        tReal d=0;
        for(auto& Sk:mSinf) {
            //Sinf=Hext
            Sk=(*iHext);
            d+=Sk*Sk;
            
            iHext++;
        }
        //Sinf is normalized
        d=sqrt(d);
        if (d>0) {
            for(auto& Sk:mSinf) {
                Sk/=d;
            }
        }
    }
  
    inline void setSpinDirectionsAtRelaxationTime(const std::vector<tReal>& Hext) {
 
        //iterator on Hex
        auto iHext=Hext.begin();
 
        //iterator on Sinf
        tReal *iSk=mSinf.data();
 
        //end iteraton on Sinf
        const tReal* eS=mSinf.data();eS+=mSinf.size();
        
        //Sinf=0
        memset(iSk,0,sizeof(tReal)*mSinf.size());
               
        tReal d=0;
        while ((iHext!=Hext.end()) && (iSk!=eS)) {
 
            //Sinf=Hext
            (*iSk)=(*iHext);
            d+=(*iSk)*(*iSk);
            
            iHext++;
            iSk++;
        }
 
        //Sinf is normalized
        d=sqrt(d);
        if (d>0) {
            for(auto& Sk:mSinf) {
                Sk/=d;
            }
        }
        
        
    }
  
    
    //implemented templated methods
    //===================
    
    inline tBoolean open(const SM_Beam& beam) {
 
        tBoolean ok=true;
        if (getIndex()==0) {//only stochastic output index index 0 creates the data
 
            //system must be a LL system to work on this stochastic output
            const SM_LandauLifschitzSystem* llSystem=dynamic_cast< const SM_LandauLifschitzSystem*>(&beam.getSystem());
            if (llSystem==null) {
                return false;
            }
            
            //number of raws values number is null
            tIndex &nRaws=getPackedRawValuesNumber();
            nRaws=0;
            
            
            //get the number of time steps for the simultaion
            tInteger nDrawns=beam.getStepsNumber();
            //maximum number of draws 
            nDrawns/=getDrawnStepsNumber();
            nDrawns+=1;
            
            //E(t)=displaystyle \frac{1}{S} . \sum_{s=0}^{s=S-1}   \f$
            //drawn value per time
            std::valarray<tReal>& E=getPackedRawValues();
            
            //time of drawn
            mTimes.resize(nDrawns);
            memset(&mTimes[0],0,mTimes.size()*sizeof(tReal));
            
            //initialize drawn values to 0
            E.resize(mTimes.size());
            memset(&E[0],0,E.size()*sizeof(tReal));
            
            //set the outputs
            std::valarray<tReal>& outputs=getOutputValues();
            outputs.resize(mTimes.size());
            memset(&outputs[0],0,outputs.size()*sizeof(tReal));
            
            
            //open the file for writting
            mPackedFileName=CORE_IO::GetAbsolutePath(this->getOutputPath()+"/"+this->getPrefix());
            mFileName=mPackedFileName;
            if (getPackedSimulationsNumber()>1) {
                mPackedFileName+="-"+std::to_string(getPackedSimulationsIndex())+"p";
            }
            mPackedFileName+=".dn";
            mFileName+=".dn";
            if (hasLogPerPackedSimulations() || (getPackedSimulationsNumber()==1)) {
                std::ofstream file(mPackedFileName.c_str(),std::ios::out);
                if (file) {
                    printHeader(file,beam);
                    
                } else {
                    std::cout<<"Fatal Error : "<<mPackedFileName<<" can not be opened \n";
                    return false;
                }
            }
        }
        
        return ok;
    }
    
    inline tBoolean open(const tIndex& s,const SM_System& system) {
        tBoolean ok=true;
        if (getIndex()==0) {
            //initialize the pointer values for the simulation
            miE=&getPackedRawValues()[0];
            miTime=&mTimes[0];
            mRegisteredTimesNumber=0;
        }
        return ok;
    }
   
 
   
 
    inline tBoolean  store(const SM_System& system) {
      
        tBoolean ok=true;
        if ((getRootIndex()==-1) || (getIndex()==0)) {
         
            //compute mean of S over the networks
            
            //std::cout<<"store for SM_StochaticOutputDecreasingNoise."<<getIndex()<<" RootIndex:"<<getRootIndex()<<"\n";
            system.computeMagneticMomentDirectionsMeanOverNetworks(0,mBarS);
            //std::cout<<"end store for SM_StochaticOutputDecreasingNoise."<<getIndex()<<"\n";
        }         
        
        
        if (getIndex()==0) {
 
            //ll system casting
            const SM_LandauLifschitzSystem* llSystem=dynamic_cast< const SM_LandauLifschitzSystem*>(&system);
          
            
            if (mRegisteredTimesNumber>=mTimes.size()) {
                throw CORE_Exception("stochmagnet/main/stochasticOutput",
                                     "SM_StochasticDecreasingNoise::store()",
                                     "error in times dimension");
            }
            //std::cout<<"compute registered data ("<<getIndex()<<")\n";
            //compute |S-S_inf|^2
            const tReal *iSinf_k=mSinf.data();
            tReal dk,d=0;
            for(auto& Sk:mBarS) {
                dk=(Sk-(*iSinf_k));
                dk*=dk;
                d+=dk;
                
                iSinf_k++;
                
            }
 
            
            //get the noise rate at t
            tReal iepsilon2=llSystem->getNoiseRate();
            //eps:=eps^2
            iepsilon2*=iepsilon2;
            //eps:=1/eps^2
            if (iepsilon2>functions_numeric::getEpsilon<tReal>()) iepsilon2=1./iepsilon2;
            else iepsilon2=1.;
            
           
         
            //register the drawn value
            //compute \frac{1}{epsilon_t^2} . \sum_p |S_\infinity-\bar S(t)|^2
            (*miE)+=d*iepsilon2;
            
            
            //time registration of the draw
            *miTime=llSystem->getTime();
            // if ((*miTime)==llSystem->getTimeStepper().getTimeStep()) {
            //     std::cout<<" store t="<<(*miTime)<<" E="<<(*miE)<<"\n";
            // }
            
            //iterators at next drawn
            miE++;
            miTime++;
            mRegisteredTimesNumber++;
        }
        
        //std::cout<<"end SM_StochDecrNoise::store("<<getIndex()<<")\n";
        return ok;
        
    }
    
    inline tBoolean close(const tIndex& s,const SM_System& system,const tBoolean& hasSucceeded) {
        tBoolean ok=true;
 
        //register the number of drawns
        //the s-th simulation has succeeded : mEs[k]+=\bar \mu_s[k]
        if (getIndex()==0) getPackedRawValuesNumber()++;
        
        return ok;
    }
 
   
    inline const std::valarray<tReal>& getTimes() const {
        return mTimes;
    }
    inline const tIndex& getRegisteredTimesNumber() const {
        return mRegisteredTimesNumber;
    }
    inline tBoolean close(const SM_Beam& beam) {
        if (getIndex()==0) {
            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(mPackedFileName.c_str(),std::ios::app);
                flushStochasticOutput(file,getPackedRawValuesNumber(),getPackedRawValues());
                file.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(mFileName,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();
            }
        }
        
    }
 
    virtual tString toString() const override {
        std::stringstream ret;
        ret<<SuperClass::toString();
        ret<<"\t Sinf="<<functions_array::toString(mSinf)<<"\n";
        return ret.str();
    }
   
    
private:
    inline void printHeader(std::ofstream& file,const SM_Beam& beam) {
        //header is the data
        //get the system
        file<<"#Decreasing Noise  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();
        functions_string::replaceAll("\n","\n#",str);
        file<<"#"<<str<<"\n";
 
        //print strochastic output
        str="StochasticOutput:"+toString();
        functions_string::replaceAll("\n","\n#",str);
        file<<"#"<<str<<"\n";
        file<<"# bar S(s,t) = frac{1}{P} sum_p S(s,t,p) \n";
        file<<"# E(t)= frac{1}{S} . sum_{s=0}^{s=S-1}  frac{1}{varepsilon_t^2}.  frac{1}{S_{inf}-bar S(s,t)|^2} \n";
        file<<"#";
        for(const auto& d:getOutputDescription()) file<<d<<"\t";
        file<<"\n";
    }
     
  
    inline void flushStochasticOutput(std::ofstream& file,
                                      const tIndex& rawValuesNumber,
                                      const std::valarray<tReal>& rawValues) {
        std::valarray<tReal>& outputs=getOutputValues();
        const tReal *iRawValues=&rawValues[0];
        const tReal *iTimes=&mTimes[0];
        tReal *iOutputs=&outputs[0];
        for(tIndex i=0;i<mRegisteredTimesNumber;i++) {
            (*iOutputs)=(*iRawValues);
            if (rawValuesNumber>0)  (*iOutputs)/=rawValuesNumber;
            file<<(*iTimes)<<"\t"<<(*iOutputs)<<"\t"<<rawValuesNumber<<"\n";
            iRawValues++;
            iTimes++;
        }
 
        
        file<<"\n";
    }
 
    
 
 
};
 
#endif