gvector.cpp

来自「一个由Mike Gashler完成的机器学习方面的includes neural」· C++ 代码 · 共 292 行

/*
	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 "GVector.h"
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "GBits.h"
#include "GMacros.h"

// static
double GVector::ComputeDotProduct(const double* pA, const double* pB, int nSize)
{
	int n;
	double dVal = 0;
	for(n = 0; n < nSize; n++)
		dVal += pA[n] * pB[n];
	return dVal;
}

// static
double GVector::ComputeDotProduct(const double* pOrigin, const double* pTarget, const double* pVector, int nSize)
{
	int n;
	double dVal = 0;
	for(n = 0; n < nSize; n++)
		dVal += (pTarget[n] - pOrigin[n]) * pVector[n];
	return dVal;
}

// static
double GVector::ComputeSquaredDistance(double* pA, double* pB, int nSize)
{
	double dist = 0;
	double d;
	int n;
	for(n = 0; n < nSize; n++)
	{
		d = pA[n] - pB[n];
		dist += (d * d);
	}
	return dist;
}

// static
double GVector::ComputeSquaredMagnitude(const double* pVector, int nSize)
{
	double dMag = 0;
	int i;
	for(i = 0; i < nSize; i++)
		dMag += (pVector[i] * pVector[i]);
	return dMag;
}

// static
double GVector::LNormComputeMagnitude(double l, const double* pVector, int nSize)
{
	double dMag = 0;
	int i;
	for(i = 0; i < nSize; i++)
		dMag += pow(ABS(pVector[i]), l);
	return pow(dMag, 1.0 / l);
}

// static
double GVector::ComputeCorrelation(const double* pA, const double* pB, int nDims)
{
	double dDotProd = ComputeDotProduct(pA, pB, nDims);
	if(dDotProd == 0)
		return 0;
	return dDotProd / (sqrt(ComputeSquaredMagnitude(pA, nDims)) * sqrt(ComputeSquaredMagnitude(pB, nDims)));
}

// static
void GVector::Normalize(double* pVector, int nSize)
{
	double dMag = sqrt(ComputeSquaredMagnitude(pVector, nSize));
	int i;
	for(i = 0; i < nSize; i++)
		pVector[i] /= dMag;
}

// static
void GVector::LNormNormalize(double l, double* pVector, int nSize)
{
	double dMag = LNormComputeMagnitude(l, pVector, nSize);
	int i;
	for(i = 0; i < nSize; i++)
		pVector[i] /= dMag;
}

// static
int GVector::FindMax(double* pVector, int nSize)
{
	int nIndex = rand() % nSize;
	int n;
	for(n = 0; n < nSize; n++)
	{
		if(pVector[n] > pVector[nIndex])
			nIndex = n;
	}
	return nIndex;
}

// static
void GVector::GetRandomVector(double* pOutVector, int nSize)
{
	int i;
	for(i = 0; i < nSize; i++)
		pOutVector[i] = GBits::GetRandomGaussian();
	Normalize(pOutVector, nSize);
}

// static
bool GVector::TestPivot(double* pCenter, double* pNormal, double* pVector, int nDims)
{
	int n;
	double dVal = 0;
	for(n = 0; n < nDims; n++)
		dVal += (pVector[n] - pCenter[n]) * pNormal[n];
	return dVal >= 0;
}

// static
void GVector::Add(double* pDest, double* pSource, int nDims)
{
	int i;
	for(i = 0; i < nDims; i++)
		pDest[i] += pSource[i];
}

// static
void GVector::Subtract(double* pDest, double* pSource, int nDims)
{
	int i;
	for(i = 0; i < nDims; i++)
		pDest[i] -= pSource[i];
}

// static
void GVector::Multiply(double* pVector, double dScalar, int nDims)
{
	int i;
	for(i = 0; i < nDims; i++)
		pVector[i] *= dScalar;
}

// static
void GVector::SetToZero(double* pVector, int nDims)
{
	int i;
	for(i = 0; i < nDims; i++)
		pVector[i] = 0;
}

void GVector::GetMostEccentricMaxs(int* pOutPoints, int nPoints, bool bMax, double* pVector, int nDims)
{
	if(nPoints < 1)
		return;

	// Handle the case where we don't have enough points
	int i;
	if(nDims <= nPoints)
	{
		for(i = 0; i < nPoints; i++)
			pOutPoints[i] = i * nDims / nPoints;
		return;
	}

	// Find the most eccentric of the eccentric points until we have fewer than we need
	int nSign = bMax ? 1 : -1;
	GTEMPBUF(int, pWorkBuffers, nDims * 2);
	int* pPrev = &pWorkBuffers[nDims];
	int* pCurrent = pWorkBuffers;
	for(i = 0; i < nDims; i++)
		pCurrent[i] = i;
	int nPointCount = nDims;
	int nPrevPointCount = 0;
	int a, b, c, nPos, maxGap;
	while(nPointCount > nPoints)
	{
		// Swap the buffers
		int* pTemp = pPrev;
		pPrev = pCurrent;
		pCurrent = pTemp;

		// Compute the max gap
		maxGap = nDims * 5 / nPointCount; // 5 is the max allowed gap ratio. Smaller values restrict to more even-ness (and increase computation time)

		// Find all the local maxes
		nPos = 0;
		a = pPrev[0];
		b = pPrev[1];
		if(b - a > maxGap || (pVector[a] - pVector[b]) * nSign > 0)
			pCurrent[nPos++] = a;
		for(i = 1; i < nPointCount - 1; i++)
		{
			a = pPrev[i - 1];
			b = pPrev[i];
			c = pPrev[i + 1];
			if(b - a > maxGap || c - b > maxGap || GBits::GetSign((pVector[b] - pVector[a]) * nSign) + GBits::GetSign((pVector[b] - pVector[c]) * nSign) > 0)
				pCurrent[nPos++] = b;
		}
		a = pPrev[i - 1];
		b = pPrev[i];
		if(b - a > maxGap || (pVector[b] - pVector[a]) * nSign > 0)
			pCurrent[nPos++] = b;

		nPrevPointCount = nPointCount;
		nPointCount = nPos;
	}

	// Use all the points in pCurrent and some of the points in pPrev
	GAssert(nPointCount <= nPoints && nPrevPointCount > nPoints, "unexpected array sizes");
	nPos = 0;
	int nPrevPos = 0;
	int nNeeded;
	for(i = 0; i < nPoints; i++)
	{
		// Skip excess points in pPrev
		nNeeded = nPoints - i;
		while(		(nPos >= nPointCount || pPrev[nPrevPos] != pCurrent[nPos]) &&
				nPrevPointCount - nPrevPos > nNeeded &&
				i * nPrevPointCount / nPoints > nPrevPos
			)
			nPrevPos++;

		// Take the next point
		pOutPoints[i] = pPrev[nPrevPos];
		if(nPos < nPointCount && pCurrent[nPos] == pPrev[nPrevPos])
			nPos++;
		nPrevPos++;
	}
}

void GVector::GetMostEccentricPoints(int* pOutPoints, int nPoints, double* pVector, int nDims)
{
	// Get a set of maxs and a set of mins
	int nMins = nPoints / 2;
	int nMaxs = nMins;
	if(nMaxs + nMins < nPoints)
		nMaxs++;
	GTEMPBUF(int, pMaxs, nPoints);
	int* pMins = &pMaxs[nMaxs];
	GVector::GetMostEccentricMaxs(pMaxs, nMaxs, true, pVector, nDims);
	GVector::GetMostEccentricMaxs(pMins, nMins, false, pVector, nDims);

	// Merge the points
	int i;
	int nMaxPos = 0;
	int nMinPos = 0;
	for(i = 0; i < nPoints; i++)
	{
		if(nMaxPos >= nMaxs)
			pOutPoints[i] = pMins[nMinPos++];
		else if(nMinPos >= nMins)
			pOutPoints[i] = pMaxs[nMaxPos++];
		else if(pMins[nMinPos] < pMaxs[nMaxPos])
			pOutPoints[i] = pMins[nMinPos++];
		else
			pOutPoints[i] = pMaxs[nMaxPos++];
	}
}

void GVector::InterpolateIndexes(int nIndexes, double* pInIndexes, double* pOutIndexes, float fRatio, int nCorrIndexes, double* pCorrIndexes1, double* pCorrIndexes2)
{
	GAssert(nCorrIndexes >= 2, "need at least two correlated indexes (at least the two extremes)");
	int nCorr = 0;
	double fInvRatio = (float)1 - fRatio;
	double fIndex, fWeight, f0, f1;
	int i;
	for(i = 0; i < nIndexes; i++)
	{
		fIndex = pInIndexes[i];
		while(nCorr < nCorrIndexes - 2 && fIndex >= pCorrIndexes1[nCorr + 1])
			nCorr++;
		fWeight = (fIndex - pCorrIndexes1[nCorr]) / (pCorrIndexes1[nCorr + 1] - pCorrIndexes1[nCorr]);
		f0 = fInvRatio * pCorrIndexes1[nCorr] + fRatio * pCorrIndexes2[nCorr];
		f1 = fInvRatio * pCorrIndexes1[nCorr + 1] + fRatio * pCorrIndexes2[nCorr + 1];
		pOutIndexes[i] = ((float)1 - fWeight) * f0 + fWeight * f1;
	}
}