manifold.cpp

来自「一个由Mike Gashler完成的机器学习方面的includes neural」· C++ 代码 · 共 1,682 行 · 第 1/5 页

			searcher.Iterate();
			d = critic.GetBestError();
			if(d < dBest)
			{
				dBest = d;
				i = 0;
			}
			else
			{
				i++;
				if(i > 1000)
					break;
			}
			//printf("Error: %f\n", d);
		}
		double* pBestVec = critic.GetBestYet();
		//printf("Scale=%f, Add=%f\n", pBestVec[0], pBestVec[1]);
		//printf("order errors: %d\n", MeasureSingleDimOrderErrors());
		return critic.GetBestError();
	}
};



double LengthOfSineFunc(double x)
{
	double d = cos(x);
	return sqrt(d * d + 1.0);
}

double LengthOfSwissRoll(double x)
{
#ifdef WIN32
	GAssert(false, "not implemented yet for Win32");
	return 0;
#else
	return (x * sqrt(x * x + 1) + asinh(x)) / 2;
#endif
}

class SwissRollModel : public ManifoldModel
{
protected:

public:
	enum ManifoldType
	{
		SWISS_ROLL,
		S_CURVE,
		SPIRALS,
	};

	SwissRollModel(ManifoldType eType, int nPoints, bool bMask, int nNeighbors, double dSquishingRate, bool bComputeIdeal, int nSupervisedPoints, SwissRollModel* pTrainedModel) : ManifoldModel()
	{
		// 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 ARFF data
		if(bComputeIdeal)
			m_pIdealResults = new double[nPoints * (eType == SPIRALS ? 1 : 2)];
		double t;
		int n;
		m_pData = new GArffData(nPoints);
		if(eType == SWISS_ROLL)
		{
			// Load the image mask (if necessary)
			GImage imageMask;
			if(bMask)
			{
				if(!imageMask.LoadPNGFile("mask.png"))
					GAssert(false, "failed to load mask");
			}

			for(n = 0; n < nPoints; n++)
			{
				t = ((double)n * 8) / nPoints;
				while(true)
				{
					double* pVector = new double[3];
					pVector[0] = (t + 2) * sin(t) + 14;
					pVector[1] = GBits::GetRandomDouble() * 12 - 6;
					pVector[2] = (t + 2) * cos(t);
					m_pData->AddVector(pVector);
					if(bComputeIdeal)
					{
						m_pIdealResults[2 * n] = pVector[1];
						m_pIdealResults[2 * n + 1] = LengthOfSwissRoll(t + 2);/* - LengthOfSwissRoll(2);*/
					}
					if(bMask)
					{
						int x = (int)(n * imageMask.GetWidth() / nPoints);
						int y = (int)((pVector[1] + 6) * imageMask.GetHeight() / 12);
						GColor c = imageMask.SafeGetPixel(x, y);
						if(gGreen(c) < 128)
							break;
					}
					else
						break;
				}
			}
		}
		else if(eType == S_CURVE)
		{
			for(n = 0; n < nPoints; n++)
			{
				t = ((double)n * 2.2 * PI - .1 * PI) / nPoints;
				double* pVector = new double[3];
				pVector[0] = 1.0 - sin(t);
				pVector[1] = t;
				pVector[2] = GBits::GetRandomDouble() * 2;
				m_pData->AddVector(pVector);
				if(bComputeIdeal)
				{
					m_pIdealResults[2 * n] = pVector[2];
					m_pIdealResults[2 * n + 1] = (n > 0 ? GMath::Integrate(LengthOfSineFunc, 0, t, n + 30) : 0);
				}
			}
		}
		else if(eType == SPIRALS)
		{
			double dHeight = 3;
			double dWraps = 1.5;
			double dSpiralLength = sqrt((dWraps * 2.0 * PI) * (dWraps * 2.0 * PI) + dHeight * dHeight);
			double dTotalLength = 2.0 * (dSpiralLength + 1); // radius = 1
			double d;
			for(n = 0; n < nPoints; n++)
			{
				t = ((double)n * dTotalLength) / nPoints;
				double* pVector = new double[3];
				if(t < dSpiralLength)
				{
					d = (dSpiralLength - t) * dWraps * 2 * PI / dSpiralLength; // d = radians
					pVector[0] = -cos(d);
					pVector[1] = dHeight * t / dSpiralLength;
					pVector[2] = -sin(d);
				}
				else if(t - 2.0 - dSpiralLength >= 0)
				{
					d = (t - 2.0 - dSpiralLength) * dWraps * 2 * PI / dSpiralLength; // d = radians
					pVector[0] = cos(d);
					pVector[1] = dHeight * (dSpiralLength - (t - 2.0 - dSpiralLength)) / dSpiralLength;
					pVector[2] = sin(d);
				}
				else
				{
					d = (t - dSpiralLength) / 2.0; // 2 = diameter
					pVector[0] = 2.0 * d - 1.0;
					pVector[1] = dHeight;
					pVector[2] = 0;
				}
				m_pData->AddVector(pVector);
				if(bComputeIdeal)
					m_pIdealResults[n] = dTotalLength * n / nPoints;
			}
		}

		// Allocate the sculpter
		m_pSculpter = new GManifoldSculpting(m_pData->GetSize(), m_pRelation->GetInputCount(), nNeighbors);
		m_pSculpter->SetData(m_pRelation, m_pData);
		m_pSculpter->SetSquishingRate(dSquishingRate);
		m_pSculpter->SetSmoothingAdvantage(10);

		// Init the sculpter
		m_pSculpter->SquishBegin(/*nTargetDimensions*/(eType == SPIRALS ? 1 : 2));

		// Set the supervised points
		if(nSupervisedPoints > 0)
		{
			int i;
			for(i = 0; i < nSupervisedPoints; i++)
			{
				int nPoint = rand() % nPoints;
				m_pSculpter->SetVector(nPoint, pTrainedModel->GetSculpter()->GetVector(nPoint), false);
			}
		}
	}

