Merge branch 'master' of dbv.inf-cv.uni-jena.de:nice/nice-core into cmake_mergebranch

core/algebra/EigValues.cpp

@@ -18,20 +18,38 @@ void
 EVArnoldi::getEigenvalues ( const GenericMatrix & data, Vector & eigenvalues,
                             Matrix & eigenvectors, uint k )
 {
+  /////////////////////////////////////  
+  ///////// check input arguments /////
+  /////////////////////////////////////
 
   if ( data.rows () != data.cols () )
   {
     throw ( "EVArnoldi: matrix has to be quadratic" );
   }
+  
   if ( k <= 0 )
   {
     throw ( "EVArnoldi: please use k>0." );
   }
+  
+  // did we specify more eigenvalues than the matrix can actually have?
+  if ( k > data.rows() )
+  {
+    throw ( "EVArnoldi: specified k is larger then dimension of matrix! Aborting..." );
+  }  
+  
+  //////////////////////////////////////
+  ///////// initialize variables ///////
+  //////////////////////////////////////
 
   if ( verbose )
-    cerr << "Initialize Matrices";
+    cerr << "Initialize Matrices" << std::endl;
 
   uint n = data.cols ();
+  
+  if ( verbose )
+    std::cerr << "EVArnoldi: n: " << n << " k: " << k << std::endl;
+  
   NICE::Matrix rmatrix ( n, k, 0 ); //=eigenvectors
   NICE::Matrix qmatrix ( n, k, 0 ); //saves data =eigenvectors
   eigenvalues.resize ( k );
@@ -52,12 +70,24 @@ EVArnoldi::getEigenvalues ( const GenericMatrix & data, Vector & eigenvalues,
 //       rmatrix ( j, i ) = 0.5;
       //TODO the random initialization might help, but it is bad for reproducibility :(
 
+  ////////////////////////////////////
+  ///////// start computation  ///////
+  ////////////////////////////////////      
+  if ( verbose )
+    std::cerr << "EVArnoldi: start main computation" << std::endl;    
+      
   //reduceddim
   double delta = 1.0;
   uint iteration = 0;
   while ( delta > mindelta && iteration < maxiterations )
   {
-    NICE::Vector rold ( rmatrix.getColumn ( k - 1 ) );
+    //replace Vector rold by matrix rold to check for convergence of all eigenvectors
+    //NICE::Vector rold ( rmatrix.getColumn ( k - 1 ) );
+    NICE::Matrix rold(rmatrix);
+    
+    if ( verbose )
+      std::cerr << "EVArnoldi: start for loop over reduced dims" << std::endl;
+    
     // meta-comment: i is an index for the iteration, j is an index for the basis
     // element (1 <= j <= k)
     for ( uint reduceddim = 0; reduceddim < k; reduceddim++ )
@@ -82,16 +112,34 @@ EVArnoldi::getEigenvalues ( const GenericMatrix & data, Vector & eigenvalues,
           ( qmatrix.getColumn ( j ).
             scalarProduct ( rmatrix.getColumn ( reduceddim ) ) );
     }
-    //convergence stuff
-    NICE::Vector diff = rold - rmatrix.getColumn ( k - 1 );
-    delta = diff.normL2 ();
+    
+    if ( verbose )
+      std::cerr << "EVArnoldi: ended for loop over reduced dims" << std::endl;
+    
+    //convergence stuff (replaced by checking all eigenvectors instead of a single one
+    //NICE::Vector diff = rold - rmatrix.getColumn ( k - 1 );
+    //delta = diff.normL2 ();
+    NICE::Vector tmpDiff;
+    double norm_tmpDiff;
+    delta = 0.0;
+    for ( uint j = 0; j < k; j++ )
+    {
+      tmpDiff = rold.getColumn(j) - rmatrix.getColumn(j);
+      norm_tmpDiff = tmpDiff.normL2();
+      if (norm_tmpDiff > delta)
+        delta = norm_tmpDiff;
+    }
     iteration++;
 
     if ( verbose )
       cerr << "EVArnoldi: [" << iteration << "] delta=" << delta << endl;
   }
+  
+  if ( verbose )
+    std::cerr << "EVArnoldi: while-loop done" << std::endl;
+  
   eigenvectors = rmatrix;
-
+  
   for ( uint i = 0; i < k; i++ )
   {
     NICE::Vector tmp;
@@ -99,4 +147,45 @@ EVArnoldi::getEigenvalues ( const GenericMatrix & data, Vector & eigenvalues,
     data.multiply ( tmp, eigenvectors.getColumn ( i ) );
     eigenvalues[i] = tmp.scalarProduct ( eigenvectors.getColumn ( i ) );
   }
+
+  
+  // post-processing: ensure that eigenvalues are in correct order!
+  
+  if ( this->b_verifyDecreasingOrder && (k > 1) )
+  {
+    NICE::VectorT< int > ewPermutation;
+    eigenvalues.sortDescend ( ewPermutation );
+    
+    NICE::Vector tmpRow;
+    int tmpIdx (-1);
+    for ( uint i = 0; i < k; i++ )
+    { 
+      if ( i == ewPermutation[i] )
+      {
+        //identity - nothing to do here
+      }
+      else
+      {
+        if ( tmpIdx == -1 )
+        {
+          tmpIdx = i;
+          tmpRow = eigenvectors.getColumn ( i );
+          eigenvectors.getColumnRef ( i ) = eigenvectors.getColumn ( ewPermutation[i] );
+        }
+        else
+        {
+          if ( tmpIdx != ewPermutation[i] )
+          {
+            eigenvectors.getColumnRef ( i ) = eigenvectors.getColumn ( ewPermutation[i] );
+          }
+          else // tmpIdx == ewPermutation[i]
+          {
+            eigenvectors.getColumnRef ( i ) = tmpRow;
+            tmpIdx = -1;
+          }
+        }
+      }
+    }
+    
+  }// sorting is only useful if we compute more then 1 ew  
 }

core/algebra/EigValues.h

@@ -1,7 +1,7 @@
 /**
 * @file EigValues.h
 * @brief Computing eigenvalues and eigenvectors
-* @author Michael Koch,Erik Rodner
+* @author Michael Koch,Erik Rodner, Alexander Freytag
 * @date 08/19/2008
 
 */
@@ -44,12 +44,14 @@ class EVArnoldi : public EigValues
     uint maxiterations;
     double mindelta;
     bool verbose;
+    bool b_verifyDecreasingOrder;
 
   public:
-    EVArnoldi ( bool verbose = false, uint _maxiterations = 100, double _mindelta = 0.01 )
+    EVArnoldi ( bool verbose = false, uint _maxiterations = 100, double _mindelta = 0.01, bool b_verifyDecreasingOrder = true)
         : maxiterations ( _maxiterations ), mindelta ( _mindelta )
     {
       this->verbose = verbose;
+      this->b_verifyDecreasingOrder = b_verifyDecreasingOrder;
     };
 
     /**

core/basics/Config.cpp

@@ -232,6 +232,7 @@ void Config::restore (istream & is, int format)
         continue;
 
       line = StringTools::chomp ( line );
+      len = line.size();
 
   #if defined DEBUGCONFIG
       DEBUGPRINT ("Config: (%d) '%s' (len = %d) / %s\n", (int)count, line.c_str(), (int)len,
@@ -274,7 +275,7 @@ void Config::restore (istream & is, int format)
       while ( (line[p]!='=') && (p<len) ) p++;
       if ( (p >= len-1) || (p<=i) ) continue;
       key = line.substr( i, p-i );
-      value = line.substr( p+1, len-p );
+      value = line.substr( p+1 ); // with only one argument, substr copies from the specified position until the end of the string
 
       StringTools::normalize_string(value);
       StringTools::normalize_string(key);

core/basics/Persistent.h

@@ -126,6 +126,42 @@ class Persistent {
     {
       save ( s, format );
     }
+    
+    virtual inline std::string createStartTag( const std::string & s_varName ) const 
+    {
+      return "<" + s_varName + ">";
+    }
+    
+    virtual inline std::string createEndTag( const std::string & s_varName ) const 
+    {
+      return "</" + s_varName + ">";
+    }
+    
+    virtual inline std::string removeStartTag( const std::string & s_varName ) const 
+    {
+      return s_varName.substr( 1, s_varName.length() -2 );
+    }
+    
+    virtual inline std::string removeEndTag( const std::string & s_varName ) const 
+    {
+      return s_varName.substr( 2, s_varName.length() -3 );
+    }
+
+    virtual inline bool isStartTag( const std::string & s_varName,  const std::string & s_matchingVarName ) const 
+    {
+      return (    ( s_varName.length() > 2 ) &&   // min length
+                  ( s_varName.substr( 0, 1 ).compare( "<" ) == 0 ) &&  // proper beginning of a start-tag
+                  ( removeStartTag(s_varName).compare(s_matchingVarName) == 0 ) // start tag matches given section name
+             );
+    }     
+    
+    virtual inline bool isEndTag( const std::string & s_varName,  const std::string & s_matchingVarName ) const 
+    {
+      return (    ( s_varName.length() > 3 ) &&   // min length
+                  ( s_varName.substr( 0, 2 ).compare( "</" ) == 0 ) &&  // proper beginning of an end-tag
+                  ( removeEndTag(s_varName).compare(s_matchingVarName) == 0 ) // end tag matches given section name
+             );
+    }    
 
 };
 

core/basics/StringTools.cpp

@@ -136,8 +136,29 @@ void StringTools::trimbounds ( std::string & value, char trimChar )
 /** @author Matti Bickel */
 std::string StringTools::trim(string s, const std::string& drop)
 {
-    std::string r=s.erase(s.find_last_not_of(drop)+1);
-    return r.erase(0,r.find_first_not_of(drop));
+    size_t findLastIndex = s.find_last_not_of(drop);
+    if ( findLastIndex < s.size() )
+    {
+      std::string r=s.erase(findLastIndex+1);
+      
+      size_t findFirstIndex = r.find_first_not_of(drop);
+      
+      if ( findFirstIndex > 0 )
+        return r.erase(0,findFirstIndex);
+      else
+        return r;
+      
+    } else
+    {
+      size_t findFirstIndex = s.find_first_not_of(drop);
+
+      if ( findFirstIndex > 0 )
+        return s.erase(0,findFirstIndex);
+      else
+        return s;
+      
+    }
+      
 }
 
 std::string StringTools::chomp(string s)

core/progs/imageDemo.cpp

@@ -26,11 +26,15 @@ int main ( int argc, char **argv )
 	ColorImage srcColor ( argv[1] );
 
   // show the image and wait for the window being closed manually
+#ifdef NICE_USELIB_GLUT
 #ifdef NICE_USELIB_QT
   showImage ( srcColor );
 #else
   cerr << "Visualization disabled: no QT library available" << endl;
 #endif
+#else 
+  cerr << "Visualization disabled: no GLUT library available" << endl;
+#endif
 
   // simple grayvalue image
   Image src;
@@ -48,9 +52,11 @@ int main ( int argc, char **argv )
   ColorImage gradientColor;
   imageToPseudoColor ( gradient, gradientColor );
 
+#ifdef NICE_USELIB_GLUT
 #ifdef NICE_USELIB_QT
   showImage ( gradientColor );
 #endif
+#endif
 
 
 }

core/vector/Algorithms.cpp

@@ -105,7 +105,21 @@ void choleskySolveMatrixLargeScale ( const Matrix & G, Matrix & B )
 	#pragma message NICE_WARNING("LinAl is not installed: choleskyInvertLargeScale will not be optimized.")
   #define CHOLESKYLINAL_WARNING
   #endif
-	choleskyInvert ( G, B );
+  if (G.rows() != G.cols())
+    fthrow(Exception, "Matrix G is not quadratic !");
+  if ( G.rows() != B.rows() )
+    fthrow(Exception, "Matrices sizes do not fit together.");
+  if (B.rows() == B.cols())
+    choleskyInvert ( G, B );
+  else
+  {
+    Vector b(B.rows(),0.0);
+    for (size_t i=0;i<B.cols();i++)
+    {
+      choleskySolveLargeScale (G, B.getColumn(i), b);
+      B.getColumnRef(i) = b;
+    }
+  }
 #endif
 }
 

