From e3a521ced44c583e07512e5ec3528ea9da460bcf Mon Sep 17 00:00:00 2001
From: Thomas Lichtenegger <thomas.lichtenegger@jku.at>
Date: Tue, 13 Mar 2018 11:12:22 +0100
Subject: [PATCH] New path model MarkovPath is going to replace
 multiIntervalPath. Allows for multiple databases, extension of path and
 adaption of weights during runtime.

---
 src/recurrence/Make/files                     |   1 +
 src/recurrence/recBase/recBase.C              |   6 +
 src/recurrence/recBase/recBase.H              |   2 +
 .../recPath/MarkovPath/MarkovPath.C           | 274 ++++++++++++++++++
 .../recPath/MarkovPath/MarkovPath.H           | 114 ++++++++
 src/recurrence/recPath/recPath/recPath.H      |   4 +-
 6 files changed, 400 insertions(+), 1 deletion(-)
 create mode 100644 src/recurrence/recPath/MarkovPath/MarkovPath.C
 create mode 100644 src/recurrence/recPath/MarkovPath/MarkovPath.H

diff --git a/src/recurrence/Make/files b/src/recurrence/Make/files
index cacbc10d..7ad4a818 100644
--- a/src/recurrence/Make/files
+++ b/src/recurrence/Make/files
@@ -10,6 +10,7 @@ recPath/recPath/recPath.C
 recPath/recPath/newRecPath.C
 recPath/simpleRandomPath/simpleRandomPath.C
 recPath/noPath/noPath.C
+recPath/MarkovPath/MarkovPath.C
 recPath/multiIntervalPath/multiIntervalPath.C
 recStatAnalysis/recStatAnalysis/recStatAnalysis.C
 recStatAnalysis/recStatAnalysis/newRecStatAnalysis.C
diff --git a/src/recurrence/recBase/recBase.C b/src/recurrence/recBase/recBase.C
index ad9c593d..91b48f4b 100644
--- a/src/recurrence/recBase/recBase.C
+++ b/src/recurrence/recBase/recBase.C
@@ -114,6 +114,12 @@ recStatAnalysis& recBase::recStatA()
    return recStatAnalysis_(); 
 }
 
