Merge pull request #191 from PauloCarvalhoRJ/Issues_20180729

Issues 20180729

GammaRay.pro

@@ -582,7 +582,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 = 4.5
+VERSION = 4.5.1
 
 # Define a preprocessor macro so we can get the application version in application code.
 DEFINES += APP_VERSION=\\\"$$VERSION\\\"

calculator/icalcproperty.cpp

@@ -14,5 +14,7 @@ QString ICalcProperty::getScriptCompatibleName()
 	compatibleName = compatibleName.replace( ')', '_' );
     compatibleName = compatibleName.replace( '-', '_' );
     compatibleName = compatibleName.replace( '+', '_' );
+    compatibleName = compatibleName.replace( '>', '_' );
+    compatibleName = compatibleName.replace( '<', '_' );
     return compatibleName;
 }

dialogs/bidistributionmodelingdialog.cpp

@@ -89,7 +89,7 @@ void BidistributionModelingDialog::onParameters()
         QString title = "Bidistr. model for " + m_atX->getContainingFile()->getName() + ": " + m_atX->getName() + " X " + m_atY->getName();
 
         //assumes that both variables are from the same data file
-        DataFile* data_file = (DataFile*)m_atX->getContainingFile();
+        DataFile* data_file = dynamic_cast<DataFile*>(m_atX->getContainingFile());
 
         //load the data
         data_file->loadData();

dialogs/machinelearningdialog.cpp

@@ -168,8 +168,8 @@ void MachineLearningDialog::runCARTClassify()
     int maxSplitsContinuous = gpf.getParameterByName<GSLibParInt*>("MaxSplitsContinuous")->_value;
 
     //Get the selected data files.
-    DataFile* trainingDataFile = (DataFile*)m_trainingFileSelector->getSelectedFile();
-    DataFile* outputDataFile = (DataFile*)m_outputFileSelector->getSelectedFile();
+    DataFile* trainingDataFile = dynamic_cast<DataFile*>(m_trainingFileSelector->getSelectedFile());
+    DataFile* outputDataFile = dynamic_cast<DataFile*>(m_outputFileSelector->getSelectedFile());
 
     //load the data from filesystem.
     trainingDataFile->loadData();
@@ -230,8 +230,8 @@ void MachineLearningDialog::runCARTRegression()
     int maxSplitsContinuous = gpf.getParameterByName<GSLibParInt*>("MaxSplitsContinuous")->_value;
 
     //Get the selected data files.
-    DataFile* trainingDataFile = (DataFile*)m_trainingFileSelector->getSelectedFile();
-    DataFile* outputDataFile = (DataFile*)m_outputFileSelector->getSelectedFile();
+    DataFile* trainingDataFile = dynamic_cast<DataFile*>( m_trainingFileSelector->getSelectedFile() );
+    DataFile* outputDataFile = dynamic_cast<DataFile*>( m_outputFileSelector->getSelectedFile() );
 
     //load the data from filesystem.
     trainingDataFile->loadData();
@@ -346,8 +346,8 @@ void MachineLearningDialog::runRandomForestClassify()
     }
 
     //Get the selected data files.
-    DataFile* trainingDataFile = (DataFile*)m_trainingFileSelector->getSelectedFile();
-    DataFile* outputDataFile = (DataFile*)m_outputFileSelector->getSelectedFile();
+    DataFile* trainingDataFile = dynamic_cast<DataFile*>( m_trainingFileSelector->getSelectedFile() );
+    DataFile* outputDataFile = dynamic_cast<DataFile*>( m_outputFileSelector->getSelectedFile() );
 
     //load the data from filesystem.
     trainingDataFile->loadData();
@@ -444,8 +444,8 @@ void MachineLearningDialog::runRandomForestRegression()
     }
 
     //Get the selected data files.
-    DataFile* trainingDataFile = (DataFile*)m_trainingFileSelector->getSelectedFile();
-    DataFile* outputDataFile = (DataFile*)m_outputFileSelector->getSelectedFile();
+    DataFile* trainingDataFile = dynamic_cast<DataFile*>(m_trainingFileSelector->getSelectedFile());
+    DataFile* outputDataFile = dynamic_cast<DataFile*>(m_outputFileSelector->getSelectedFile());
 
     //load the data from filesystem.
     trainingDataFile->loadData();