core/vector/VectorT.h

@@ -390,7 +390,7 @@ public:
     inline void sortDescend();
 	 
     /**
-    * @brief sort elements in an descending order.
+    * @brief sort elements in an descending order. Permutation is only correct if all elements are different.
     */
     inline void sortDescend(VectorT<int> & permutation);
 

core/vector/VectorT.tcc

@@ -742,21 +742,32 @@ void VectorT<ElementType>::sortDescend() {
 }
 
 /**
- * @brief sort elements in an descending order.
+ * @brief sort elements in an descending order. Permutation is only correct if all elements are different.
  */
 template<typename ElementType>
 void VectorT<ElementType>::sortDescend(VectorT<int> & permutation) {
-//	VectorT<ElementType> tmp_cp(*this);
-	ippsSortDescend_I(getDataPointer(), this->dataSize);
-//	permutation.resize(this->size());
-//	for (int i = 0; i < this->size(); i++)
-//	{
-//		ElementTyp entry((*this)[i]);
-//		for (int i = 0; i < this->size(); i++)
-//		{
-//			
-//		}
-//	}
+  //copy elements to extract ordering information lateron
+  VectorT<ElementType> tmp_cp(*this);
+
+  // sort the elements
+  ippsSortDescend_I(getDataPointer(), this->dataSize);
+  
+  //compute permutation
+  permutation.resize(this->size());
+
+  int idxSelf ( 0 );
+  for (VectorT<ElementType>::const_iterator itSelf = (*this).begin(); itSelf != (*this).end(); itSelf++, idxSelf++)
+  {
+    int idxOrig ( 0 );
+    for ( VectorT<ElementType>::const_iterator itOrig = tmp_cp.begin(); itOrig != tmp_cp.end(); itOrig++, idxOrig++)
+    {
+      if ( *itOrig == *itSelf )
+      {
+        permutation[idxOrig] = idxSelf;
+        break;
+      }
+    }
+  }
 }
 
 template<typename ElementType>