SM_System.h

#ifndef SM_System_H
#define SM_System_H
 
//base class
#include "SM_Object.h"
 
//network
#include "SM_Network.h"
 
//material
#include "SM_Material.h"
 
//operator
#include "SM_Operator.h"
 
//array manipulating
#include <valarray>
 
class SM_System : public SM_Object {
   
    //attributes
private :
 
 
    //associations
    
 
    //network
    CORE_UniquePointer<SM_Network> mNetwork;
 
    //material of the network
    SM_Material mMaterial;
   
 
    //operator
    std::map<tString,CORE_UniquePointer<SM_Operator>> mOperators;   
   
    //initial direction of the magnetic spin moment  at time 0 of size mParticlesNumber x SM_Constants::DIM 
    SM_RealField mS0;
    
    // direction of the magnetic spin moment  at time  t of size mParticlesNumber x SM_Constants::DIM 
    SM_RealField mSt;
    
    //number of steppers
    tIndex mStepIndex;
 
    //noise rate temperature
    tReal mNoiseRateTemperature;
 
    //stochastic noise
    tReal mStochasticNoise;
    
    
protected:
    // CONSTRUCTORS 
    SM_System(void) {
        
       
        mS0.setElementsNumber(1);
        mS0={0,1,0};
        mNoiseRateTemperature=0;
        mStochasticNoise=0;
        mStepIndex=0;
    }
    
    
    // DESTRUCTORS 
    virtual ~SM_System(void) {
    }
 
 
public : 
 
     
    //MEMORY
    
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+getContentsMemorySize();
    }
 
    virtual tMemSize getContentsMemorySize() const {
        tMemSize mem=SM_Object::getContentsMemorySize();
        mem+=(mNetwork.get()==null)?0:mNetwork->getContentsMemorySize();
        mem+=mMaterial.getContentsMemorySize();
        mem+=mS0.getContentsMemorySize();
        mem+=mSt.getContentsMemorySize();
        std::for_each(mOperators.cbegin(),mOperators.cend(),
                      [&](const auto& op) {
                          mem+=op.first.size();
                          mem+=op.second->getMemorySize();
                      });
        return mem;
    }
    
    //SET & GET methods
 
    
    //Network data
    //=============
public:
    inline void setNetwork(CORE_UniquePointer<SM_Network> net) {
        mNetwork=std::move(net);
    }
    
    inline const SM_Network& getNetwork() const {
        return *mNetwork.get();
    }
    inline SM_Network& getNetwork() {
        return *mNetwork.get();
    }
 
    inline tBoolean hasNetwork() const {
        return (mNetwork.get()!=null);
    };
    
    //material
    //========
    
    inline const SM_Material& getMaterial() const {
        return mMaterial;
    }
    
    inline SM_Material& getMaterial()  {
        return mMaterial;
    }
    
   
    
    //magnetic moment data
    //====================
 
public:
    
    inline void setInitialMagneticMomentDirections(const std::valarray<tReal>& S0)  {
        mS0=S0;
        mS0.normalize();
    }
  
    inline void setInitialMagneticMomentDirections(const SM_RealField& S0)  {
        mS0=S0;
        mS0.normalize();
    }
    
    inline void setInitialMagneticMomentDirections(const tIndex& N,
                                                   const std::array<tReal,SM_Constants::DIM>& S0)  {
        mS0.setElementsNumber(N);
        mS0.initialize(S0);
        mS0.normalize();
    }
    virtual tBoolean updateInitialMagneticMomentDirections(const SM_RealField& S0)  {
      
        if (mS0.getElementsNumber()==S0.getElementsNumber()) {
            mS0=S0;
            return true;
        }
        return false;
    }
   
    inline const SM_RealField& getInitialMagneticMomentDirections()  const {
        return mS0;
    }
     
    inline SM_RealField& getInitialMagneticMomentDirections()   {
        return mS0;
    }
 
    inline const SM_RealField& getMagneticMomentDirections()  const {
        return mSt;
    }
  
    inline SM_RealField& getMagneticMomentDirections()  {
        return mSt;
    }
 
    virtual tBoolean initializeMagneticMomentDirections() {
        mStepIndex=0;
        if (mSt.getElementsNumber()==mS0.getElementsNumber()) {
            mSt=mS0;
            return true;
        }
        return false;
    }
   
    //operators
    //==========
    
    inline void addOperator(CORE_UniquePointer<SM_Operator> op) {
        if (op.get()!=null) mOperators[op->getName()]=std::move(op);
    }
    
    inline const SM_SliceOperator* getOperator(const tString& name) const {
        auto op=mOperators.find(name);
        if (op!=mOperators.end()) {
            return &(op->second->getSliceOperator());
        }
        return null;
    }
    inline SM_SliceOperator* getOperator(const tString& name)  {
        auto op=mOperators.find(name);
        if (op!=mOperators.end()) {
            return &(op->second->getSliceOperator());
        }
        return null;
    }
 
    inline void copyOperators(const SM_System& system) {
        const std::map<tString,CORE_UniquePointer<SM_Operator>>&  operators=system.getOperators();
        CORE_UniquePointer<SM_Operator> cop;
        for_each(operators.begin(),operators.end(),
                 [&](const auto& op) {
                     cop=op.second->NewInstance();
                     cop->copy(*op.second.get());
                     mOperators[op.first]=std::move(cop);
                 });
    }
    inline const std::map<tString,CORE_UniquePointer<SM_Operator>>& getOperators() const {
        return mOperators;
    }
 
    inline std::map<tString,CORE_UniquePointer<SM_Operator>>& getOperators()  {
        return mOperators;
    }
    inline tIndex getOperatorsNumber() const {
        return mOperators.size();
    }
    inline void getOperatorsName(std::vector<tString>& names) const {
        names.clear();
        for(const auto& op : mOperators) {
            names.push_back(op.first);
        }
        
    }
 
 
 
 