@@ -501,7 +501,7 @@ void MachineLearningDialog::runRandomForestRegression()
 
 bool MachineLearningDialog::isClassification()
 {
-    DataFile* trainingFile = (DataFile*)m_trainingFileSelector->getSelectedFile();
+    DataFile* trainingFile = dynamic_cast<DataFile*>(m_trainingFileSelector->getSelectedFile());
     if( trainingFile ){
         Attribute* at = trainingFile->getAttributeFromGEOEASIndex(
                     m_trainingDependentVariableSelector->getSelectedVariableGEOEASIndex() );

dialogs/nscoredialog.cpp

@@ -123,7 +123,7 @@ void NScoreDialog::onHistogram()
     }
 
     //get the original data file.
-    DataFile* original_data_file = (DataFile*)m_attribute->getContainingFile();
+    DataFile* original_data_file = dynamic_cast<DataFile*>(m_attribute->getContainingFile());
 
     //determine whether it is a point set or grid
     bool isPointSet = ( original_data_file->getFileType() == "POINTSET" );
@@ -144,7 +144,7 @@ void NScoreDialog::onSave()
     }
 
     //get the original data file before n-score
-    DataFile* original_data_file = (DataFile*)m_attribute->getContainingFile();
+    DataFile* original_data_file = dynamic_cast<DataFile*>(m_attribute->getContainingFile());
 
     //determine whether the original data file is a point set or grid
     bool isPointSet = ( original_data_file->getFileType() == "POINTSET" );

dialogs/postikdialog.cpp

@@ -124,7 +124,7 @@ void PostikDialog::onConfigureAndRun()
     }
 
     //if the user set a data file to better characterize the cumulative distribution (tabulated quantiles)...
