& v){
+ return to_str(v.begin(), v.end(),"map");
+}
+
+
+
+#define COMPILER_ASSERT(expr) enum { CompilerAssertAtLine##__LINE__ = sizeof( char[(expr) ? +1 : -1] ) }
+
+
+// ----------------------------
+// Platform Compatability Stuff
+// ----------------------------
+
+#ifdef WINDOWS
+typedef UINT_PTR uintptr_t;
+// typedef INT_PTR ptrdiff_t;
+#else
+int _stricmp(const char* szA, const char* szB);
+int _strnicmp(const char* szA, const char* szB, int len);
+long filelength(int filedes);
+#endif
+
+///\brief Verify that \a expected and \a got are equal for test code. Unlike Assert, this check does not disappear in optimized builds.
+///
+///If expected==got then does nothing. Otherwise prints to stderr:
+///
+///+///Test for equality failed: ---------test_descr goes here ---------------
+///
+///Expected: ------------expected goes here---------------
+///Got : ------------got goes here ---------------
+///

+///
+///Then it throws an exception using ThrowError
+///
+///Calls operator==(const T1&,const T2&) to determine equality.
+///
+///Calls GClasses::to_str to form the string representation of \a expected
+///and \a got
+///
+///\param expected The value expected from specifications
+///
+///\param got The value actually produced by the code
+///
+///\param test_descr A short test description to allow a human to
+/// easily find the failing test in the code and
+/// understand why it was written and have some help
+/// in diagnosing the bug.
+template
+void TestEqual(const T1& expected, const T2& got, std::string test_descr){
+ using std::endl;
+ if(!(expected == got)){
+ std::cerr
+ << endl
+ << "Test for equality failed: " << test_descr << endl
+ << endl
+ << "Expected: " << GClasses::to_str(expected) << endl
+ << "Got : " << GClasses::to_str(got) << endl
+ ;
+ throw Ex("Test for equality failed: ", test_descr);
+ }
+}
+
+///"Specialization" of TestEqual for c-strings done using overloading
+void TestEqual(char const* expected, char const* got, std::string desc);
+
+///"Specialization" of TestEqual for c-strings done using overloading
+void TestEqual(char const* expected, char* got, std::string desc);
+
+///"Specialization" of TestEqual for c-strings done using overloading
+void TestEqual(char* expected, char* got, std::string desc);
+
+///\brief Verify that \a expectedSubstring is a substring of \a got for test code. Unlike Assert, this check does not disappear in optimized builds.
+///
+///If \a got contains \a expectedSubstring then does nothing.
+///Otherwise prints to stderr:
+///
+///+///Substring match failed: ---------test_descr goes here ---------------
+///
+///Expected substring: ------------expectedSubstring goes here---------------
+///Got : ------------got goes here ---------------
+///