public:
 
 
    //noise function
    //===============
 
    
    virtual void setStochasticNoise(const tReal& eps) {
        mStochasticNoise=eps;
    }
    inline const tReal& getStochasticNoise() const {
        return mStochasticNoise;
    }
 
    virtual void setNoiseTemperature(const tReal& T) {
        mNoiseRateTemperature=T;
    }
 
    inline const tReal& getNoiseTemperature() const {
        return mNoiseRateTemperature;
    }
 
    virtual tULLInt getGeneratedRandomNumbers() const=0;
 
 
    //index of the step
    //=================
    
    inline const tIndex& getStepIndex() const {
        return mStepIndex;
    }
    inline tIndex& getStepIndex() {
        return mStepIndex;
    }
 
    //computed functions
    //==================
    
    virtual void normalize(const tBoolean& isAdimensionized) {
        
      
        //compute characteristic values of material
        getMaterial().computeCharacteristicValues(isAdimensionized);
      
        if (isAdimensionized) {
           
            //adimensionize operators
            for(auto& op:mOperators) {
                op.second->adimensionize(getMaterial().getCharacteristicField());
            }
 
        }
    }
 
    //Esperence computations
    //=======================
 
    virtual void sumOverNetworks(const int& rootNetwork,tIndex& n,const tIndex& nValues,tReal* values) const {
        
    }
    inline void computeMagneticMomentDirectionsSumOverNetworks(const int& rootNetwork,std::array<tReal,SM_Constants::DIM>& sumS)const {
        tIndex p=getNetwork().getParticlesNumber();
 
        //compute sum_p S(p) for all p in current network
        ComputeFieldSum(0,p,mSt,sumS);
 
        //sum over all networks
        sumOverNetworks(rootNetwork,p,SM_Constants::DIM,sumS.data());
    }
    inline void computeMagneticMomentDirectionsMeanOverNetworks(const int& rootNetwork,
                                                                std::array<tReal,SM_Constants::DIM>& meanS)const {
 
        //compute sum_p S(p) for all p in current network
        tIndex p=getNetwork().getParticlesNumber();
        ComputeFieldSum(0,p,mSt,meanS);
        
        //sum over all networks
        //for root netwok the mean is the mean of all netwoks for other network it is the mean of its netwok
     
        sumOverNetworks(rootNetwork,p,SM_Constants::DIM,meanS.data());
 
        //mean of S
        for(auto& meanSk:meanS) {
            meanSk/=p;
        }
    }
    
    inline void computeMagneticMomentDirectionsSum(const tInteger& start,
                                                   const tInteger& end,
                                                   std::array<tReal,SM_Constants::DIM>& sumS) const {
        ComputeFieldSum(start,end,mSt,sumS);
     }
    
    inline static void ComputeFieldSum(const tInteger& start,
                                       const tInteger& end,
                                       const SM_RealField& F,
                                       std::array<tReal,SM_Constants::DIM>& sumF) {
       
        //iterators on sumF
        tReal *iSumF=sumF.data();
        const tReal *eSumF=iSumF;eSumF+=SM_Constants::DIM;
        
        //sumF=0
        memset(iSumF,0,sizeof(tReal)*SM_Constants::DIM);
        
        //iterators on F
        if (end>start) {
            const tReal *iF=&F[0];iF+=start*SM_Constants::DIM;
            const tReal *eF=&F[0];eF+=end*SM_Constants::DIM;
            //loop on all elements in [start,end[
            while (iF!=eF) {
                iSumF=sumF.data();
                while (iSumF!=eSumF) {
                    (*iSumF)+=(*iF);
                    //itertaor at next coordinate sum
                    iF++;
                    iSumF++;
                }
            }
        }
    }
    
    
    virtual void discretize() {
 
        //get the network data
        const SM_Network& network=getNetwork();
 
        //get the material
        const SM_Material& material=getMaterial();
        
        //discretize magnetic moments direction at t
        mSt.setElementsNumber(network.getParticlesNumber()+network.getHaloParticlesNumber());
        
        //discretize the operators on the network
        for(auto& op : mOperators) {
            op.second->discretize(network,material);
        }
 
     
 
    }
 
  
    inline void updateOperatorsState(const tIndex& timeIndex,
                                     const SM_Network& network,
                                     const SM_Material& material,
                                     const SM_RealField& S) {
        std::map<tString,CORE_UniquePointer<SM_Operator>>& operators=getOperators();
        for(auto& op : operators) {
            op.second->updateState(timeIndex,network,material,S);
        }
    }
    
    //Magnetic fields methods
    //=========================
    
    inline void computeMagneticField(const tIndex& timeIndex,
                                     const SM_Network& network,
                                     const SM_Material& material,
                                     const SM_RealField& S,
                                     SM_RealField&  H) const {
        tString routineName="SM_System::computeMagneticField()";
        CORE_Profiler::StartCallRoutine(routineName);
        
        
        //H=0
        if (H.getSize()>0) {
            H=0.;
            
            //B+=Bop
            const std::map<tString,CORE_UniquePointer<SM_Operator>>& operators=getOperators();
            for(const auto& op : operators) {
                op.second->computeMagneticField(timeIndex,network,material,S,1,H);
            }
        }
        CORE_Profiler::EndCallRoutine(routineName);
    }
  
 
 
 
 
    //spin energy methods
    //====================
    
    inline tReal computeSpinEnergy(const tIndex& i,
                                   const SM_Network& network,
                                   const SM_Material& material,
                                   const SM_RealField& S) const {
 
        return computeSpinEnergy(i,0,network,material,S);
        
    }
    inline tReal computeSpinEnergy(const tIndex& i,
                                   const tIndex& timeIndex,
                                   const SM_Network& network,
                                   const SM_Material& material,
                                   const SM_RealField& S) const {
        
      
        //number of real values : s= P x D
        tReal E=0;
        const std::map<tString,CORE_UniquePointer<SM_Operator>>& operators=getOperators();
        for(const auto& op : operators) {
            E+=op.second->computeSpinEnergy(i,timeIndex,network,material,S);
        }
        return E;
    }
   
 
    //energy methods
    //===============
   
 