-    DataFile *cg_dataForDist = (DataFile*)m_fileForDistSelector->getSelectedFile();
+    DataFile *cg_dataForDist = dynamic_cast<DataFile*>(m_fileForDistSelector->getSelectedFile());
     if( cg_dataForDist ){
         //the optional data to better inform the cumulative distribution
         m_gpf_postik->getParameter<GSLibParFile*>(6)->_path = cg_dataForDist->getPath();

dialogs/variogramanalysisdialog.cpp

@@ -133,7 +133,7 @@ void VariogramAnalysisDialog::onOpenVarMapParameters()
         //Construct an object composition based on the parameter file template for the gamv program.
         m_gpf_varmap = new GSLibParameterFile( "varmap" );
         //Get the data file (can be either point set or grid).
-        DataFile* input_data_file = (DataFile*)m_head->getContainingFile();
+        DataFile* input_data_file = dynamic_cast<DataFile*>(m_head->getContainingFile());
         //loads data in file.
         input_data_file->loadData();
         //get the variables indexes in parent data file
@@ -697,7 +697,7 @@ void VariogramAnalysisDialog::onSaveVariogramModel()
 void VariogramAnalysisDialog::onVarNReals()
 {
     //Get the data file (can be either point set or grid).
-    DataFile* input_data_file = (DataFile*)m_head->getContainingFile();
+    DataFile* input_data_file = dynamic_cast<DataFile*>(m_head->getContainingFile());
 
     //determine whether the data is regular (grid) or irregular (point cloud)
     bool is_irregular = (input_data_file->getFileType().compare("POINTSET") == 0);

domain/cartesiangrid.cpp

@@ -380,7 +380,7 @@ void CartesianGrid::getAllVariables(std::vector<IJAbstractVariable *> &result)
     getAllObjects( all_contained_objects );
     std::vector<ProjectComponent*>::iterator it = all_contained_objects.begin();
     for(; it != all_contained_objects.end(); ++it){
-        ProjectComponent* pc = (ProjectComponent*)(*it);
+        ProjectComponent* pc = *it;
 		if( pc->isAttribute() ){
             result.push_back( dynamic_cast<Attribute*>(pc) );
         }

domain/project.cpp

@@ -582,7 +582,7 @@ void Project::freeLoadedData()
             File * fileAspect = (File*)(*it);
             //... if it's a data file ...
             if( fileAspect->isDataFile() ){
-                DataFile* dataFileAspect = (DataFile*) fileAspect;
+                DataFile* dataFileAspect = dynamic_cast<DataFile*>( fileAspect );
                 //... frees its loaded data (if any).
                 dataFileAspect->freeLoadedData();
             }

geostats/geostatsutils.cpp

@@ -188,7 +188,7 @@ MatrixNXM<double> GeostatsUtils::makeCovMatrix(std::multiset<GridCell> &samples,
         int dim = samples.size();
         for( int i = 0; i < dim; ++i ){
             covMatrix( dim, i ) = 1.0; //last row with ones
-            covMatrix( i, dim ) = 1.0; //last columns with ones
+			covMatrix( i, dim ) = 1.0; //last column with ones
         }
         covMatrix( dim, dim ) = 0.0; //last element is zero
     }
@@ -262,8 +262,8 @@ void GeostatsUtils::getValuedNeighborsTopoOrdered(GridCell &cell,
         //build the list as a set to get free ordering
         std::set<IJKDelta> *deltas = new std::set<IJKDelta>();
         for( int dk = 0; dk <= nSlicesAround/2; ++dk){
-            for( int dj = 0; dj <= nColsAround/2; ++dj ){
-                for( int di = 0; di <= nRowsAround/2; ++di){
+			for( int dj = 0; dj <= nRowsAround/2; ++dj ){
+				for( int di = 0; di <= nColsAround/2; ++di){
                     deltas->insert( IJKDelta( di, dj, dk) );
                 }
             }
@@ -278,8 +278,8 @@ void GeostatsUtils::getValuedNeighborsTopoOrdered(GridCell &cell,
         std::copy(deltas->begin(), deltas->end(), std::back_inserter(*deltasV));
         //store the list in the cache for later reuse.
         IJKDeltasCache::cache.emplace( IJKDeltasCacheKey( nColsAround,
-                                                         nRowsAround,
-                                                         nSlicesAround ), deltasV );
+														  nRowsAround,
+														  nSlicesAround ), deltasV );
         //we don't need the std::set anymore
         delete deltas;
     }
@@ -303,8 +303,8 @@ void GeostatsUtils::getValuedNeighborsTopoOrdered(GridCell &cell,
             int jj = indexes[iIndex]._j;
             int kk = indexes[iIndex]._k;
             //...if the index is within the grid limits...
-            if( ii >= 0 && ii < row_limit &&
-                jj >= 0 && jj < column_limit &&
+			if( ii >= 0 && ii < column_limit &&
+				jj >= 0 && jj < row_limit &&
                 kk >= 0 && kk < slice_limit ){
                 //...get the value corresponding to the cell index.
                 double value = cg->dataIJK( cell._dataIndex, ii, jj, kk );

geostats/matrixmxn.h

@@ -4,6 +4,8 @@
 #include <vector>
 #include <cmath>
 #include "domain/application.h"
+#include "spectral/spectral.h"
+#include "spectral/svd.h"
 
 /** A generic N lines x M columns matrix template. */
 
@@ -14,11 +16,14 @@ class MatrixNXM
     /** Constructor that initializes the matrix elements with a value. */
     MatrixNXM(unsigned int n, unsigned int m, T initValue = 0.0 );
 
+	/** Constructor that initializes the matrix from the values in a spectral::array object. */
+	MatrixNXM( const spectral::array& array );
+
     /** Returns the number of rows. */
-    T getN(){ return _n; }
+	T getN() const { return _n; }
 
     /** Returns the number of columns. */
-    T getM(){ return _m; }
+	T getM() const { return _m; }
 
     /** Operator () for l-value element access: e.g.: a(1,2) = 20.0; .
      * Due to performance concern, no range check is performed.
@@ -40,7 +45,17 @@ class MatrixNXM
     /** Inverts this matrix. It is assumed that the matrix is square, no check in this regard is performed, though
      * there is a check for singularity (prints a message and aborts calculation.)
      */
-    void invert();
+	void invertWithGaussJordan();
+
+	/** Inverts this matrix with Singular Values Decomposition (works with non-square matrices).
+	 */
+	void invertWithSVD();
+
+	/** Returns a new matrix that is a transpose of this matrix.*/
+	MatrixNXM<T> getTranspose( ) const;
+
+	/** Returns a spectral::array object equivalent to this matrix. */
+	spectral::array toSpectralArray() const;
 
 private:
     /** Number of rows. */
@@ -61,7 +76,88 @@ MatrixNXM<T>::MatrixNXM(unsigned int n, unsigned int m, T initValue ) :
 {}
 
 template <typename T>
-void MatrixNXM<T>::invert(){
+MatrixNXM<T>::MatrixNXM( const spectral::array& array ) :
+	_n( array.M() ), //yes, N is swapped with M indeed.
+	_m( array.N() ),
+	_values( _n*_m, 0.0 )
+{
+	for( int i = 0; i < _n; ++i )
+		for( int j = 0; j < _m; ++j )
+			(*this)(i,j) = array(i,j);
+}
+
+template <typename T>
+void MatrixNXM<T>::invertWithSVD(){
+	MatrixNXM<T> &a = *this;
+
+	//Make a spectral-compatible copy of this matrix.
+	spectral::array A( (spectral::index)a.getN(), (spectral::index)a.getM() );
+	for( int i = 0; i < a.getN(); ++i)
+		for( int j = 0; j < a.getM(); ++j)
+			A(i,j) = a(i,j);
+
+	//Get the U, Sigma and V* matrices from SVD on A.
+	spectral::SVD svd = spectral::svd( A );
+	spectral::array U = svd.U();
+	spectral::array Sigma = svd.S();
+	spectral::array V = svd.V(); //SVD yields V* already transposed, that is V, but to check, you must transpose V
+								 //to get A = U.Sigma.V*
+
+	//Make a full Sigma matrix (to be compatible with multiplication with the other matrices)
+	{
+		spectral::array SigmaTmp( (spectral::index)a.getN(), (spectral::index)a.getM() );
+		for( int i = 0; i < a.getN(); ++i)
+			SigmaTmp(i, i) = Sigma.d_[i];
+		Sigma = SigmaTmp;
+	}
+
+	//Make U*
+	//U contains only real numbers, thus U's transpose conjugate is its transpose.
+	spectral::array Ustar( U );
+	//transpose U to get U*
+	{
+		Eigen::MatrixXd tmp = spectral::to_2d( Ustar );
+		tmp.transposeInPlace();
+		Ustar = spectral::to_array( tmp );
+	}
+
+	//Make Sigmadagger (pseudoinverse of Sigma)
+	spectral::array SigmaDagger( Sigma );
+	{
+		//compute reciprocals of the non-zero elements in the main diagonal.
+		for( int i = 0; i < a.getN(); ++i){
+			double value = SigmaDagger(i, i);
+			//only absolute values greater than the machine epsilon are considered non-zero.
+			if( std::abs( value ) > std::numeric_limits<double>::epsilon() )
+				SigmaDagger( i, i ) = 1.0 / value;
+			else
+				SigmaDagger( i, i ) = 0.0;
+		}
+		//transpose
+		Eigen::MatrixXd tmp = spectral::to_2d( SigmaDagger );
+		tmp.transposeInPlace();
+		SigmaDagger = spectral::to_array( tmp );
+	}
+
+	//Make Adagger (pseudoinverse of A) by "reversing" the transform U.Sigma.V*,
+	//hence, Adagger = V.Sigmadagger.Ustar .
+	spectral::array Adagger;
+	{
+		Eigen::MatrixXd eigenV = spectral::to_2d( V );
+		Eigen::MatrixXd eigenSigmadagger = spectral::to_2d( SigmaDagger );
+		Eigen::MatrixXd eigenUstar = spectral::to_2d( Ustar );
+		Eigen::MatrixXd eigenAdagger = eigenV * eigenSigmadagger * eigenUstar;
+		Adagger = spectral::to_array( eigenAdagger );
+	}
+
+	//Return the result.
+	for( int i = 0; i < a.getN(); ++i)
+		for( int j = 0; j < a.getM(); ++j)
+			a(i,j) = Adagger(i,j);
+}
+
+template <typename T>
+void MatrixNXM<T>::invertWithGaussJordan(){
     /* This is a Gauss-Jordan method implemented from the code in Numerical Recipes, 3rd edition.
        It is intended to solve a linear system, but it was modified just perform invertion.
     */
@@ -117,7 +213,7 @@ void MatrixNXM<T>::invert(){
     }
 }
 
-//TODO: naive matrix multiplication, improve performance (e.g. parallel)
+//TODO: naive matrix multiplication, improve performance (e.g. parallel) or use spectral::'s classes/methods
 template <typename T>
 MatrixNXM<T> MatrixNXM<T>::operator*(const MatrixNXM<T>& b) {
    MatrixNXM<T>& a = *this;
@@ -130,4 +226,25 @@ MatrixNXM<T> MatrixNXM<T>::operator*(const MatrixNXM<T>& b) {
 }
 
 
+template <typename T>
+spectral::array MatrixNXM<T>::toSpectralArray() const{
+	const MatrixNXM<T> &a = *this;
+	//Make a spectral-compatible copy of this matrix.
+	spectral::array A( (spectral::index)a.getN(), (spectral::index)a.getM() );
+	for( int i = 0; i < a.getN(); ++i)
+		for( int j = 0; j < a.getM(); ++j)
+			A(i,j) = a(i,j);
+	return A;
+}
+
+template <typename T>
+MatrixNXM<T> MatrixNXM<T>::getTranspose( ) const {
+   spectral::array a = spectral::transpose( (*this).toSpectralArray() );
+   MatrixNXM<T> result( a.M_, a.N_);
+   for(uint i = 0; i < a.M_; ++i)
+	   for(uint j = 0; j < a.N_; ++j)
+			result(i,j) = a(i,j);
+   return result;
+}
+
 #endif // MATRIXMXN_H