manifold.cpp

一个非常有用的开源代码

// ********************************************************
// This is demo code. You may derive from, use, modify, and
// distribute it without limitation for any purpose.
// Obviously you don't get a warranty or an assurance of
// fitness for a particular purpose with this code. Your
// welcome to remove this header and claim original
// authorship. I really don't care.
// ********************************************************

#include "Manifold.h"
#include "../GClasses/GManifold.h"
#include "../GClasses/GArff.h"
#include "../GClasses/GBits.h"
#include "../GClasses/GFile.h"
#include "../GClasses/GTime.h"
#include "../GClasses/GThread.h"
#include "../GClasses/GVector.h"
#include <math.h>
/*
void ReadFileLine(char* pBuf, int nSize, FILE* pFile)
{
	while(nSize > 0)
	{
		int n = fgetc(pFile);
		*pBuf = n;
		pBuf++;
		if(n == '\n')
			break;
		nSize--;
	}
}

// Load points from LLE output file
FILE* pLLEFile = fopen("hlleoutput.dat", "r");
char* pLine = new char[nDataPoints * 20];
Holder<char*> hLine(pLine);
GAssert(pLLEFile, "failed to open lle.dat file");
int d;
for(d = 0; d < 2; d++)
{
	ReadFileLine(pLine, nDataPoints * 20 - 3, pLLEFile);
	GAssert(pLine, "failed to read line");
	char* pData = pLine;
	int i;
	for(i = 0; i < nDataPoints; i++)
	{
		while(*pData > '\0' && *pData <= ' ')
			pData++;
		double val = 50 * atof(pData);
		//printf("%f\n", val);
		m_pSculpter->GetDataPoint(i)[1-d] = val;
		while(*pData > ' ')
			pData++;
	}
}
int i;
for(i = 0; i < nDataPoints; i++)
{
	m_pSculpter->GetDataPoint(i)[2] = 0;
}
*/

class ManifoldModel
{
protected:
	GArffRelation* m_pRelation;
	GArffData* m_pData;
	GManifoldSculpting* m_pSculpter;

public:
	ManifoldModel(bool bMask)
	{
		// Make the relation
		m_pRelation = new GArffRelation();
		m_pRelation->AddAttribute(new GArffAttribute(true, 0, NULL)); // x
		m_pRelation->AddAttribute(new GArffAttribute(true, 0, NULL)); // y
		m_pRelation->AddAttribute(new GArffAttribute(true, 0, NULL)); // z
	
		// Make the swiss roll data set
		srand(30); // so we always get consistent results

		// Load the image mask (if necessary)
		GImage imageMask;
bMask = false;
		if(bMask)
		{
			if(!imageMask.LoadPNGFile("mask.png"))
				GAssert(false, "failed to load mask");
		}

		// Make the ARFF data
		double t;
		int n;
		int nDataPoints = 2000;
		m_pData = new GArffData(nDataPoints);
		for(n = 0; n < nDataPoints; n++)
		{
			t = ((double)n * 8) / nDataPoints;
			while(true)
			{
				double* pVector = new double[3];
				pVector[0] = (t + 2) * sin(t) + 14;
				pVector[1] = ((double)GBits::GetRandomUint() / 0xffffffff) * 12 - 6;
				pVector[2] = (t + 2) * cos(t);
				m_pData->AddVector(pVector);
				if(bMask)
				{
					int x = (int)(n * imageMask.GetWidth() / nDataPoints);
					int y = (int)((pVector[1] + 6) * imageMask.GetHeight() / 12);
					GColor c = imageMask.SafeGetPixel(x, y);
					if(gGreen(c) < 128)
						break;
				}
				else
					break;
			}
		}
	
		// Allocate the squisher
		int nNeighbors = 40;
		m_pSculpter = new GManifoldSculpting(m_pData->GetSize(), m_pRelation->GetAttributeCount(), nNeighbors);
		m_pSculpter->SetData(m_pRelation, m_pData);
		m_pSculpter->SetSquishingRate(.99);
		m_pSculpter->SetSmoothingAdvantage(10);
	
		// Init the squisher
		m_pSculpter->SquishBegin(/*nTargetDimensions*/2);	
	}

	~ManifoldModel()
	{
		delete(m_pSculpter);
		delete(m_pData);
		delete(m_pRelation);
	}

	double* GetPoint(int n)
	{
		return m_pSculpter->GetDataPoint(n);
	}

	GManifoldSculpting* GetSculpter()
	{
		return m_pSculpter;
	}

