Merge pull request #245 from PauloCarvalhoRJ/release_6.3

Release 6.3
PauloCarvalhoRJ · Apr 8, 2020 · 85d4a47 · 85d4a47
2 parents 7372508 + 007b39c
commit 85d4a47
Show file tree

Hide file tree

Showing 16 changed files with 671 additions and 85 deletions.
diff --git a/GammaRay.pro b/GammaRay.pro
@@ -775,7 +775,7 @@ win32 {
 # The application version
 # Don't forget to update the Util::importSettingsFromPreviousVersion() method to
 # enable the import of registry/user settings of previous versions.
-VERSION = 6.2
+VERSION = 6.3
 
 # Define a preprocessor macro so we can get the application version in application code.
 DEFINES += APP_VERSION=\\\"$$VERSION\\\"

diff --git a/README.md b/README.md
@@ -8,6 +8,7 @@ GammaRay is a graphical user interface (GUI) that automates geostatistical workf
 If you enjoyed this project, you might also enjoy GeostatsPy: https://github.com/GeostatsGuy/GeostatsPy<br>
 
 VERSION HISTORY:<br>
+&nbsp;&nbsp;&nbsp;Version 6.3   - improvements to data imputation with MCMC and parameter experiments for Automatic Variogram Fitting.<br>
 &nbsp;&nbsp;&nbsp;Version 6.2   - data imputation with Markov Chains-Monte Carlo simulation and filtering of data sets.<br>
 &nbsp;&nbsp;&nbsp;Version 6.1   - create facies transition matrices and 3D picking & probing.<br>
 <strong>Version 6.0   - Markov Chain Random Field Simulation and numerous new features and enhancements<br></strong>

diff --git a/dialogs/automaticvarfitdialog.cpp b/dialogs/automaticvarfitdialog.cpp
diff --git a/dialogs/automaticvarfitdialog.h b/dialogs/automaticvarfitdialog.h
@@ -59,12 +59,34 @@ private Q_SLOTS:
             double finTempI, double finTempF, int finTempSteps,
             double hopFactI, double hopFactF, int hopFactSteps
             );
-    void runExperimentsWithLSRS();
-    void runExperimentsWithPSO();
-    void runExperimentsWithGenetic();
+    void runExperimentsWithLSRS( const AutomaticVarFitExperimentsDialog& expParDiag );
+    void runExperimentsWithLSRS(
+            int seedI,      int seedF,      int seedSteps,
+            double nLinesI, double nLinesF, int nLinesSteps
+            );
+    void runExperimentsWithPSO(const AutomaticVarFitExperimentsDialog &expParDiag);
+    void runExperimentsWithPSO(
+            int seedI,             int seedF,             int seedSteps,
+            double nParticlesI,    double nParticlesF,    int nParticlesSteps,
+            double inertiaI,       double inertiaF,       int inertiaSteps,
+            double acceleration1I, double acceleration1F, int acceleration1Steps,
+            double acceleration2I, double acceleration2F, int acceleration2Steps
+            );
+    void runExperimentsWithGenetic(const AutomaticVarFitExperimentsDialog& expParDiag);
+    void runExperimentsWithGenetic(
+            int    seedI,         int seedF,            int seedSteps,
+            double popSizeI,      double popSizeF,      int popSizeSteps,
+            double selSizeI,      double selSizeF,      int selSizeSteps,
+            double xOverProbI,    double xOverProbF,    int xOverProbSteps,
+            double pointOfXOverI, double pointOfXOverF, int pointOfXOverSteps,
+            double mutRateI,      double mutRateF,      int mutRateSteps
+            );
 
+    // First element in each pair: curve caption.
+    // Second element in each pair: curve values.
     void showConvergenceCurves(
-            const std::vector< std::vector< double > >& curves ) const;
+            QString chartTitle,
+            const std::vector< std::pair< QString, std::vector< double > > >& curves ) const;
 };
 
 #endif // AUTOMATICVARFITDIALOG_H
diff --git a/dialogs/automaticvarfitdialog.ui b/dialogs/automaticvarfitdialog.ui
diff --git a/dialogs/mcmcdataimputationdialog.cpp b/dialogs/mcmcdataimputationdialog.cpp
@@ -54,6 +54,9 @@ MCMCDataImputationDialog::MCMCDataImputationDialog(QWidget *parent) :
     m_ftmSelector = new FileSelectorWidget( FileSelectorType::FaciesTransitionMatrices );
     ui->frmCmbFTM->layout()->addWidget( m_ftmSelector );
 
+    m_enforceFtmSelector = new FileSelectorWidget( FileSelectorType::FaciesTransitionMatrices, true );
+    ui->frmCmbEnforceFTM->layout()->addWidget( m_enforceFtmSelector );
+
     onRemakeDistributionCombos();
 }
 