+///
+///Then it throws an exception using ThrowError
+///
+///\param expectedSubstring The value expected as a substring from
+/// specifications
+///
+///\param got The value actually produced by the code
+///
+///\param test_descr A short test description to allow a human to
+/// easily find the failing test in the code and
+/// understand why it was written and have some help
+/// in diagnosing the bug.
+void TestContains(std::string expectedSubstring, std::string got,
+ std::string desc);
+
+} // namespace GClasses
+
+#endif // __GERROR_H__
diff --git a/src/ai/waffles/GHeap.cpp b/src/ai/waffles/GHeap.cpp
new file mode 100644
index 0000000..65e18e7
--- /dev/null
+++ b/src/ai/waffles/GHeap.cpp
@@ -0,0 +1,32 @@
+/*
+ Copyright (C) 2006, Mike Gashler
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ see http://www.gnu.org/copyleft/lesser.html
+*/
+
+#include "GHeap.h"
+#include "GError.h"
+
+using namespace GClasses;
+
+// virtual
+GHeap::~GHeap()
+{
+ clear();
+}
+
+void GHeap::clear()
+{
+ while(m_pCurrentBlock)
+ {
+ char* pNext = *(char**)m_pCurrentBlock;
+ delete[] m_pCurrentBlock;
+ m_pCurrentBlock = pNext;
+ }
+ m_nCurrentPos = m_nMinBlockSize;
+}
\ No newline at end of file
diff --git a/src/ai/waffles/GHeap.h b/src/ai/waffles/GHeap.h
new file mode 100644
index 0000000..bed8370
--- /dev/null
+++ b/src/ai/waffles/GHeap.h
@@ -0,0 +1,104 @@
+/*
+ Copyright (C) 2006, Mike Gashler
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ see http://www.gnu.org/copyleft/lesser.html
+*/
+
+#ifndef __GHEAP_H__
+#define __GHEAP_H__
+
+#include
+#include
+#include "GError.h"
+
+namespace GClasses {
+
+#define BITS_PER_POINTER (sizeof(void*) * 8)
+#define ALIGN_DOWN(p) (((p) / BITS_PER_POINTER) * BITS_PER_POINTER)
+#define ALIGN_UP(p) ALIGN_DOWN((p) + BITS_PER_POINTER - 1)
+
+/// Provides a heap in which to put strings or whatever
+/// you need to store. If you need to allocate space for
+/// a lot of small objects, it's much more efficient to
+/// use this class than the C++ heap. Plus, you can
+/// delete them all by simply deleting the heap. You can't,
+/// however, reuse the space for individual objects in
+/// this heap.
+class GHeap
+{
+protected:
+ char* m_pCurrentBlock;
+ size_t m_nMinBlockSize;
+ size_t m_nCurrentPos;
+
+public:
+ GHeap(size_t nMinBlockSize)
+ {
+ m_pCurrentBlock = NULL;
+ m_nMinBlockSize = nMinBlockSize;
+ m_nCurrentPos = nMinBlockSize;
+ }
+
+ virtual ~GHeap();
+
+ /// Deletes all the blocks and frees up memory
+ void clear();
+
+ /// Allocate space in the heap and copy a string to it. Returns
+ /// a pointer to the string
+ char* add(const char* szString)
+ {
+ return add(szString, (int)strlen(szString));
+ }
+
+ /// Allocate space in the heap and copy a string to it. Returns
+ /// a pointer to the string
+ char* add(const char* pString, size_t nLength)
+ {
+ char* pNewString = allocate(nLength + 1);
+ memcpy(pNewString, pString, nLength);
+ pNewString[nLength] = '\0';
+ return pNewString;
+ }
+
+ /// Allocate space in the heap and return a pointer to it
+ char* allocate(size_t nLength)
+ {
+ if(m_nCurrentPos + nLength > m_nMinBlockSize)
+ {
+ char* pNewBlock = new char[sizeof(char*) + std::max(nLength, m_nMinBlockSize)];
+ *(char**)pNewBlock = m_pCurrentBlock;
+ m_pCurrentBlock = pNewBlock;
+ m_nCurrentPos = 0;
+ }
+ char* pNewBytes = m_pCurrentBlock + sizeof(char*) + m_nCurrentPos;
+ m_nCurrentPos += nLength;
+ return pNewBytes;
+ }
+
+ /// Allocate space in the heap and return a pointer to it
+ char* allocAligned(size_t nLength)
+ {
+ size_t nAlignedCurPos = ALIGN_UP(m_nCurrentPos);
+ if(nAlignedCurPos + nLength > m_nMinBlockSize)
+ {
+ char* pNewBlock = new char[sizeof(char*) + std::max(nLength, m_nMinBlockSize)];
+ *(char**)pNewBlock = m_pCurrentBlock;
+ m_pCurrentBlock = pNewBlock;
+ m_nCurrentPos = 0;
+ nAlignedCurPos = 0;
+ }
+ char* pNewBytes = m_pCurrentBlock + sizeof(char*) + nAlignedCurPos;
+ m_nCurrentPos = nAlignedCurPos + nLength;
+ return pNewBytes;
+ }
+};
+
+} // namespace GClasses
+
+#endif // __GHEAP_H__
diff --git a/src/ai/waffles/GHolders.cpp b/src/ai/waffles/GHolders.cpp
new file mode 100644
index 0000000..ce9c044
--- /dev/null
+++ b/src/ai/waffles/GHolders.cpp
@@ -0,0 +1,41 @@
+#include "GHolders.h"
+#include "GError.h"
+
+namespace GClasses {
+
+#ifdef _DEBUG
+GTempBufSentinel::GTempBufSentinel(void* pBuf)
+: m_pBuf(pBuf)
+{
+ *(char*)pBuf = 'S';
+}
+
+GTempBufSentinel::~GTempBufSentinel()
+{
+ GAssert(*(char*)m_pBuf == 'S'); // buffer overrun!
+}
+#endif // _DEBUG
+
+void GOverrunSentinel::Check()
+{
+ if(m_sentinel != 0x5e47143a)
+ throw Ex("buffer overrun!");
+}
+
+
+void verboten()
+{
+ throw Ex("Tried to copy a holder. (The method that facilitates this should have been private to catch this at compile time.)");
+}
+
+void FileHolder::reset(FILE* pFile)
+{
+ if(m_pFile && pFile != m_pFile)
+ {
+ if(fclose(m_pFile) != 0)
+ GAssert(false);
+ }
+ m_pFile = pFile;
+}
+
+} // namespace GClasses
diff --git a/src/ai/waffles/GHolders.h b/src/ai/waffles/GHolders.h
new file mode 100644
index 0000000..fde5618
--- /dev/null
+++ b/src/ai/waffles/GHolders.h
@@ -0,0 +1,445 @@
+/*
+ Copyright (C) 2006, Mike Gashler
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ see http://www.gnu.org/copyleft/lesser.html
+*/
+
+#ifndef __GHOLDERS_H__
+#define __GHOLDERS_H__
+
+#include
+#include
+#ifdef WINDOWS
+# include
+# pragma warning(disable: 4996)
+#else
+# ifdef __FreeBSD__
+# include
+# else
+# include
+# endif
+#endif
+#include
+
+namespace GClasses {
+
+/// The threshold over which a temporary buffer will be placed on the heap instead of the stack.
+#define MAX_STACK_TEMP_BUFFER 1024
+
+/// A helper class used by the GTEMPBUF macro
+class GTempBufHelper
+{
+public:
+ char* m_pBuf;
+ GTempBufHelper(size_t nSize)
+ {
+ m_pBuf = ((nSize > MAX_STACK_TEMP_BUFFER) ? new char[nSize] : NULL);
+ }
+
+ ~GTempBufHelper()
+ {
+ delete[] m_pBuf;
+ }
+};
+
+
+// Macro for allocating a temporary buffer
+#ifdef _DEBUG
+/// This is a helper class used by the debug version of the GTEMPBUF macro to help detect buffer overruns.
+class GTempBufSentinel
+{
+ void* m_pBuf;
+
+public:
+ GTempBufSentinel(void* pBuf);
+ ~GTempBufSentinel();
+};
+
+/// A macro for allocating a temporary buffer. If the buffer is small, it will use alloca to put it on
+/// the stack. If the buffer is big, it will allocate it on the heap and use a holder to ensure that it
+/// is properly deleted.
+# define GTEMPBUF(typ, var, cnt)\
+ GTempBufHelper var##__(sizeof(typ) * (cnt) + 1);\
+ typ* var = (((sizeof(typ) * (cnt)) <= MAX_STACK_TEMP_BUFFER) ? (typ*)alloca(sizeof(typ) * (cnt) + 1) : (typ*)var##__.m_pBuf);\
+ GTempBufSentinel var##_sentinel(&var[cnt]);
+#else
+
+/// A macro for allocating a temporary buffer. If the buffer is small, it will use alloca to put it on
+/// the stack. If the buffer is big, it will allocate it on the heap and use a holder to ensure that it
+/// is properly deleted.
+# define GTEMPBUF(typ, var, cnt)\
+ GTempBufHelper var##__(sizeof(typ) * (cnt));\
+ typ* var = (((sizeof(typ) * (cnt)) <= MAX_STACK_TEMP_BUFFER) ? (typ*)alloca(sizeof(typ) * (cnt)) : (typ*)var##__.m_pBuf);
+#endif
+
+
+void verboten();
+
+/// This class is very similar to the standard C++ class auto_ptr,
+/// except it throws an exception if you try to make a copy of it.
+/// This way, it will fail early if you use it in a manner that
+/// could result in non-deterministic behavior. (For example, if you
+/// create a vector of auto_ptrs, wierd things happen if an oom exception
+/// is thrown while resizing the buffer--part of the data will be lost
+/// when it reverts back to the original buffer. But if you make a
+/// vector of these, it will fail quickly, thus alerting you to the
+/// issue.)
+template
+class Holder
+{
+private:
+ T* m_p;
+
+public:
+ Holder(T* p = NULL)
+ {
+ //COMPILER_ASSERT(sizeof(T) > sizeof(double));
+ m_p = p;
+ }
+
+private:
+ /// Private copy constructor so that the attempts to copy a
+ /// holder are caught at compile time rather than at run time
+ Holder(const Holder& other)
+ {
+ verboten();
+ }
+public:
+ /// Deletes the object that is being held
+ ~Holder()
+ {
+ delete(m_p);
+ }
+private:
+ /// Private operator= so that the attempts to copy a
+ /// holder are caught at compile time rather than at run time
+ const Holder& operator=(const Holder& other)
+ {
+ verboten();
+ return *this;
+ }
+public:
+ /// Deletes the object that is being held, and sets the holder
+ /// to hold p. Will not delete the held pointer if the new
+ /// pointer is the same.
+ void reset(T* p = NULL)
+ {
+ if(p != m_p)
+ {
+ delete(m_p);
+ m_p = p;
+ }
+ }
+
+ /// Returns a pointer to the object being held
+ T* get()
+ {
+ return m_p;
+ }
+
+ /// Releases the object. (After calling this method, it is your job to delete the object.)
+ T* release()
+ {
+ T* pTmp = m_p;
+ m_p = NULL;
+ return pTmp;
+ }
+
+ T& operator*() const
+ {
+ return *m_p;
+ }
+
+ T* operator->() const
+ {
+ return m_p;
+ }
+};
+
+/// Just like Holder, except for arrays
+template
+class ArrayHolder
+{
+private:
+ T* m_p;
+
+public:
+ ArrayHolder(T* p = NULL)
+ {
+ m_p = p;
+ }
+private:
+ ArrayHolder(const ArrayHolder& other)
+ {
+ verboten();
+ }
+public:
+ /// Deletes the array of objects being held
+ ~ArrayHolder()
+ {
+ delete[] m_p;
+ }
+private:
+ const ArrayHolder& operator=(const ArrayHolder& other)
+ {
+ verboten();
+ return *this;
+ }
+public:
+ /// Deletes the array of objects being held and sets this holder to hold NULL
+ void reset(T* p = NULL)
+ {
+ if(p != m_p)
+ {
+ delete[] m_p;
+ m_p = p;
+ }
+ }
+
+ /// Returns a pointer to the first element of the array being held
+ T* get()
+ {
+ return m_p;
+ }
+
+ /// Releases the array. (After calling this method, it is your job to delete the array.)
+ T* release()
+ {
+ T* pTmp = m_p;
+ m_p = NULL;
+ return pTmp;
+ }
+
+ T& operator[](size_t n)
+ {
+ return m_p[n];
+ }
+};
+
+/// Closes a file when this object goes out of scope
+class FileHolder
+{
+private:
+ FILE* m_pFile;
+
+public:
+ FileHolder()
+ {
+ m_pFile = NULL;
+ }
+
+ FileHolder(FILE* pFile)
+ {
+ m_pFile = pFile;
+ }
+
+ /// Close the file
+ ~FileHolder()
+ {
+ if(m_pFile)
+ fclose(m_pFile);
+ }
+
+ /// Close the file and set this holder to hold NULL
+ void reset(FILE* pFile = NULL);
+
+ /// Returns a pointer to the FILE being held
+ FILE* get()
+ {
+ return m_pFile;
+ }
+
+ /// Releases the FILE (it is now your job to close it) and sets this holder to hold NULL
+ FILE* release()
+ {
+ FILE* pFile = m_pFile;
+ m_pFile = NULL;
+ return pFile;
+ }
+};
+
+/// Deletes all of the pointers in a vector when this object goes out of scope.
+template
+class VectorOfPointersHolder
+{
+protected:
+ std::vector& m_vec;
+
+public:
+ VectorOfPointersHolder(std::vector& vec)
+ : m_vec(vec)
+ {
+ }
+
+ /// Deletes all of the pointers in the vector
+ ~VectorOfPointersHolder()
+ {
+ for(typename std::vector::iterator it = m_vec.begin(); it != m_vec.end(); it++)
+ delete(*it);
+ }
+};
+
+/// A helper class used by the smart_ptr class.
+template
+class smart_ptr_ref_counter
+{
+public:
+ T* m_p;
+ size_t m_refCount;
+
+ smart_ptr_ref_counter(T* p)
+ : m_p(p), m_refCount(1)
+ {
+ }
+
+ ~smart_ptr_ref_counter()
+ {
+ delete(m_p);
+ }
+};
+
+/// A reference-counting smart-pointer.
+template
+class smart_ptr
+{
+protected:
+ smart_ptr_ref_counter* m_pRefCounter;
+
+ template
+ friend class smart_ptr;
+public:
+ smart_ptr()
+ : m_pRefCounter(NULL)
+ {
+ }
+
+ smart_ptr(const smart_ptr& other)
+ : m_pRefCounter(other.m_pRefCounter)
+ {
+ if(m_pRefCounter)
+ m_pRefCounter->m_refCount++;
+ }
+
+ ///\brief Make a smart pointer of a different type (T) that shares
+ ///referent and ownership with \a other
+ ///
+ ///This is mainly for creating smart_pointer from
+ ///smart_pointer.
+ ///
+ ///The basic use is:\code
+ ///smart_ptr original_ptr(new int);
+ ///int const* type_tag = NULL;
+ ///smart_ptr new_ptr(original_ptr, type_tag);
+ ///\endcode
+ ///
+ ///\warning Using this code for anything other than creating a
+ /// smart_pointer from smart_pointer is
+ /// dangerous and may do bad things due to pointer
+ /// aliasing. You do it at your own risk.
+ ///
+ ///\post get() == (T*)other.get() and when \a other goes out of
+ /// scope, it decreases the reference count for *this and
+ /// vice-versa
+ template
+ smart_ptr(const smart_ptr~~& other, T*):
+ m_pRefCounter((smart_ptr_ref_counter*)other.m_pRefCounter)
+ {
+ if(m_pRefCounter)
+ m_pRefCounter->m_refCount++;
+ }
+
+
+
+ template
+ smart_ptr(S* pThat)
+ {
+ m_pRefCounter = new smart_ptr_ref_counter(pThat);
+ }
+
+ ~smart_ptr()
+ {
+ reset();
+ }
+
+ smart_ptr& operator=(const smart_ptr& that)
+ {
+ if(that.m_pRefCounter == m_pRefCounter)
+ return *this;
+ reset();
+ if(that.m_pRefCounter)
+ {
+ m_pRefCounter = that.m_pRefCounter;
+ m_pRefCounter->m_refCount++;
+ }
+ return *this;
+ }
+
+ template
+ smart_ptr& operator=(S* pThat)
+ {
+ reset();
+ if(pThat)
+ m_pRefCounter = new smart_ptr_ref_counter(pThat);
+ return *this;
+ }
+
+ T& operator*() const
+ {
+ return *m_pRefCounter->m_p;
+ }
+
+ T* operator->() const
+ {
+ return m_pRefCounter->m_p;
+ }
+
+ T* get() const
+ {
+ return m_pRefCounter ? m_pRefCounter->m_p : NULL;
+ }
+
+ void reset()
+ {
+ if(m_pRefCounter)
+ {
+ if(--m_pRefCounter->m_refCount == 0)
+ delete(m_pRefCounter);
+ m_pRefCounter = NULL;
+ }
+ }
+
+ size_t refCount()
+ {
+ return m_pRefCounter->m_refCount;
+ }
+};
+
+/// Placing these on the stack can help catch buffer overruns
+class GOverrunSentinel
+{
+protected:
+ unsigned int m_sentinel;
+
+public:
+ GOverrunSentinel() : m_sentinel(0x5e47143a)
+ {
+ }
+
+ ~GOverrunSentinel()
+ {
+ Check();
+ }
+
+ void Check();
+};
+
+
+
+} // namespace GClasses
+
+#endif // __GHOLDERS_H__
diff --git a/src/ai/waffles/GLearner.cpp b/src/ai/waffles/GLearner.cpp
new file mode 100644
index 0000000..a5001ed
--- /dev/null
+++ b/src/ai/waffles/GLearner.cpp
@@ -0,0 +1,248 @@
+/*
+ Copyright (C) 2006, Mike Gashler
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ see http://www.gnu.org/copyleft/lesser.html
+*/
+
+#include "GLearner.h"
+#include
+#include
+#include "GError.h"
+#include "GVec.h"
+#include "GHeap.h"
+#include "GDom.h"
+#include "GNeuralNet.h"
+#include "GTransform.h"
+#include "GRand.h"
+#include
+#include
+
+using std::vector;
+
+namespace GClasses {
+
+// ---------------------------------------------------------------
+
+GTransducer::GTransducer(GRand& rand)
+: m_rand(rand)
+{
+}
+
+GTransducer::~GTransducer()
+{
+}
+
+// ---------------------------------------------------------------
+
+GSupervisedLearner::GSupervisedLearner(GRand& rand)
+: GTransducer(rand), m_pFilterFeatures(NULL), m_pFilterLabels(NULL), m_pCalibrations(NULL)
+{
+}
+
+GSupervisedLearner::GSupervisedLearner(GDomNode* pNode, GLearnerLoader& ll)
+: GTransducer(ll.rand()), m_pFilterFeatures(NULL), m_pFilterLabels(NULL)
+{
+ GDomNode* pFeatureFilter = pNode->fieldIfExists("_ff");
+ if(pFeatureFilter)
+ m_pFilterFeatures = ll.loadIncrementalTransform(pFeatureFilter);
+ GDomNode* pLabelFilter = pNode->fieldIfExists("_fl");
+ if(pLabelFilter)
+ m_pFilterLabels = ll.loadIncrementalTransform(pLabelFilter);
+ m_pRelFeatures = GRelation::deserialize(pNode->field("_rf"));
+ m_pRelLabels = GRelation::deserialize(pNode->field("_rl"));
+ m_pCalibrations = NULL;
+ GDomNode* pCalibs = pNode->fieldIfExists("cal");
+ if(pCalibs)
+ {
+ GDomListIterator it(pCalibs);
+ size_t labelDims = m_pRelLabels->size();
+ if(it.remaining() != labelDims)
+ throw Ex("The number of calibrations does not match the number of labels");
+ m_pCalibrations = new GNeuralNet*[labelDims];
+ for(size_t i = 0; i < labelDims; i++)
+ {
+ m_pCalibrations[i] = new GNeuralNet(it.current(), ll);
+ it.advance();
+ }
+ }
+}
+
+GSupervisedLearner::~GSupervisedLearner()
+{
+ if(m_pCalibrations)
+ {
+ size_t labelDims = m_pRelLabels->size();
+ for(size_t i = 0; i < labelDims; i++)
+ delete(m_pCalibrations[i]);
+ delete[] m_pCalibrations;
+ }
+ delete(m_pFilterFeatures);
+ delete(m_pFilterLabels);
+}
+
+
+// virtual
+void GSupervisedLearner::clearFeatureFilter()
+{
+ delete(m_pFilterFeatures);
+ m_pFilterFeatures = NULL;
+}
+
+void GSupervisedLearner::wrapFeatures(GIncrementalTransform* pFilter)
+{
+ if(m_pFilterFeatures)
+ m_pFilterFeatures = new GIncrementalTransformChainer(pFilter, m_pFilterFeatures);
+ else
+ m_pFilterFeatures = pFilter;
+}
+
+// virtual
+void GSupervisedLearner::clearLabelFilter()
+{
+ delete(m_pFilterLabels);
+ m_pFilterLabels = NULL;
+}
+
+
+void GSupervisedLearner::predict(const double* pIn, double* pOut)
+{
+ if(m_pFilterFeatures)
+ {
+ double* pInnerFeatures = m_pFilterFeatures->innerBuf();
+ m_pFilterFeatures->transform(pIn, pInnerFeatures);
+ if(m_pFilterLabels)
+ {
+ double* pInnerLabels = m_pFilterLabels->innerBuf();
+ predictInner(pInnerFeatures, pInnerLabels);
+ m_pFilterLabels->untransform(pInnerLabels, pOut);
+ }
+ else
+ predictInner(pInnerFeatures, pOut);
+ }
+ else
+ {
+ if(m_pFilterLabels)
+ {
+ double* pInnerLabels = m_pFilterLabels->innerBuf();
+ predictInner(pIn, pInnerLabels);
+ m_pFilterLabels->untransform(pInnerLabels, pOut);
+ }
+ else
+ predictInner(pIn, pOut);
+ }
+}
+
+
+// ---------------------------------------------------------------
+
+void GIncrementalLearner::beginIncrementalLearning(sp_relation& pFeatureRel, sp_relation& pLabelRel)
+{
+ m_pRelFeatures = pFeatureRel;
+ m_pRelLabels = pLabelRel;
+ if(m_pFilterFeatures)
+ m_pFilterFeatures->train(pFeatureRel);
+ if(m_pFilterLabels)
+ m_pFilterLabels->train(pLabelRel);
+ beginIncrementalLearningInner(m_pFilterFeatures ? m_pFilterFeatures->after() : pFeatureRel, m_pFilterLabels ? m_pFilterLabels->after() : pLabelRel);
+}
+
+void GIncrementalLearner::trainIncremental(const double* pIn, const double* pOut)
+{
+ const double* pInInner;
+ const double* pOutInner;
+ if(m_pFilterFeatures)
+ {
+ pInInner = m_pFilterFeatures->innerBuf();
+ m_pFilterFeatures->transform(pIn, m_pFilterFeatures->innerBuf());
+ }
+ else
+ pInInner = pIn;
+ if(m_pFilterLabels)
+ {
+ pOutInner = m_pFilterLabels->innerBuf();
+ m_pFilterLabels->transform(pOut, m_pFilterLabels->innerBuf());
+ }
+ else
+ pOutInner = pOut;
+ trainIncrementalInner(pInInner, pOutInner);
+}
+
+// ---------------------------------------------------------------
+
+// virtual
+GIncrementalTransform* GLearnerLoader::loadIncrementalTransform(GDomNode* pNode)
+{
+ const char* szClass = pNode->field("class")->asString();
+ if(szClass[0] == 'G')
+ {
+ if(szClass[1] < 'N')
+ {
+ if(szClass[1] < 'I')
+ {
+ if(strcmp(szClass, "GDataAugmenter") == 0)
+ return new GDataAugmenter(pNode, *this);
+ else if(strcmp(szClass, "GDiscretize") == 0)
+ return new GDiscretize(pNode, *this);
+ }
+ else
+ {
+ if(strcmp(szClass, "GIncrementalTransformChainer") == 0)
+ return new GIncrementalTransformChainer(pNode, *this);
+ }
+ }
+ else
+ {
+ if(szClass[1] < 'P')
+ {
+ if(strcmp(szClass, "GNoiseGenerator") == 0)
+ return new GNoiseGenerator(pNode, *this);
+ else if(strcmp(szClass, "GNominalToCat") == 0)
+ return new GNominalToCat(pNode, *this);
+ else if(strcmp(szClass, "GNormalize") == 0)
+ return new GNormalize(pNode, *this);
+ }
+ else
+ {
+ if(strcmp(szClass, "GPairProduct") == 0)
+ return new GPairProduct(pNode, *this);
+ else if(strcmp(szClass, "GPCA") == 0)
+ return new GPCA(pNode, *this);
+ else if(strcmp(szClass, "GReservoir") == 0)
+ return new GReservoir(pNode, *this);
+ }
+ }
+ }
+ if(m_throwIfClassNotFound)
+ throw Ex("Unrecognized class: ", szClass);
+ return NULL;
+}
+
+// virtual
+GSupervisedLearner* GLearnerLoader::loadSupervisedLearner(GDomNode* pNode)
+{
+ return loadIncrementalLearner(pNode);
+}
+
+// virtual
+GIncrementalLearner* GLearnerLoader::loadIncrementalLearner(GDomNode* pNode)
+{
+ const char* szClass = pNode->field("class")->asString();
+ if(szClass[0] == 'G')
+ {
+ if(strcmp(szClass, "GNeuralNet") == 0)
+ return new GNeuralNet(pNode, *this);
+ else if(strcmp(szClass, "GReservoirNet") == 0)
+ return new GReservoirNet(pNode, *this);
+ }
+ if(m_throwIfClassNotFound)
+ throw Ex("Unrecognized class: ", szClass);
+ return NULL;
+}
+
+
+} // namespace GClasses
diff --git a/src/ai/waffles/GLearner.h b/src/ai/waffles/GLearner.h
new file mode 100644
index 0000000..1bab3ae
--- /dev/null
+++ b/src/ai/waffles/GLearner.h
@@ -0,0 +1,362 @@
+/*
+ Copyright (C) 2006, Mike Gashler
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ see http://www.gnu.org/copyleft/lesser.html
+*/
+
+#ifndef __GLEARNER_H__
+#define __GLEARNER_H__
+
+#include "GMatrix.h"
+
+namespace GClasses {
+
+class GNeuralNet;
+
+
+class GDom;
+class GDomNode;
+class GDistribution;
+class GCategoricalDistribution;
+class GNormalDistribution;
+class GUniformDistribution;
+class GUnivariateDistribution;
+class GIncrementalTransform;
+class GSparseMatrix;
+class GCollaborativeFilter;
+class GNeuralNet;
+class GLearnerLoader;
+
+
+
+
+// nRep and nFold are zero-indexed
+typedef void (*RepValidateCallback)(void* pThis, size_t nRep, size_t nFold, size_t labelDims, double* pFoldResults);
+
+
+/// This is the base class of supervised learning algorithms (that may or may not
+/// have an internal model allowing them to generalize rows that were not available
+/// at training time). Note that the literature typically refers to supervised learning
+/// algorithms that can't generalize (because they lack an internal hypothesis model)
+/// as "Semi-supervised". (You cannot generalize with a semi-supervised algorithm--you have to
+/// train again with the new rows.)
+class GTransducer
+{
+protected:
+ GRand& m_rand;
+
+public:
+ GTransducer(GRand& rand);
+ virtual ~GTransducer();
+
+ /// Returns a reference to the random number generator associated with this object.
+ GRand& rand() { return m_rand; }
+
+protected:
+
+ /// Returns true iff this algorithm can implicitly handle nominal features. If it
+ /// cannot, then the GNominalToCat transform will be used to convert nominal
+ /// features to continuous values before passing them to it.
+ virtual bool canImplicitlyHandleNominalFeatures() { return true; }
+
+ /// Returns true iff this algorithm can implicitly handle continuous features. If it
+ /// cannot, then the GDiscretize transform will be used to convert continuous
+ /// features to nominal values before passing them to it.
+ virtual bool canImplicitlyHandleContinuousFeatures() { return true; }
+
+ /// Returns true if this algorithm supports any feature value, or if it does not
+ /// implicitly handle continuous features. If a limited range of continuous values is
+ /// supported, returns false and sets pOutMin and pOutMax to specify the range.
+ virtual bool supportedFeatureRange(double* pOutMin, double* pOutMax) { return true; }
+
+ /// Returns true iff this algorithm supports missing feature values. If it cannot,
+ /// then an imputation filter will be used to predict missing values before any
+ /// feature-vectors are passed to the algorithm.
+ virtual bool canImplicitlyHandleMissingFeatures() { return true; }
+
+ /// Returns true iff this algorithm can implicitly handle nominal labels (a.k.a.
+ /// classification). If it cannot, then the GNominalToCat transform will be
+ /// used during training to convert nominal labels to continuous values, and
+ /// to convert categorical predictions back to nominal labels.
+ virtual bool canImplicitlyHandleNominalLabels() { return true; }
+
+ /// Returns true iff this algorithm can implicitly handle continuous labels (a.k.a.
+ /// regression). If it cannot, then the GDiscretize transform will be
+ /// used during training to convert nominal labels to continuous values, and
+ /// to convert nominal predictions back to continuous labels.
+ virtual bool canImplicitlyHandleContinuousLabels() { return true; }
+
+ /// Returns true if this algorithm supports any label value, or if it does not
+ /// implicitly handle continuous labels. If a limited range of continuous values is
+ /// supported, returns false and sets pOutMin and pOutMax to specify the range.
+ virtual bool supportedLabelRange(double* pOutMin, double* pOutMax) { return true; }
+};
+
+
+/// This is the base class of algorithms that learn with supervision and
+/// have an internal hypothesis model that allows them to generalize
+/// rows that were not available at training time.
+class GSupervisedLearner : public GTransducer
+{
+protected:
+ GIncrementalTransform* m_pFilterFeatures;
+ GIncrementalTransform* m_pFilterLabels;
+ sp_relation m_pRelFeatures;
+ sp_relation m_pRelLabels;
+ GNeuralNet** m_pCalibrations;
+
+public:
+ /// General-purpose constructor
+ GSupervisedLearner(GRand& rand);
+
+ /// Deserialization constructor
+ GSupervisedLearner(GDomNode* pNode, GLearnerLoader& ll);
+
+ /// Destructor
+ virtual ~GSupervisedLearner();
+
+
+ /// Returns true because fully supervised learners have an internal
+ /// model that allows them to generalize previously unseen rows.
+ virtual bool canGeneralize() { return true; }
+
+ /// Returns a smart-pointer to the feature relation (meta-data about the input attributes).
+ /// (Note that this relation describes outer data, and may contain types that are not
+ /// supported by the inner algorithm.)
+ sp_relation relFeatures() { return m_pRelFeatures; }
+
+ /// Returns a smart-pointer to the label relation (meta-data about the output attributes).
+ /// (Note that this relation describes outer data, and may contain types that are not
+ /// supported by the inner algorithm.)
+ sp_relation relLabels() { return m_pRelLabels; }
+
+ /// Returns the current feature filter (or NULL if none has been set).
+ GIncrementalTransform* featureFilter() { return m_pFilterFeatures; }
+
+ /// Returns the current label filter (or NULL if none has been set).
+ GIncrementalTransform* labelFilter() { return m_pFilterLabels; }
+
+ /// Clears the filter for features.
+ virtual void clearFeatureFilter();
+
+ /// Wrap whatever feature filter is currently set with the specified filter
+ void wrapFeatures(GIncrementalTransform* pFilter);
+
+ /// Clears the filter for labels.
+ virtual void clearLabelFilter();
+
+
+ /// Evaluate pIn to compute a prediction for pOut. The model must be trained
+ /// (by calling train) before the first time that this method is called.
+ /// pIn and pOut should point to arrays of doubles of the same size as the
+ /// number of columns in the training matrices that were passed to the train
+ /// method.
+ void predict(const double* pIn, double* pOut);
+
+
+ /// Discards all training for the purpose of freeing memory.
+ /// If you call this method, you must train before making any predictions.
+ /// No settings or options are discarded, so you should be able to
+ /// train again without specifying any other parameters and still get
+ /// a comparable model.
+ virtual void clear() = 0;
+
+
+protected:
+ /// This is the implementation of the model's training algorithm. (This method is called by train).
+ virtual void trainInner(GMatrix& features, GMatrix& labels) = 0;
+
+ /// This is the implementation of the model's prediction algorithm. (This method is called by predict).
+ virtual void predictInner(const double* pIn, double* pOut) = 0;
+
+
+ /// This method determines which data filters (normalize, discretize,
+ /// and/or nominal-to-cat) are needed and trains them.
+ void setupFilters(GMatrix& features, GMatrix& labels);
+};
+
+///\brief Converts a GSupervisedLearner to a string
+std::string to_str(const GSupervisedLearner& learner);
+
+
+
+
+
+/// This is the base class of supervised learning algorithms that
+/// can learn one row at a time.
+class GIncrementalLearner : public GSupervisedLearner
+{
+public:
+ /// General-purpose constructor
+ GIncrementalLearner(GRand& rand)
+ : GSupervisedLearner(rand)
+ {
+ }
+
+ /// Deserialization constructor
+ GIncrementalLearner(GDomNode* pNode, GLearnerLoader& ll)
+ : GSupervisedLearner(pNode, ll)
+ {
+ }
+
+ /// Destructor
+ virtual ~GIncrementalLearner()
+ {
+ }
+
+ /// Only the GFilter class should return true to this method
+ virtual bool isFilter() { return false; }
+
+ /// Returns true
+ virtual bool canTrainIncrementally() { return true; }
+
+ /// You must call this method before you call trainIncremental.
+ /// Unlike "train", this method does not automatically set up any filters (even
+ /// if you have automatic filter setup enabled). Rather,
+ /// it assumes that you have already set up any filters that you wish to use.
+ /// Behavior is undefined if you change the filters after this method is called.
+ void beginIncrementalLearning(sp_relation& pFeatureRel, sp_relation& pLabelRel);
+
+ /// Pass a single input row and the corresponding label to
+ /// incrementally train this model
+ void trainIncremental(const double* pIn, const double* pOut);
+
+
+protected:
+ /// Prepare the model for incremental learning.
+ virtual void beginIncrementalLearningInner(sp_relation& pFeatureRel, sp_relation& pLabelRel) = 0;
+
+ /// Refine the model with the specified pattern.
+ virtual void trainIncrementalInner(const double* pIn, const double* pOut) = 0;
+};
+
+
+
+/// This class is for loading various learning algorithms from a DOM. When any
+/// learning algorithm is saved, it calls baseDomNode, which creates (among other
+/// things) a field named "class" which specifies the class name of the algorithm.
+/// This class contains methods that will recognize any of the classes in this library
+/// and load them. If it doesn't recognize a class, it will either return NULL or throw
+/// and exception, depending on the flags you pass to the constructor.
+/// Obviously this loader won't recognize any classes that you make. Therefore, you should
+/// overload the corresponding method in this class with a new method that will first
+/// recognize and load your classes, and then call these methods to handle other types.
+class GLearnerLoader
+{
+protected:
+ bool m_throwIfClassNotFound;
+ GRand& m_rand;
+
+public:
+ /// Constructor. If throwIfClassNotFound is true, then all of the methods in this
+ /// class will throw an exception of the DOM refers to an unrecognized class.
+ /// If throwIfClassNotFound is false, then NULL will be returned if the class
+ /// is not recognized.
+ GLearnerLoader(GRand& rand, bool throwIfClassNotFound = true)
+ : m_throwIfClassNotFound(throwIfClassNotFound), m_rand(rand)
+ {
+ }
+
+ virtual ~GLearnerLoader() {}
+
+ /// Loads an incremental transform (or a two-way incremental transform) from a DOM.
+ virtual GIncrementalTransform* loadIncrementalTransform(GDomNode* pNode);
+
+ /// Loads a supervised learning algorithm (or an incremental learner) from a DOM.
+ virtual GSupervisedLearner* loadSupervisedLearner(GDomNode* pNode);
+
+ /// Loads an incremental learner from a DOM.
+ virtual GIncrementalLearner* loadIncrementalLearner(GDomNode* pNode);
+
+
+ /// Returns the random number generator associated with this object.
+ GRand& rand() { return m_rand; }
+};
+
+
+
+/// Always outputs the label mean (for continuous labels) and the most common
+/// class (for nominal labels).
+class GBaselineLearner : public GSupervisedLearner
+{
+protected:
+ std::vector m_prediction;
+
+public:
+ /// General-purpose constructor
+ GBaselineLearner(GRand& rand);
+
+ /// Deserialization constructor
+ GBaselineLearner(GDomNode* pNode, GLearnerLoader& ll);
+
+ /// Destructor
+ virtual ~GBaselineLearner();
+
+
+ /// Marshal this object into a DOM, which can then be converted to a variety of serial formats.
+ virtual GDomNode* serialize(GDom* pDoc) const;
+
+ /// See the comment for GSupervisedLearner::clear
+ virtual void clear();
+
+ /// This model has no parameters to tune, so this method is a noop.
+ void autoTune(GMatrix& features, GMatrix& labels);
+
+protected:
+ /// See the comment for GSupervisedLearner::trainInner
+ virtual void trainInner(GMatrix& features, GMatrix& labels);
+
+ /// See the comment for GSupervisedLearner::predictInner
+ virtual void predictInner(const double* pIn, double* pOut);
+
+ /// See the comment for GSupervisedLearner::predictDistributionInner
+ virtual void predictDistributionInner(const double* pIn, GPrediction* pOut);
+};
+
+
+/// This is an implementation of the identity function. It might be
+/// useful, for example, as the observation function in a GRecurrentModel
+/// if you want to create a Jordan network.
+class GIdentityFunction : public GSupervisedLearner
+{
+protected:
+ size_t m_labelDims;
+ size_t m_featureDims;
+
+public:
+ /// General-purpose constructor
+ GIdentityFunction(GRand& rand);
+
+ /// Deserialization constructor
+ GIdentityFunction(GDomNode* pNode, GLearnerLoader& ll);
+
+ /// Destructor
+ virtual ~GIdentityFunction();
+
+ /// Marshal this object into a DOM, which can then be converted to a variety of serial formats.
+ virtual GDomNode* serialize(GDom* pDoc) const;
+
+ /// See the comment for GSupervisedLearner::clear
+ virtual void clear();
+
+protected:
+ /// See the comment for GSupervisedLearner::trainInner
+ virtual void trainInner(GMatrix& features, GMatrix& labels);
+
+ /// See the comment for GSupervisedLearner::predictInner
+ virtual void predictInner(const double* pIn, double* pOut);
+
+ /// See the comment for GSupervisedLearner::predictDistributionInner
+ virtual void predictDistributionInner(const double* pIn, GPrediction* pOut);
+};
+
+
+} // namespace GClasses
+
+#endif // __GLEARNER_H__
+
diff --git a/src/ai/waffles/GMatrix.cpp b/src/ai/waffles/GMatrix.cpp
new file mode 100644
index 0000000..37ef267
--- /dev/null
+++ b/src/ai/waffles/GMatrix.cpp
@@ -0,0 +1,3646 @@
+/*
+ Copyright (C) 2006, Mike Gashler
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ see http://www.gnu.org/copyleft/lesser.html
+*/
+
+#include "GMatrix.h"
+#include "GError.h"
+#include "GVec.h"
+#include "GHeap.h"
+#include "GDom.h"
+#include
+#include "GLearner.h"
+#include "GRand.h"
+#include "GTokenizer.h"
+#include
+#include
+#include
+#include
+#include
+
+using std::vector;
+using std::string;
+using std::ostream;
+using std::ostringstream;
+
+namespace GClasses {
+
+// static
+smart_ptr GRelation::deserialize(GDomNode* pNode)
+{
+ sp_relation sp;
+ if(pNode->fieldIfExists("name"))
+ sp = new GArffRelation(pNode);
+ else if(pNode->fieldIfExists("attrs"))
+ sp = new GUniformRelation(pNode);
+ else
+ sp = new GMixedRelation(pNode);
+ return sp;
+}
+
+void GRelation::print(ostream& stream, const GMatrix* pData, size_t precision) const
+{
+ stream.precision(precision);
+
+ // Write the relation title
+ stream << "@RELATION ";
+ if(type() == ARFF)
+ stream << ((GArffRelation*)this)->name();
+ else
+ stream << "Untitled";
+ stream << "\n\n";
+
+ // Write the attributes
+ for(size_t i = 0; i < size(); i++)
+ {
+ stream << "@ATTRIBUTE ";
+ printAttrName(stream, i);
+ stream << "\t";
+ if(valueCount(i) == 0)
+ stream << "real";
+ else
+ {
+ stream << "{";
+ for(size_t j = 0; j < valueCount(i); j++)
+ {
+ if(j > 0)
+ stream << ",";
+ printAttrValue(stream, i, (double)j);
+ }
+ stream << "}";
+ }
+ stream << "\n";
+ }
+
+ // Write the data
+ stream << "\n@DATA\n";
+ if(!pData)
+ return;
+ for(size_t i = 0; i < pData->rows(); i++)
+ printRow(stream, pData->row(i), ",");
+}
+
+// virtual
+void GRelation::printAttrName(std::ostream& stream, size_t column) const
+{
+ stream << "attr_" << column;
+}
+
+// virtual
+void GRelation::printAttrValue(ostream& stream, size_t column, double value) const
+{
+ size_t valCount = valueCount(column);
+ if(valCount == 0)
+ {
+ if(value == UNKNOWN_REAL_VALUE)
+ stream << "?";
+ else
+ stream << value;
+ }
+ else
+ {
+ int val = (int)value;
+ if(val < 0)
+ stream << "?";
+ else if(val >= (int)valCount)
+ throw Ex("value out of range");
+ else if(val < 26)
+ {
+ char tmp[2];
+ tmp[0] = 'a' + val;
+ tmp[1] = '\0';
+ stream << tmp;
+ }
+ else
+ stream << "_" << val;
+ }
+}
+
+// virtual
+bool GRelation::isCompatible(const GRelation& that) const
+{
+ if(this == &that)
+ return true;
+ if(size() != that.size())
+ return false;
+ for(size_t i = 0; i < size(); i++)
+ {
+ if(valueCount(i) != that.valueCount(i))
+ return false;
+ }
+ return true;
+}
+
+void GRelation::printRow(ostream& stream, const double* pRow, const char* separator) const
+{
+ size_t j = 0;
+ if(j < size())
+ {
+ printAttrValue(stream, j, *pRow);
+ pRow++;
+ j++;
+ }
+ for(; j < size(); j++)
+ {
+ stream << separator;
+ printAttrValue(stream, j, *pRow);
+ pRow++;
+ }
+ stream << "\n";
+}
+
+size_t GRelation::countRealSpaceDims(size_t nFirstAttr, size_t nAttrCount) const
+{
+ size_t nDims = 0;
+ for(size_t i = 0; i < nAttrCount; i++)
+ {
+ size_t nValues = valueCount(nFirstAttr + i);
+ if(nValues == 0)
+ nDims += 2;
+ else
+ nDims += nValues;
+ }
+ return nDims;
+}
+
+void GRelation::toRealSpace(const double* pIn, double* pOut, size_t nFirstAttr, size_t nAttrCount) const
+{
+ size_t nDims = 0;
+ size_t i, j, k, nValues;
+ for(i = 0; i < nAttrCount; i++)
+ {
+ nValues = valueCount(nFirstAttr + i);
+ if(nValues == 0)
+ {
+ pOut[nDims++] = pIn[i];
+ if(pIn[i] == UNKNOWN_REAL_VALUE)
+ pOut[nDims++] = UNKNOWN_REAL_VALUE;
+ else
+ pOut[nDims++] = pIn[i] * pIn[i];
+ }
+ else
+ {
+ k = nDims;
+ for(j = 0; j < nValues; j++)
+ pOut[nDims++] = 0;
+ if(pIn[i] >= 0) // For unknown discrete values, set to all zeros.
+ {
+ GAssert(pIn[i] >= 0 && pIn[i] < nValues);
+ pOut[k + (size_t)pIn[i]] = 1;
+ }
+ }
+ }
+}
+
+
+
+
+
+
+
+
+
+
+
+GUniformRelation::GUniformRelation(GDomNode* pNode)
+{
+ m_attrCount = (size_t)pNode->field("attrs")->asInt();
+ m_valueCount = (size_t)pNode->field("vals")->asInt();
+}
+
+// virtual
+GDomNode* GUniformRelation::serialize(GDom* pDoc) const
+{
+ GDomNode* pNode = pDoc->newObj();
+ pNode->addField(pDoc, "attrs", pDoc->newInt(m_attrCount));
+ pNode->addField(pDoc, "vals", pDoc->newInt(m_valueCount));
+ return pNode;
+}
+
+// virtual
+void GUniformRelation::deleteAttribute(size_t index)
+{
+ if(index >= m_attrCount)
+ throw Ex("Index out of range");
+ m_attrCount--;
+}
+
+// virtual
+bool GUniformRelation::isCompatible(const GRelation& that) const
+{
+ if(that.type() == GRelation::UNIFORM)
+ {
+ if(m_attrCount == ((GUniformRelation*)&that)->m_attrCount && m_valueCount == ((GUniformRelation*)&that)->m_valueCount)
+ return true;
+ else
+ return false;
+ }
+ else
+ return GRelation::isCompatible(that);
+}
+
+
+
+
+GMixedRelation::GMixedRelation()
+{
+}
+
+GMixedRelation::GMixedRelation(vector& attrValues)
+{
+ m_valueCounts.reserve(attrValues.size());
+ for(vector::iterator it = attrValues.begin(); it != attrValues.end(); it++)
+ addAttr(*it);
+}
+
+GMixedRelation::GMixedRelation(GDomNode* pNode)
+{
+ m_valueCounts.clear();
+ GDomNode* pValueCounts = pNode->field("vals");
+ GDomListIterator it(pValueCounts);
+ m_valueCounts.reserve(it.remaining());
+ for( ; it.current(); it.advance())
+ m_valueCounts.push_back((size_t)it.current()->asInt());
+}
+
+GMixedRelation::GMixedRelation(const GRelation* pCopyMe)
+{
+ copy(pCopyMe);
+}
+
+GMixedRelation::GMixedRelation(const GRelation* pCopyMe, size_t firstAttr, size_t attrCount)
+{
+ addAttrs(pCopyMe, firstAttr, attrCount);
+}
+
+// virtual
+GMixedRelation::~GMixedRelation()
+{
+}
+
+// virtual
+GDomNode* GMixedRelation::serialize(GDom* pDoc) const
+{
+ GDomNode* pNode = pDoc->newObj();
+ GDomNode* pValueCounts = pNode->addField(pDoc, "vals", pDoc->newList());
+ for(size_t i = 0; i < m_valueCounts.size(); i++)
+ pValueCounts->addItem(pDoc, pDoc->newInt(m_valueCounts[i]));
+ return pNode;
+}
+
+// virtual
+GRelation* GMixedRelation::clone() const
+{
+ GMixedRelation* pNewRelation = new GMixedRelation();
+ pNewRelation->addAttrs(this, 0, size());
+ return pNewRelation;
+}
+
+// virtual
+GRelation* GMixedRelation::cloneSub(size_t start, size_t count) const
+{
+ GMixedRelation* pNewRelation = new GMixedRelation();
+ pNewRelation->addAttrs(this, start, count);
+ return pNewRelation;
+}
+
+void GMixedRelation::addAttrs(const GRelation* pCopyMe, size_t firstAttr, size_t attrCount)
+{
+ if(firstAttr + attrCount > pCopyMe->size())
+ {
+ if(attrCount == (size_t)-1)
+ attrCount = pCopyMe->size() - firstAttr;
+ else
+ throw Ex("out of range");
+ }
+ for(size_t i = 0; i < attrCount; i++)
+ copyAttr(pCopyMe, firstAttr + i);
+}
+
+void GMixedRelation::addAttrs(size_t attrCount, size_t valueCount)
+{
+ for(size_t i = 0; i < attrCount; i++)
+ addAttr(valueCount);
+}
+
+void GMixedRelation::copy(const GRelation* pCopyMe)
+{
+ flush();
+ if(pCopyMe)
+ addAttrs(pCopyMe);
+}
+
+// virtual
+void GMixedRelation::flush()
+{
+ m_valueCounts.clear();
+}
+
+void GMixedRelation::addAttr(size_t nValues)
+{
+ if(type() == ARFF)
+ {
+ ostringstream oss;
+ oss << "attr_";
+ oss << size();
+ string s = oss.str();
+ ((GArffRelation*)this)->addAttribute(s.c_str(), nValues, NULL);
+ }
+ else
+ m_valueCounts.push_back(nValues);
+}
+
+// virtual
+void GMixedRelation::copyAttr(const GRelation* pThat, size_t nAttr)
+{
+ if(nAttr >= pThat->size())
+ throw Ex("attribute index out of range");
+ addAttr(pThat->valueCount(nAttr));
+}
+
+// virtual
+bool GMixedRelation::areContinuous(size_t first, size_t count) const
+{
+ size_t c = first;
+ for(size_t i = 0; i < count && i + c < size(); i++)
+ {
+ if(valueCount(c) != 0)
+ return false;
+ c++;
+ }
+ return true;
+}
+
+// virtual
+bool GMixedRelation::areNominal(size_t first, size_t count) const
+{
+ size_t c = first;
+ for(size_t i = 0; i < count && i + c < size(); i++)
+ {
+ if(valueCount(c) == 0)
+ return false;
+ c++;
+ }
+ return true;
+}
+
+// virtual
+void GMixedRelation::swapAttributes(size_t nAttr1, size_t nAttr2)
+{
+ std::swap(m_valueCounts[nAttr1], m_valueCounts[nAttr2]);
+}
+
+// virtual
+void GMixedRelation::deleteAttribute(size_t nAttr)
+{
+ m_valueCounts.erase(m_valueCounts.begin() + nAttr);
+}
+
+// virtual
+void GMixedRelation::setAttrValueCount(size_t nAttr, size_t nValues)
+{
+ m_valueCounts[nAttr] = nValues;
+}
+
+// ------------------------------------------------------------------
+
+GDomNode* GArffAttribute::serialize(GDom* pDoc, size_t valCount) const
+{
+ GDomNode* pNode = pDoc->newObj();
+ pNode->addField(pDoc, "name", pDoc->newString(m_name.c_str()));
+ if(valCount > 0 && m_values.size() >= valCount)
+ {
+ GDomNode* pVals = pNode->addField(pDoc, "vals", pDoc->newList());
+ for(size_t i = 0; i < valCount; i++)
+ pVals->addItem(pDoc, pDoc->newString(m_values[i].c_str()));
+ }
+ else
+ pNode->addField(pDoc, "vc", pDoc->newInt(valCount));
+ return pNode;
+}
+
+class GArffTokenizer : public GTokenizer
+{
+public:
+ GCharSet m_whitespace, m_spaces, m_space, m_valEnd, m_valEnder, m_valHardEnder, m_argEnd, m_newline, m_commaNewlineTab;
+
+ GArffTokenizer(const char* szFilename) : GTokenizer(szFilename),
+ m_whitespace("\t\n\r "), m_spaces(" \t"), m_space(" "), m_valEnd(",}\n"), m_valEnder(" ,\t}\n"), m_valHardEnder(",}\t\n"), m_argEnd(" \t\n{\r"), m_newline("\n"), m_commaNewlineTab(",\n\t") {}
+ GArffTokenizer(const char* pFile, size_t len) : GTokenizer(pFile, len),
+ m_whitespace("\t\n\r "), m_spaces(" \t"), m_space(" "), m_valEnd(",}\n"), m_valEnder(" ,\t}\n"), m_valHardEnder(",}\t\n"), m_argEnd(" \t\n{\r"), m_newline("\n"), m_commaNewlineTab(",\n\t") {}
+ virtual ~GArffTokenizer() {}
+};
+
+GArffRelation::GArffRelation()
+: m_name("Untitled")
+{
+}
+
+GArffRelation::GArffRelation(GDomNode* pNode)
+{
+ m_name = pNode->field("name")->asString();
+ GDomListIterator it(pNode->field("attrs"));
+ while(it.current())
+ {
+ GDomNode* pVals = it.current()->fieldIfExists("vals");
+ const char* valName = it.current()->field("name")->asString();
+ if(pVals)
+ {
+ vector valNames;
+ GDomListIterator it2(pVals);
+ while(it2.current())
+ {
+ valNames.push_back(it2.current()->asString());
+ it2.advance();
+ }
+ addAttribute(valName, valNames.size(), &valNames);
+ }
+ else
+ addAttribute(valName, it.current()->field("vc")->asInt(), NULL);
+ it.advance();
+ }
+}
+
+GArffRelation::~GArffRelation()
+{
+}
+
+// virtual
+GDomNode* GArffRelation::serialize(GDom* pDoc) const
+{
+ GDomNode* pNode = pDoc->newObj();
+ pNode->addField(pDoc, "name", pDoc->newString(m_name.c_str()));
+ GDomNode* pAttrs = pNode->addField(pDoc, "attrs", pDoc->newList());
+ for(size_t i = 0; i < m_attrs.size(); i++)
+ pAttrs->addItem(pDoc, m_attrs[i].serialize(pDoc, m_valueCounts[i]));
+ return pNode;
+}
+
+// virtual
+void GArffRelation::flush()
+{
+ GAssert(m_attrs.size() == m_valueCounts.size());
+ m_attrs.clear();
+ GMixedRelation::flush();
+}
+
+// virtual
+GRelation* GArffRelation::clone() const
+{
+ GArffRelation* pNewRelation = new GArffRelation();
+ pNewRelation->addAttrs(this);
+ pNewRelation->setName(name());
+ return pNewRelation;
+}
+
+// virtual
+GRelation* GArffRelation::cloneSub(size_t start, size_t count) const
+{
+ GArffRelation* pNewRelation = new GArffRelation();
+ pNewRelation->addAttrs(this, start, count);
+ pNewRelation->setName(name());
+ return pNewRelation;
+}
+
+void GArffRelation::addAttribute(const char* szName, size_t nValues, vector* pValues)
+{
+ GAssert(m_attrs.size() == m_valueCounts.size());
+ size_t index = m_valueCounts.size();
+ m_valueCounts.push_back(nValues);
+ m_attrs.resize(index + 1);
+ if(szName)
+ m_attrs[index].m_name = szName;
+ else
+ {
+ m_attrs[index].m_name = "attr_";
+ std::ostringstream oss;
+ oss << index;
+ m_attrs[index].m_name += oss.str();
+ }
+ if(pValues)
+ {
+ if(nValues != pValues->size())
+ throw Ex("mismatching value counts");
+ for(size_t i = 0; i < nValues; i++)
+ m_attrs[index].m_values.push_back((*pValues)[i]);
+ }
+}
+
+// virtual
+void GArffRelation::copyAttr(const GRelation* pThat, size_t nAttr)
+{
+ if(nAttr >= pThat->size())
+ throw Ex("attribute index out of range");
+ if(pThat->type() == ARFF)
+ {
+ size_t index = m_valueCounts.size();
+ GArffRelation* pOther = (GArffRelation*)pThat;
+ addAttribute(pOther->m_attrs[nAttr].m_name.c_str(), pOther->m_valueCounts[nAttr], NULL);
+ for(size_t i = 0; i < pOther->m_attrs[nAttr].m_values.size(); i++)
+ m_attrs[index].m_values.push_back(pOther->m_attrs[nAttr].m_values[i]);
+ }
+ else
+ addAttribute(NULL, pThat->valueCount(nAttr), NULL);
+}
+
+void GArffRelation::setName(const char* szName)
+{
+ m_name = szName;
+}
+
+void GArffRelation::parseAttribute(GArffTokenizer& tok)
+{
+ tok.skip(tok.m_spaces);
+ string name = tok.nextArg(tok.m_argEnd);
+ //std::cerr << "Attr:" << name << "\n"; //DEBUG
+ tok.skip(tok.m_spaces);
+ char c = tok.peek();
+ if(c == '{')
+ {
+ tok.advance(1);
+ GAssert(m_attrs.size() == m_valueCounts.size());
+ size_t index = m_valueCounts.size();
+ m_attrs.resize(index + 1);
+ while(true)
+ {
+ tok.nextArg(tok.m_valEnd);
+ char* szVal = tok.trim(tok.m_whitespace);
+ if(*szVal == '\0')
+ throw Ex("Empty value specified on line ", to_str(tok.line()));
+ if(*szVal == '\'')
+ {
+ size_t len = strlen(szVal);
+ if(len > 1 && szVal[len - 1] == '\'')
+ {
+ szVal[len - 1] = '\0';
+ szVal++;
+ }
+ }
+ else if(*szVal == '"')
+ {
+ size_t len = strlen(szVal);
+ if(len > 1 && szVal[len - 1] == '"')
+ {
+ szVal[len - 1] = '\0';
+ szVal++;
+ }
+ }
+ m_attrs[index].m_values.push_back(szVal);
+ char c = tok.peek();
+ if(c == ',')
+ tok.advance(1);
+ else if(c == '}')
+ break;
+ else if(c == '\n')
+ throw Ex("Expected a '}' but got new-line on line ", to_str(tok.line()));
+ else
+ throw Ex("inconsistency");
+ }
+ m_valueCounts.push_back(m_attrs[index].m_values.size());
+ if(name.length() > 0)
+ m_attrs[index].m_name = name;
+ else
+ {
+ m_attrs[index].m_name = "attr_";
+ std::ostringstream oss;
+ oss << index;
+ m_attrs[index].m_name += oss.str();
+ }
+ }
+ else
+ {
+ const char* szType = tok.nextUntil(tok.m_whitespace);
+ if( _stricmp(szType, "CONTINUOUS") == 0 ||
+ _stricmp(szType, "REAL") == 0 ||
+ _stricmp(szType, "NUMERIC") == 0 ||
+ _stricmp(szType, "INTEGER") == 0)
+ {
+ addAttribute(name.c_str(), 0, NULL);
+ }
+ else if(_stricmp(szType, "STRING") == 0)
+ addAttribute(name.c_str(), -1, NULL);
+ else
+ throw Ex("Unsupported attribute type: (", szType, "), at line ", to_str(tok.line()));
+ }
+ tok.skipTo(tok.m_newline);
+ tok.advance(1);
+}
+
+// virtual
+void GArffRelation::printAttrName(std::ostream& stream, size_t column) const
+{
+ stream << GRelation::quote(attrName(column));
+}
+
+// static
+std::string GRelation::quote(const std::string aString){
+ typedef std::string::const_iterator iter;
+
+ //If the string has no bad characters, just return a copy
+ std::size_t firstBad = aString.find_first_of(",' %\\\"");
+ std::string ret(aString);
+ if(firstBad == string::npos){
+ return ret;
+ }
+
+ //The string has bad characters, start over
+ ret.clear();
+
+ //If the string has no apostrophes, just quote it with single quotes
+ std::size_t firstApostrophe = aString.find_first_of('\'');
+ if(firstApostrophe == string::npos){
+ ret.push_back('\'');
+ ret.append(aString);
+ ret.push_back('\'');
+ return ret;
+ }
+
+
+ //Otherwise, use backslash to quote every character
+ ret.reserve(2*aString.size());
+ for(iter c=aString.begin();c != aString.end(); ++c)
+ {
+ if(*c == ',' || *c == '\'' ||
+ *c == ' ' || *c == '%' ||
+ *c == '\\' || *c == '"')
+ {
+ ret.push_back('\\');
+ }
+ ret.push_back(*c);
+ }
+ return ret;
+}
+
+
+// virtual
+void GArffRelation::printAttrValue(ostream& stream, size_t column, double value) const
+{
+ size_t valCount = valueCount(column);
+ if(valCount == 0)
+ {
+ if(value == UNKNOWN_REAL_VALUE)
+ stream << "?";
+ else
+ stream << value;
+ }
+ else
+ {
+ int val = (int)value;
+ if(val < 0)
+ stream << "?";
+ else if(val >= (int)valCount)
+ throw Ex("value out of range");
+ else if(m_attrs[column].m_values.size() > 0)
+ stream << GRelation::quote(m_attrs[column].m_values[val]);
+ else if(val < 26)
+ {
+ char tmp[2];
+ tmp[0] = 'a' + val;
+ tmp[1] = '\0';
+ stream << tmp;
+ }
+ else
+ stream << "_" << val;
+ }
+}
+
+// virtual
+bool GArffRelation::isCompatible(const GRelation& that) const
+{
+ if(that.type() == GRelation::ARFF)
+ {
+ if(this == &that)
+ return true;
+ if(!GRelation::isCompatible(that))
+ return false;
+ for(size_t i = 0; i < size() ; i++)
+ {
+ if(((GArffRelation*)this)->attrName(i)[0] != '\0' && ((GArffRelation*)&that)->attrName(i)[0] != '\0' && strcmp(((GArffRelation*)this)->attrName(i), ((GArffRelation*)&that)->attrName(i)) != 0)
+ return false;
+ size_t vals = valueCount(i);
+ if(vals != 0)
+ {
+ const GArffAttribute& attrThis = m_attrs[i];
+ const GArffAttribute& attrThat = ((const GArffRelation*)&that)->m_attrs[i];
+ size_t named = std::min(attrThis.m_values.size(), attrThat.m_values.size());
+ for(size_t j = 0; j < named; j++)
+ {
+ if( attrThis.m_values[j].length() != 0 &&
+ attrThat.m_values[j].length() != 0 &&
+ strcmp(attrThis.m_values[j].c_str(), attrThat.m_values[j].c_str()) != 0)
+ return false;
+ }
+ }
+ }
+ return true;
+ }
+ else
+ return GRelation::isCompatible(that);
+}
+
+int GArffRelation::findEnumeratedValue(size_t nAttr, const char* szValue) const
+{
+ size_t nValueCount = valueCount(nAttr);
+ size_t actualValCount = m_attrs[nAttr].m_values.size();
+ if(nValueCount > actualValCount)
+ throw Ex("some values have no names");
+ size_t i;
+ for(i = 0; i < nValueCount; i++)
+ {
+ if(_stricmp(m_attrs[nAttr].m_values[i].c_str(), szValue) == 0)
+ return (int)i;
+ }
+ return UNKNOWN_DISCRETE_VALUE;
+}
+
+const char* GArffRelation::attrName(size_t nAttr) const
+{
+ return m_attrs[nAttr].m_name.c_str();
+}
+
+
+int GArffRelation::addAttrValue(size_t nAttr, const char* szValue)
+{
+ int val = (int)m_valueCounts[nAttr]++;
+ GAssert(m_attrs[nAttr].m_values.size() == (size_t)val);
+ m_attrs[nAttr].m_values.push_back(szValue);
+ return val;
+}
+
+// virtual
+void GArffRelation::setAttrValueCount(size_t nAttr, size_t nValues)
+{
+ m_attrs[nAttr].m_values.clear();
+ GMixedRelation::setAttrValueCount(nAttr, nValues);
+}
+
+// virtual
+void GArffRelation::swapAttributes(size_t nAttr1, size_t nAttr2)
+{
+ GMixedRelation::swapAttributes(nAttr1, nAttr2);
+ std::swap(m_attrs[nAttr1], m_attrs[nAttr2]);
+}
+
+// virtual
+void GArffRelation::deleteAttribute(size_t nAttr)
+{
+ m_attrs.erase(m_attrs.begin() + nAttr);
+ GMixedRelation::deleteAttribute(nAttr);
+}
+
+double GArffRelation::parseValue(size_t attr, const char* val)
+{
+ size_t values = valueCount(attr);
+ if(values == 0)
+ {
+ if(strcmp(val, "?") == 0)
+ return UNKNOWN_REAL_VALUE;
+ else
+ {
+ if((*val >= '0' && *val <= '9') || *val == '-' || *val == '.')
+ {
+ }
+ else
+ throw Ex("Invalid real value, ", val, ". Expected it to start with one of {0-9,.,-}.");
+ return atof(val);
+ }
+ }
+ else
+ {
+ if(strcmp(val, "?") == 0)
+ return UNKNOWN_DISCRETE_VALUE;
+ else
+ {
+ size_t v = (size_t)-1;
+ for(size_t j = 0; j < values; j++)
+ {
+ if(_stricmp(val, m_attrs[attr].m_values[j].c_str()) == 0)
+ {
+ v = j;
+ break;
+ }
+ }
+ if(v == (size_t)-1)
+ {
+ if(*val >= '0' && *val <= '9')
+ v = atoi(val);
+ else
+ {
+ string sChoices;
+ for(size_t j = 0; j < values; j++)
+ {
+ if(j != 0)
+ sChoices += ',';
+ sChoices += m_attrs[attr].m_values[j].c_str();
+ }
+ throw Ex("Invalid categorical value, ", val, ". Expected one of {", sChoices, "}");
+ }
+ }
+ return double(v);
+ }
+ }
+}
+
+
+// ------------------------------------------------------------------
+
+GMatrix::GMatrix(sp_relation& pRelation, GHeap* pHeap)
+: m_pRelation(pRelation), m_pHeap(pHeap)
+{
+}
+
+GMatrix::GMatrix(size_t rows, size_t cols, GHeap* pHeap)
+: m_pHeap(pHeap)
+{
+ m_pRelation = new GUniformRelation(cols, 0);
+ newRows(rows);
+}
+
+GMatrix::GMatrix(vector& attrValues, GHeap* pHeap)
+: m_pHeap(pHeap)
+{
+ m_pRelation = new GMixedRelation(attrValues);
+}
+
+GMatrix::GMatrix(const GMatrix& orig):m_pHeap(orig.m_pHeap){
+ copy(&orig);
+ m_pRelation = m_pRelation->clone();
+}
+
+GMatrix& GMatrix::operator=(const GMatrix& orig){
+ copy(&orig);
+ m_pRelation = m_pRelation->clone();
+ return *this;
+}
+
+GMatrix::GMatrix(GDomNode* pNode, GHeap* pHeap)
+: m_pHeap(pHeap)
+{
+ m_pRelation = GRelation::deserialize(pNode->field("rel"));
+ GDomNode* pRows = pNode->field("pats");
+ GDomListIterator it(pRows);
+ reserve(it.remaining());
+ size_t dims = (size_t)m_pRelation->size();
+ double* pPat;
+ for(size_t i = 0; it.current(); it.advance())
+ {
+ GDomNode* pRow = it.current();
+ GDomListIterator it2(pRow);
+ if(it2.remaining() != dims)
+ throw Ex("Row ", to_str(i), " has an unexpected number of values");
+ pPat = newRow();
+ for( ; it2.current(); it2.advance())
+ {
+ *pPat = it2.current()->asDouble();
+ pPat++;
+ }
+ i++;
+ }
+}
+
+GMatrix::~GMatrix()
+{
+ flush();
+}
+
+bool GMatrix::operator==(const GMatrix& other) const{
+ //Check if relation is compatible
+ if(!relation()->isCompatible(*other.relation())){
+ return false;
+ }
+ //Check if same size
+ if(!(rows()==other.rows() && cols() == other.cols())){
+ return false;
+ }
+ //Check if have same entries
+ const std::size_t c = cols();
+ std::vector::const_iterator rThis, rOther;
+ for(rThis = m_rows.begin(), rOther = other.m_rows.begin();
+ rThis != m_rows.end(); ++rThis, ++rOther){
+ if(!std::equal(*rThis, c+*rThis, *rOther)){
+ return false;
+ }
+ }
+ return true;
+}
+
+
+void GMatrix::flush()
+{
+ if(!m_pHeap)
+ {
+ for(vector::iterator it = m_rows.begin(); it != m_rows.end(); it++)
+ delete[] (*it);
+ }
+ m_rows.clear();
+}
+
+inline bool IsRealValue(const char* szValue)
+{
+ if(*szValue == '-')
+ szValue++;
+ if(*szValue == '.')
+ szValue++;
+ if(*szValue >= '0' && *szValue <= '9')
+ return true;
+ return false;
+}
+
+double GMatrix_parseValue(GArffRelation* pRelation, size_t col, const char* szVal, GTokenizer& tok)
+{
+ size_t vals = pRelation->valueCount(col);
+ if(vals == 0)
+ {
+ // Continuous attribute
+ if(*szVal == '\0' || (*szVal == '?' && szVal[1] == '\0'))
+ return UNKNOWN_REAL_VALUE;
+ else
+ {
+ if(!IsRealValue(szVal))
+ throw Ex("Expected a numeric value at line ", to_str(tok.line()), ", col ", to_str(tok.col()));
+ return atof(szVal);
+ }
+ }
+ else if(vals != (size_t)-1)
+ {
+ // Nominal attribute
+ if(*szVal == '\0' || (*szVal == '?' && szVal[1] == '\0'))
+ return UNKNOWN_DISCRETE_VALUE;
+ else
+ {
+ int nVal = pRelation->findEnumeratedValue(col, szVal);
+ if(nVal == UNKNOWN_DISCRETE_VALUE)
+ throw Ex("Unrecognized enumeration value '", szVal, "' for attribute ", to_str(col), " at line ", to_str(tok.line()), ", col ", to_str(tok.col()));
+ return (double)nVal;
+ }
+ }
+ else
+ return 0.0;
+}
+
+
+void GMatrix::col(size_t index, double* pOutVector)
+{
+ for(size_t i = 0; i < rows(); i++)
+ *(pOutVector++) = row(i)[index];
+}
+
+void GMatrix::setCol(size_t index, const double* pVector)
+{
+ for(size_t i = 0; i < rows(); i++)
+ row(i)[index] = *(pVector++);
+}
+
+void GMatrix::add(GMatrix* pThat, bool transpose)
+{
+ if(transpose)
+ {
+ size_t c = (size_t)cols();
+ if(rows() != (size_t)pThat->cols() || c != pThat->rows())
+ throw Ex("expected matrices of same size");
+ for(size_t i = 0; i < rows(); i++)
+ {
+ double* pRow = row(i);
+ for(size_t j = 0; j < c; j++)
+ *(pRow++) += pThat->row(j)[i];
+ }
+ }
+ else
+ {
+ size_t c = cols();
+ if(rows() != pThat->rows() || c != pThat->cols())
+ throw Ex("expected matrices of same size");
+ for(size_t i = 0; i < rows(); i++)
+ GVec::add(row(i), pThat->row(i), c);
+ }
+}
+
+
+void GMatrix::subtract(GMatrix* pThat, bool transpose)
+{
+ if(transpose)
+ {
+ size_t c = (size_t)cols();
+ if(rows() != (size_t)pThat->cols() || c != pThat->rows())
+ throw Ex("expected matrices of same size");
+ for(size_t i = 0; i < rows(); i++)
+ {
+ double* pRow = row(i);
+ for(size_t j = 0; j < c; j++)
+ *(pRow++) -= pThat->row(j)[i];
+ }
+ }
+ else
+ {
+ size_t c = cols();
+ if(rows() != pThat->rows() || c != pThat->cols())
+ throw Ex("expected matrices of same size");
+ for(size_t i = 0; i < rows(); i++)
+ GVec::subtract(row(i), pThat->row(i), c);
+ }
+}
+
+void GMatrix::multiply(double scalar)
+{
+ size_t c = cols();
+ for(size_t i = 0; i < rows(); i++)
+ GVec::multiply(row(i), scalar, c);
+}
+
+void GMatrix::multiply(const double* pVectorIn, double* pVectorOut, bool transpose)
+{
+ size_t rowCount = rows();
+ size_t colCount = cols();
+ if(transpose)
+ {
+ GVec::setAll(pVectorOut, 0.0, colCount);
+ for(size_t i = 0; i < rowCount; i++)
+ GVec::addScaled(pVectorOut, *(pVectorIn++), row(i), colCount);
+ }
+ else
+ {
+ for(size_t i = 0; i < rowCount; i++)
+ *(pVectorOut++) = GVec::dotProduct(row(i), pVectorIn, colCount);
+ }
+}
+
+// static
+GMatrix* GMatrix::multiply(GMatrix& a, GMatrix& b, bool transposeA, bool transposeB)
+{
+ if(transposeA)
+ {
+ if(transposeB)
+ {
+ size_t dims = a.rows();
+ if((size_t)b.cols() != dims)
+ throw Ex("dimension mismatch");
+ size_t w = b.rows();
+ size_t h = a.cols();
+ GMatrix* pOut = new GMatrix(h, w);
+ for(size_t y = 0; y < h; y++)
+ {
+ double* pRow = pOut->row(y);
+ for(size_t x = 0; x < w; x++)
+ {
+ double* pB = b[x];
+ double sum = 0;
+ for(size_t i = 0; i < dims; i++)
+ sum += a[i][y] * pB[i];
+ *(pRow++) = sum;
+ }
+ }
+ return pOut;
+ }
+ else
+ {
+ size_t dims = a.rows();
+ if(b.rows() != dims)
+ throw Ex("dimension mismatch");
+ size_t w = b.cols();
+ size_t h = a.cols();
+ GMatrix* pOut = new GMatrix(h, w);
+ for(size_t y = 0; y < h; y++)
+ {
+ double* pRow = pOut->row(y);
+ for(size_t x = 0; x < w; x++)
+ {
+ double sum = 0;
+ for(size_t i = 0; i < dims; i++)
+ sum += a[i][y] * b[i][x];
+ *(pRow++) = sum;
+ }
+ }
+ return pOut;
+ }
+ }
+ else
+ {
+ if(transposeB)
+ {
+ size_t dims = (size_t)a.cols();
+ if((size_t)b.cols() != dims)
+ throw Ex("dimension mismatch");
+ size_t w = b.rows();
+ size_t h = a.rows();
+ GMatrix* pOut = new GMatrix(h, w);
+ for(size_t y = 0; y < h; y++)
+ {
+ double* pRow = pOut->row(y);
+ double* pA = a[y];
+ for(size_t x = 0; x < w; x++)
+ *(pRow++) = GVec::dotProduct(pA, b[x], dims);
+ }
+ return pOut;
+ }
+ else
+ {
+ size_t dims = (size_t)a.cols();
+ if(b.rows() != dims)
+ throw Ex("dimension mismatch");
+ size_t w = b.cols();
+ size_t h = a.rows();
+ GMatrix* pOut = new GMatrix(h, w);
+ for(size_t y = 0; y < h; y++)
+ {
+ double* pRow = pOut->row(y);
+ double* pA = a[y];
+ for(size_t x = 0; x < w; x++)
+ {
+ double sum = 0;
+ for(size_t i = 0; i < dims; i++)
+ sum += pA[i] * b[i][x];
+ *(pRow++) = sum;
+ }
+ }
+ return pOut;
+ }
+ }
+}
+
+GMatrix* GMatrix::transpose()
+{
+ size_t r = rows();
+ size_t c = (size_t)cols();
+ GMatrix* pTarget = new GMatrix(c, r);
+ for(size_t i = 0; i < c; i++)
+ {
+ double* pRow = pTarget->row(i);
+ for(size_t j = 0; j < r; j++)
+ *(pRow++) = row(j)[i];
+ }
+ return pTarget;
+}
+
+double GMatrix::trace()
+{
+ size_t min = std::min((size_t)cols(), rows());
+ double sum = 0;
+ for(size_t n = 0; n < min; n++)
+ sum += row(n)[n];
+ return sum;
+}
+
+size_t GMatrix::toReducedRowEchelonForm()
+{
+ size_t nLead = 0;
+ double* pRow;
+ size_t rowCount = rows();
+ size_t colCount = cols();
+ for(size_t nRow = 0; nRow < rowCount; nRow++)
+ {
+ // Find the next pivot (swapping rows as necessary)
+ size_t i = nRow;
+ while(std::abs(row(i)[nLead]) < 1e-9)
+ {
+ if(++i >= rowCount)
+ {
+ i = nRow;
+ if(++nLead >= colCount)
+ return nRow;
+ }
+ }
+ if(i > nRow)
+ swapRows(i, nRow);
+
+ // Scale the pivot to 1
+ pRow = row(nRow);
+ GVec::multiply(pRow + nLead, 1.0 / pRow[nLead], colCount - nLead);
+
+ // Elliminate all values above and below the pivot
+ for(i = 0; i < rowCount; i++)
+ {
+ if(i != nRow)
+ GVec::addScaled(row(i) + nLead, -row(i)[nLead], pRow + nLead, colCount - nLead);
+ }
+
+ nLead++;
+ }
+ return rowCount;
+}
+
+bool GMatrix::gaussianElimination(double* pVector)
+{
+ if(rows() != (size_t)cols())
+ throw Ex("Expected a square matrix");
+ double d;
+ double* pRow;
+ size_t rowCount = rows();
+ size_t colCount = cols();
+ for(size_t nRow = 0; nRow < rowCount; nRow++)
+ {
+ size_t i;
+ for(i = nRow; i < rowCount && std::abs(row(i)[nRow]) < 1e-4; i++)
+ {
+ }
+ if(i >= rowCount)
+ continue;
+ if(i > nRow)
+ {
+ swapRows(i, nRow);
+ d = pVector[i];
+ pVector[i] = pVector[nRow];
+ pVector[nRow] = d;
+ }
+
+ // Scale the pivot to 1
+ pRow = row(nRow);
+ d = 1.0 / pRow[nRow];
+ GVec::multiply(pRow + nRow, d, colCount - nRow);
+ pVector[nRow] *= d;
+
+ // Elliminate all values above and below the pivot
+ for(i = 0; i < rowCount; i++)
+ {
+ if(i != nRow)
+ {
+ d = -row(i)[nRow];
+ GVec::addScaled(row(i) + nRow, d, pRow + nRow, colCount - nRow);
+ pVector[i] += d * pVector[nRow];
+ }
+ }
+ }
+
+ // Arbitrarily assign null-space values to 1
+ for(size_t nRow = 0; nRow < rowCount; nRow++)
+ {
+ if(row(nRow)[nRow] < 0.5)
+ {
+ if(std::abs(pVector[nRow]) >= 1e-4)
+ return false;
+ for(size_t i = 0; i < rowCount; i++)
+ {
+ if(i == nRow)
+ {
+ pVector[nRow] = 1;
+ row(nRow)[nRow] = 1;
+ }
+ else
+ {
+ pVector[i] -= row(i)[nRow];
+ row(i)[nRow] = 0;
+ }
+ }
+ }
+ }
+ return true;
+}
+
+GMatrix* GMatrix::cholesky(bool tolerant)
+{
+ size_t rowCount = rows();
+ size_t colCount = (size_t)cols();
+ GMatrix* pOut = new GMatrix(m_pRelation);
+ pOut->newRows(rowCount);
+ double d;
+ for(size_t j = 0; j < rowCount; j++)
+ {
+ size_t i;
+ for(i = 0; i < j; i++)
+ {
+ d = 0;
+ for(size_t k = 0; k < i; k++)
+ d += (pOut->row(i)[k] * pOut->row(j)[k]);
+ if(std::abs(pOut->row(i)[i]) < 1e-12)
+ pOut->row(i)[i] = 1e-10;
+ pOut->row(j)[i] = (1.0 / pOut->row(i)[i]) * (row(i)[j] - d);
+ }
+ d = 0;
+ for(size_t k = 0; k < i; k++)
+ d += (pOut->row(i)[k] * pOut->row(j)[k]);
+ d = row(j)[i] - d;
+ if(d < 0)
+ {
+ if(d > -1e-12)
+ d = 0; // it's probably just rounding error
+ else if(tolerant)
+ d = -d;
+ else
+ throw Ex("not positive definite");
+ }
+ pOut->row(j)[i] = sqrt(d);
+ for(i++; i < colCount; i++)
+ pOut->row(j)[i] = 0;
+ }
+ return pOut;
+}
+
+void GMatrix::LUDecomposition()
+{
+ size_t colCount = cols();
+ double* pRow = row(0);
+ for(size_t i = 1; i < colCount; i++)
+ pRow[i] /= pRow[0];
+ for(size_t i = 1; i < colCount; i++)
+ {
+ for(size_t j = i; j < colCount; j++)
+ { // do a column of L
+ double sum = 0.0;
+ for(size_t k = 0; k < i; k++)
+ sum += row(j)[k] * row(k)[i];
+ row(j)[i] -= sum;
+ }
+ if(i == colCount - 1)
+ continue;
+ for(size_t j = i + 1; j < colCount; j++)
+ { // do a row of U
+ double sum = 0.0;
+ for(size_t k = 0; k < i; k++)
+ sum += row(i)[k] * row(k)[j];
+ row(i)[j] = (row(i)[j] - sum) / row(i)[i];
+ }
+ }
+}
+
+/*
+void GMatrix::invert()
+{
+ if(rows() != (size_t)cols())
+ throw Ex("only square matrices supported");
+ if(rows() == 1)
+ {
+ row(0)[0] = 1.0 / row(0)[0];
+ return;
+ }
+
+ // Do LU decomposition (I think this is the Doolittle algorithm)
+ int colCount = cols();
+ double* pRow = row(0);
+ for(int i = 1; i < colCount; i++)
+ pRow[i] /= pRow[0];
+ for(int i = 1; i < colCount; i++)
+ {
+ for(int j = i; j < colCount; j++)
+ { // do a column of L
+ double sum = 0.0;
+ for(int k = 0; k < i; k++)
+ sum += row(j)[k] * row(k)[i];
+ row(j)[i] -= sum;
+ }
+ if(i == colCount - 1)
+ continue;
+ for(int j = i + 1; j < colCount; j++)
+ { // do a row of U
+ double sum = 0.0;
+ for(int k = 0; k < i; k++)
+ sum += row(i)[k] * row(k)[j];
+ row(i)[j] = (row(i)[j] - sum) / row(i)[i];
+ }
+ }
+
+ // Invert L
+ for(int i = 0; i < colCount; i++)
+ {
+ for(int j = i; j < colCount; j++ )
+ {
+ double x = 1.0;
+ if ( i != j )
+ {
+ x = 0.0;
+ for(int k = i; k < j; k++ )
+ x -= row(j)[k] * row(k)[i];
+ }
+ row(j)[i] = x / row(j)[j];
+ }
+ }
+
+ // Invert U
+ for(int i = 0; i < colCount; i++)
+ {
+ for(int j = i; j < colCount; j++ )
+ {
+ if( i == j )
+ continue;
+ double sum = 0.0;
+ for (int k = i; k < j; k++ )
+ sum += row(k)[j] * ((i == k) ? 1.0 : row(i)[k]);
+ row(i)[j] = -sum;
+ }
+ }
+
+ // A^-1 = U^-1 x L^-1
+ for(int i = 0; i < colCount; i++ )
+ {
+ for(int j = 0; j < colCount; j++ )
+ {
+ double sum = 0.0;
+ for(int k = ((i > j) ? i : j); k < colCount; k++)
+ sum += ((j == k) ? 1.0 : row(j)[k]) * row(k)[i];
+ row(j)[i] = sum;
+ }
+ }
+}
+*/
+void GMatrix::inPlaceSquareTranspose()
+{
+ size_t size = rows();
+ if(size != (size_t)cols())
+ throw Ex("Expected a square matrix");
+ for(size_t a = 0; a < size; a++)
+ {
+ for(size_t b = a + 1; b < size; b++)
+ std::swap(row(a)[b], row(b)[a]);
+ }
+}
+
+double GMatrix_pythag(double a, double b)
+{
+ double at = std::abs(a);
+ double bt = std::abs(b);
+ if(at > bt)
+ {
+ double ct = bt / at;
+ return at * sqrt(1.0 + ct * ct);
+ }
+ else if(bt > 0.0)
+ {
+ double ct = at / bt;
+ return bt * sqrt(1.0 + ct * ct);
+ }
+ else
+ return 0.0;
+}
+
+double GMatrix_takeSign(double a, double b)
+{
+ return (b >= 0.0 ? std::abs(a) : -std::abs(a));
+}
+
+void GMatrix::singularValueDecomposition(GMatrix** ppU, double** ppDiag, GMatrix** ppV, bool throwIfNoConverge, size_t maxIters)
+{
+ if(rows() >= (size_t)cols())
+ singularValueDecompositionHelper(ppU, ppDiag, ppV, throwIfNoConverge, maxIters);
+ else
+ {
+ GMatrix* pTemp = transpose();
+ Holder hTemp(pTemp);
+ pTemp->singularValueDecompositionHelper(ppV, ppDiag, ppU, throwIfNoConverge, maxIters);
+ (*ppV)->inPlaceSquareTranspose();
+ (*ppU)->inPlaceSquareTranspose();
+ }
+}
+
+double GMatrix_safeDivide(double n, double d)
+{
+ if(d == 0.0 && n == 0.0)
+ return 0.0;
+ else
+ {
+ double t = n / d;
+ //GAssert(t > -1e200, "prob");
+ return t;
+ }
+}
+
+void GMatrix::fixNans()
+{
+ size_t colCount = cols();
+ for(size_t i = 0; i < rows(); i++)
+ {
+ double* pRow = row(i);
+ for(size_t j = 0; j < colCount; j++)
+ {
+ if(*pRow >= -1e308 && *pRow < 1e308)
+ {
+ }
+ else
+ *pRow = (i == (size_t)j ? 1.0 : 0.0);
+ pRow++;
+ }
+ }
+}
+
+void GMatrix::singularValueDecompositionHelper(GMatrix** ppU, double** ppDiag, GMatrix** ppV, bool throwIfNoConverge, size_t maxIters)
+{
+ int m = (int)rows();
+ int n = (int)cols();
+ if(m < n)
+ throw Ex("Expected at least as many rows as columns");
+ int i, j, k;
+ int l = 0;
+ int p, q;
+ double c, f, h, s, x, y, z;
+ double norm = 0.0;
+ double g = 0.0;
+ double scale = 0.0;
+ GMatrix* pU = new GMatrix(m, m);
+ Holder hU(pU);
+ pU->setAll(0.0);
+ pU->copyColumns(0, this, 0, n);
+ double* pSigma = new double[n];
+ ArrayHolder hSigma(pSigma);
+ GMatrix* pV = new GMatrix(n, n);
+ Holder hV(pV);
+ pV->setAll(0.0);
+ GTEMPBUF(double, temp, n);
+
+ // Householder reduction to bidiagonal form
+ for(int i = 0; i < n; i++)
+ {
+ // Left-hand reduction
+ temp[i] = scale * g;
+ l = i + 1;
+ g = 0.0;
+ s = 0.0;
+ scale = 0.0;
+ if(i < m)
+ {
+ for(k = i; k < m; k++)
+ scale += std::abs(pU->row(k)[i]);
+ if(scale != 0.0)
+ {
+ for(k = i; k < m; k++)
+ {
+ pU->row(k)[i] = GMatrix_safeDivide(pU->row(k)[i], scale);
+ double t = pU->row(k)[i];
+ s += t * t;
+ }
+ f = pU->row(i)[i];
+ g = -GMatrix_takeSign(sqrt(s), f);
+ h = f * g - s;
+ pU->row(i)[i] = f - g;
+ if(i != n - 1)
+ {
+ for(j = l; j < n; j++)
+ {
+ s = 0.0;
+ for(k = i; k < m; k++)
+ s += pU->row(k)[i] * pU->row(k)[j];
+ f = GMatrix_safeDivide(s, h);
+ for(k = i; k < m; k++)
+ pU->row(k)[j] += f * pU->row(k)[i];
+ }
+ }
+ for(k = i; k < m; k++)
+ pU->row(k)[i] *= scale;
+ }
+ }
+ pSigma[i] = scale * g;
+
+ // Right-hand reduction
+ g = 0.0;
+ s = 0.0;
+ scale = 0.0;
+ if(i < m && i != n - 1)
+ {
+ for(k = l; k < n; k++)
+ scale += std::abs(pU->row(i)[k]);
+ if(scale != 0.0)
+ {
+ for(k = l; k < n; k++)
+ {
+ pU->row(i)[k] = GMatrix_safeDivide(pU->row(i)[k], scale);
+ double t = pU->row(i)[k];
+ s += t * t;
+ }
+ f = pU->row(i)[l];
+ g = -GMatrix_takeSign(sqrt(s), f);
+ h = f * g - s;
+ pU->row(i)[l] = f - g;
+ for(k = l; k < n; k++)
+ temp[k] = GMatrix_safeDivide(pU->row(i)[k], h);
+ if(i != m - 1)
+ {
+ for(j = l; j < m; j++)
+ {
+ s = 0.0;
+ for(k = l; k < n; k++)
+ s += pU->row(j)[k] * pU->row(i)[k];
+ for(k = l; k < n; k++)
+ pU->row(j)[k] += s * temp[k];
+ }
+ }
+ for(k = l; k < n; k++)
+ pU->row(i)[k] *= scale;
+ }
+ }
+ norm = std::max(norm, std::abs(pSigma[i]) + std::abs(temp[i]));
+ }
+
+ // Accumulate right-hand transform
+ for(int i = n - 1; i >= 0; i--)
+ {
+ if(i < n - 1)
+ {
+ if(g != 0.0)
+ {
+ for(j = l; j < n; j++)
+ pV->row(i)[j] = GMatrix_safeDivide(GMatrix_safeDivide(pU->row(i)[j], pU->row(i)[l]), g); // (double-division to avoid underflow)
+ for(j = l; j < n; j++)
+ {
+ s = 0.0;
+ for(k = l; k < n; k++)
+ s += pU->row(i)[k] * pV->row(j)[k];
+ for(k = l; k < n; k++)
+ pV->row(j)[k] += s * pV->row(i)[k];
+ }
+ }
+ for(j = l; j < n; j++)
+ {
+ pV->row(i)[j] = 0.0;
+ pV->row(j)[i] = 0.0;
+ }
+ }
+ pV->row(i)[i] = 1.0;
+ g = temp[i];
+ l = i;
+ }
+
+ // Accumulate left-hand transform
+ for(i = n - 1; i >= 0; i--)
+ {
+ l = i + 1;
+ g = pSigma[i];
+ if(i < n - 1)
+ {
+ for(j = l; j < n; j++)
+ pU->row(i)[j] = 0.0;
+ }
+ if(g != 0.0)
+ {
+ g = GMatrix_safeDivide(1.0, g);
+ if(i != n - 1)
+ {
+ for(j = l; j < n; j++)
+ {
+ s = 0.0;
+ for(k = l; k < m; k++)
+ s += pU->row(k)[i] * pU->row(k)[j];
+ f = GMatrix_safeDivide(s, pU->row(i)[i]) * g;
+ for(k = i; k < m; k++)
+ pU->row(k)[j] += f * pU->row(k)[i];
+ }
+ }
+ for(j = i; j < m; j++)
+ pU->row(j)[i] *= g;
+ }
+ else
+ {
+ for(j = i; j < m; j++)
+ pU->row(j)[i] = 0.0;
+ }
+ pU->row(i)[i] += 1.0;
+ }
+
+ // Diagonalize the bidiagonal matrix
+ for(k = n - 1; k >= 0; k--) // For each singular value
+ {
+ for(size_t iter = 1; iter <= maxIters; iter++)
+ {
+ // Test for splitting
+ bool flag = true;
+ for(l = k; l >= 0; l--)
+ {
+ q = l - 1;
+ if(std::abs(temp[l]) + norm == norm)
+ {
+ flag = false;
+ break;
+ }
+ if(std::abs(pSigma[q]) + norm == norm)
+ break;
+ }
+
+ if(flag)
+ {
+ c = 0.0;
+ s = 1.0;
+ for(i = l; i <= k; i++)
+ {
+ f = s * temp[i];
+ temp[i] *= c;
+ if(std::abs(f) + norm == norm)
+ break;
+ g = pSigma[i];
+ h = GMatrix_pythag(f, g);
+ pSigma[i] = h;
+ h = GMatrix_safeDivide(1.0, h);
+ c = g * h;
+ s = -f * h;
+ for(j = 0; j < m; j++)
+ {
+ y = pU->row(j)[q];
+ z = pU->row(j)[i];
+ pU->row(j)[q] = y * c + z * s;
+ pU->row(j)[i] = z * c - y * s;
+ }
+ }
+ }
+
+ z = pSigma[k];
+ if(l == k)
+ {
+ // Detect convergence
+ if(z < 0.0)
+ {
+ // Singular value should be positive
+ pSigma[k] = -z;
+ for(j = 0; j < n; j++)
+ pV->row(k)[j] *= -1.0;
+ }
+ break;
+ }
+ if(throwIfNoConverge && iter >= maxIters)
+ throw Ex("failed to converge");
+
+ // Shift from bottom 2x2 minor
+ x = pSigma[l];
+ q = k - 1;
+ y = pSigma[q];
+ g = temp[q];
+ h = temp[k];
+ f = GMatrix_safeDivide(((y - z) * (y + z) + (g - h) * (g + h)), (2.0 * h * y));
+ g = GMatrix_pythag(f, 1.0);
+ f = GMatrix_safeDivide(((x - z) * (x + z) + h * (GMatrix_safeDivide(y, (f + GMatrix_takeSign(g, f))) - h)), x);
+
+ // QR transform
+ c = 1.0;
+ s = 1.0;
+ for(j = l; j <= q; j++)
+ {
+ i = j + 1;
+ g = temp[i];
+ y = pSigma[i];
+ h = s * g;
+ g = c * g;
+ z = GMatrix_pythag(f, h);
+ temp[j] = z;
+ c = GMatrix_safeDivide(f, z);
+ s = GMatrix_safeDivide(h, z);
+ f = x * c + g * s;
+ g = g * c - x * s;
+ h = y * s;
+ y = y * c;
+ for(p = 0; p < n; p++)
+ {
+ x = pV->row(j)[p];
+ z = pV->row(i)[p];
+ pV->row(j)[p] = x * c + z * s;
+ pV->row(i)[p] = z * c - x * s;
+ }
+ z = GMatrix_pythag(f, h);
+ pSigma[j] = z;
+ if(z != 0.0)
+ {
+ z = GMatrix_safeDivide(1.0, z);
+ c = f * z;
+ s = h * z;
+ }
+ f = c * g + s * y;
+ x = c * y - s * g;
+ for(p = 0; p < m; p++)
+ {
+ y = pU->row(p)[j];
+ z = pU->row(p)[i];
+ pU->row(p)[j] = y * c + z * s;
+ pU->row(p)[i] = z * c - y * s;
+ }
+ }
+ temp[l] = 0.0;
+ temp[k] = f;
+ pSigma[k] = x;
+ }
+ }
+
+ // Sort the singular values from largest to smallest
+ for(i = 1; i < n; i++)
+ {
+ for(j = i; j > 0; j--)
+ {
+ if(pSigma[j - 1] >= pSigma[j])
+ break;
+ pU->swapColumns(j - 1, j);
+ pV->swapRows(j - 1, j);
+ std::swap(pSigma[j - 1], pSigma[j]);
+ }
+ }
+
+ // Return results
+ pU->fixNans();
+ pV->fixNans();
+ *ppU = hU.release();
+ *ppDiag = hSigma.release();
+ *ppV = hV.release();
+}
+
+GMatrix* GMatrix::pseudoInverse()
+{
+ GMatrix* pU;
+ double* pDiag;
+ GMatrix* pV;
+ size_t colCount = cols();
+ size_t rowCount = rows();
+ if(rowCount < (size_t)colCount)
+ {
+ GMatrix* pTranspose = transpose();
+ Holder hTranspose(pTranspose);
+ pTranspose->singularValueDecompositionHelper(&pU, &pDiag, &pV, false, 80);
+ }
+ else
+ singularValueDecompositionHelper(&pU, &pDiag, &pV, false, 80);
+ Holder hU(pU);
+ ArrayHolder hDiag(pDiag);
+ Holder hV(pV);
+ GMatrix sigma(rowCount < (size_t)colCount ? colCount : rowCount, rowCount < (size_t)colCount ? rowCount : colCount);
+ sigma.setAll(0.0);
+ size_t m = std::min(rowCount, colCount);
+ for(size_t i = 0; i < m; i++)
+ {
+ if(std::abs(pDiag[i]) > 1e-9)
+ sigma[i][i] = GMatrix_safeDivide(1.0, pDiag[i]);
+ else
+ sigma[i][i] = 0.0;
+ }
+ GMatrix* pT = GMatrix::multiply(*pU, sigma, false, false);
+ Holder hT(pT);
+ if(rowCount < (size_t)colCount)
+ return GMatrix::multiply(*pT, *pV, false, false);
+ else
+ return GMatrix::multiply(*pV, *pT, true, true);
+}
+
+// static
+GMatrix* GMatrix::kabsch(GMatrix* pA, GMatrix* pB)
+{
+ GMatrix* pCovariance = GMatrix::multiply(*pA, *pB, true, false);
+ Holder hCov(pCovariance);
+ GMatrix* pU;
+ double* pDiag;
+ GMatrix* pV;
+ pCovariance->singularValueDecomposition(&pU, &pDiag, &pV);
+ Holder hU(pU);
+ delete[] pDiag;
+ Holder hV(pV);
+ GMatrix* pK = GMatrix::multiply(*pV, *pU, true, true);
+ return pK;
+}
+
+// static
+GMatrix* GMatrix::align(GMatrix* pA, GMatrix* pB)
+{
+ size_t columns = pA->cols();
+ GTEMPBUF(double, mean, columns);
+ pA->centroid(mean);
+ GMatrix* pAA = pA->clone();
+ Holder hAA(pAA);
+ pAA->centerMeanAtOrigin();
+ GMatrix* pBB = pB->clone();
+ Holder hBB(pBB);
+ pBB->centerMeanAtOrigin();
+ GMatrix* pK = GMatrix::kabsch(pBB, pAA);
+ Holder hK(pK);
+ hAA.reset(NULL);
+ GMatrix* pAligned = GMatrix::multiply(*pBB, *pK, false, true);
+ Holder hAligned(pAligned);
+ hBB.reset(NULL);
+ for(vector::iterator it = pAligned->m_rows.begin(); it != pAligned->m_rows.end(); it++)
+ GVec::add(*it, mean, columns);
+ return hAligned.release();
+}
+
+double GMatrix::determinant()
+{
+ // Check size
+ size_t n = rows();
+ if(n != cols())
+ throw Ex("Only square matrices are supported");
+
+ // Convert to a triangular matrix
+ double epsilon = 1e-10;
+ GMatrix* pC = this->clone();
+ Holder hC(pC);
+ GMatrix& C = *pC;
+ GTEMPBUF(size_t, Kp, 2 * n);
+ size_t* Lp = Kp + n;
+ size_t l, ko, lo;
+ double po,t0;
+ bool nonSingular = true;
+ size_t k = 0;
+ while(nonSingular && k < n)
+ {
+ po = C[k][k];
+ lo = k;
+ ko = k;
+ for(size_t i = k; i < n; i++)
+ for(size_t j = k; j < n; j++)
+ if(std::abs(C[i][j]) > std::abs(po))
+ {
+ po = C[i][j];
+ lo = i;
+ ko = j;
+ }
+ Lp[k] = lo;
+ Kp[k] = ko;
+ if(std::abs(po) < epsilon)
+ {
+ nonSingular = false;
+ //throw Ex("Failed to compute determinant. Pivot too small.");
+ }
+ else
+ {
+ if(lo != k)
+ {
+ for(size_t j = k; j < n; j++)
+ {
+ t0 = C[k][j];
+ C[k][j] = C[lo][j];
+ C[lo][j] = t0;
+ }
+ }
+ if(ko != k)
+ {
+ for(size_t i = 0; i < n; i++)
+ {
+ t0 = C[i][k];
+ C[i][k] = C[i][ko];
+ C[i][ko] = t0;
+ }
+ }
+ for(size_t i = k + 1; i < n; i++)
+ {
+ C[i][k] /= po;
+ for(size_t j = k + 1; j < n; j++)
+ C[i][j] -= C[i][k] * C[k][j];
+ }
+ k++;
+ }
+ }
+ if(nonSingular && std::abs(C[n - 1][n - 1]) < epsilon)
+ nonSingular = false;
+
+ // Compute determinant
+ if(!nonSingular)
+ return 0.0;
+ else
+ {
+ double det = 1.0;
+ for(k = 0; k < n; k++)
+ det *= C[k][k];
+ l = 0;
+ for(k = 0; k < n - 1; k++)
+ {
+ if(Lp[k] != k)
+ l++;
+ if(Kp[k] != k)
+ l++;
+ }
+ if((l % 2) != 0)
+ det = -det;
+ return det;
+ }
+}
+
+void GMatrix::makeIdentity()
+{
+ size_t rowCount = rows();
+ size_t colCount = cols();
+ for(size_t nRow = 0; nRow < rowCount; nRow++)
+ GVec::setAll(row(nRow), 0.0, colCount);
+ size_t nMin = std::min((size_t)colCount, rowCount);
+ for(size_t i = 0; i < nMin; i++)
+ row(i)[i] = 1.0;
+}
+
+void GMatrix::mirrorTriangle(bool upperToLower)
+{
+ size_t n = std::min(rows(), (size_t)cols());
+ if(upperToLower)
+ {
+ for(size_t i = 0; i < n; i++)
+ {
+ for(size_t j = i + 1; j < n; j++)
+ row(j)[i] = row(i)[j];
+ }
+ }
+ else
+ {
+ for(size_t i = 0; i < n; i++)
+ {
+ for(size_t j = i + 1; j < n; j++)
+ row(i)[j] = row(j)[i];
+ }
+ }
+}
+
+double GMatrix::eigenValue(const double* pEigenVector)
+{
+ // Find the element with the largest magnitude
+ size_t nEl = 0;
+ size_t colCount = cols();
+ for(size_t i = 1; i < colCount; i++)
+ {
+ if(std::abs(pEigenVector[i]) > std::abs(pEigenVector[nEl]))
+ nEl = i;
+ }
+ return GVec::dotProduct(row(nEl), pEigenVector, colCount) / pEigenVector[nEl];
+}
+
+void GMatrix::eigenVector(double eigenvalue, double* pOutVector)
+{
+ GAssert(rows() == (size_t)cols()); // Expected a square matrix
+ size_t rowCount = rows();
+ for(size_t i = 0; i < rowCount; i++)
+ row(i)[i] = row(i)[i] - eigenvalue;
+ GVec::setAll(pOutVector, 0.0, rowCount);
+ if(!gaussianElimination(pOutVector))
+ throw Ex("no solution");
+ GVec::normalize(pOutVector, rowCount);
+}
+
+GMatrix* GMatrix::eigs(size_t nCount, double* pEigenVals, GRand* pRand, bool mostSignificant)
+{
+ size_t dims = cols();
+ if(nCount > dims)
+ throw Ex("Can't have more eigenvectors than columns");
+ if(rows() != (size_t)dims)
+ throw Ex("expected a square matrix");
+
+/*
+ // The principle components of the Cholesky (square-root) matrix are the same as
+ // the eigenvectors of this matrix.
+ GMatrix* pDeviation = cholesky();
+ Holder hDeviation(pDeviation);
+ GMatrix* pData = pDeviation->transpose();
+ Holder hData(pData);
+ size_t s = pData->rows();
+ for(size_t i = 0; i < s; i++)
+ {
+ double* pRow = pData->newRow();
+ GVec::copy(pRow, pData->row(i), dims);
+ GVec::multiply(pRow, -1, dims);
+ }
+
+ // Extract the principle components
+ GMatrix* pOut = new GMatrix(m_pRelation);
+ pOut->newRows(nCount);
+ for(size_t i = 0; i < nCount; i++)
+ {
+ pData->principalComponentAboutOrigin(pOut->row(i), dims, pRand);
+ pData->removeComponentAboutOrigin(pOut->row(i), dims);
+ }
+*/
+
+ // Use the power method to compute the first few eigenvectors. todo: we really should use the Lanczos method instead
+ GMatrix* pOut = new GMatrix(m_pRelation);
+ pOut->newRows(nCount);
+ GMatrix* pA;
+ if(mostSignificant)
+ pA = clone();
+ else
+ pA = pseudoInverse();
+ Holder hA(pA);
+ GTEMPBUF(double, pTemp, dims);
+ for(size_t i = 0; i < nCount; i++)
+ {
+ // Use the power method to compute the next eigenvector
+ double* pX = pOut->row(i);
+ pRand->spherical(pX, dims);
+ for(size_t j = 0; j < 100; j++) // todo: is there a better way to detect convergence?
+ {
+ pA->multiply(pX, pTemp);
+ GVec::copy(pX, pTemp, dims);
+ GVec::safeNormalize(pX, dims, pRand);
+ }
+
+ // Compute the corresponding eigenvalue
+ double lambda = pA->eigenValue(pX);
+ if(pEigenVals)
+ pEigenVals[i] = lambda;
+
+ // Deflate (subtract out the eigenvector)
+ for(size_t j = 0; j < dims; j++)
+ {
+ double* pRow = pA->row(j);
+ for(size_t k = 0; k < dims; k++)
+ {
+ *pRow = *pRow - lambda * pX[j] * pX[k];
+ pRow++;
+ }
+ }
+ }
+
+ return pOut;
+}
+/*
+GMatrix* GMatrix::leastSignificantEigenVectors(size_t nCount, GRand* pRand)
+{
+ GMatrix* pInv = clone();
+ Holder hInv(pInv);
+ pInv->invert();
+ GMatrix* pOut = pInv->mostSignificantEigenVectors(nCount, pRand);
+ double eigenvalue;
+ for(size_t i = 0; i < nCount; i++)
+ {
+ eigenvalue = 1.0 / pInv->eigenValue(pOut->row(i));
+ GMatrix* cp = clone();
+ Holder hCp(cp);
+ cp->eigenVector(eigenvalue, pOut->row(i));
+ }
+ return pOut;
+}
+*/
+double* GMatrix::newRow()
+{
+ size_t nAttributes = m_pRelation->size();
+ double* pNewRow;
+ if(m_pHeap)
+ pNewRow = (double*)m_pHeap->allocate(sizeof(double) * nAttributes);
+ else
+ pNewRow = new double[nAttributes];
+ m_rows.push_back(pNewRow);
+ return pNewRow;
+}
+
+void GMatrix::takeRow(double* pRow)
+{
+ m_rows.push_back(pRow);
+}
+
+void GMatrix::newRows(size_t nRows)
+{
+ reserve(m_rows.size() + nRows);
+ for(size_t i = 0; i < nRows; i++)
+ newRow();
+}
+
+void GMatrix::fromVector(const double* pVec, size_t nRows)
+{
+ if(rows() < nRows)
+ newRows(nRows - rows());
+ else
+ {
+ while(rows() > nRows)
+ deleteRow(0);
+ }
+ size_t nCols = m_pRelation->size();
+ for(size_t r = 0; r < nRows; r++)
+ {
+ double* pRow = row(r);
+ GVec::copy(pRow, pVec, nCols);
+ pVec += nCols;
+ }
+}
+
+void GMatrix::toVector(double* pVec)
+{
+ size_t nCols = cols();
+ for(size_t i = 0; i < rows(); i++)
+ {
+ GVec::copy(pVec, row(i), nCols);
+ pVec += nCols;
+ }
+}
+
+void GMatrix::setAll(double val)
+{
+ size_t colCount = cols();
+ for(size_t i = 0; i < rows(); i++)
+ GVec::setAll(row(i), val, colCount);
+}
+
+void GMatrix::copy(const GMatrix* pThat)
+{
+ m_pRelation = pThat->m_pRelation;
+ flush();
+ newRows(pThat->rows());
+ copyColumns(0, pThat, 0, m_pRelation->size());
+}
+
+GMatrix* GMatrix::clone()
+{
+ GMatrix* pOther = new GMatrix(relation());
+ pOther->newRows(rows());
+ pOther->copyColumns(0, this, 0, cols());
+ return pOther;
+}
+
+GMatrix* GMatrix::cloneSub(size_t rowStart, size_t colStart, size_t rowCount, size_t colCount)
+{
+ if(rowStart + rowCount > rows())
+ throw Ex("row index out of range");
+ sp_relation pSubRel = (colCount == cols() ? m_pRelation : m_pRelation->cloneSub(colStart, colCount));
+ GMatrix* pThat = new GMatrix(pSubRel);
+ pThat->newRows(rowCount);
+ for(size_t i = 0; i < rowCount; i++)
+ GVec::copy(pThat->row(i), row(rowStart + i) + colStart, colCount);
+ return pThat;
+}
+
+void GMatrix::copyRow(const double* pRow)
+{
+ double* pNewRow = newRow();
+ GVec::copy(pNewRow, pRow, m_pRelation->size());
+}
+
+void GMatrix::copyColumns(size_t nDestStartColumn, const GMatrix* pSource, size_t nSourceStartColumn, size_t nColumnCount)
+{
+ if(rows() != pSource->rows())
+ throw Ex("expected datasets to have the same number of rows");
+ size_t count = rows();
+ for(size_t i = 0; i < count; i++)
+ GVec::copy(row(i) + nDestStartColumn, pSource->row(i) + nSourceStartColumn, nColumnCount);
+}
+
+GMatrix* GMatrix::attrSubset(size_t firstAttr, size_t attrCount)
+{
+ if(firstAttr + attrCount > m_pRelation->size())
+ throw Ex("index out of range");
+ sp_relation relNew;
+ if(relation()->type() == GRelation::UNIFORM)
+ {
+ GUniformRelation* pNewRelation = new GUniformRelation(attrCount, relation()->valueCount(firstAttr));
+ relNew = pNewRelation;
+ }
+ else
+ {
+ GMixedRelation* pNewRelation = new GMixedRelation();
+ pNewRelation->addAttrs(m_pRelation.get(), firstAttr, attrCount);
+ relNew = pNewRelation;
+ }
+ GMatrix* pNewData = new GMatrix(relNew);
+ pNewData->newRows(rows());
+ pNewData->copyColumns(0, this, firstAttr, attrCount);
+ return pNewData;
+}
+
+void GMatrix::swapRows(size_t a, size_t b)
+{
+ std::swap(m_rows[a], m_rows[b]);
+}
+
+void GMatrix::swapColumns(size_t nAttr1, size_t nAttr2)
+{
+ if(nAttr1 == nAttr2)
+ return;
+ m_pRelation = m_pRelation->clone();
+ m_pRelation->swapAttributes(nAttr1, nAttr2);
+ size_t nCount = rows();
+ double tmp;
+ double* pRow;
+ for(size_t i = 0; i < nCount; i++)
+ {
+ pRow = row(i);
+ tmp = pRow[nAttr1];
+ pRow[nAttr1] = pRow[nAttr2];
+ pRow[nAttr2] = tmp;
+ }
+}
+
+void GMatrix::deleteColumn(size_t index)
+{
+ m_pRelation = m_pRelation->clone();
+ m_pRelation->deleteAttribute(index);
+ size_t nCount = rows();
+ double* pRow;
+ size_t nAttrCountBefore = m_pRelation->size();
+ for(size_t i = 0; i < nCount; i++)
+ {
+ pRow = row(i);
+ for(size_t j = index; j < nAttrCountBefore; j++)
+ pRow[j] = pRow[j + 1];
+ }
+}
+
+double* GMatrix::releaseRow(size_t index)
+{
+ size_t last = m_rows.size() - 1;
+ double* pRow = m_rows[index];
+ m_rows[index] = m_rows[last];
+ m_rows.pop_back();
+ return pRow;
+}
+
+void GMatrix::deleteRow(size_t index)
+{
+ double* pRow = releaseRow(index);
+ if(!m_pHeap)
+ delete[] pRow;
+}
+
+double* GMatrix::releaseRowPreserveOrder(size_t index)
+{
+ double* pRow = m_rows[index];
+ m_rows.erase(m_rows.begin() + index);
+ return pRow;
+}
+
+void GMatrix::deleteRowPreserveOrder(size_t index)
+{
+ double* pRow = releaseRowPreserveOrder(index);
+ if(!m_pHeap)
+ delete[] pRow;
+}
+
+void GMatrix::releaseAllRows()
+{
+ m_rows.clear();
+}
+
+// static
+GMatrix* GMatrix::mergeHoriz(GMatrix* pSetA, GMatrix* pSetB)
+{
+ if(pSetA->rows() != pSetB->rows())
+ throw Ex("Expected same number of rows");
+ GArffRelation* pRel = new GArffRelation();
+ sp_relation spRel;
+ spRel = pRel;
+ GRelation* pRelA = pSetA->relation().get();
+ GRelation* pRelB = pSetB->relation().get();
+ size_t nSetADims = pRelA->size();
+ size_t nSetBDims = pRelB->size();
+ pRel->addAttrs(pRelA);
+ pRel->addAttrs(pRelB);
+ GMatrix* pNewSet = new GMatrix(spRel);
+ Holder hNewSet(pNewSet);
+ pNewSet->reserve(pSetA->rows());
+ double* pNewRow;
+ for(size_t i = 0; i < pSetA->rows(); i++)
+ {
+ pNewRow = pNewSet->newRow();
+ GVec::copy(pNewRow, pSetA->row(i), nSetADims);
+ GVec::copy(&pNewRow[nSetADims], pSetB->row(i), nSetBDims);
+ }
+ return hNewSet.release();
+}
+
+void GMatrix::shuffle(GRand& rand, GMatrix* pExtension)
+{
+ if(pExtension)
+ {
+ if(pExtension->rows() != rows())
+ throw Ex("Expected pExtension to have the same number of rows");
+ for(size_t n = m_rows.size(); n > 0; n--)
+ {
+ size_t r = (size_t)rand.next(n);
+ std::swap(m_rows[r], m_rows[n - 1]);
+ std::swap(pExtension->m_rows[r], pExtension->m_rows[n - 1]);
+ }
+ }
+ else
+ {
+ for(size_t n = m_rows.size(); n > 0; n--)
+ std::swap(m_rows[(size_t)rand.next(n)], m_rows[n - 1]);
+ }
+}
+
+void GMatrix::shuffle2(GRand& rand, GMatrix& other)
+{
+ for(size_t n = m_rows.size(); n > 0; n--)
+ {
+ size_t r = (size_t)rand.next(n);
+ std::swap(m_rows[r], m_rows[n - 1]);
+ std::swap(other.m_rows[r], other.m_rows[n - 1]);
+ }
+}
+
+void GMatrix::shuffleLikeCards()
+{
+ for(size_t i = 0; i < rows(); i++)
+ {
+ size_t n = i;
+ while(n & 1)
+ n = (n >> 1);
+ n = (n >> 1) + rows() / 2;
+ std::swap(m_rows[i], m_rows[n]);
+ }
+}
+
+double GMatrix::entropy(size_t nColumn)
+{
+ // Count the number of occurrences of each value
+ GAssert(m_pRelation->valueCount(nColumn) > 0); // continuous attributes are not supported
+ size_t nPossibleValues = m_pRelation->valueCount(nColumn);
+ GTEMPBUF(size_t, pnCounts, nPossibleValues);
+ size_t nTotalCount = 0;
+ memset(pnCounts, '\0', m_pRelation->valueCount(nColumn) * sizeof(size_t));
+ size_t nRows = rows();
+ for(size_t n = 0; n < nRows; n++)
+ {
+ int nValue = (int)row(n)[nColumn];
+ if(nValue < 0)
+ {
+ GAssert(nValue == UNKNOWN_DISCRETE_VALUE);
+ continue;
+ }
+ GAssert(nValue < (int)nPossibleValues);
+ pnCounts[nValue]++;
+ nTotalCount++;
+ }
+ if(nTotalCount == 0)
+ return 0;
+
+ // Total up the entropy
+ double dEntropy = 0;
+ double dRatio;
+ for(size_t n = 0; n < nPossibleValues; n++)
+ {
+ if(pnCounts[n] > 0)
+ {
+ dRatio = (double)pnCounts[n] / nTotalCount;
+ dEntropy -= dRatio * log(dRatio);
+ }
+ }
+ return M_LOG2E * dEntropy;
+}
+
+void GMatrix::splitByPivot(GMatrix* pGreaterOrEqual, size_t nAttribute, double dPivot, GMatrix* pExtensionA, GMatrix* pExtensionB)
+{
+ if(pExtensionA && pExtensionA->rows() != rows())
+ throw Ex("Expected pExtensionA to have the same number of rows as this dataset");
+ GAssert(pGreaterOrEqual->m_pHeap == m_pHeap);
+ size_t nUnknowns = 0;
+ double* pRow;
+ size_t n;
+ for(n = rows() - 1; n >= nUnknowns && n < rows(); n--)
+ {
+ pRow = row(n);
+ if(pRow[nAttribute] == UNKNOWN_REAL_VALUE)
+ {
+ std::swap(m_rows[nUnknowns], m_rows[n]);
+ if(pExtensionA)
+ std::swap(pExtensionA->m_rows[nUnknowns], pExtensionA->m_rows[n]);
+ nUnknowns++;
+ n++;
+ }
+ else if(pRow[nAttribute] >= dPivot)
+ {
+ pGreaterOrEqual->takeRow(releaseRow(n));
+ if(pExtensionA)
+ pExtensionB->takeRow(pExtensionA->releaseRow(n));
+ }
+ }
+
+ // Send all the unknowns to the side with more rows
+ if(pGreaterOrEqual->rows() > rows() - nUnknowns)
+ {
+ for(; n < rows(); n--)
+ {
+ pGreaterOrEqual->takeRow(releaseRow(n));
+ if(pExtensionA)
+ pExtensionB->takeRow(pExtensionA->releaseRow(n));
+ }
+ }
+}
+
+void GMatrix::splitByNominalValue(GMatrix* pSingleClass, size_t nAttr, int nValue, GMatrix* pExtensionA, GMatrix* pExtensionB)
+{
+ GAssert(pSingleClass->m_pHeap == m_pHeap);
+ for(size_t i = rows() - 1; i < rows(); i--)
+ {
+ double* pVec = row(i);
+ if((int)pVec[nAttr] == nValue)
+ {
+ pSingleClass->takeRow(releaseRow(i));
+ if(pExtensionA)
+ pExtensionB->takeRow(pExtensionA->releaseRow(i));
+ }
+ }
+}
+
+void GMatrix::splitBySize(GMatrix* pOtherData, size_t nOtherRows)
+{
+ GAssert(pOtherData->m_pHeap == m_pHeap);
+ if(nOtherRows > rows())
+ throw Ex("row count out of range");
+ size_t targetSize = pOtherData->rows() + nOtherRows;
+ pOtherData->reserve(targetSize);
+ while(pOtherData->rows() < targetSize)
+ pOtherData->takeRow(releaseRow(rows() - 1));
+}
+
+void GMatrix::mergeVert(GMatrix* pData)
+{
+ if(relation()->type() == GRelation::ARFF && pData->relation()->type() == GRelation::ARFF && relation().get() != pData->relation().get())
+ {
+ // Make an value mapping for pData
+ GArffRelation* pThis = (GArffRelation*)relation().get();
+ GArffRelation* pThat = (GArffRelation*)pData->relation().get();
+ if(pThis->size() != pThat->size())
+ throw Ex("Mismatching number of columns");
+ vector< vector > valueMap;
+ valueMap.resize(pThis->size());
+ for(size_t i = 0; i < pThis->size(); i++)
+ {
+ if(strcmp(pThis->attrName(i), pThat->attrName(i)) != 0)
+ throw Ex("The name of attribute ", to_str(i), " does not match");
+ if(pThis->valueCount(i) == 0 && pThat->valueCount(i) != 0)
+ throw Ex("Attribute ", to_str(i), " is continuous in one matrix and nominal in the other");
+ if(pThis->valueCount(i) != 0 && pThat->valueCount(i) == 0)
+ throw Ex("Attribute ", to_str(i), " is continuous in one matrix and nominal in the other");
+ vector& vm = valueMap[i];
+ GArffAttribute& attrThis = pThis->m_attrs[i];
+ GArffAttribute& attrThat = pThat->m_attrs[i];
+ for(size_t j = 0; j < pThat->valueCount(i); j++)
+ {
+ if(attrThis.m_values.size() >= pThis->valueCount(i) && attrThat.m_values.size() >= j && attrThat.m_values[j].length() > 0)
+ {
+ int newVal = pThis->findEnumeratedValue(i, attrThat.m_values[j].c_str());
+ if(newVal == UNKNOWN_DISCRETE_VALUE)
+ newVal = pThis->addAttrValue(i, attrThat.m_values[j].c_str());
+ vm.push_back(newVal);
+ }
+ else
+ vm.push_back(j);
+ }
+ }
+
+ // Merge the data and map the values in pData to match those in this Matrix with the same name
+ for(size_t j = 0; j < pData->rows(); j++)
+ {
+ double* pRow = pData->row(j);
+ takeRow(pRow);
+ for(size_t i = 0; i < pThis->size(); i++)
+ {
+ if(pThis->valueCount(i) != 0 && *pRow != UNKNOWN_DISCRETE_VALUE)
+ {
+ vector& vm = valueMap[i];
+ int oldVal = (int)*pRow;
+ GAssert(oldVal >= 0 && (size_t)oldVal < vm.size());
+ *pRow = (double)vm[oldVal];
+ }
+ pRow++;
+ }
+ }
+ pData->releaseAllRows();
+ }
+ else
+ {
+ if(!relation()->isCompatible(*pData->relation().get()))
+ throw Ex("The two matrices have incompatible relations");
+ for(size_t i = 0; i < pData->rows(); i++)
+ takeRow(pData->row(i));
+ pData->releaseAllRows();
+ }
+}
+
+double GMatrix::mean(size_t nAttribute)
+{
+ if(nAttribute >= cols())
+ throw Ex("attribute index out of range");
+ double sum = 0;
+ size_t missing = 0;
+ for(vector::iterator it = m_rows.begin(); it != m_rows.end(); it++)
+ {
+ if((*it)[nAttribute] == UNKNOWN_REAL_VALUE)
+ missing++;
+ else
+ sum += (*it)[nAttribute];
+ }
+ size_t count = m_rows.size() - missing;
+ if(count > 0)
+ return sum / count;
+ else
+ {
+ throw Ex("at least one value is required to compute a mean");
+ return 0.0;
+ }
+}
+
+
+void GMatrix::centroid(double* pOutMeans)
+{
+ size_t c = cols();
+ for(size_t n = 0; n < c; n++)
+ pOutMeans[n] = mean(n);
+}
+
+double GMatrix::variance(size_t nAttr, double mean)
+{
+ double d;
+ double dSum = 0;
+ size_t nMissing = 0;
+ for(vector::iterator it = m_rows.begin(); it != m_rows.end(); it++)
+ {
+ if((*it)[nAttr] == UNKNOWN_REAL_VALUE)
+ {
+ nMissing++;
+ continue;
+ }
+ d = (*it)[nAttr] - mean;
+ dSum += (d * d);
+ }
+ size_t nCount = m_rows.size() - nMissing;
+ if(nCount > 1)
+ return dSum / (nCount - 1);
+ else
+ return 0; // todo: wouldn't UNKNOWN_REAL_VALUE be better here?
+}
+
+void GMatrix::minAndRange(size_t nAttribute, double* pMin, double* pRange)
+{
+ double dMin = 1e300;
+ double dMax = -1e300;
+ for(vector::iterator it = m_rows.begin(); it != m_rows.end(); it++)
+ {
+ if((*it)[nAttribute] == UNKNOWN_REAL_VALUE)
+ continue;
+ if((*it)[nAttribute] < dMin)
+ dMin = (*it)[nAttribute];
+ if((*it)[nAttribute] > dMax)
+ dMax = (*it)[nAttribute];
+ }
+ if(dMax >= dMin)
+ {
+ *pMin = dMin;
+ *pRange = dMax - dMin;
+ }
+ else
+ {
+ *pMin = UNKNOWN_REAL_VALUE;
+ *pRange = UNKNOWN_REAL_VALUE;
+ }
+}
+
+void GMatrix::minAndRangeUnbiased(size_t nAttribute, double* pMin, double* pRange)
+{
+ double min, range, d;
+ minAndRange(nAttribute, &min, &range);
+ d = .5 * (range * (rows() + 1) / (rows() - 1) - range);
+ *pMin = (min - d);
+ *pRange = (range + d);
+}
+
+void GMatrix::normalize(size_t nAttribute, double dInputMin, double dInputRange, double dOutputMin, double dOutputRange)
+{
+ GAssert(dInputRange > 0);
+ double dScale = dOutputRange / dInputRange;
+ for(vector::iterator it = m_rows.begin(); it != m_rows.end(); it++)
+ {
+ (*it)[nAttribute] -= dInputMin;
+ (*it)[nAttribute] *= dScale;
+ (*it)[nAttribute] += dOutputMin;
+ }
+}
+
+/*static*/ double GMatrix::normalize(double dVal, double dInputMin, double dInputRange, double dOutputMin, double dOutputRange)
+{
+ GAssert(dInputRange > 0);
+ dVal -= dInputMin;
+ dVal /= dInputRange;
+ dVal *= dOutputRange;
+ dVal += dOutputMin;
+ return dVal;
+}
+
+double GMatrix::baselineValue(size_t nAttribute)
+{
+ if(m_pRelation->valueCount(nAttribute) == 0)
+ return mean(nAttribute);
+ int j;
+ int val;
+ int nValues = (int)m_pRelation->valueCount(nAttribute);
+ GTEMPBUF(size_t, counts, nValues + 1);
+ memset(counts, '\0', sizeof(size_t) * (nValues + 1));
+ for(vector::iterator it = m_rows.begin(); it != m_rows.end(); it++)
+ {
+ val = (int)(*it)[nAttribute] + 1;
+ counts[val]++;
+ }
+ val = 1;
+ for(j = 2; j <= nValues; j++)
+ {
+ if(counts[j] > counts[val])
+ val = j;
+ }
+ return (double)(val - 1);
+}
+
+bool GMatrix::isAttrHomogenous(size_t col)
+{
+ if(m_pRelation->valueCount(col) > 0)
+ {
+ int d;
+ vector::iterator it = m_rows.begin();
+ for( ; it != m_rows.end(); it++)
+ {
+ d = (int)(*it)[col];
+ if(d != UNKNOWN_DISCRETE_VALUE)
+ {
+ it++;
+ break;
+ }
+ }
+ for( ; it != m_rows.end(); it++)
+ {
+ int t = (int)(*it)[col];
+ if(t != d && t != UNKNOWN_DISCRETE_VALUE)
+ return false;
+ }
+ }
+ else
+ {
+ double d;
+ vector::iterator it = m_rows.begin();
+ for( ; it != m_rows.end(); it++)
+ {
+ d = (*it)[col];
+ if(d != UNKNOWN_REAL_VALUE)
+ {
+ it++;
+ break;
+ }
+ }
+ for( ; it != m_rows.end(); it++)
+ {
+ double t = (*it)[col];
+ if(t != d && t != UNKNOWN_REAL_VALUE)
+ return false;
+ }
+ }
+ return true;
+}
+
+bool GMatrix::isHomogenous()
+{
+ for(size_t i = 0; i < cols(); i++)
+ {
+ if(!isAttrHomogenous(i))
+ return false;
+ }
+ return true;
+}
+
+void GMatrix::replaceMissingValuesWithBaseline(size_t nAttr)
+{
+ double bl = baselineValue(nAttr);
+ size_t count = rows();
+ if(m_pRelation->valueCount(nAttr) == 0)
+ {
+ for(size_t i = 0; i < count; i++)
+ {
+ if(row(i)[nAttr] == UNKNOWN_REAL_VALUE)
+ row(i)[nAttr] = bl;
+ }
+ }
+ else
+ {
+ for(size_t i = 0; i < count; i++)
+ {
+ if(row(i)[nAttr] == UNKNOWN_DISCRETE_VALUE)
+ row(i)[nAttr] = bl;
+ }
+ }
+}
+
+void GMatrix::replaceMissingValuesRandomly(size_t nAttr, GRand* pRand)
+{
+ GTEMPBUF(size_t, indexes, rows());
+
+ // Find the rows that are not missing values in this attribute
+ size_t* pCur = indexes;
+ double dOk = m_pRelation->valueCount(nAttr) == 0 ? -1e300 : 0;
+ for(size_t i = 0; i < rows(); i++)
+ {
+ if(row(i)[nAttr] >= dOk)
+ {
+ *pCur = i;
+ pCur++;
+ }
+ }
+
+ // Replace missing values
+ size_t nonMissing = pCur - indexes;
+ for(size_t i = 0; i < rows(); i++)
+ {
+ if(row(i)[nAttr] < dOk)
+ row(i)[nAttr] = row(indexes[(size_t)pRand->next(nonMissing)])[nAttr];
+ }
+}
+
+void GMatrix::principalComponent(double* pOutVector, const double* pMean, GRand* pRand)
+{
+ // Initialize the out-vector to a random direction
+ size_t dims = cols();
+ pRand->spherical(pOutVector, dims);
+
+ // Iterate
+ double* pVector;
+ size_t nCount = rows();
+ GTEMPBUF(double, pAccumulator, dims);
+ double d;
+ double mag = 0;
+ for(size_t iters = 0; iters < 200; iters++)
+ {
+ GVec::setAll(pAccumulator, 0.0, dims);
+ for(size_t n = 0; n < nCount; n++)
+ {
+ pVector = row(n);
+ d = GVec::dotProduct(pMean, pVector, pOutVector, dims);
+ double* pAcc = pAccumulator;
+ const double* pM = pMean;
+ for(size_t j = 0; j < dims; j++)
+ *(pAcc++) += d * (*(pVector++) - *(pM++));
+ }
+ GVec::copy(pOutVector, pAccumulator, dims);
+ GVec::safeNormalize(pOutVector, dims, pRand);
+ d = GVec::squaredMagnitude(pAccumulator, dims);
+ if(iters < 6 || d - mag > 1e-8)
+ mag = d;
+ else
+ break;
+ }
+}
+
+void GMatrix::principalComponentAboutOrigin(double* pOutVector, GRand* pRand)
+{
+ // Initialize the out-vector to a random direction
+ size_t dims = cols();
+ pRand->spherical(pOutVector, dims);
+
+ // Iterate
+ double* pVector;
+ size_t nCount = rows();
+ GTEMPBUF(double, pAccumulator, dims);
+ double d;
+ double mag = 0;
+ for(size_t iters = 0; iters < 200; iters++)
+ {
+ GVec::setAll(pAccumulator, 0.0, dims);
+ for(size_t n = 0; n < nCount; n++)
+ {
+ pVector = row(n);
+ d = GVec::dotProduct(pVector, pOutVector, dims);
+ double* pAcc = pAccumulator;
+ for(size_t j = 0; j < dims; j++)
+ *(pAcc++) += d * *(pVector++);
+ }
+ GVec::copy(pOutVector, pAccumulator, dims);
+ GVec::safeNormalize(pOutVector, dims, pRand);
+ d = GVec::squaredMagnitude(pAccumulator, dims);
+ if(iters < 6 || d - mag > 1e-8)
+ mag = d;
+ else
+ break;
+ }
+}
+
+void GMatrix::principalComponentIgnoreUnknowns(double* pOutVector, const double* pMean, GRand* pRand)
+{
+ if(!doesHaveAnyMissingValues())
+ {
+ principalComponent(pOutVector, pMean, pRand);
+ return;
+ }
+
+ // Initialize the out-vector to a random direction
+ size_t dims = cols();
+ pRand->spherical(pOutVector, dims);
+
+ // Iterate
+ double* pVector;
+ size_t nCount = rows();
+ GTEMPBUF(double, pAccumulator, dims);
+ double d;
+ double mag = 0;
+ for(size_t iters = 0; iters < 200; iters++)
+ {
+ GVec::setAll(pAccumulator, 0.0, dims);
+ for(size_t n = 0; n < nCount; n++)
+ {
+ pVector = row(n);
+ d = GVec::dotProductIgnoringUnknowns(pMean, pVector, pOutVector, dims);
+ double* pAcc = pAccumulator;
+ const double* pM = pMean;
+ for(size_t j = 0; j < dims; j++)
+ {
+ if(*pVector != UNKNOWN_REAL_VALUE)
+ (*pAcc) += d * (*pVector - *pM);
+ pVector++;
+ pAcc++;
+ pM++;
+ }
+ }
+ GVec::copy(pOutVector, pAccumulator, dims);
+ GVec::safeNormalize(pOutVector, dims, pRand);
+ d = GVec::squaredMagnitude(pAccumulator, dims);
+ if(iters < 6 || d - mag > 1e-8)
+ mag = d;
+ else
+ break;
+ }
+}
+
+void GMatrix::weightedPrincipalComponent(double* pOutVector, const double* pMean, const double* pWeights, GRand* pRand)
+{
+ // Initialize the out-vector to a random direction
+ size_t dims = cols();
+ pRand->spherical(pOutVector, dims);
+
+ // Iterate
+ double* pVector;
+ size_t nCount = rows();
+ GTEMPBUF(double, pAccumulator, dims);
+ double d;
+ double mag = 0;
+ for(size_t iters = 0; iters < 200; iters++)
+ {
+ GVec::setAll(pAccumulator, 0.0, dims);
+ const double* pW = pWeights;
+ for(size_t n = 0; n < nCount; n++)
+ {
+ pVector = row(n);
+ d = GVec::dotProduct(pMean, pVector, pOutVector, dims);
+ double* pAcc = pAccumulator;
+ const double* pM = pMean;
+ for(size_t j = 0; j < dims; j++)
+ *(pAcc++) += (*pW) * d * (*(pVector++) - *(pM++));
+ pW++;
+ }
+ GVec::copy(pOutVector, pAccumulator, dims);
+ GVec::safeNormalize(pOutVector, dims, pRand);
+ d = GVec::squaredMagnitude(pAccumulator, dims);
+ if(iters < 6 || d - mag > 1e-8)
+ mag = d;
+ else
+ break;
+ }
+}
+
+double GMatrix::eigenValue(const double* pMean, const double* pEigenVector, GRand* pRand)
+{
+ // Use the element of the eigenvector with the largest magnitude,
+ // because that will give us the least rounding error when we compute the eigenvalue.
+ size_t dims = cols();
+ size_t index = GVec::indexOfMaxMagnitude(pEigenVector, dims, pRand);
+
+ // The eigenvalue is the factor by which the eigenvector is scaled by the covariance matrix,
+ // so we compute just the part of the covariance matrix that we need to see how much the
+ // max-magnitude element of the eigenvector is scaled.
+ double d = 0;
+ for(size_t i = 0; i < dims; i++)
+ d += covariance(index, pMean[index], i, pMean[i]) * pEigenVector[i];
+ return d / pEigenVector[index];
+}
+
+void GMatrix::removeComponent(const double* pMean, const double* pComponent)
+{
+ size_t dims = cols();
+ size_t nCount = rows();
+ for(size_t i = 0; i < nCount; i++)
+ {
+ double* pVector = row(i);
+ double d = GVec::dotProductIgnoringUnknowns(pMean, pVector, pComponent, dims);
+ for(size_t j = 0; j < dims; j++)
+ {
+ if(*pVector != UNKNOWN_REAL_VALUE)
+ (*pVector) -= d * pComponent[j];
+ pVector++;
+ }
+ }
+}
+
+void GMatrix::removeComponentAboutOrigin(const double* pComponent)
+{
+ size_t dims = cols();
+ size_t nCount = rows();
+ for(size_t i = 0; i < nCount; i++)
+ {
+ double* pVector = row(i);
+ double d = GVec::dotProduct(pVector, pComponent, dims);
+ for(size_t j = 0; j < dims; j++)
+ {
+ (*pVector) -= d * pComponent[j];
+ pVector++;
+ }
+ }
+}
+
+void GMatrix::centerMeanAtOrigin()
+{
+ //Calculate mean
+ size_t dims = cols();
+ GTEMPBUF(double, mean, dims);
+ centroid(mean);
+ //Skip non-continuous rows by setting their mean to 0
+ for(unsigned i = 0; i < dims; ++i){
+ if(relation()->valueCount(i) != 0){ mean[i] = 0; }
+ }
+ //Subtract the new mean from all rows
+ for(vector::iterator it = m_rows.begin();
+ it != m_rows.end(); it++){
+ GVec::subtract(*it, mean, dims);
+ }
+}
+
+size_t GMatrix::countPrincipalComponents(double d, GRand* pRand)
+{
+ size_t dims = cols();
+ GMatrix tmpData(relation(), heap());
+ tmpData.copy(this);
+ tmpData.centerMeanAtOrigin();
+ GTEMPBUF(double, vec, dims);
+ double thresh = d * d * tmpData.sumSquaredDistance(NULL);
+ size_t i;
+ for(i = 1; i < dims; i++)
+ {
+ tmpData.principalComponentAboutOrigin(vec, pRand);
+ tmpData.removeComponentAboutOrigin(vec);
+ if(tmpData.sumSquaredDistance(NULL) < thresh)
+ break;
+ }
+ return i;
+}
+
+double GMatrix::sumSquaredDistance(const double* pPoint)
+{
+ size_t dims = relation()->size();
+ double err = 0;
+ if(pPoint)
+ {
+ for(size_t i = 0; i < rows(); i++)
+ err += GVec::squaredDistance(pPoint, row(i), dims);
+ }
+ else
+ {
+ for(size_t i = 0; i < rows(); i++)
+ err += GVec::squaredMagnitude(row(i), dims);
+ }
+ return err;
+}
+
+double GMatrix::columnSumSquaredDifference(GMatrix& that, size_t col)
+{
+ if(that.rows() != rows())
+ throw Ex("Mismatching number of rows");
+ if(col >= cols() || col >= that.cols())
+ throw Ex("column index out of range");
+ double sse = 0.0;
+ if(relation()->valueCount(col) == 0)
+ {
+ for(size_t i = 0; i < rows(); i++)
+ {
+ double d = row(i)[col] - that.row(i)[col];
+ sse += (d * d);
+ }
+ }
+ else
+ {
+ for(size_t i = 0; i < rows(); i++)
+ {
+ if((int)row(i)[col] != (int)that.row(i)[col])
+ sse++;
+ }
+ }
+ return sse;
+}
+
+double GMatrix::sumSquaredDifference(GMatrix& that, bool transpose)
+{
+ if(transpose)
+ {
+ size_t colCount = (size_t)cols();
+ if(rows() != (size_t)that.cols() || colCount != that.rows())
+ throw Ex("expected matrices of same size");
+ double err = 0;
+ for(size_t i = 0; i < rows(); i++)
+ {
+ double* pRow = row(i);
+ for(size_t j = 0; j < colCount; j++)
+ {
+ double d = *(pRow++) - that[j][i];
+ err += (d * d);
+ }
+ }
+ return err;
+ }
+ else
+ {
+ if(this->rows() != that.rows() || this->cols() != that.cols())
+ throw Ex("mismatching sizes");
+ size_t colCount = cols();
+ double d = 0;
+ for(size_t i = 0; i < rows(); i++)
+ d += GVec::squaredDistance(this->row(i), that[i], colCount);
+ return d;
+ }
+}
+
+double GMatrix::linearCorrelationCoefficient(size_t attr1, double attr1Origin, size_t attr2, double attr2Origin)
+{
+ double sx = 0;
+ double sy = 0;
+ double sxy = 0;
+ double mx, my;
+ double* pPat;
+ size_t count = rows();
+ size_t i;
+ for(i = 0; i < count; i++)
+ {
+ pPat = row(i);
+ mx = pPat[attr1] - attr1Origin;
+ my = pPat[attr2] - attr2Origin;
+ if(pPat[attr1] == UNKNOWN_REAL_VALUE || pPat[attr2] == UNKNOWN_REAL_VALUE)
+ continue;
+ break;
+ }
+ if(i >= count)
+ return 0;
+ double d, x, y;
+ size_t j = 1;
+ for(i++; i < count; i++)
+ {
+ pPat = row(i);
+ if(pPat[attr1] == UNKNOWN_REAL_VALUE || pPat[attr2] == UNKNOWN_REAL_VALUE)
+ continue;
+ x = pPat[attr1] - attr1Origin;
+ y = pPat[attr2] - attr2Origin;
+ d = (double)j / (j + 1);
+ sx += (x - mx) * (x - mx) * d;
+ sy += (y - my) * (y - my) * d;
+ sxy += (x - mx) * (y - my) * d;
+ mx += (x - mx) / (j + 1);
+ my += (y - my) / (j + 1);
+ j++;
+ }
+ if(sx == 0 || sy == 0 || sxy == 0)
+ return 0;
+ return (sxy / j) / (sqrt(sx / j) * sqrt(sy / j));
+}
+
+double GMatrix::covariance(size_t nAttr1, double dMean1, size_t nAttr2, double dMean2)
+{
+ size_t nRowCount = rows();
+ double* pVector;
+ double dSum = 0;
+ for(size_t i = 0; i < nRowCount; i++)
+ {
+ pVector = row(i);
+ dSum += ((pVector[nAttr1] - dMean1) * (pVector[nAttr2] - dMean2));
+ }
+ return dSum / (nRowCount - 1);
+}
+
+GMatrix* GMatrix::covarianceMatrix()
+{
+ size_t colCount = cols();
+ GMatrix* pOut = new GMatrix(colCount, colCount);
+
+ // Compute the deviations
+ GTEMPBUF(double, pMeans, colCount);
+ for(size_t i = 0; i < colCount; i++)
+ pMeans[i] = mean(i);
+
+ // Compute the covariances for half the matrix
+ for(size_t i = 0; i < colCount; i++)
+ {
+ double* pRow = pOut->row(i);
+ pRow += i;
+ for(size_t n = i; n < colCount; n++)
+ *(pRow++) = covariance(i, pMeans[i], n, pMeans[n]);
+ }
+
+ // Fill out the other half of the matrix
+ for(size_t i = 1; i < colCount; i++)
+ {
+ double* pRow = pOut->row(i);
+ for(size_t n = 0; n < i; n++)
+ *(pRow++) = pOut->row(n)[i];
+ }
+ return pOut;
+}
+
+class Row_Binary_Predicate_Functor
+{
+protected:
+ size_t m_dim;
+
+public:
+ Row_Binary_Predicate_Functor(size_t dim) : m_dim(dim)
+ {
+ }
+
+ bool operator() (double* pA, double* pB) const
+ {
+ return pA[m_dim] < pB[m_dim];
+ }
+};
+
+void GMatrix::sort(size_t dim)
+{
+ Row_Binary_Predicate_Functor comparer(dim);
+ std::sort(m_rows.begin(), m_rows.end(), comparer);
+}
+
+void GMatrix::sortPartial(size_t row, size_t col)
+{
+ Row_Binary_Predicate_Functor comparer(col);
+ vector::iterator targ = m_rows.begin() + row;
+ std::nth_element(m_rows.begin(), targ, m_rows.end(), comparer);
+}
+
+void GMatrix::reverseRows()
+{
+ std::reverse(m_rows.begin(), m_rows.end());
+}
+
+double GMatrix_PairedTTestHelper(void* pThis, double x)
+{
+ double v = *(double*)pThis;
+ return pow(1.0 + x * x / v, -(v + 1) / 2);
+}
+
+void GMatrix::pairedTTest(size_t* pOutV, double* pOutT, size_t attr1, size_t attr2, bool normalize)
+{
+ double* pPat;
+ double a, b, m;
+ double asum = 0;
+ double asumOfSquares = 0;
+ double bsum = 0;
+ double bsumOfSquares = 0;
+ size_t rowCount = rows();
+ for(size_t i = 0; i < rowCount; i++)
+ {
+ pPat = row(i);
+ a = pPat[attr1];
+ b = pPat[attr2];
+ if(normalize)
+ {
+ m = (a + b) / 2;
+ a /= m;
+ b /= m;
+ }
+ asum += a;
+ asumOfSquares += (a * a);
+ bsum += b;
+ bsumOfSquares += (b * b);
+ }
+ double amean = asum / rowCount;
+ double avariance = (asumOfSquares / rowCount - amean * amean) * rowCount / (rowCount - 1);
+ double bmean = bsum / rowCount;
+ double bvariance = (bsumOfSquares / rowCount - bmean * bmean) * rowCount / (rowCount - 1);
+ double grand = sqrt((avariance + bvariance) / rowCount);
+ *pOutV = 2 * rowCount - 2;
+ *pOutT = std::abs(bmean - amean) / grand;
+}
+
+void GMatrix::wilcoxonSignedRanksTest(size_t attr1, size_t attr2, double tolerance, int* pNum, double* pWMinus, double* pWPlus)
+{
+ // Make sorted list of differences
+ GHeap heap(1024);
+ GMatrix tmp(0, 2, &heap); // col 0 holds the absolute difference. col 1 holds the sign.
+ for(size_t i = 0; i < rows(); i++)
+ {
+ double* pPat = row(i);
+ double absdiff = std::abs(pPat[attr2] - pPat[attr1]);
+ if(absdiff > tolerance)
+ {
+ double* pStat = tmp.newRow();
+ pStat[0] = absdiff;
+ if(pPat[attr1] < pPat[attr2])
+ pStat[1] = -1;
+ else
+ pStat[1] = 1;
+ }
+ }
+
+ // Convert column 0 to ranks
+ tmp.sort(0);
+ double prev = UNKNOWN_REAL_VALUE;
+ size_t index = 0;
+ size_t j;
+ double ave;
+ for(size_t i = 0; i < tmp.rows(); i++)
+ {
+ double* pPat = tmp[i];
+ if(std::abs(pPat[0] - prev) >= tolerance)
+ {
+ ave = (double)(index + 1 + i) / 2;
+ for(j = index; j < i; j++)
+ {
+ double* pStat = tmp[j];
+ pStat[0] = ave;
+ }
+ prev = pPat[0];
+ index = i;
+ }
+ }
+ ave = (double)(index + 1 + tmp.rows()) / 2;
+ for(j = index; j < tmp.rows(); j++)
+ {
+ double* pStat = tmp[j];
+ pStat[0] = ave;
+ }
+
+ // Sum up the scores
+ double a = 0;
+ double b = 0;
+ for(size_t i = 0; i < tmp.rows(); i++)
+ {
+ double* pStat = tmp[i];
+ if(pStat[1] > 0)
+ a += pStat[0];
+ else if(pStat[1] < 0)
+ b += pStat[0];
+ else
+ {
+ a += 0.5 * pStat[0];
+ b += 0.5 * pStat[0];
+ }
+ }
+ *pNum = tmp.rows();
+ *pWMinus = b;
+ *pWPlus = a;
+}
+
+size_t GMatrix::countValue(size_t attribute, double value)
+{
+ size_t count = 0;
+ for(size_t i = 0; i < rows(); i++)
+ {
+ if(row(i)[attribute] == value)
+ count++;
+ }
+ return count;
+}
+
+bool GMatrix::doesHaveAnyMissingValues()
+{
+ size_t dims = m_pRelation->size();
+ for(size_t j = 0; j < dims; j++)
+ {
+ if(m_pRelation->valueCount(j) == 0)
+ {
+ for(size_t i = 0; i < rows(); i++)
+ {
+ if(row(i)[j] == UNKNOWN_REAL_VALUE)
+ return true;
+ }
+ }
+ else
+ {
+ for(size_t i = 0; i < rows(); i++)
+ {
+ if(row(i)[j] == UNKNOWN_DISCRETE_VALUE)
+ return true;
+ }
+ }
+ }
+ return false;
+}
+
+void GMatrix::ensureDataHasNoMissingReals()
+{
+ size_t dims = m_pRelation->size();
+ for(size_t i = 0; i < rows(); i++)
+ {
+ double* pPat = row(i);
+ for(size_t j = 0; j < dims; j++)
+ {
+ if(m_pRelation->valueCount(j) != 0)
+ continue;
+ if(pPat[i] == UNKNOWN_REAL_VALUE)
+ throw Ex("Missing values in continuous attributes are not supported");
+ }
+ }
+}
+
+void GMatrix::ensureDataHasNoMissingNominals()
+{
+ size_t dims = m_pRelation->size();
+ for(size_t i = 0; i < rows(); i++)
+ {
+ double* pPat = row(i);
+ for(size_t j = 0; j < dims; j++)
+ {
+ if(m_pRelation->valueCount(j) == 0)
+ continue;
+ if((int)pPat[i] == UNKNOWN_DISCRETE_VALUE)
+ throw Ex("Missing values in nominal attributes are not supported");
+ }
+ }
+}
+
+void GMatrix::print(ostream& stream)
+{
+ m_pRelation->print(stream, this, 14);
+}
+
+double GMatrix::measureInfo()
+{
+ size_t c = cols();
+ double dInfo = 0;
+ for(size_t n = 0; n < c; n++)
+ {
+ if(m_pRelation->valueCount(n) == 0)
+ {
+ if(rows() > 1)
+ {
+ double m = mean(n);
+ dInfo += variance(n, m);
+ }
+ }
+ else
+ dInfo += entropy(n);
+ }
+ return dInfo;
+}
+
+double GMatrix::sumSquaredDiffWithIdentity()
+{
+ size_t m = std::min(rows(), (size_t)cols());
+ double err = 0;
+ double d;
+ for(size_t i = 0; i < m; i++)
+ {
+ double* pRow = row(i);
+ for(size_t j = 0; j < m; j++)
+ {
+ d = *(pRow++);
+ if(i == j)
+ d -= 1;
+ err += (d * d);
+ }
+ }
+ return err;
+}
+/*
+bool GMatrix::leastCorrelatedVector(double* pOut, GMatrix* pThat, GRand* pRand)
+{
+ if(rows() != pThat->rows() || cols() != pThat->cols())
+ throw Ex("Expected matrices with the same dimensions");
+ GMatrix* pC = GMatrix::multiply(*pThat, *this, false, true);
+ Holder hC(pC);
+ GMatrix* pE = GMatrix::multiply(*pC, *pC, true, false);
+ Holder hE(pE);
+ double d = pE->sumSquaredDiffWithIdentity();
+ if(d < 0.001)
+ return false;
+ GMatrix* pF = pE->mostSignificantEigenVectors(rows(), pRand);
+ Holder hF(pF);
+ GVec::copy(pOut, pF->row(rows() - 1), rows());
+ return true;
+}
+*/
+bool GMatrix::leastCorrelatedVector(double* pOut, GMatrix* pThat, GRand* pRand)
+{
+ if(rows() != pThat->rows() || cols() != pThat->cols())
+ throw Ex("Expected matrices with the same dimensions");
+ GMatrix* pC = GMatrix::multiply(*pThat, *this, false, true);
+ Holder hC(pC);
+ GMatrix* pD = GMatrix::multiply(*pThat, *pC, true, false);
+ Holder hD(pD);
+ double d = pD->sumSquaredDifference(*this, true);
+ if(d < 1e-9)
+ return false;
+ pD->subtract(this, true);
+ pD->principalComponentAboutOrigin(pOut, pRand);
+ return true;
+/*
+ GMatrix* pE = GMatrix::multiply(*pD, *pD, true, false);
+ Holder hE(pE);
+ GMatrix* pF = pE->mostSignificantEigenVectors(1, pRand);
+ Holder hF(pF);
+ GVec::copy(pOut, pF->row(0), rows());
+ return true;
+*/
+}
+
+double GMatrix::dihedralCorrelation(GMatrix* pThat, GRand* pRand)
+{
+ size_t colCount = cols();
+ if(rows() == 1)
+ return std::abs(GVec::dotProduct(row(0), pThat->row(0), colCount));
+ GTEMPBUF(double, pBuf, rows() + 2 * colCount);
+ double* pA = pBuf + rows();
+ double* pB = pA + colCount;
+ if(!leastCorrelatedVector(pBuf, pThat, pRand))
+ return 1.0;
+ multiply(pBuf, pA, true);
+ if(!pThat->leastCorrelatedVector(pBuf, this, pRand))
+ return 1.0;
+ pThat->multiply(pBuf, pB, true);
+ return std::abs(GVec::correlation(pA, pB, colCount));
+}
+
+void GMatrix::project(double* pDest, const double* pPoint)
+{
+ size_t dims = cols();
+ GVec::setAll(pDest, 0.0, dims);
+ for(size_t i = 0; i < rows(); i++)
+ {
+ double* pBasis = row(i);
+ GVec::addScaled(pDest, GVec::dotProduct(pPoint, pBasis, dims), pBasis, dims);
+ }
+}
+
+void GMatrix::project(double* pDest, const double* pPoint, const double* pOrigin)
+{
+ size_t dims = cols();
+ GVec::copy(pDest, pOrigin, dims);
+ for(size_t i = 0; i < rows(); i++)
+ {
+ double* pBasis = row(i);
+ GVec::addScaled(pDest, GVec::dotProduct(pOrigin, pPoint, pBasis, dims), pBasis, dims);
+ }
+}
+
+
+
+
+std::string to_str(const GMatrix& m){
+ std::stringstream out;
+ out.precision(14);
+ if(!m.relation()->areContinuous(0,m.cols())){
+ out << '[';
+ for(unsigned j = 0; j < m.cols(); ++j){
+ if(j != 0){ out << ", "; }
+ if(m.relation()->areContinuous(j,1)){
+ out << "continuous";
+ }else{
+ out << "discrete:" << m.relation()->valueCount(j);
+ }
+ }
+ out << "]\n";
+ }else{
+ out << "Continuous matrix:\n";
+ }
+
+ out << '[';
+ for(unsigned i=0; i < m.rows(); ++i){
+ for(unsigned j=0; j < m.cols(); ++j){
+ if(j != 0){ out << ", "; }
+ out << m[i][j];
+ }
+ if(i+1 < m.rows()){
+ out << "\n";
+ }
+ }
+ out << "]\n";
+ return out.str();
+}
+
+
+
+
+GMatrixArray::GMatrixArray(sp_relation& pRelation)
+: m_pRelation(pRelation)
+{
+}
+
+GMatrixArray::GMatrixArray(size_t cols)
+{
+ m_pRelation = new GUniformRelation(cols, 0);
+}
+
+GMatrixArray::~GMatrixArray()
+{
+ flush();
+}
+
+void GMatrixArray::flush()
+{
+ for(vector::iterator it = m_sets.begin(); it != m_sets.end(); it++)
+ delete(*it);
+ m_sets.clear();
+}
+
+GMatrix* GMatrixArray::newSet(size_t rows)
+{
+ GMatrix* pData = new GMatrix(m_pRelation);
+ m_sets.push_back(pData);
+ pData->newRows(rows);
+ return pData;
+}
+
+void GMatrixArray::newSets(size_t count, size_t rows)
+{
+ m_sets.reserve(m_sets.size() + count);
+ for(size_t i = 0; i < count; i++)
+ newSet(rows);
+}
+
+size_t GMatrixArray::largestSet()
+{
+ vector::iterator it = m_sets.begin();
+ size_t biggestRows = (*it)->rows();
+ size_t biggestIndex = 0;
+ size_t i = 1;
+ for(it++; it != m_sets.end(); it++)
+ {
+ if((*it)->rows() > biggestRows)
+ {
+ biggestRows = (*it)->rows();
+ biggestIndex = i;
+ }
+ i++;
+ }
+ return biggestIndex;
+}
+
+size_t GMatrixArray::countEmptySets()
+{
+ size_t count = 0;
+ for(vector::iterator it = m_sets.begin(); it != m_sets.end(); it++)
+ {
+ if((*it)->rows() == 0)
+ count++;
+ }
+ return count;
+}
+
+} // namespace GClasses
diff --git a/src/ai/waffles/GMatrix.h b/src/ai/waffles/GMatrix.h
new file mode 100644
index 0000000..25cebc5
--- /dev/null
+++ b/src/ai/waffles/GMatrix.h
@@ -0,0 +1,1211 @@
+/*
+ Copyright (C) 2006, Mike Gashler
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ see http://www.gnu.org/copyleft/lesser.html
+*/
+
+#ifndef __GMATRIX_H__
+#define __GMATRIX_H__
+
+#include
+#include
+#include
+#include
+
+#include "GError.h"
+#include "GHolders.h"
+
+
+namespace GClasses {
+
+#define UNKNOWN_REAL_VALUE -1e308
+
+// Why do we need a different value for unknown discrete values? Because it's
+// common to store discrete values in an integer. Integers can't store -1e308,
+// and we can't use -1 for unknown reals b/c it's not an insane value.
+#define UNKNOWN_DISCRETE_VALUE -1
+
+
+class GMatrix;
+class GPrediction;
+class GRand;
+class GHeap;
+class GDom;
+class GDomNode;
+class GArffTokenizer;
+class GDistanceMetric;
+class GSimpleAssignment;
+
+/// \brief Holds the metadata for a dataset.
+///
+/// The metadata includes which attributes are continuous or nominal,
+/// and how many values each nominal attribute supports.
+class GRelation
+{
+public:
+ enum RelationType
+ {
+ UNIFORM,
+ MIXED,
+ ARFF
+ };
+
+ GRelation() {}
+ virtual ~GRelation() {}
+
+ /// \brief Returns the type of relation
+ virtual RelationType type() const = 0;
+
+ /// \brief Marshal this object into a DOM, which can then be
+ /// converted to a variety of serial formats.
+ virtual GDomNode* serialize(GDom* pDoc) const = 0;
+
+ /// \brief Returns the number of attributes (columns)
+ virtual size_t size() const = 0;
+
+ /// \brief Returns the number of values in the specified
+ /// attribute. (Returns 0 for continuous attributes.)
+ virtual size_t valueCount(size_t nAttr) const = 0;
+
+ /// \brief Returns true of all of the attributes in the specified
+ /// range are continuous
+ virtual bool areContinuous(size_t first = 0, size_t count = (size_t)-1) const = 0;
+
+ /// \brief Returns true of all of the attributes in the specified
+ /// range are nominal
+ virtual bool areNominal(size_t first, size_t count) const = 0;
+
+ /// \brief Makes a deep copy of this relation
+ virtual GRelation* clone() const = 0;
+
+ /// \brief Makes a deep copy of the specified subset of this relation
+ virtual GRelation* cloneSub(size_t start, size_t count) const = 0;
+
+ /// \brief Deletes the specified attribute
+ virtual void deleteAttribute(size_t index) = 0;
+
+ /// \brief Swaps two attributes
+ virtual void swapAttributes(size_t nAttr1, size_t nAttr2) = 0;
+
+ /// \brief Prints as an ARFF file to the specified stream. (pData
+ /// can be NULL if data is not available)
+ void print(std::ostream& stream, const GMatrix* pData, size_t precision) const;
+
+ /// \brief Prints the specified attribute name to a stream
+ virtual void printAttrName(std::ostream& stream, size_t column) const;
+
+ /// \brief Prints the specified value to a stream
+ virtual void printAttrValue(std::ostream& stream, size_t column, double value) const;
+
+ /// \brief Returns the name of the attribute with index \a nAttr as
+ /// a standard string object or "" if the atribute has no name
+ ///
+ /// \param nAttr the index of the attribute whose name is returned
+ ///
+ /// \return the name of the attribute with index \a nAttr as a
+ /// standard string object or "" if the atribute has no name
+ virtual std::string attrNameStr(std::size_t nAttr) const{ return ""; }
+
+ /// \brief Returns true iff this and that have the same number of
+ /// values for each attribute
+ virtual bool isCompatible(const GRelation& that) const;
+
+ /// \brief Print a single row of data in ARFF format
+ void printRow(std::ostream& stream, const double* pRow, const char* separator) const;
+
+ /// \brief Counts the size of the corresponding real-space vector
+ size_t countRealSpaceDims(size_t nFirstAttr, size_t nAttrCount) const;
+
+ /// \brief Converts a row (pIn) to a real-space vector (pOut) (pIn
+ /// should point to the nFirstAttr'th element, not the first
+ /// element)
+ void toRealSpace(const double* pIn, double* pOut, size_t nFirstAttr, size_t nAttrCount) const;
+
+ /// \brief Converts a real-space vector (pIn) to a row (pOut)
+ ///
+ /// nFirstAttr and nAttrCount refer to the row indexes
+ void fromRealSpace(const double* pIn, double* pOut, size_t nFirstAttr, size_t nAttrCount, GRand* pRand) const;
+
+ /// \brief Converts a real-space vector (pIn) to an array of
+ /// predictions (pOut)
+ ///
+ /// nFirstAttr and nAttrCount refer to the prediction indexes
+ void fromRealSpace(const double* pIn, GPrediction* pOut, size_t nFirstAttr, size_t nAttrCount) const;
+
+ /// \brief Load from a DOM.
+ static smart_ptr deserialize(GDomNode* pNode);
+
+ /// \brief Saves to a file
+ void save(const GMatrix* pData, const char* szFilename, size_t precision) const;
+
+
+ protected:
+ /// \brief Returns a copy of aString modified to escape internal
+ /// instances of comma, apostrophe, space, percent, back-slash, and
+ /// double-quote
+ static std::string quote(const std::string aString);
+};
+
+typedef smart_ptr sp_relation;
+
+
+/// \brief A relation with a minimal memory footprint that assumes all
+/// attributes are continuous, or all of them are nominal and have the
+/// same number of possible values.
+class GUniformRelation : public GRelation
+{
+protected:
+ size_t m_attrCount;
+ size_t m_valueCount;
+
+public:
+ GUniformRelation(size_t attrCount, size_t valueCount = 0)
+ : m_attrCount(attrCount), m_valueCount(valueCount)
+ {
+ }
+
+ GUniformRelation(GDomNode* pNode);
+
+ virtual RelationType type() const { return UNIFORM; }
+
+ /// \brief Serializes this object
+ virtual GDomNode* serialize(GDom* pDoc) const;
+
+ /// \brief Returns the number of attributes (columns)
+ virtual size_t size() const { return m_attrCount; }
+
+ /// \brief Returns the number of values in each nominal attribute
+ /// (or 0 if the attributes are continuous)
+ virtual size_t valueCount(size_t) const { return m_valueCount; }
+
+ /// \brief See the comment for GRelation::areContinuous
+ virtual bool areContinuous(size_t, size_t) const { return m_valueCount == 0; }
+
+ /// \brief See the comment for GRelation::areNominal
+ virtual bool areNominal(size_t, size_t) const { return m_valueCount != 0; }
+
+ /// \brief Returns a copy of this object
+ virtual GRelation* clone() const { return new GUniformRelation(m_attrCount, m_valueCount); }
+
+ /// \brief Returns a deep copy of the specified subset of this
+ /// relation
+ virtual GRelation* cloneSub(size_t, size_t count) const { return new GUniformRelation(count, m_valueCount); }
+
+ /// \brief Drop the specified attribute
+ virtual void deleteAttribute(size_t index);
+
+ /// \brief Swap two attributes (since all attributes are identical, does nothing)
+ virtual void swapAttributes(size_t, size_t) {}
+
+ /// \brief Returns true iff this and that have the same number of
+ /// values for each attribute
+ virtual bool isCompatible(const GRelation& that) const;
+};
+
+
+
+class GMixedRelation : public GRelation
+{
+protected:
+ std::vector m_valueCounts;
+
+public:
+ /// \brief Makes an empty relation
+ GMixedRelation();
+
+ /// \brief Construct a mixed relation. attrValues specifies the
+ /// number of nominal values in each attribute (column), or 0 for
+ /// continuous attributes.
+ GMixedRelation(std::vector& attrValues);
+
+ /// \brief Loads from a DOM.
+ GMixedRelation(GDomNode* pNode);
+
+ /// \brief Makes a copy of pCopyMe
+ GMixedRelation(const GRelation* pCopyMe);
+
+ /// \brief Makes a copy of the specified range of pCopyMe
+ GMixedRelation(const GRelation* pCopyMe, size_t firstAttr, size_t attrCount);
+
+ virtual ~GMixedRelation();
+
+ virtual RelationType type() const { return MIXED; }
+
+ /// \brief Marshalls this object to a DOM, which can be saved to a
+ /// variety of serial formats.
+ virtual GDomNode* serialize(GDom* pDoc) const;
+
+ /// \brief Makes a deep copy of this relation
+ virtual GRelation* clone() const;
+
+ /// \brief Makes a deep copy of the specified subset of this
+ /// relation
+ virtual GRelation* cloneSub(size_t start, size_t count) const;
+
+ /// \brief Deletes all the attributes
+ virtual void flush();
+
+ /// \brief If nValues is zero, adds a real attribute. If nValues is
+ /// > 0, adds an attribute with "nValues" nominal values
+ void addAttr(size_t nValues);
+
+ /// \brief Adds "attrCount" new attributes, each with "valueCount"
+ /// values. (Use valueCount=0 for continuous attributes.)
+ void addAttrs(size_t attrCount, size_t valueCount);
+
+ /// \brief Copies the specified attributes and adds them to this
+ /// relation. If attrCount < 0, then it will copy all attributes
+ /// from firstAttr to the end.
+ void addAttrs(const GRelation* pCopyMe, size_t firstAttr = 0, size_t attrCount = (size_t)-1);
+
+ /// \brief Flushes this relation and then copies all of the
+ /// attributes from pCopyMe
+ void copy(const GRelation* pCopyMe);
+
+ /// \brief Adds a copy of the specified attribute to this relation
+ virtual void copyAttr(const GRelation* pThat, size_t nAttr);
+
+ /// \brief Returns the total number of attributes in this relation
+ virtual size_t size() const
+ {
+ return m_valueCounts.size();
+ }
+
+ /// \brief Returns the number of nominal values in the specified
+ /// attribute
+ virtual size_t valueCount(size_t nAttr) const
+ {
+ return m_valueCounts[nAttr];
+ }
+
+ /// \brief Sets the number of values for this attribute
+ virtual void setAttrValueCount(size_t nAttr, size_t nValues);
+
+ /// \brief Returns true iff all attributes in the specified range are continuous
+ virtual bool areContinuous(size_t first = 0, size_t count = (size_t)-1) const;
+
+ /// \brief Returns true iff all attributes in the specified range are nominal
+ virtual bool areNominal(size_t first = 0, size_t count = (size_t)-1) const;
+
+ /// \brief Swaps two columns
+ virtual void swapAttributes(size_t nAttr1, size_t nAttr2);
+
+ /// \brief Deletes an attribute.
+ virtual void deleteAttribute(size_t nAttr);
+};
+
+
+class GArffAttribute
+{
+public:
+ std::string m_name;
+ std::vector m_values;
+
+ GDomNode* serialize(GDom* pDoc, size_t valCount) const;
+};
+
+
+/// \brief ARFF = Attribute-Relation File Format. This stores richer
+/// information than GRelation. This includes a name, a name for each
+/// attribute, and names for each supported nominal value.
+class GArffRelation : public GMixedRelation
+{
+friend class GMatrix;
+protected:
+ std::string m_name;
+ std::vector m_attrs;
+
+public:
+ /// General-purpose constructor
+ GArffRelation();
+
+ /// Deserializing constructor
+ GArffRelation(GDomNode* pNode);
+
+ virtual ~GArffRelation();
+
+ virtual RelationType type() const { return ARFF; }
+
+ /// \brief Marshalls this object to a DOM, which can be saved to a
+ /// variety of serial formats.
+ virtual GDomNode* serialize(GDom* pDoc) const;
+
+ /// \brief Returns a deep copy of this object
+ virtual GRelation* clone() const;
+
+ /// \brief Makes a deep copy of the specified subset of this relation
+ virtual GRelation* cloneSub(size_t start, size_t count) const;
+
+ /// \brief Deletes all the attributes
+ virtual void flush();
+
+ /// \brief Prints the specified attribute name to a stream
+ virtual void printAttrName(std::ostream& stream, size_t column) const;
+
+ /// \brief Prints the specified value to a stream
+ virtual void printAttrValue(std::ostream& stream, size_t column, double value) const;
+
+ /// \brief Returns true iff the attributes in both relations have
+ /// the same names, the same number of values, and the names of
+ /// those values all match. (Empty strings are considered to match
+ /// everything.)
+ virtual bool isCompatible(const GRelation& that) const;
+
+ /// \brief Adds a new attribute (column) to the relation
+ void addAttribute(const char* szName, size_t nValues, std::vector* pValues);
+
+ /// \brief Adds a copy of the specified attribute to this relation
+ virtual void copyAttr(const GRelation* pThat, size_t nAttr);
+
+ /// \brief Returns the name of the relation
+ const char* name() const { return m_name.c_str(); }
+
+ /// \brief Sets the name of this relation
+ void setName(const char* szName);
+
+ /// \brief Returns the name of the attribute with index \a nAttr
+ ///
+ /// \param nAttr the index of the attribute whose name is returned
+ ///
+ /// \return the name of the attribute with index \a nAttr as a
+ /// standard string object or "" if the atribute has no name
+ const char* attrName(size_t nAttr) const;
+
+ /// \brief Returns the name of the attribute with index \a nAttr as
+ /// a standard string object or "" if the atribute has no name
+ ///
+ /// \param nAttr the index of the attribute whose name is returned
+ ///
+ /// \return the name of the attribute with index \a nAttr as a
+ /// standard string object or "" if the atribute has no name
+ virtual std::string attrNameStr(std::size_t nAttr) const {
+ return attrName(nAttr); }
+
+
+ /// \brief Adds a new possible value to a nominal attribute. Returns
+ /// the numerical form of the new value.
+ int addAttrValue(size_t nAttr, const char* szValue);
+
+ /// \brief Sets the number of values for the specified attribute
+ virtual void setAttrValueCount(size_t nAttr, size_t nValues);
+
+ /// \brief Swaps two columns
+ virtual void swapAttributes(size_t nAttr1, size_t nAttr2);
+
+ /// \brief Deletes an attribute
+ virtual void deleteAttribute(size_t nAttr);
+
+ /// \brief Returns the nominal index for the specified attribute
+ /// with the given value
+ int findEnumeratedValue(size_t nAttr, const char* szValue) const;
+
+ /// \brief Parses a value
+ double parseValue(size_t attr, const char* val);
+
+ /// \brief Parses the meta-data for an attribute
+ void parseAttribute(GArffTokenizer& tok);
+
+};
+
+/// \brief Represents a matrix or a database table.
+///
+/// Elements can be discrete or continuous.
+///
+/// References a GRelation object, which stores the meta-information about each column.
+class GMatrix
+{
+protected:
+ sp_relation m_pRelation;
+ GHeap* m_pHeap;
+ std::vector m_rows;
+
+public:
+ /// \brief Construct a rows x cols matrix with all elements of the
+ /// matrix assumed to be continuous.
+ ///
+ /// It is okay to initially set rows to 0 and later call newRow to
+ /// add each row. Adding columns later, however, is not very
+ /// computationally efficient.)
+ GMatrix(size_t rows, size_t cols, GHeap* pHeap = NULL);
+
+ /// \brief Construct a matrix with a mixed relation. That is, one
+ /// with some continuous attributes (columns), and some nominal
+ /// attributes (columns).
+ ///
+ /// attrValues specifies the number of nominal values suppored in
+ /// each attribute (column), or 0 for a continuous attribute.
+ ///
+ /// Initially, this matrix will have 0 rows, but you can add more
+ /// rows by calling newRow or newRows.
+ GMatrix(std::vector& attrValues, GHeap* pHeap = NULL);
+
+ /// \brief Create an empty matrix whose attributes/column types are
+ /// specified by pRelation
+ ///
+ /// pRelation is a smart-pointer to a relation, which specifies the
+ /// type of each attribute (column) in the data set.
+ ///
+ /// Initially, this matrix will have 0 rows, but you can add more
+ /// rows by calling newRow or newRows.
+ GMatrix(sp_relation& pRelation, GHeap* pHeap = NULL);
+
+ ///\brief Copy-constructor
+ ///
+ ///Copies \a orig, making a new relation object and new storage for
+ ///the rows (with the same content), but uses the same GHeap object
+ ///as \a orig
+ ///
+ ///\param orig the GMatrix object to copy
+ GMatrix(const GMatrix& orig);
+
+
+ //I put the operator= right after the copy constructor because you
+ //should always have both or neither in a class, this makes that
+ //easy to verify
+
+ ///\brief Make *this into a copy of orig
+ ///
+ ///Copies \a orig, making a new relation object and new storage for
+ ///the rows (with the same content), but uses the same GHeap object
+ ///as \a orig
+ ///
+ ///\param orig the GMatrix object to copy
+ ///
+ ///\return a reference to this GMatrix object
+ GMatrix& operator=(const GMatrix& orig);
+
+
+ /// \brief Load from a DOM.
+ GMatrix(GDomNode* pNode, GHeap* pHeap = NULL);
+
+ ~GMatrix();
+
+ /// \brief Returns true iff all the entries in *this and \a
+ /// other are identical and their relations are compatible, and they
+ /// are the same size
+ ///
+ /// \return true iff all the entries in *this and \a other are
+ /// identical, their relations are compatible, and they are the same
+ /// size
+ bool operator==(const GMatrix& other) const;
+
+ /// \brief Adds a new row to the dataset. (The values in the row are
+ /// not initialized)
+ double* newRow();
+
+ /// \brief Adds "nRows" uninitialized rows to the data set
+ void newRows(size_t nRows);
+
+ /// \brief Matrix add
+ ///
+ /// Adds the values in pThat to this. (If transpose is true, adds
+ /// the transpose of pThat to this.) Both datasets must have the
+ /// same dimensions. Behavior is undefined for nominal columns.
+ void add(GMatrix* pThat, bool transpose);
+
+ /// \brief Returns a new dataset that contains a subset of the
+ /// attributes in this dataset
+ GMatrix* attrSubset(size_t firstAttr, size_t attrCount);
+
+ /// \brief This computes the square root of this matrix. (If you
+ /// take the matrix that this returns and multiply it by its
+ /// transpose, you should get the original dataset again.)
+ /// (Returns a lower-triangular matrix.)
+ ///
+ /// Behavior is undefined if there are nominal attributes. If
+ /// tolerant is true, it will return even if it cannot compute
+ /// accurate results. If tolerant is false (the default) and this
+ /// matrix is not positive definate, it will throw an exception.
+ GMatrix* cholesky(bool tolerant = false);
+
+ /// \brief Makes a deep copy of this dataset
+ GMatrix* clone();
+
+ /// \brief Makes a deep copy of the specified rectangular region of
+ /// this matrix
+ GMatrix* cloneSub(size_t rowStart, size_t colStart, size_t rowCount, size_t colCount);
+
+ /// \brief Copies the specified column into pOutVector
+ void col(size_t index, double* pOutVector);
+
+ /// \brief Returns the number of columns in the dataset
+ size_t cols() const { return m_pRelation->size(); }
+
+ /// \brief Copies all the data from pThat. (Just references the same
+ /// relation)
+ void copy(const GMatrix* pThat);
+
+ /// \brief Copies the specified block of columns from pSource to
+ /// this dataset.
+ ///
+ /// pSource must have the same number of rows as this dataset.
+ void copyColumns(size_t nDestStartColumn, const GMatrix* pSource, size_t nSourceStartColumn, size_t nColumnCount);
+
+ /// \brief Adds a copy of the row to the data set
+ void copyRow(const double* pRow);
+
+ /// \brief Computes the determinant of this matrix
+ double determinant();
+
+
+ /// \brief Computes the eigenvalue that corresponds to the specified
+ /// eigenvector of this matrix
+ double eigenValue(const double* pEigenVector);
+
+ /// \brief Computes the eigenvector that corresponds to the
+ /// specified eigenvalue of this matrix. Note that this method
+ /// trashes this matrix, so make a copy first if you care.
+ void eigenVector(double eigenvalue, double* pOutVector);
+
+ /// \brief Computes y in the equation M*y=x (or y=M^(-1)x), where M
+ /// is this dataset, which must be a square matrix, and x is pVector
+ /// as passed in, and y is pVector after the call.
+ ///
+ /// If there are multiple solutions, it finds the one for which all
+ /// the variables in the null-space have a value of 1. If there are
+ /// no solutions, it returns false. Note that this method trashes
+ /// this dataset (so make a copy first if you care).
+ bool gaussianElimination(double* pVector);
+
+ /// \brief Returns the heap used to allocate rows for this dataset
+ GHeap* heap() { return m_pHeap; }
+
+ /// \brief Performs an in-place LU-decomposition, such that the
+ /// lower triangle of this matrix (including the diagonal) specifies
+ /// L, and the uppoer triangle of this matrix (not including the
+ /// diagonal) specifies U, and all values of U along the diagonal
+ /// are ones. (The upper triangle of L and the lower triangle of U
+ /// are all zeros.)
+ void LUDecomposition();
+
+ /// \brief This computes K=kabsch(A,B), such that K is an n-by-n
+ /// matrix, where n is pA->cols(). K is the optimal orthonormal
+ /// rotation matrix to align A and B, such that A(K^T) minimizes
+ /// sum-squared error with B, and BK minimizes sum-squared error
+ /// with A. (This rotates around the origin, so typically you will
+ /// want to subtract the centroid from both pA and pB before calling
+ /// this.)
+ static GMatrix* kabsch(GMatrix* pA, GMatrix* pB);
+
+ /// \brief This uses the Kabsch algorithm to rotate and translate pB
+ /// in order to minimize RMS with pA. (pA and pB must have the same
+ /// number of rows and columns.)
+ static GMatrix* align(GMatrix* pA, GMatrix* pB);
+
+
+
+ /// \brief Sets this dataset to an identity matrix. (It doesn't
+ /// change the number of columns or rows. It just stomps over
+ /// existing values.)
+ void makeIdentity();
+
+ /// \brief copies one of the triangular submatrices over the other,
+ /// making a symmetric matrix.
+ ///
+ /// \param upperToLower If true, copies the upper triangle of this
+ /// matrix over the lower triangle. Otherwise,
+ /// copies the lower triangle of this matrix
+ /// over the upper triangle
+ void mirrorTriangle(bool upperToLower);
+
+ /// \brief Merges two datasets side-by-side. The resulting dataset
+ /// will contain the attributes of both datasets. Both pSetA and
+ /// pSetB (and the resulting dataset) must have the same number of
+ /// rows
+ static GMatrix* mergeHoriz(GMatrix* pSetA, GMatrix* pSetB);
+
+ /// \brief Steals all the rows from pData and adds them to this set.
+ /// (You still have to delete pData.) Both datasets must have the
+ /// same number of columns.
+ void mergeVert(GMatrix* pData);
+
+ /// \brief Computes nCount eigenvectors and the corresponding
+ /// eigenvalues using the power method (which is only accurate if a
+ /// small number of eigenvalues/vectors are needed.)
+ ///
+ /// If mostSignificant is true, the largest eigenvalues are
+ /// found. If mostSignificant is false, the smallest eigenvalues are
+ /// found.
+ GMatrix* eigs(size_t nCount, double* pEigenVals, GRand* pRand, bool mostSignificant);
+
+ /// \brief Multiplies every element in the dataset by scalar.
+ /// Behavior is undefined for nominal columns.
+ void multiply(double scalar);
+
+ /// \brief Multiplies this matrix by the column vector pVectorIn to
+ /// get pVectorOut.
+ ///
+ /// (If transpose is true, then it multiplies the transpose of this
+ /// matrix by pVectorIn to get pVectorOut.)
+ ///
+ /// pVectorIn should have
+ /// the same number of elements as columns (or rows if transpose is
+ /// true)
+ ///
+ /// pVectorOut should have the same number of elements as
+ /// rows (or cols, if transpose is true.)
+ ///
+ /// \note if transpose is true, then pVectorIn is treated as a
+ /// row vector and is multiplied by this matrix to get pVectorOut.
+ void multiply(const double* pVectorIn, double* pVectorOut, bool transpose = false);
+
+ /// \brief Matrix multiply.
+ ///
+ /// For convenience, you can also specify that neither, one, or both
+ /// of the inputs are virtually transposed prior to the
+ /// multiplication. (If you want the results to come out transposed,
+ /// you can use the equality (AB)^T=(B^T)(A^T) to figure out how to
+ /// specify the parameters.)
+ static GMatrix* multiply(GMatrix& a, GMatrix& b, bool transposeA, bool transposeB);
+
+
+ /// \brief Computes the Moore-Penrose pseudoinverse of this matrix
+ /// (using the SVD method). You are responsible to delete the
+ /// matrix this returns.
+ GMatrix* pseudoInverse();
+
+ /// \brief Returns a relation object, which holds meta-data about
+ /// the attributes (columns)
+ sp_relation& relation() { return m_pRelation; }
+
+ /// \brief Returns a relation object, which holds meta-data about
+ /// the attributes (columns) (const version)
+ smart_ptr relation() const {
+ return smart_ptr (m_pRelation, NULL);
+ }
+
+ /// \brief Allocates space for the specified number of patterns (to
+ /// avoid superfluous resizing)
+ void reserve(size_t n) { m_rows.reserve(n); }
+
+ /// \brief Returns the number of rows in the dataset
+ size_t rows() const { return m_rows.size(); }
+
+
+ /// \brief Sets the relation for this dataset
+ void setRelation(sp_relation& pRelation) { m_pRelation = pRelation; }
+
+ /// \brief Performs SVD on A, where A is this m-by-n matrix.
+ ///
+ /// You are responsible to delete(*ppU), delete(*ppV), and delete[]
+ /// *ppDiag.
+ ///
+ /// \param ppU *ppU will be set to an m-by-m matrix where the
+ /// columns are the *eigenvectors of A(A^T).
+ ///
+ /// \param ppDiag *ppDiag will be set to an array of n doubles
+ /// holding the square roots of the corresponding
+ /// eigenvalues.
+ ///
+ /// \param ppV *ppV will be set to an n-by-n matrix where the rows
+ /// are the eigenvectors of (A^T)A.
+ ///
+ /// \param throwIfNoConverge if true, throws an exception if the SVD
+ /// solver does not converge. does nothing
+ /// otherwise
+ ///
+ /// \param maxIters the maximum number of iterations to perform in
+ /// the SVD solver
+ void singularValueDecomposition(GMatrix** ppU, double** ppDiag, GMatrix** ppV, bool throwIfNoConverge = false, size_t maxIters = 80);
+
+ /// \brief Matrix subtract. Subtracts the values in *pThat from *this.
+ ///
+ /// (If transpose is true, subtracts the transpose of *pThat from
+ /// this.) Both datasets must have the same dimensions. Behavior is
+ /// undefined for nominal columns.
+ ///
+ /// \param pThat pointer to the matrix to subtract from *this
+ ///
+ /// \param transpose If true, the transpose of *pThat is subtracted.
+ /// If false, *pThat is subtracted
+ void subtract(GMatrix* pThat, bool transpose);
+
+ /// \brief Returns the sum squared difference between this matrix
+ /// and an identity matrix
+ double sumSquaredDiffWithIdentity();
+
+ /// \brief Adds an already-allocated row to this dataset. The row must
+ /// be allocated in the same heap that this dataset uses. (There is no way
+ /// for this method to verify that, so be careful.)
+ void takeRow(double* pRow);
+
+ /// \brief Converts the matrix to reduced row echelon form
+ size_t toReducedRowEchelonForm();
+
+ /// \brief Copies all the data from this dataset into pVector.
+ ///
+ /// pVector must be big enough to hold rows() x cols() doubles.
+ void toVector(double* pVector);
+
+
+ /// \brief Returns the sum of the diagonal elements
+ double trace();
+
+ /// \brief Returns a pointer to a new dataset that is this dataset
+ /// transposed. (All columns in the returned dataset will be
+ /// continuous.)
+ ///
+ /// The returned matrix is newly allocated on the system heap with
+ /// operator new and must be deleted by the caller.
+ ///
+ /// \return A pointer to a new dataset that is this dataset
+ /// transposed. All columns in the returned dataset will be
+ /// continuous. The caller is responsible for deleting the
+ /// returned dataset.
+ GMatrix* transpose();
+
+ /// \brief Copies the data from pVector over this dataset.
+ ///
+ /// nRows specifies the number of rows of data in pVector.
+ void fromVector(const double* pVector, size_t nRows);
+
+ /// \brief Returns a pointer to the specified row
+ inline double* row(size_t index) { return m_rows[index]; }
+
+ /// \brief Returns a pointer to the specified row
+ inline double* operator [](size_t index) { return m_rows[index]; }
+
+ /// \brief Returns a const pointer to the specified row
+ inline const double* row(size_t index) const { return m_rows[index]; }
+
+ /// \brief Returns a const pointer to the specified row
+ inline const double* operator [](size_t index) const {
+ return m_rows[index]; }
+
+ /// \brief Sets all elements in this dataset to the specified value
+ void setAll(double val);
+
+ /// \brief Copies pVector over the specified column
+ void setCol(size_t index, const double* pVector);
+
+ /// \brief Swaps the two specified rows
+ void swapRows(size_t a, size_t b);
+
+ /// \brief Swaps two columns
+ void swapColumns(size_t nAttr1, size_t nAttr2);
+
+ /// \brief Deletes a column
+ void deleteColumn(size_t index);
+
+ /// \brief Swaps the specified row with the last row, and then
+ /// releases it from the dataset.
+ ///
+ /// If this dataset does not have its own heap, then you must delete
+ /// the row this returns
+ double* releaseRow(size_t index);
+
+ /// \brief Swaps the specified row with the last row, and then deletes it.
+ void deleteRow(size_t index);
+
+ /// \brief Releases the specified row from the dataset and shifts
+ /// everything after it up one slot.
+ ///
+ /// If this dataset does not have its own heap, then you must delete
+ /// the row this returns
+ double* releaseRowPreserveOrder(size_t index);
+
+ /// \brief Deletes the specified row and shifts everything after it
+ /// up one slot
+ void deleteRowPreserveOrder(size_t index);
+
+ /// \brief Replaces any occurrences of NAN in the matrix with the
+ /// corresponding values from an identity matrix.
+ void fixNans();
+
+ /// \brief Deletes all the data
+ void flush();
+
+ /// \brief Abandons (leaks) all the rows of data
+ void releaseAllRows();
+
+ /// \brief Randomizes the order of the rows.
+ ///
+ /// If pExtension is non-NULL, then it will also be shuffled such
+ /// that corresponding rows are preserved.
+ void shuffle(GRand& rand, GMatrix* pExtension = NULL);
+
+ /// \brief Shuffles the order of the rows. Also shuffles the rows in
+ /// "other" in the same way, such that corresponding rows are
+ /// preserved.
+ void shuffle2(GRand& rand, GMatrix& other);
+
+ /// \brief This is an inferior way to shuffle the data
+ void shuffleLikeCards();
+
+ /// \brief Sorts the data from smallest to largest in the specified
+ /// dimension
+ void sort(size_t nDimension);
+
+ /// This partially sorts the specified column, such that the specified row
+ /// will contain the same row as if it were fully sorted, and previous
+ /// rows will contain a value <= to it in that column, and later rows
+ /// will contain a value >= to it in that column. Unlike sort, which
+ /// has O(m*log(m)) complexity, this method has O(m) complexity. This might
+ /// be useful, for example, for efficiently finding the row with a median
+ /// value in some attribute, or for separating data by a threshold in
+ /// some value.
+ void sortPartial(size_t row, size_t col);
+
+ /// \brief Reverses the row order
+ void reverseRows();
+
+ /// \brief Sorts rows according to the specified compare
+ /// function. (Return true to indicate that the first row comes
+ /// before the second row.)
+ template
+ void sort(CompareFunc& pComparator)
+ {
+ std::sort(m_rows.begin(), m_rows.end(), pComparator);
+ }
+
+ /// \brief Splits this set of data into two sets. Values
+ /// greater-than-or-equal-to dPivot stay in this data set. Values
+ /// less than dPivot go into pLessThanPivot
+ ///
+ /// If pExtensionA is non-NULL, then it will also split pExtensionA
+ /// such that corresponding rows are preserved.
+ void splitByPivot(GMatrix* pGreaterOrEqual, size_t nAttribute, double dPivot, GMatrix* pExtensionA = NULL, GMatrix* pExtensionB = NULL);
+
+ /// \brief Moves all rows with the specified value in the specified
+ /// attribute into pSingleClass
+ ///
+ /// If pExtensionA is non-NULL, then it will also split pExtensionA
+ /// such that corresponding rows are preserved.
+ void splitByNominalValue(GMatrix* pSingleClass, size_t nAttr, int nValue, GMatrix* pExtensionA = NULL, GMatrix* pExtensionB = NULL);
+
+ /// \brief Removes the last nOtherRows rows from this data set and
+ /// puts them in pOtherData
+ void splitBySize(GMatrix* pOtherData, size_t nOtherRows);
+
+ /// \brief Measures the entropy of the specified attribute
+ double entropy(size_t nColumn);
+
+ /// \brief Finds the min and the range of the values of the
+ /// specified attribute
+ void minAndRange(size_t nAttribute, double* pMin, double* pRange);
+
+ /// \brief Estimates the actual min and range based on a random sample
+ void minAndRangeUnbiased(size_t nAttribute, double* pMin, double* pRange);
+
+ /// \brief Shifts the data such that the mean occurs at the origin.
+ /// Only continuous values are affected. Nominal values are left
+ /// unchanged.
+ void centerMeanAtOrigin();
+
+ /// \brief Computes the arithmetic mean of the values in the
+ /// specified column
+ double mean(size_t nAttribute);
+
+
+ /// \brief Computes the arithmetic means of all attributes
+ void centroid(double* pOutCentroid);
+
+ /// \brief Computes the average variance of a single attribute
+ double variance(size_t nAttr, double mean);
+
+ /// \brief Normalizes the specified attribute values
+ void normalize(size_t nAttribute, double dInputMin, double dInputRange, double dOutputMin, double dOutputRange);
+
+ /// \brief Normalize a value from the input min and range to the
+ /// output min and range
+ static double normalize(double dVal, double dInputMin, double dInputRange, double dOutputMin, double dOutputRange);
+
+ /// \brief Returns the mean if the specified attribute is
+ /// continuous, otherwise returns the most common nominal value in
+ /// the attribute.
+ double baselineValue(size_t nAttribute);
+
+ /// \brief Returns true iff the specified attribute contains
+ /// homogenous values. (Unknowns are counted as homogenous with
+ /// anything)
+ bool isAttrHomogenous(size_t col);
+
+ /// \brief Returns true iff each of the last labelDims columns in
+ /// the data are homogenous
+ bool isHomogenous();
+
+ /// \brief Replace missing values with the appropriate measure of
+ /// central tendency.
+ ///
+ /// If the specified attribute is continuous, replaces all
+ /// missing values in that attribute with the mean. If the
+ /// specified attribute is nominal, replaces all missing values in
+ /// that attribute with the most common value.
+ void replaceMissingValuesWithBaseline(size_t nAttr);
+
+ /// \brief Replaces all missing values by copying a randomly
+ /// selected non-missing value in the same attribute.
+ void replaceMissingValuesRandomly(size_t nAttr, GRand* pRand);
+
+ /// \brief This is an efficient algorithm for iteratively computing
+ /// the principal component vector (the eigenvector of the
+ /// covariance matrix) of the data.
+ ///
+ /// See "EM Algorithms for PCA and SPCA"
+ /// by Sam Roweis, 1998 NIPS.
+ ///
+ /// The size of pOutVector will be the number of columns in this matrix.
+ /// (To compute the next principal component, call RemoveComponent,
+ /// then call this method again.)
+ void principalComponent(double* pOutVector, const double* pMean, GRand* pRand);
+
+ /// \brief Computes the first principal component assuming the mean
+ /// is already subtracted out of the data
+ void principalComponentAboutOrigin(double* pOutVector, GRand* pRand);
+
+ /// \brief Computes principal components, while ignoring missing
+ /// values
+ void principalComponentIgnoreUnknowns(double* pOutVector, const double* pMean, GRand* pRand);
+
+ /// \brief Computes the first principal component of the data with
+ /// each row weighted according to the vector pWeights. (pWeights
+ /// must have an element for each row.)
+ void weightedPrincipalComponent(double* pOutVector, const double* pMean, const double* pWeights, GRand* pRand);
+
+ /// \brief Computes the eigenvalue that corresponds to \a *pEigenvector.
+ ///
+ /// After you compute the principal component, you can call this to
+ /// obtain the eigenvalue that corresponds to that principal
+ /// component vector (eigenvector).
+ double eigenValue(const double* pMean, const double* pEigenVector, GRand* pRand);
+
+ /// \brief Removes the component specified by pComponent from the
+ /// data. (pComponent should already be normalized.)
+ ///
+ /// This might be useful, for example, to remove the first principal
+ /// component from the data so you can then proceed to compute the
+ /// second principal component, and so forth.
+ void removeComponent(const double* pMean, const double* pComponent);
+
+ /// \brief Removes the specified component assuming the mean is zero.
+ void removeComponentAboutOrigin(const double* pComponent);
+
+ /// \brief Computes the minimum number of principal components
+ /// necessary so that less than the specified portion of the
+ /// deviation in the data is unaccounted for.
+ ///
+ /// For example, if the data projected onto the first 3 principal
+ /// components contains 90 percent of the deviation that the
+ /// original data contains, then if you pass the value 0.1 to this
+ /// method, it will return 3.
+ size_t countPrincipalComponents(double d, GRand* pRand);
+
+ /// \brief Computes the sum-squared distance between pPoint and all
+ /// of the points in the dataset.
+ ///
+ /// If pPoint is NULL, it computes the sum-squared distance with the origin.
+ ///
+ /// \note that this is equal to the sum of all the eigenvalues times
+ /// the number of dimensions, so you can efficiently compute
+ /// eigenvalues as the difference in sumSquaredDistance with
+ /// the mean after removing the corresponding component, and
+ /// then dividing by the number of dimensions. This is more
+ /// efficient than calling eigenValue.
+ double sumSquaredDistance(const double* pPoint);
+
+ /// \brief Computes the sum-squared distance between the specified
+ /// column of this and that. If the column is a nominal attribute,
+ /// then Hamming distance is used.
+ double columnSumSquaredDifference(GMatrix& that, size_t col);
+
+ /// \brief Computes the squared distance between this and that.
+ ///
+ /// If transpose is true, computes the difference between this and
+ /// the transpose of that.
+ double sumSquaredDifference(GMatrix& that, bool transpose = false);
+
+ /// \brief Computes the linear coefficient between the two specified
+ /// attributes.
+ ///
+ /// Usually you will want to pass the mean values for attr1Origin
+ /// and attr2Origin.
+ double linearCorrelationCoefficient(size_t attr1, double attr1Origin, size_t attr2, double attr2Origin);
+
+ /// \brief Computes the covariance between two attributes
+ double covariance(size_t nAttr1, double dMean1, size_t nAttr2, double dMean2);
+
+ /// \brief Computes the covariance matrix of the data
+ GMatrix* covarianceMatrix();
+
+ /// \brief Performs a paired T-Test with data from the two specified
+ /// attributes.
+ ///
+ /// pOutV will hold the degrees of freedom. pOutT will hold the T-value.
+ /// You can use GMath::tTestAlphaValue to convert these to a P-value.
+ void pairedTTest(size_t* pOutV, double* pOutT, size_t attr1, size_t attr2, bool normalize);
+
+ /// \brief Performs the Wilcoxon signed ranks test from the two
+ /// specified attributes.
+ ///
+ /// If two values are closer than tolerance, they are considered to
+ /// be equal.
+ void wilcoxonSignedRanksTest(size_t attr1, size_t attr2, double tolerance, int* pNum, double* pWMinus, double* pWPlus);
+
+ /// \brief Prints the data to the specified stream
+ void print(std::ostream& stream);
+
+ /// \brief Returns the number of ocurrences of the specified value
+ /// in the specified attribute
+ size_t countValue(size_t attribute, double value);
+
+ /// \brief Returns true iff this matrix is missing any values.
+ bool doesHaveAnyMissingValues();
+
+ /// \brief Throws an exception if this data contains any missing
+ /// values in a continuous attribute
+ void ensureDataHasNoMissingReals();
+
+ /// \brief Throws an exception if this data contains any missing
+ /// values in a nominal attribute
+ void ensureDataHasNoMissingNominals();
+
+ /// \brief Computes the sum entropy of the data (or the sum variance
+ /// for continuous attributes)
+ double measureInfo();
+
+ /// \brief Computes the vector in this subspace that has the
+ /// greatest distance from its projection into pThat subspace.
+ ///
+ /// Returns true if the results are computed.
+ ///
+ /// Returns false if the subspaces are so nearly parallel that pOut
+ /// cannot be computed with accuracy.
+ bool leastCorrelatedVector(double* pOut, GMatrix* pThat, GRand* pRand);
+
+ /// \brief Computes the cosine of the dihedral angle between this
+ /// subspace and pThat subspace
+ double dihedralCorrelation(GMatrix* pThat, GRand* pRand);
+
+ /// \brief Projects pPoint onto this hyperplane (where each row
+ /// defines one of the orthonormal basis vectors of this hyperplane)
+ ///
+ /// This computes (A^T)Ap, where A is this matrix, and p is pPoint.
+ void project(double* pDest, const double* pPoint);
+
+ /// \brief Projects pPoint onto this hyperplane (where each row
+ /// defines one of the orthonormal basis vectors of this hyperplane)
+ void project(double* pDest, const double* pPoint, const double* pOrigin);
+
+ /// \brief Performs a bipartite matching between the rows of \a a
+ /// and \a b using the Linear Assignment Problem (LAP) optimizer
+ ///
+ /// Treats the rows of the matrices \a a and \a b as vectors and
+ /// calculates the distances between these vectors using \a cost.
+ /// Returns an optimal assignment from rows of \a a to rows of \a b
+ /// that minimizes sum of the costs of the assignments.
+ ///
+ /// Each row is considered to be a vector in multidimensional space.
+ /// The cost is the distance given by \a cost when called on each
+ /// row of \a a and row of \a b in turn. The cost must not be
+ /// \f$-\infty\f$ for any pair of rows. Other than that, there are no
+ /// limitations on the cost function.
+ ///
+ /// Because of the limitations of GDistanceMetric, \a a and \a b
+ /// must have the same number of columns.
+ ///
+ /// If \f$m\f$ is \f$\max(rows(a), rows(b))\f$ then this routine
+ /// requires \f$\Theta(rows(a) \cdot rows(b))\f$ memory and
+ /// \f$O(m^3)\f$ time.
+ ///
+ /// \param a the matrix containing the vectors of set a. Must have
+ /// the same number of columns as the matrix containing the
+ /// vectors of set b. Each row is considered to be a
+ /// vector in multidimensional space.
+ ///
+ /// \param b the matrix containing the vectors of set b. Must have
+ /// the same number of columns as the matrix containing the
+ /// vectors of set a. Each row is considered to be a
+ /// vector in multidimensional space.
+ ///
+ /// \param metric given a row of \a a and a row of \a b, returns the
+ /// cost of assigning a to b.
+ ///
+ /// \return the optimal assignment in which each of the rows of \a a
+ /// or \a b (whichever has fewer rows) is assigned to a row
+ /// of the other matrix
+ static GSimpleAssignment bipartiteMatching(GMatrix& a, GMatrix& b, GDistanceMetric& metric);
+
+protected:
+ double determinantHelper(size_t nEndRow, size_t* pColumnList);
+ void inPlaceSquareTranspose();
+ void singularValueDecompositionHelper(GMatrix** ppU, double** ppDiag, GMatrix** ppV, bool throwIfNoConverge, size_t maxIters);
+};
+
+
+
+/// \brief This is a special holder that guarantees the data set will
+/// release all of its data before it is deleted
+class GReleaseDataHolder
+{
+protected:
+ GMatrix* m_pData;
+
+public:
+ GReleaseDataHolder(GMatrix* pData)
+ {
+ m_pData = pData;
+ }
+
+ ~GReleaseDataHolder()
+ {
+ m_pData->releaseAllRows();
+ }
+};
+
+/// \brief This class guarantees that the rows in b are merged
+/// vertically back into a when this object goes out of scope.
+class GMergeDataHolder
+{
+protected:
+ GMatrix& m_a;
+ GMatrix& m_b;
+
+public:
+ GMergeDataHolder(GMatrix& a, GMatrix& b) : m_a(a), m_b(b) {}
+ ~GMergeDataHolder()
+ {
+ m_a.mergeVert(&m_b);
+ }
+};
+
+
+/// \brief Represents an array of matrices or datasets that all have
+/// the same number of columns.
+class GMatrixArray
+{
+protected:
+ sp_relation m_pRelation;
+ std::vector m_sets;
+
+public:
+ GMatrixArray(sp_relation& pRelation);
+ GMatrixArray(size_t cols);
+ ~GMatrixArray();
+ std::vector& sets() { return m_sets; }
+
+ /// \brief Adds a new dataset to the array and preallocates the
+ /// specified number of rows
+ GMatrix* newSet(size_t rows);
+
+ /// \brief Adds count new datasets to the array, and preallocates
+ /// the specified number of rows in each one
+ void newSets(size_t count, size_t rows);
+
+ /// \brief Deletes all the datasets
+ void flush();
+
+ /// \brief Returns the index of the largest data set
+ size_t largestSet();
+
+ /// \brief Returns the number of data sets with zero rows
+ size_t countEmptySets();
+};
+
+} // namespace GClasses
+
+#endif // __GMATRIX_H__
diff --git a/src/ai/waffles/GNeuralNet.cpp b/src/ai/waffles/GNeuralNet.cpp
new file mode 100644
index 0000000..641aace
--- /dev/null
+++ b/src/ai/waffles/GNeuralNet.cpp
@@ -0,0 +1,1266 @@
+/*
+ Copyright (C) 2006, Mike Gashler
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ see http://www.gnu.org/copyleft/lesser.html
+*/
+
+#include "GNeuralNet.h"
+#include "GActivation.h"
+#include "GError.h"
+#include "GRand.h"
+#include "GVec.h"
+#include "GDom.h"
+#include "GTransform.h"
+
+using std::vector;
+
+namespace GClasses {
+
+void GNeuron::resetWeights(GRand* pRand, double inputCenter)
+{
+ for(vector::iterator weight = m_weights.begin(); weight != m_weights.end(); weight++)
+ *weight = pRand->normal() * 0.1;
+
+ // Remove all bias (todo: this has very little effect since the weights are small already--why even bother?)
+ double& bias = m_weights[0];
+ for(vector::iterator weight = m_weights.begin() + 1; weight != m_weights.end(); weight++)
+ bias -= inputCenter * (*weight);
+}
+
+// ----------------------------------------------------------------------
+
+void GNeuralNetLayer::resetWeights(GRand* pRand, double inputCenter)
+{
+ for(vector::iterator neuron = m_neurons.begin(); neuron != m_neurons.end(); neuron++)
+ neuron->resetWeights(pRand, inputCenter);
+}
+
+// ----------------------------------------------------------------------
+
+GBackProp::GBackProp(GNeuralNet* pNN)
+: m_pNN(pNN)
+{
+ // Initialize structures to mirror the neural network
+ m_layers.resize(m_pNN->m_layers.size());
+ for(size_t i = 0; i < m_layers.size(); i++)
+ {
+ GBackPropLayer& layer = m_layers[i];
+ layer.m_neurons.resize(pNN->m_layers[i].m_neurons.size());
+ for(size_t j = 0; j < layer.m_neurons.size(); j++)
+ {
+ GBackPropNeuron& neuron = layer.m_neurons[j];
+ neuron.m_weights.resize(pNN->m_layers[i].m_neurons[j].m_weights.size());
+ }
+ }
+}
+
+void GBackProp::backPropLayer(GNeuralNetLayer* pNNFromLayer, GNeuralNetLayer* pNNToLayer, GBackPropLayer* pBPFromLayer, GBackPropLayer* pBPToLayer, size_t fromBegin)
+{
+ vector::iterator nnFrom, nnCur;
+ vector::iterator bpFrom, bpCur;
+ vector::iterator nn_w;
+
+ // Sum the error times weight for all the children
+ nnFrom = pNNFromLayer->m_neurons.begin();
+ bpFrom = pBPFromLayer->m_neurons.begin();
+ nn_w = nnFrom->m_weights.begin() + 1 + fromBegin;
+ bpCur = pBPToLayer->m_neurons.begin();
+ while(bpCur != pBPToLayer->m_neurons.end())
+ {
+ bpCur->m_error = (*nn_w) * bpFrom->m_error; // use "=" for first pass
+ nn_w++;
+ bpCur++;
+ }
+ nnFrom++;
+ bpFrom++;
+ while(bpFrom != pBPFromLayer->m_neurons.end())
+ {
+ nn_w = nnFrom->m_weights.begin() + 1 + fromBegin;
+ bpCur = pBPToLayer->m_neurons.begin();
+ while(bpCur != pBPToLayer->m_neurons.end())
+ {
+ bpCur->m_error += (*nn_w) * bpFrom->m_error; // use "+=" for subsequent passes
+ nn_w++;
+ bpCur++;
+ }
+ nnFrom++;
+ bpFrom++;
+ }
+
+ // Multiply by the derivative of the activation function
+ nnCur = pNNToLayer->m_neurons.begin();
+ bpCur = pBPToLayer->m_neurons.begin();
+ while(bpCur != pBPToLayer->m_neurons.end())
+ {
+ bpCur->m_error *= pNNToLayer->m_pActivationFunction->derivativeOfNet(nnCur->m_net, nnCur->m_activation);
+ nnCur++;
+ bpCur++;
+ }
+}
+
+void GBackProp::backPropFromSingleNode(GNeuron& nnFrom, GBackPropNeuron& bpFrom, GNeuralNetLayer* pNNToLayer, GBackPropLayer* pBPToLayer)
+{
+ vector::iterator nnCur;
+ vector::iterator bpCur;
+ vector::iterator nn_w;
+
+ // Sum the error times weight for all the children
+ nn_w = nnFrom.m_weights.begin() + 1;
+ bpCur = pBPToLayer->m_neurons.begin();
+ while(bpCur != pBPToLayer->m_neurons.end())
+ {
+ bpCur->m_error = (*nn_w) * bpFrom.m_error;
+ nn_w++;
+ bpCur++;
+ }
+
+ // Multiply by the derivative of the activation function
+ nnCur = pNNToLayer->m_neurons.begin();
+ bpCur = pBPToLayer->m_neurons.begin();
+ while(bpCur != pBPToLayer->m_neurons.end())
+ {
+ bpCur->m_error *= pNNToLayer->m_pActivationFunction->derivativeOfNet(nnCur->m_net, nnCur->m_activation);
+ nnCur++;
+ bpCur++;
+ }
+}
+
+void GBackProp::backPropLayer2(GNeuralNetLayer* pNNFromLayer1, GNeuralNetLayer* pNNFromLayer2, GNeuralNetLayer* pNNToLayer, GBackPropLayer* pBPFromLayer1, GBackPropLayer* pBPFromLayer2, GBackPropLayer* pBPToLayer, size_t pass)
+{
+ vector::iterator nnFrom, nnCur;
+ vector::iterator bpFrom, bpCur;
+ vector::iterator nn_w;
+
+ double sum, alpha;
+ int w = 1;
+ nnCur = pNNToLayer->m_neurons.begin();
+ bpCur = pBPToLayer->m_neurons.begin();
+ while(bpCur != pBPToLayer->m_neurons.end())
+ {
+ // Sum error from the first previous layer
+ sum = 0;
+ nnFrom = pNNFromLayer1->m_neurons.begin();
+ bpFrom = pBPFromLayer1->m_neurons.begin();
+ while(bpFrom != pBPFromLayer1->m_neurons.end())
+ {
+ sum += nnFrom->m_weights[w] * bpFrom->m_error;
+ nnFrom++;
+ bpFrom++;
+ }
+
+ // Sum error from the second previous layer
+ if(pNNFromLayer2)
+ {
+ nnFrom = pNNFromLayer2->m_neurons.begin();
+ bpFrom = pBPFromLayer2->m_neurons.begin();
+ while(bpFrom != pBPFromLayer2->m_neurons.end())
+ {
+ sum += nnFrom->m_weights[w] * bpFrom->m_error;
+ nnFrom++;
+ bpFrom++;
+ }
+ }
+
+ // Multiply by derivative of activation function
+ sum *= pNNToLayer->m_pActivationFunction->derivativeOfNet(nnCur->m_net, nnCur->m_activation);
+
+ // Average with error computed from previous passes
+ alpha = 1.0 / pass;
+ bpCur->m_error *= (1.0 - alpha);
+ bpCur->m_error += (alpha * sum);
+
+ nnCur++;
+ bpCur++;
+ w++;
+ }
+}
+
+void GBackProp::adjustWeights(GNeuralNetLayer* pNNFromLayer, GNeuralNetLayer* pNNToLayer, GBackPropLayer* pBPFromLayer, double learningRate, double momentum)
+{
+ vector::iterator nnFrom;
+ vector::iterator bpFrom;
+ nnFrom = pNNFromLayer->m_neurons.begin();
+ bpFrom = pBPFromLayer->m_neurons.begin();
+ while(bpFrom != pBPFromLayer->m_neurons.end())
+ {
+ vector~~