	virtual ~SwissRollModel()
	{
	}

	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");
			GArffData* pPreprocessedData = GPCA::DoPCA(m_pRelation, m_pData);
			Holder<GArffData*> hPreprocessedData(pPreprocessedData);
			m_pSculpter->SetData(m_pRelation, pPreprocessedData);
		}
*/
	}

	void DoSemiSupervisedThing()
	{
		// Make another swiss roll
		int nNeighbors = 24;
		GManifoldSculpting* pSculpter = new GManifoldSculpting(SWISS_ROLL_POINTS, 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 < SWISS_ROLL_POINTS; n++)
		{
			t = ((double)n * 8) / SWISS_ROLL_POINTS;
			values[0] = (t + 2) * sin(t) + 14;
			values[1] = ((double)GBits::GetRandomUint() / 0xffffffff) * 12 - 6;
			values[2] = (t + 2) * cos(t);
			pSculpter->SetVector(n, values, true);
		}

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

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

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

// virtual
double Compute2DErrorCritic::ComputeError(double* pVector)
{
	return m_pModel->Compute2DError(pVector[0], pVector[1], pVector[2], pVector[3], pVector[4]);
}

// virtual
double Compute1DErrorCritic::ComputeError(double* pVector)
{
	return m_pModel->Compute1DError(pVector[0], pVector[1]);
}


char* MakeEightDigitInt(int n, char* szBuf)
{
	sprintf(szBuf, "%d", n);
	int len = strlen(szBuf);
	int i;
	for(i = 7 - len; i >= 0; i--)
		szBuf[i] = '0';
	sprintf(szBuf + 8 - len, "%d", n);
	return szBuf;
}

class ImageModel : public ManifoldModel
{
protected:
	GImage* m_pImages;

public:
	ImageModel(int nSkip, int nImageCount, int nImageWidth, int nImageHeight, int nTargetDims, int nNeighbors, double dSquishingRate, const char* szFilenamePrefix, int nSupervisedPoints, ImageModel* pTrainedModel, bool bComputeIdeal) : ManifoldModel()
	{
		printf("Image Count=%d, wid=%d, hgt=%d, target dims=%d, neighbors=%d, squishing rate=%f, supervised points=%d\n", nImageCount, nImageWidth, nImageHeight, nTargetDims, nNeighbors, dSquishingRate, nSupervisedPoints);
		int nAttrs;
		m_pImages = new GImage[nImageCount];
		nAttrs = nImageWidth * nImageHeight;

		// Make the relation
		m_pRelation = new GArffRelation();
		int i, x, y;
		for(i = 0; i < nAttrs; i++)
			m_pRelation->AddAttribute(new GArffAttribute(true, 0, NULL));

		// Load the images into the ARFF data
		m_pData = new GArffData(nImageCount);
		const char* szAppPath = ControllerBase::GetAppPath();
		int nLen = strlen(szAppPath);
		GTEMPBUF(char, szFilename, nLen + 32);
		strcpy(szFilename, szAppPath);
		strcat(szFilename, szFilenamePrefix);
		nLen += strlen(szFilenamePrefix);
		double* pVector;
		int nGrayScaleValue;
		GColor col;
		int nSmoothingAdvantage = 10;
		nNeighbors = FACE_NEIGHBORS;
		char szTmp[9];
		for(i = 0; i < nImageCount; i++)
		{
			MakeEightDigitInt(i + 1 + nSkip, szTmp);
			strcpy(szFilename + nLen, szTmp);
			strcpy(szFilename + nLen + 8, ".png");
			if(!m_pImages[i].LoadPNGFile(szFilename))
				throw "Failed to load image";
			GAssert((int)m_pImages[i].GetWidth() == nImageWidth, "unexpected size");
			GAssert((int)m_pImages[i].GetHeight() == nImageHeight, "unexpected size");
			pVector = new double[nImageWidth * nImageHeight];
			m_pData->AddVector(pVector);