+void recBase::extendPath()
+{
+    recPath_         ->  getRecPath();
+    recModel_        ->  writeRecPath();
+}
+
 void recBase::updateRecFields()
 {
     recModel_->updateRecFields(); 
diff --git a/src/recurrence/recBase/recBase.H b/src/recurrence/recBase/recBase.H
index 3a1d2061..4c45b087 100644
--- a/src/recurrence/recBase/recBase.H
+++ b/src/recurrence/recBase/recBase.H
@@ -85,6 +85,8 @@ public:
     recModel& recM();
     
     recStatAnalysis& recStatA();
+    
+    void extendPath();
  
     void updateRecFields();
     
diff --git a/src/recurrence/recPath/MarkovPath/MarkovPath.C b/src/recurrence/recPath/MarkovPath/MarkovPath.C
new file mode 100644
index 00000000..1c8e773f
--- /dev/null
+++ b/src/recurrence/recPath/MarkovPath/MarkovPath.C
@@ -0,0 +1,274 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+    
+    Contributing authors:
+    Thomas Lichtenegger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    CFDEMcoupling academic is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "MarkovPath.H"
+#include "Random.H"
+#include "recModel.H"
+#include "addToRunTimeSelectionTable.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(MarkovPath, 0);
+
+addToRunTimeSelectionTable
+(
+    recPath,
+    MarkovPath,
+    dictionary
+);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+MarkovPath::MarkovPath
+(
+    const dictionary& dict,
+    recBase& base
+)
+:
+    recPath(dict, base),
+    propsDict_(dict.subDict(typeName + "Props")),
+    meanIntervalSteps_(propsDict_.lookupOrDefault<label>("meanIntervalSteps",-1)),
+    numIntervals_(base.recM().numIntervals()),
+    recSteps_(0),
+    intervalSizes_(numIntervals_),
+    intervalSizesCumulative_(numIntervals_),
+    Pjump_(0.0),
+    intervalWeights_(propsDict_.lookupOrDefault<scalarList>("intervalWeights",scalarList(numIntervals_,1.0))),
+    intervalWeightsCumulative_(intervalWeights_)
+{    
+    if(meanIntervalSteps_<0)
+    {
+        // if no mean interval length for consecutive steps is specified, use 1/5 from first interval
+        meanIntervalSteps_ = (label) (0.2 * intervalSizes_[0]);
+    }
+    
+    for(int i=0;i<numIntervals_;i++)
+    {
+        intervalSizes_[i] = base.recM().numRecFields(i);
+    }
+    
+    // normalize weights and create cumulative distribution
+    weightsNormalization();
+    weightsCumulation();
+    
+    for(int i=0;i<numIntervals_;i++)
+    {
+        scalar sum1 = 0.0;
+        for(int j=0;j<=i;j++)
+        {
+            sum1 += intervalSizes_[j];
+        }
+        intervalSizesCumulative_[i] = sum1;
+    }
+    
+    // given a jump probability of P, the probability of finding a chain of length N is
+    // P(N) = (1 - P)^N * P, and the mean length E(N) = (1 - P) / P
+    Pjump_ = 1.0 / (1 + meanIntervalSteps_);
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+MarkovPath::~MarkovPath()
+{}
+
+// * * * * * * * * * * * * * protected Member Functions  * * * * * * * * * * * * //
+
+void MarkovPath::getRecPath()
+{
+    label numRecIntervals = 0;
+    
+    if (Pstream::master())
+    {
+        computeRecPath();
+        numRecIntervals=virtualTimeIndexList_.size();
+    }
+    
+    Pstream::scatter(numRecIntervals);
+    Pstream::scatter(recSteps_);
+    
+    if (not Pstream::master())
+    {
+        virtualTimeIndexList_.setSize(numRecIntervals);
+    }
+    
+    Pstream::scatter(virtualTimeIndexList_);
+    
+    if (verbose_)
+    {
+        Info << "\nRecurrence path communicated to all processors.\n" << endl;
+    }
+}
+
+
+void MarkovPath::computeRecPath()
+{
+    if(base_.recM().totRecSteps() == 1)
+    {
+        Info<< "\nPrimitive recurrence path with one element.\n" << endl;
+        return;
+    }
+   
+    // extend path if
+    // a) it does already exist --> add a single interval
+    // b) it does not exist yet --> until number of initial steps is reached
+    if(recSteps_ > base_.recM().totRecSteps() )
+    {
+        extendPath();
+        Info << "\nExtending recurrence path done\n" << endl;
+        return;
+    }
+
+    while(recSteps_ <= base_.recM().totRecSteps() )
+    {
+        extendPath();
+    }
+    Info<< "\nComputing recurrence path done\n" << endl;
+}
+
+void MarkovPath::extendPath()
+{
+    Random ranGen(osRandomInteger());
+    
+    SymmetricSquareMatrix<scalar>& recurrenceMatrix( base_.recM().recurrenceMatrix() );
+    
+    label seqStart=0;
+    label seqEnd=0;
+    label virtualTimeIndex=0;
+    
+    // if previous intervals exist, perform a jump, otherwise start with step 0
+    if(virtualTimeIndexList_.size() > 0 )
+    {
+        virtualTimeIndex = virtualTimeIndexList_.last().second();
+
+        // jump to similar state in same or other database
+        scalar randInterval = ranGen.scalar01();
+    
+        label interval = numIntervals_-1;
+        for(int i = numIntervals_-2 ;i >= 0; i--)
+        {
+           if (randInterval < intervalWeightsCumulative_[i]) interval=i;
+        }
+
+        label startLoop = 0;
+        if (interval > 0) startLoop = intervalSizesCumulative_[interval-1];
+        label endLoop = intervalSizesCumulative_[interval] - meanIntervalSteps_;
+    
+        scalar nextMinimum(GREAT);
+        for (label j = startLoop; j <= endLoop; j++)
+        {
+            if(abs(j - virtualTimeIndex) < meanIntervalSteps_) continue;
+            if (recurrenceMatrix[j][virtualTimeIndex] < nextMinimum)
+            {
+                nextMinimum = recurrenceMatrix[j][virtualTimeIndex];
+                seqStart = j+1;
+            }
+        }
+    
+        virtualTimeIndex = seqStart;
+    }
+    
+    // take a series of consecutive steps
+    bool takeAnotherStep = true;
+    while(takeAnotherStep)
+    {
+        virtualTimeIndex++;
+
+        // interval border?
+        for(int i = 0;i < numIntervals_; i++)
+        {
+            if (intervalSizesCumulative_[i] - 1 == virtualTimeIndex) takeAnotherStep=false;
+        }
+        
+        // take another step according to jump probability?
+        scalar randJump = ranGen.scalar01();
+        if (randJump < Pjump_) takeAnotherStep=false;
+    }
+    
+    seqEnd = virtualTimeIndex;
+    
+    // add interval to recurrence path
+    labelPair seqStartEnd(seqStart,seqEnd);
+    virtualTimeIndexList_.append(seqStartEnd);
+    recSteps_ += seqEnd - seqStart + 1;
+}
+
+
+void MarkovPath::weightsCumulation()
+{
+    for(int i=0;i<numIntervals_;i++)
+    {
+        scalar sum1 = 0.0;
+        for(int j=0;j<=i;j++)
+        {
+            sum1 += intervalWeights_[j];
+        }
+        intervalWeightsCumulative_[i] = sum1;
+    }
+}
+
+
+void MarkovPath::weightsNormalization()
+{
+    scalar wsum = 0.0;
+    for(int i=0;i<numIntervals_;i++)
+    {
+        wsum += intervalWeights_[i];
+    }
+    
+    for(int i=0;i<numIntervals_;i++)
+    {
+        intervalWeights_[i] /= wsum;
+    }
+}
+
+
+void MarkovPath::updateIntervalWeights(scalarList newWeights)
+{
+    // check if number of weights is correct
+    if(newWeights.size() != numIntervals_)
+    {
+        FatalError <<"number of new weights not equal number of intervals\n" << abort(FatalError);
+    }
+    
+    for(int i=0; i<numIntervals_; i++)
+    {
+        intervalWeights_[i] = newWeights[i];
+    }
+    
+    weightsNormalization();
+    weightsCumulation();
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/recurrence/recPath/MarkovPath/MarkovPath.H b/src/recurrence/recPath/MarkovPath/MarkovPath.H
new file mode 100644
index 00000000..98337a63
--- /dev/null
+++ b/src/recurrence/recPath/MarkovPath/MarkovPath.H
@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+    
+    Contributing authors:
+    Thomas Lichtenegger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    CFDEMcoupling academic is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+    
+Description
+    A recurrence database consisting of N separate intervals is assumed with separation
+    times t0 (start time), t1, ... tN (end time) and weights w0, ... wN-1.
+\*---------------------------------------------------------------------------*/
+
+#ifndef MarkovPath_H
+#define MarkovPath_H
+
+#include "recPath.H"
+#include "scalarList.H"
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class MarkovPath Declaration
+\*---------------------------------------------------------------------------*/
+
+class MarkovPath
+:
+    public recPath
+{
+protected:
+
+    // Protected data
+
+    dictionary propsDict_;
+    
+    void extendPath();
+    
+    void getRecPath();
+
+    void computeRecPath();
+    
+    void weightsCumulation();
+    
+    void weightsNormalization();
+
+    label meanIntervalSteps_;
+    
+    label numIntervals_;
+    
+    label recSteps_;
+    
+    labelList intervalSizes_;
+    
+    labelList intervalSizesCumulative_;
+    
+    scalar Pjump_;
+    
+    scalarList intervalWeights_;
+    
+    scalarList intervalWeightsCumulative_;
+    
+
+public:
+
+    //- Runtime type information
+    TypeName("MarkovPath");
+
+    // Constructors
+
+        //- Construct from components
+        MarkovPath
+        (
+            const dictionary& dict,
+            recBase& base
+        );
+
+
+    // Destructor
+
+        virtual ~MarkovPath();
+
+
+    // Member Functions
+
+    void updateIntervalWeights(scalarList);
+  
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/recurrence/recPath/recPath/recPath.H b/src/recurrence/recPath/recPath/recPath.H
index f023e0ff..60fc0abb 100644
--- a/src/recurrence/recPath/recPath/recPath.H
+++ b/src/recurrence/recPath/recPath/recPath.H
@@ -102,8 +102,10 @@ public:
             const dictionary& dict,
             recBase& base
         );
-    
   
+
+    virtual void updateIntervalWeights(scalarList)
+    {}
 };