	void PreProcess(bool bLLE)
	{
		if(bLLE)
		{
			// Compute data range
			double dTmp, dInputRangeX, dInputRangeY, dLLERangeX, dLLERangeY;
			m_pData->GetMinAndRange(0, &dTmp, &dInputRangeX);
			m_pData->GetMinAndRange(1, &dTmp, &dInputRangeY);

			printf("Pre-processing with LLE...\n");
			GArffData* pPreprocessedData = NULL;
			int nLen;
			char* pFile = GFile::LoadFileToBuffer("lledata.arff", &nLen);
			if(pFile)
			{
				// Load the LLE data from a file (since my LLE implementation is so slow)
				Holder<char*> hFile(pFile);
				GArffRelation* pTmpRelation = GArffRelation::ParseFile(&pPreprocessedData, pFile, nLen);
				Holder<GArffRelation*> hTmpRelation(pTmpRelation);
				GAssert(pTmpRelation && pPreprocessedData, "failed to parse lle file");
			}
			else
			{
				// Compute the LLE transformation of the data
				pPreprocessedData = GLLE::DoLLE(m_pRelation, m_pData, 14);
				m_pRelation->SaveArffFile(pPreprocessedData, "lledata.arff");
			}
			Holder<GArffData*> hPreprocessedData(pPreprocessedData);

			// Normalize the data
			pPreprocessedData->GetMinAndRange(0, &dTmp, &dLLERangeX);
			pPreprocessedData->GetMinAndRange(1, &dTmp, &dLLERangeY);
			pPreprocessedData->Normalize(0, 0, dLLERangeX, 0, dInputRangeX);
			pPreprocessedData->Normalize(1, 0, dLLERangeY, 0, dInputRangeY);

			// Set the new data
			m_pSculpter->SetData(m_pRelation, pPreprocessedData);
			printf("Done\n");
		}
		else
		{
			printf("Done\n");
			GVector eigenValues;
			GArffData* pPreprocessedData = GPCA::DoPCA(m_pRelation, m_pData, &eigenValues);
			Holder<GArffData*> hPreprocessedData(pPreprocessedData);
			m_pSculpter->SetData(m_pRelation, pPreprocessedData);
		}
	}

	void DoSemiSupervisedThing()
	{
		// Make another swiss roll
		int nDataPoints = 2000;
		int nNeighbors = 24;
		GManifoldSculpting* pSculpter = new GManifoldSculpting(nDataPoints, 3, nNeighbors);
		pSculpter->SetSquishingRate(.99);
		pSculpter->SetSmoothingAdvantage(10);
		srand(0); // Make sure we always get consistent results
		double values[3];
		double t;
		int n;
		for(n = 0; n < nDataPoints; n++)
		{
			t = ((double)n * 8) / nDataPoints;
			values[0] = (t + 2) * sin(t) + 14;
			values[1] = ((double)GBits::GetRandomUint() / 0xffffffff) * 12 - 6;
			values[2] = (t + 2) * cos(t);
			pSculpter->SetDataPoint(n, values, true);
		}

		// Init the squisher
		pSculpter->SquishBegin(2);

		// Set some supervised points
		for(n = 0; n < nDataPoints; n++)
		{
			if(rand() % 20 == 0)
				pSculpter->SetDataPoint(n, m_pSculpter->GetDataPoint(n), false);
		}

		// Swap in the new squisher
		delete(m_pSculpter);
		m_pSculpter = pSculpter;
	}

	void ToMatrix(const char* szFilename)
	{
		FILE* pFile = fopen(szFilename, "w");
		GAssert(pFile, "failed to create file");
		int nDataPoints = m_pData->GetSize();
		int dim, n;
		for(dim = 0; dim < 3; dim++)
		{
			for(n = 0; n < nDataPoints; n++)
			{
				if(n > 0)
					fprintf(pFile, "\t");
				fprintf(pFile, "%f", m_pSculpter->GetDataPoint(n)[dim]);
			}
			fprintf(pFile, "\n");
		}
		fclose(pFile);
	}