@@ -93,6 +96,8 @@ void MCMCDataImputationDialog::onRunMCMC()
     mcmcSim.m_sequenceDirection = static_cast<SequenceDirection>( ui->cmbOrder->currentIndex() );
     mcmcSim.m_atVariableGroupBy = dynamic_cast<Attribute*>( m_groupByVariableSelector->getSelectedVariable() );
     mcmcSim.setNumberOfRealizations( ui->spinRealizations->value() );
+    mcmcSim.m_enforceFTM = dynamic_cast<FaciesTransitionMatrix*>( m_enforceFtmSelector->getSelectedFile() );
+    mcmcSim.m_enforceThreshold = ui->dblSpinEnforceThreshold->value();
     //----------------------------------------------------------------------------------------------------------------------------------------
 
     Application::instance()->logInfo("Commencing MCMC simulation...");

diff --git a/dialogs/mcmcdataimputationdialog.h b/dialogs/mcmcdataimputationdialog.h
@@ -29,6 +29,7 @@ class MCMCDataImputationDialog : public QDialog
     std::vector< UnivariateDistributionSelector* > m_distributionSelectors;
     VariableSelector* m_groupByVariableSelector;
     FileSelectorWidget* m_PDFSelector;
+    FileSelectorWidget* m_enforceFtmSelector;
     //----------------------------------------------
 
 private Q_SLOTS:

diff --git a/dialogs/mcmcdataimputationdialog.ui b/dialogs/mcmcdataimputationdialog.ui
@@ -372,6 +372,99 @@
           </layout>
          </widget>
         </item>
+        <item>
+         <widget class="QFrame" name="frame_4">
+          <property name="frameShape">
+           <enum>QFrame::StyledPanel</enum>
+          </property>
+          <property name="frameShadow">
+           <enum>QFrame::Raised</enum>
+          </property>
+          <layout class="QHBoxLayout" name="horizontalLayout">
+           <property name="spacing">
+            <number>0</number>
+           </property>
+           <property name="leftMargin">
+            <number>0</number>
+           </property>
+           <property name="topMargin">
+            <number>0</number>
+           </property>
+           <property name="rightMargin">
+            <number>0</number>
+           </property>
+           <property name="bottomMargin">
+            <number>0</number>
+           </property>
+           <item>
+            <widget class="QLabel" name="label_14">
+             <property name="text">
+              <string>Allow transitions only in:</string>
+             </property>
+            </widget>
+           </item>
+           <item>
+            <widget class="QFrame" name="frmCmbEnforceFTM">
+             <property name="sizePolicy">
+              <sizepolicy hsizetype="Preferred" vsizetype="Preferred">
+               <horstretch>0</horstretch>
+               <verstretch>0</verstretch>
+              </sizepolicy>
+             </property>
+             <property name="frameShape">
+              <enum>QFrame::StyledPanel</enum>
+             </property>
+             <property name="frameShadow">
+              <enum>QFrame::Raised</enum>
+             </property>
+             <layout class="QHBoxLayout" name="horizontalLayout_4">
+              <property name="spacing">
+               <number>0</number>
+              </property>
+              <property name="leftMargin">
+               <number>0</number>
+              </property>
+              <property name="topMargin">
+               <number>0</number>
+              </property>
+              <property name="rightMargin">
+               <number>0</number>
+              </property>
+              <property name="bottomMargin">
+               <number>0</number>
+              </property>
+             </layout>
+            </widget>
+           </item>
+           <item>
+            <widget class="QLabel" name="label_15">
+             <property name="sizePolicy">
+              <sizepolicy hsizetype="Minimum" vsizetype="Preferred">
+               <horstretch>0</horstretch>
+               <verstretch>0</verstretch>
+              </sizepolicy>
+             </property>
+             <property name="text">
+              <string>threshold:</string>
+             </property>
+            </widget>
+           </item>
+           <item>
+            <widget class="QDoubleSpinBox" name="dblSpinEnforceThreshold">
+             <property name="maximum">
+              <double>1.000000000000000</double>
+             </property>
+             <property name="singleStep">
+              <double>0.050000000000000</double>
+             </property>
+             <property name="value">
+              <double>0.050000000000000</double>
+             </property>
+            </widget>
+           </item>
+          </layout>
+         </widget>
+        </item>
         <item>
          <spacer name="verticalSpacer">
           <property name="orientation">

diff --git a/docs/GammaRayManual.docx b/docs/GammaRayManual.docx
diff --git a/domain/faciestransitionmatrix.cpp b/domain/faciestransitionmatrix.cpp
@@ -441,7 +441,7 @@ int FaciesTransitionMatrix::getRowIndexOfCategory(const QString &faciesName) con
     return -1;
 }
 