protected:
    inline tReal computeEnergy(const SM_Network& network,
                               const SM_Material& material,
                               const SM_RealField& S) const {
        return computeEnergy(0,network,material,S);
    }
    
    inline tReal computeEnergy(const tIndex& timeIndex,
                               const SM_Network& network,
                               const SM_Material& material,
                               const SM_RealField& S) const {
      
        //number of real values : s= P x D
        tReal E=0;
        const std::map<tString,CORE_UniquePointer<SM_Operator>>& operators=getOperators();
        for(const auto& op : operators) {
            E+=op.second->computeEnergy(timeIndex,network,material,S);
        }
        return E;
    }
 
public:
    inline const tReal& computeEnergies(std::valarray<tReal>& energies) const {
        return computeEnergies(0,energies);
    }
    
    inline const tReal& computeEnergies(const tIndex& timeIndex,std::valarray<tReal>& energies) const {
        if (energies.size()<mOperators.size()+1) energies.resize(mOperators.size()+1);
        const SM_Network& network=this->getNetwork();
        const SM_Material& material=this->getMaterial();
 
        //iterator on energy
        tReal *iE=&energies[0];
 
        //total energy
        tReal &E=*iE;
        //initialize E
        E=0;
        
        //iterator on energy of operator
        iE++;
        
        const std::map<tString,CORE_UniquePointer<SM_Operator>>& operators=getOperators();
        for(const auto& op : operators) {
            //store the energy of operator
            (*iE)=op.second->computeEnergy(timeIndex,network,material,mSt);
            //total energy contribution
            E+=(*iE);
            //iterator on energy of next operator
            iE++;
        }
        //return the total energy
        return E;
        
    }
  
    
public:
    virtual tString toString() const override {
        std::stringstream ret;
        
        //print id
        ret<<SM_Object::toString()<<"\n";
        
        //print materials
        ret<<"Material:"<<mMaterial.toString()<<"\n";
 
        //print network:
        ret<<"Network:"<<mNetwork->toString()<<"\n";
 
        //preint noise rate temperature
        ret<<"Noise Temperature:"<<getNoiseTemperature()<<"\n";
        
        //print operators
        ret<<"operators:\n";
        for(const auto& op:getOperators()) {
            ret<<"\t"<<op.first<<":"<<op.second->toString();
            ret<<"\n";
        }
        //print S(t=0)
        tUInt nPrintedValues=10;
        tDimension d;
        tIndex i,N=mS0.getElementsNumber();
        const tReal *Si_d;//iterator on S at d-coordinate of particle i
        ret<<"S0:\n";
        tInteger mod=N/nPrintedValues;//
        //mod=max(mod,1)
        mod=(mod>1)?mod:1;
        
        for (i=0;i<N;i+=mod) {
            Si_d=mS0(i);
            ret<<"\t S0["<<i<<"]=[";
            for (d=0;d<SM_Constants::DIM;d++) {
                ret<<(*Si_d)<<",";
                Si_d++;
            }
            ret.seekp(-1,ret.cur);
            ret<<"]\n";
        }
        
        
        return ret.str();
    }
 
};
 
#endif