	void FromMatrix(const char* szFilename)
	{
		FILE* pFile = fopen(szFilename, "r");
		GAssert(pFile, "failed to open the file");
		FileHolder hFile(pFile);
		int nDataPoints = m_pData->GetSize();
		char* pBuf = new char[nDataPoints * 20];
		Holder<char*> hLine(pBuf);
		char* pData;
		int d;
		for(d = 0; d < 2; d++)
		{
			pData = fgets(pBuf, nDataPoints * 20, pFile);
			GAssert(pData, "failed to read line");
			if(!pData)
				break;
			int i;
			for(i = 0; i < nDataPoints; i++)
			{
				while(*pData > '\0' && *pData <= ' ')
				pData++;
				double val = atof(pData);
				//printf("%f\n", val);
				m_pSculpter->GetDataPoint(i)[d] = val;
				while(*pData > ' ')
					pData++;
			}
		}
		int i;
		for(i = 0; i < nDataPoints; i++)
		{
			m_pSculpter->GetDataPoint(i)[2] = 0;
		}
	}
};



// -------------------------------------------------------------------------------

class ManifoldView : public ViewBase
{
protected:
	GRect m_viewRect, m_nextRect;
	ManifoldModel* m_pModel;

public:
	ManifoldView(ManifoldModel* pModel);
	virtual ~ManifoldView();

	virtual void OnChar(char c);
	virtual void OnMouseDown(int x, int y);
	virtual void OnMouseUp(int x, int y);
	virtual bool OnMousePos(int x, int y);

	void DoSemiSupervisedThing();
	void PreProcess(bool bLLE);

protected:
	virtual void Draw(SDL_Surface *pScreen);
	void DrawPoint(SDL_Surface *pScreen, int n);
};

ManifoldView::ManifoldView(ManifoldModel* pModel)
: ViewBase()
{
	m_pModel = pModel;

	// Set the view rect
	m_nextRect.Set(-30, -25, 60, 50);
	//m_nextRect.Set(-2, -6, 5, 12);
	m_viewRect = m_nextRect;
}

ManifoldView::~ManifoldView()
{
}

void ManifoldView::DrawPoint(SDL_Surface *pScreen, int n)
{
	double* pValues = m_pModel->GetPoint(n);
	int x = (int)((pValues[0] - m_viewRect.x) * m_screenRect.w / m_viewRect.w) + m_screenRect.x;
	int y = (int)((pValues[1] + .5 * pValues[2] - m_viewRect.y) * m_screenRect.h / m_viewRect.h) + m_screenRect.y;
	unsigned int col = (unsigned int)GetSpectrumColor((float)n / 2000);
	DrawDot(pScreen, x, y, col, 5);
	if(pValues[0] < m_nextRect.x)
	{
		m_nextRect.w += (m_nextRect.x - (int)pValues[0]);
		m_nextRect.x = (int)pValues[0];
	}
	else if(pValues[0] > m_nextRect.x + m_nextRect.w)
		m_nextRect.w = (int)pValues[0] - m_nextRect.x + 1;
	if(pValues[1] < m_nextRect.y)
	{
		m_nextRect.h += (m_nextRect.y - (int)pValues[1]);
		m_nextRect.y = (int)pValues[1];
	}
	else if(pValues[1] > m_nextRect.y + m_nextRect.h)
		m_nextRect.h = (int)pValues[1] - m_nextRect.y + 1;
}

/*virtual*/ void ManifoldView::Draw(SDL_Surface *pScreen)
{
	// Clear the screen
	SDL_FillRect(pScreen, NULL/*&r*/, 0x000000);
	//m_viewRect = m_nextRect;
	//m_viewRect.Set((m_viewRect.x * 9 + m_nextRect.x) / 10, (m_viewRect.y * 9 + m_nextRect.y) / 10, (m_viewRect.w * 9 + m_nextRect.w) / 10, (m_viewRect.h * 9 + m_nextRect.h) / 10);
	m_viewRect.x += GBits::GetSign(m_nextRect.x - m_viewRect.x);
	m_viewRect.y += GBits::GetSign(m_nextRect.y - m_viewRect.y);
	m_viewRect.w += GBits::GetSign(m_nextRect.w - m_viewRect.w);
	m_viewRect.h += GBits::GetSign(m_nextRect.h - m_viewRect.h);

	// Reset the next rect
	double* pValues = m_pModel->GetPoint(0);
	m_nextRect.Set((int)pValues[0], (int)pValues[1], 1, 1);

	// Plot the points in approximate order from back to front so it looks like we're
	// clipping properly--this is a cheap hacky way to do clipping, but who cares?
	int n;
	for(n = 450; n < 1000; n++)
		DrawPoint(pScreen, n);
	for(n = 0; n < 450; n++)
		DrawPoint(pScreen, n);
	for(n = 1000; n < 2000; n++)