-int FaciesTransitionMatrix::getColumnIndexOfCategory(const QString &faciesName)
+int FaciesTransitionMatrix::getColumnIndexOfCategory(const QString &faciesName) const
 {
     int count = 0;
     for( const QString& header : m_columnHeadersFaciesNames ){
@@ -485,6 +485,24 @@ int FaciesTransitionMatrix::getUpwardNextFaciesFromCumulativeFrequency
     return -1;
 }
 
+bool FaciesTransitionMatrix::hasInexistentTransitions(const FaciesTransitionMatrix &ftmToQuery, double threshold) const
+{
+    for( int i = 0; i < m_lineHeadersFaciesNames.size(); ++i ){
+        QString fromFaciesName = m_lineHeadersFaciesNames[i];
+        for( int j = 0; j < m_columnHeadersFaciesNames.size(); ++j ){
+            QString toFaciesName = m_columnHeadersFaciesNames[j];
+            if( getUpwardTransitionProbability( i, j ) == 0.0 ){
+                int iQuery = ftmToQuery.getRowIndexOfCategory( fromFaciesName );
+                int jQuery = ftmToQuery.getColumnIndexOfCategory( toFaciesName );
+                assert( iQuery>=0 && jQuery>=0 && "FaciesTransitionMatrix::hasInexistentTransitions(): queried FTM does not have a facies of this FTM.");
+                if( ftmToQuery.getUpwardTransitionProbability( iQuery, jQuery ) > threshold )
+                    return true;
+            }
+        }
+    }
+    return false;
+}
+
 QIcon FaciesTransitionMatrix::getIcon()
 {
     return QIcon(":icons32/transmat32");

diff --git a/domain/faciestransitionmatrix.h b/domain/faciestransitionmatrix.h
@@ -242,7 +242,7 @@ class FaciesTransitionMatrix : public File
     /** Returns the index of the column whose header equals the given name.
      * Returns -1 if the name is not found.
      */
-    int getColumnIndexOfCategory( const QString& faciesName );
+    int getColumnIndexOfCategory( const QString& faciesName ) const;
 
     /** Informs whether the main diagonal is zeroed.  That is: no auto-transition is present.
      * An FTM without auto-transitions is a requirement for Embedded Markov Chains (data is not
@@ -267,6 +267,17 @@ class FaciesTransitionMatrix : public File
      */
     int getUpwardNextFaciesFromCumulativeFrequency(int fromFaciesCode, double cumulativeFrequency ) const;
 
+
+    /**
+     * Returns whether the passed FTM has transitions that do not exist in this FTM.
+     * In other words, returns true if the passed FTM has a transition frequency greater than the threshold parameter
+     * where this FTM has zero frequency. Both matrices can be different, but all facies names in
+     * this FTM must exist in the queried FTM.  The criterion used is the upward transition frequency,
+     * or the value returned by getUpwardTransitionProbability().
+     * @param threshold Values below this are considered zero.
+     */
+    bool hasInexistentTransitions( const FaciesTransitionMatrix& ftmToQuery, double threshold = 0.0 ) const;
+
     // ProjectComponent interface
 public:
     virtual QIcon getIcon();

diff --git a/geostats/mcmcdataimputation.cpp b/geostats/mcmcdataimputation.cpp
@@ -6,6 +6,8 @@
 #include "domain/categorypdf.h"
 #include "domain/univariatedistribution.h"
 #include "domain/application.h"
+#include "domain/auxiliary/faciestransitionmatrixmaker.h"
+#include "domain/project.h"
 
 #include <random>
 
@@ -51,7 +53,10 @@ bool MCMCDataImputation::run()
     }
 
     //for each realization
-    for( std::vector< std::vector<double> >& imputedDataRealization : m_imputedData ){
+    for( int iReal = 0; iReal < m_imputedData.size(); ++iReal ){
+
+        //get the empty realization dataframe
+        std::vector< std::vector<double> >& imputedDataRealization = m_imputedData[iReal];
 
         //get the data set as either grouped by some variable or as is.
         std::vector< std::vector< std::vector<double> > > dataFrame;
@@ -109,7 +114,7 @@ bool MCMCDataImputation::run()
                     bool imputing = true;
 
                     //initialize the total thickness to imput with the total Z variation of the current segment
-                    double remainingUninformedThickness = m_dataSet->getSegmentHeight( currentDataRow );
+                    double remainingUninformedThickness = std::abs( currentTailZ - currentHeadZ );
 
                     //initialize the imputed segment's head coordinate with the base coordinate of the uninformed segment
                     double x0, y0, z0;
@@ -253,6 +258,70 @@ bool MCMCDataImputation::run()
             m_lastError = "MCMCDataImputation::run(): somehow the simulation completed with an empty data set.";
             return false;
         }
+
+        //check wheather the realization has allowed transitions, if user wants this check.
+        if( m_enforceFTM ){
+
+            //Creates a new segment set object to house the imputed data.
+            SegmentSet imputed_ss( "" );
+
+            //Set the same metadata of the original data set.
+            imputed_ss.setInfoFromAnotherSegmentSet( m_dataSet );
+
+            //causes a population of child Attribute objects matching the ones from the original imput data set
+            imputed_ss.setPath( m_dataSet->getPath() );
+            imputed_ss.updateChildObjectsCollection();
+
+            //loads the original data into the imputed data set
+            imputed_ss.loadData();
+
+            //adds a new Attribute corresponding to the imputed=1/0 flag, along with an extra data column
+            imputed_ss.addEmptyDataColumn( "imputed", imputed_ss.getDataLineCount() );
+
+            //replaces original data with the imputed data
+            imputed_ss.replaceDataFrame( imputedDataRealization );
+
+            //save the realization as temporary dataset
+            QString tmpFileName = Application::instance()->getProject()->generateUniqueTmpFilePath("dat");
+            imputed_ss.setPath( tmpFileName );
+            imputed_ss.writeToFS();
+
+//            //save its metadata file
+//            imputed_ss->updateMetaDataFile();
+
+//            //causes an update to the child objects in the project tree
+//            imputed_ss->setInfoFromMetadataFile();
+
+            //crate an FTM make auxiliary object.
+            FaciesTransitionMatrixMaker<DataFile> ftmMaker( &imputed_ss, m_atVariable->getAttributeGEOEASgivenIndex()-1 );
+
+            //if data file is a point or segment set...
+            if( ! imputed_ss.isRegular() ){
+                // Index == -1 means to not group by (treat entire data set as a single sequence).
+                if( m_atVariableGroupBy )
+                    ftmMaker.setGroupByColumn( m_atVariableGroupBy->getAttributeGEOEASgivenIndex()-1 );
+                else
+                    ftmMaker.setGroupByColumn( -1 );
+            }
+
+            // Compute FTM from alternating facies in data from -Z to +Z.
+            FaciesTransitionMatrix ftmOfRealization = ftmMaker.makeSimple( DataSetOrderForFaciesString::FROM_BOTTOM_TO_TOP );
+
+            // If the FTM of the realization has illegal transitions...
+            if( m_enforceFTM->hasInexistentTransitions( ftmOfRealization, m_enforceThreshold ) ){
+                //...warn user
+                Application::instance()->logWarn("Realization have forbidden transitions.  Simulating again...");
+                //...clear realization data
+                imputedDataRealization.clear();
+                //...simulate again
+                --iReal;
+            }
+
+            // delete the temporary file
+            imputed_ss.deleteFromFS();
+
+        } //discard realizations with forbidden transitions (if user wants it)
+
     } //for each realization
 
     return true;
@@ -344,5 +413,26 @@ bool MCMCDataImputation::isOKtoRun()
         return false;
     }
 
