From b3fd82e81bd42a44755a5ed3e13dbd8b0183ff87 Mon Sep 17 00:00:00 2001
From: Thomas Lichtenegger <k3b1677@relise.jku.austriangrid.at>
Date: Wed, 3 Nov 2021 16:49:16 +0100
Subject: [PATCH] Let MarkovPath provide more information how much time is
 spent in which database. Also allow user to specify minimal time spent in one
 database before switching.

---
 .../recPath/MarkovPath/MarkovPath.C           | 55 +++++++++++++++----
 .../recPath/MarkovPath/MarkovPath.H           | 10 ++++
 2 files changed, 55 insertions(+), 10 deletions(-)

diff --git a/src/recurrence/recPath/MarkovPath/MarkovPath.C b/src/recurrence/recPath/MarkovPath/MarkovPath.C
index 2ad4d38b..c17ff953 100644
--- a/src/recurrence/recPath/MarkovPath/MarkovPath.C
+++ b/src/recurrence/recPath/MarkovPath/MarkovPath.C
@@ -57,22 +57,29 @@ MarkovPath::MarkovPath
     recPath(dict, base),
     propsDict_(dict.subDict(typeName + "Props")),
     searchMinimum_(propsDict_.lookupOrDefault<bool>("searchMinimum",true)),
+    currentDatabase_(0),
     correlationSteps_(readLabel(propsDict_.lookup("correlationSteps"))),
     meanIntervalSteps_(propsDict_.lookupOrDefault<label>("meanIntervalSteps",-1)),
     minIntervalSteps_(propsDict_.lookupOrDefault<label>("minIntervalSteps",0)),
+    minStepsWithinDatabase_(propsDict_.lookupOrDefault<label>("minStepsWithinDatabase",0)),
     numIntervals_(base.recM().numIntervals()),
     recSteps_(0),
+    stepsInCurrentDatabase_(0),
     startIndex_(propsDict_.lookupOrDefault<label>("startIndex",0)),
     intervalSizes_(numIntervals_),
     intervalSizesCumulative_(numIntervals_),
+    numberOfIntervalsInEachDatabase_(numIntervals_),
     Pjump_(0.0),
     intervalWeights_(propsDict_.lookupOrDefault<scalarList>("intervalWeights",scalarList(numIntervals_,1.0))),
     intervalWeightsCumulative_(intervalWeights_),
+    timeSpentInEachDatabase_(intervalWeights_),
     ranGen(clock::getTime()+pid())
 {
     for(int i=0;i<numIntervals_;i++)
     {
         intervalSizes_[i] = base.recM().numRecFields(i);
+        numberOfIntervalsInEachDatabase_[i] = 0;
+        timeSpentInEachDatabase_[i] = 0.0;
     }
 
     if(meanIntervalSteps_<0)
@@ -113,7 +120,13 @@ MarkovPath::MarkovPath
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 MarkovPath::~MarkovPath()
-{}
+{
+    for(int i=0;i<numIntervals_;i++)
+    {
+        Info << "recurrence Path: total steps spent in database " << i << " = " << timeSpentInEachDatabase_[i]
+        << " (distributed over " << numberOfIntervalsInEachDatabase_[i] << " episodes)" << endl;
+    }
+}
 
 // * * * * * * * * * * * * * protected Member Functions  * * * * * * * * * * * * //
 
@@ -147,6 +160,13 @@ void MarkovPath::computeRecPath()
     {
         extendPath();
     }
+
+    for(int i=0;i<numIntervals_;i++)
+    {
+        Info << "recurrence Path: total steps spent so far in database " << i << " = " << timeSpentInEachDatabase_[i]
+        << " (distributed over " << numberOfIntervalsInEachDatabase_[i] << " episodes)" << endl;
+    }
+
     Info << "\nComputing recurrence path done\n" << endl;
 }
 
@@ -164,18 +184,29 @@ void MarkovPath::extendPath()
     {
         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 (stepsInCurrentDatabase_ >= minStepsWithinDatabase_)
         {
-           if (randInterval < intervalWeightsCumulative_[i]) interval=i;
+            label prevDatabase = currentDatabase_;
+            // 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;
+            }
+
+            currentDatabase_ = interval;
+            if (prevDatabase != currentDatabase_)
+            {
+                stepsInCurrentDatabase_ = 0;
+                numberOfIntervalsInEachDatabase_[prevDatabase] += 1;
+            }
         }
 
         label startLoop = 0;
-        if (interval > 0) startLoop = intervalSizesCumulative_[interval-1];
-        label endLoop = intervalSizesCumulative_[interval] - meanIntervalSteps_;
+        if (currentDatabase_ > 0) startLoop = intervalSizesCumulative_[currentDatabase_-1];
+        label endLoop = intervalSizesCumulative_[currentDatabase_] - meanIntervalSteps_;
 
         if (searchMinimum_)
         {
@@ -205,6 +236,7 @@ void MarkovPath::extendPath()
 
     // take a series of consecutive steps
     bool takeAnotherStep = true;
+    virtualTimeIndex--; // decrement because in the following at least one step is taken
     while(takeAnotherStep)
     {
         virtualTimeIndex++;
@@ -225,7 +257,10 @@ void MarkovPath::extendPath()
     // add interval to recurrence path
     labelPair seqStartEnd(seqStart,seqEnd);
     virtualTimeIndexList_.append(seqStartEnd);
-    recSteps_ += seqEnd - seqStart + 1;
+    label addSteps = seqEnd - seqStart + 1;
+    recSteps_ += addSteps;
+    stepsInCurrentDatabase_ += addSteps;
+    timeSpentInEachDatabase_[currentDatabase_] += addSteps;
 }
 
 
diff --git a/src/recurrence/recPath/MarkovPath/MarkovPath.H b/src/recurrence/recPath/MarkovPath/MarkovPath.H
index 64565296..153574e6 100644
--- a/src/recurrence/recPath/MarkovPath/MarkovPath.H
+++ b/src/recurrence/recPath/MarkovPath/MarkovPath.H
@@ -64,26 +64,36 @@ protected:
     // correlation steps determine which temporal neighborhood is excluded from nearest-neighbor search
     label correlationSteps_;
 
+    label currentDatabase_;
+
     label meanIntervalSteps_;
 
     label minIntervalSteps_;
 
+    label minStepsWithinDatabase_;
+
     label numIntervals_;
 
     label recSteps_;
 
+    label stepsInCurrentDatabase_;
+
     label startIndex_;
 
     labelList intervalSizes_;
 
     labelList intervalSizesCumulative_;
 
+    labelList numberOfIntervalsInEachDatabase_;
+
     scalar Pjump_;
 
     scalarList intervalWeights_;
 
     scalarList intervalWeightsCumulative_;
 
+    scalarList timeSpentInEachDatabase_;
+
     Random ranGen;