+    if( m_enforceFTM ){
+        if( ! m_enforceFTM->isUsable() ) {
+            m_lastError = "Facies Transition Matrix to limit transitions is not usable.  It likely contains a facies name not present the associated categorical definition.";
+            return false;
+        }
+        {
+            m_enforceFTM->readFromFS();
+
+            CategoryDefinition* cdOfVariable = m_dataSet->getCategoryDefinition( m_atVariable ); //this is verified previously
+            CategoryDefinition* cdOfFTM = m_enforceFTM->getAssociatedCategoryDefinition();
+            if( ! cdOfFTM ){
+                m_lastError = "Category definition of FTM to enforce transitions not found (nullptr).";
+                return false;
+            }
+            if( cdOfFTM != cdOfVariable ){
+                m_lastError = "Category definition of input variable must be the same object as that of the FTM used to enforce transitions.";
+                return false;
+            }
+        }
+    }
+
     return true;
 }
diff --git a/geostats/mcmcdataimputation.h b/geostats/mcmcdataimputation.h
@@ -68,6 +68,12 @@ class MCMCDataImputation
      * locations exist and the user does not provide a PDF.
      */
     CategoryPDF* m_pdfForImputationWithPreviousUnavailable;
+    /** The Facies Transition Matrix use to forbid transitions. Only transitions with count/probability - m_enforceThreshold
+     * greater than zero will be allowed in the realizations.  Set nullptr to allow any transition.
+     */
+    FaciesTransitionMatrix* m_enforceFTM;
+    /** @see m_enforceFTM */
+    double m_enforceThreshold;
     /*@}*/
 
     /** Runs the algorithm.  If false is returned, the simulation failed.  Call getLastError() to get the reasons