C++ GMatrix类代码示例

void lle(GArgReader& args)
{
	// Load the file and params
	GMatrix* pData = loadData(args.pop_string());
	Holder<GMatrix> hData(pData);
	unsigned int nSeed = getpid() * (unsigned int)time(NULL);
	GRand prng(nSeed);
	GNeighborFinder* pNF = instantiateNeighborFinder(pData, &prng, args);
	Holder<GNeighborFinder> hNF(pNF);
	int targetDims = args.pop_uint();

	// Parse Options
	while(args.size() > 0)
	{
		if(args.if_pop("-seed"))
			prng.setSeed(args.pop_uint());
		else
			throw Ex("Invalid option: ", args.peek());
	}

	// Transform the data
	GLLE transform(pNF->neighborCount(), targetDims, &prng);
	transform.setNeighborFinder(pNF);
	GMatrix* pDataAfter = transform.doit(*pData);
	Holder<GMatrix> hDataAfter(pDataAfter);
	pDataAfter->print(cout);
}

void plot_it(const char* filename, GNeuralNet& nn, GMatrix& trainFeat, GMatrix& trainLab, GMatrix& testFeat, GMatrix& testLab)
{
	GSVG svg(1000, 500);
	double xmin = trainFeat[0][0];
	double xmax = testFeat[testFeat.rows() - 1][0];
	svg.newChart(xmin, std::min(trainLab.columnMin(0), testLab.columnMin(0)), xmax, std::max(trainLab.columnMax(0), testLab.columnMax(0)));
	svg.horizMarks(20);
	svg.vertMarks(20);
	double prevx = xmin;
	double prevy = 0.0;
	double step = (xmax - xmin) / 500.0;
	GVec x(1);
	GVec y(1);
	for(x[0] = prevx; x[0] < xmax; x[0] += step)
	{
		nn.predict(x, y);
		if(prevx != x[0])
			svg.line(prevx, prevy, x[0], y[0], 0.3);
		prevx = x[0];
		prevy = y[0];
	}
	for(size_t i = 0; i < trainLab.rows(); i++)
		svg.dot(trainFeat[i][0], trainLab[i][0], 0.4, 0xff000080);
	for(size_t i = 0; i < testLab.rows(); i++)
		svg.dot(testFeat[i][0], testLab[i][0], 0.4, 0xff800000);

	std::ofstream ofs;
	ofs.open(filename);
	svg.print(ofs);
}

void addNoise(GArgReader& args)
{
	GMatrix* pData = loadData(args.pop_string());
	Holder<GMatrix> hData(pData);
	double dev = args.pop_double();

	// Parse the options
	unsigned int seed = getpid() * (unsigned int)time(NULL);
	int excludeLast = 0;
	while(args.next_is_flag())
	{
		if(args.if_pop("-seed"))
			seed = args.pop_uint();
		else if(args.if_pop("-excludelast"))
			excludeLast = args.pop_uint();
		else
			ThrowError("Invalid neighbor finder option: ", args.peek());
	}

	GRand prng(seed);
	size_t cols = pData->cols() - excludeLast;
	for(size_t r = 0; r < pData->rows(); r++)
	{
		double* pRow = pData->row(r);
		for(size_t c = 0; c < cols; c++)
			*(pRow++) += dev * prng.normal();
	}
	pData->print(cout);
}

// virtual
void GLinearRegressor::trainInner(GMatrix& features, GMatrix& labels)
{
	// Use a fast, but not-very-numerically-stable technique to compute an initial approximation for beta and epsilon
	clear();
	GMatrix* pAll = GMatrix::mergeHoriz(&features, &labels);
	Holder<GMatrix> hAll(pAll);
	GPCA pca(features.cols(), &m_rand);
	pca.train(*pAll);
	size_t inputs = features.cols();
	size_t outputs = labels.cols();
	GMatrix f(inputs, inputs);
	GMatrix l(inputs, outputs);
	for(size_t i = 0; i < inputs; i++)
	{
		GVec::copy(f[i], pca.basis(i), inputs);
		double sqmag = GVec::squaredMagnitude(f[i], inputs);
		if(sqmag > 1e-10)
			GVec::multiply(f[i], 1.0 / sqmag, inputs);
		GVec::copy(l[i], pca.basis(i) + inputs, outputs);
	}
	m_pBeta = GMatrix::multiply(l, f, true, false);
	m_pEpsilon = new double[outputs];
	m_pBeta->multiply(pca.mean(), m_pEpsilon, false);
	GVec::multiply(m_pEpsilon, -1.0, outputs);
	GVec::add(m_pEpsilon, pca.mean() + inputs, outputs);

	// Refine the results using gradient descent
	refine(features, labels, 0.06, 20, 0.75);
}

void precisionRecall(GArgReader& args)
{
	// Parse options
	unsigned int seed = getpid() * (unsigned int)time(NULL);
	bool ideal = false;
	while(args.next_is_flag())
	{
		if(args.if_pop("-seed"))
			seed = args.pop_uint();
		else if(args.if_pop("-ideal"))
			ideal = true;
		else
			ThrowError("Invalid option: ", args.peek());
	}

	// Load the data
	if(args.size() < 1)
		ThrowError("No dataset specified.");
	GMatrix* pData = loadData(args.pop_string());
	Holder<GMatrix> hData(pData);

	// Instantiate the recommender
	GRand prng(seed);
	GCollaborativeFilter* pModel = InstantiateAlgorithm(prng, args);
	Holder<GCollaborativeFilter> hModel(pModel);
	if(args.size() > 0)
		ThrowError("Superfluous argument: ", args.peek());

	// Generate precision-recall data
	GMatrix* pResults = pModel->precisionRecall(*pData, ideal);
	Holder<GMatrix> hResults(pResults);
	pResults->deleteColumn(2); // we don't need the false-positive rate column
	pResults->print(cout);
}

void GRecommenderLib::loadData(GMatrix& data, const char* szFilename)
{
	PathData pd;
	GFile::parsePath(szFilename, &pd);
	if(_stricmp(szFilename + pd.extStart, ".sparse") == 0)
	{
		GDom doc;
		doc.loadJson(szFilename);
		GSparseMatrix sm(doc.root());
		data.resize(0, 3);
		for(size_t i = 0; i < sm.rows(); i++)
		{
			GSparseMatrix::Iter rowEnd = sm.rowEnd(i);
			for(GSparseMatrix::Iter it = sm.rowBegin(i); it != rowEnd; it++)
			{
				GVec& vec = data.newRow();
				vec[0] = (double)i;
				vec[1] = (double)it->first;
				vec[2] = it->second;
			}
		}
	}
	else if(_stricmp(szFilename + pd.extStart, ".arff") == 0)
		data.loadArff(szFilename);
	else
		throw Ex("Unsupported file format: ", szFilename + pd.extStart);
}

void OdeImplicitEuler<Real>::Update (Real fTIn, Real* afXIn, Real& rfTOut,
    Real* afXOut)
{
    m_oFunction(fTIn,afXIn,m_pvData,m_kF);
    m_oDFunction(fTIn,afXIn,m_pvData,m_kDF);
    GMatrix<Real> kDG = m_kIdentity - m_fStep*m_kDF;
    GMatrix<Real> kDGInverse(m_iDim,m_iDim);
    bool bInvertible = kDG.GetInverse(kDGInverse);

    if (bInvertible)
    {
        m_kF = kDGInverse*m_kF;
        for (int i = 0; i < m_iDim; i++)
        {
            afXOut[i] = afXIn[i] + m_fStep*m_kF[i];
        }
    }
    else
    {
        size_t uiSize = m_iDim*sizeof(Real);
        System::Memcpy(afXOut,uiSize,afXIn,uiSize);
    }

    rfTOut = fTIn + m_fStep;
}

void kmeans(GArgReader& args)
{
	// Load the file and params
	GMatrix data;
	loadData(data, args.pop_string());
	int clusters = args.pop_uint();

	// Parse Options
	unsigned int nSeed = getpid() * (unsigned int)time(NULL);
	size_t reps = 1;
	while(args.size() > 0)
	{
		if(args.if_pop("-seed"))
			nSeed = args.pop_uint();
		else if(args.if_pop("-reps"))
			reps = args.pop_uint();
		else
			throw Ex("Invalid option: ", args.peek());
	}

	// Do the clustering
	GRand prng(nSeed);
	GKMeans clusterer(clusters, &prng);
	clusterer.setReps(reps);
	GMatrix* pOut = clusterer.reduce(data);
	std::unique_ptr<GMatrix> hOut(pOut);
	pOut->print(cout);
}

// virtual
void GLinearRegressor::trainInner(const GMatrix& features, const GMatrix& labels)
{
	if(!features.relation().areContinuous())
		throw Ex("GLinearRegressor only supports continuous features. Perhaps you should wrap it in a GAutoFilter.");
	if(!labels.relation().areContinuous())
		throw Ex("GLinearRegressor only supports continuous labels. Perhaps you should wrap it in a GAutoFilter.");

	// Use a fast, but not-very-numerically-stable technique to compute an initial approximation for beta and epsilon
	clear();
	GMatrix* pAll = GMatrix::mergeHoriz(&features, &labels);
	Holder<GMatrix> hAll(pAll);
	GPCA pca(features.cols());
	pca.train(*pAll);
	size_t inputs = features.cols();
	size_t outputs = labels.cols();
	GMatrix f(inputs, inputs);
	GMatrix l(inputs, outputs);
	for(size_t i = 0; i < inputs; i++)
	{
		GVec::copy(f[i].data(), pca.basis()->row(i).data(), inputs);
		double sqmag = f[i].squaredMagnitude();
		if(sqmag > 1e-10)
			f[i] *= 1.0 / sqmag;
		l[i].set(pca.basis()->row(i).data() + inputs, outputs);
	}
	m_pBeta = GMatrix::multiply(l, f, true, false);
	m_epsilon.resize(outputs);
	GVecWrapper vw(pca.centroid().data(), m_pBeta->cols());
	m_pBeta->multiply(vw.vec(), m_epsilon, false);
	m_epsilon *= -1.0;
	GVec::add(m_epsilon.data(), pca.centroid().data() + inputs, outputs);

	// Refine the results using gradient descent
	refine(features, labels, 0.06, 20, 0.75);
}

void Loader::loadTest04(GMatrix &trainFeat, GMatrix &trainLab, GMatrix &testFeat, GMatrix &testLab)
{
	size_t train_size = 128;
	size_t test_size = 256;
	
	trainFeat.resize(train_size, 1);
	trainLab.resize(train_size, 1);
	testFeat.resize(test_size, 1);
	testLab.resize(test_size, 1);
	
	size_t sines = 2;
	
	double *x, *y;
	for(size_t i = 0; i < train_size + test_size; i++)
	{
		if(i < train_size)
		{
			x = trainFeat[i];
			y = trainLab[i];
		}
		else
		{
			x = testFeat[i - train_size];
			y = testLab[i - train_size];
		}
		
		*x = i / double(train_size);
		
		*y = 16.0 - 5 * *x;
		for(size_t k = 0; k < sines; k++)
		{
			*y -= sin((4.25 * M_PI * (k + 1)) * *x);
		}
	}
}

static GMatrix matrix_rotate(int w, int h) {
    GMatrix matrix;
    matrix.setTranslate(w/2, h/2);
    matrix.preRotate(M_PI / 6);
    matrix.preTranslate(-w/2, -h/2);
    return matrix;
}

void nominalToCat(GArgReader& args)
{
	// Load the file
	GMatrix* pData = loadData(args.pop_string());
	Holder<GMatrix> hData(pData);

	// Parse Options
	int maxValues = 12;
	while(args.size() > 0)
	{
		if(args.if_pop("-maxvalues"))
			maxValues = args.pop_uint();
		else
			ThrowError("Invalid option: ", args.peek());
	}

	// Transform the data
	GNominalToCat transform(maxValues);
	transform.train(*pData);
	GMatrix* pDataNew = transform.transformBatch(*pData);
	Holder<GMatrix> hDataNew(pDataNew);

	// Print results
	pDataNew->print(cout);
}

void reducedRowEchelonForm(GArgReader& args)
{
	GMatrix* pA = loadData(args.pop_string());
	Holder<GMatrix> hA(pA);
	pA->toReducedRowEchelonForm();
	pA->print(cout);
}

/**
*  @brief
*    Returns position, rotation and scale from a matrix
*/
void PLTools::GetPosRotScale(const GMatrix &mTransform, Point3 &vPos, Quat &qRot, Point3 &vScale, bool bFlip)
{
	// Get the position
	vPos = mTransform.Translation();

	// Get the rotation of the node as quaternion
	qRot = mTransform.Rotation();

	// Flip 180 degree around the y-axis? (true for camera and spot lights)
	if (bFlip) {
		Quat qRotationOffset;

		// We have to add a rotation about the x-axis of -90 degree... is this a IGame transform bug or something other odd??
		float fAngles[] = {static_cast<float>(HALFPI), 0.0f, 0.0f};
		EulerToQuat(fAngles, qRotationOffset);
		qRot = qRotationOffset*qRot;

		// We have to 'invert the z-axis', this is no PixelLight bug or so, we decided to do so to make things more universal
		fAngles[0] = 0.0f;
		fAngles[2] = static_cast<float>(PI);
		EulerToQuat(fAngles, qRotationOffset);
		qRot = qRotationOffset*qRot;
	}

	// Look out! We REALLY need to take the parity of the transform matrix into account!
	vScale = mTransform.Scaling()*static_cast<float>(mTransform.Parity());
}

void neighbors(GArgReader& args)
{
	// Load the data
	GMatrix* pData = loadData(args.pop_string());
	Holder<GMatrix> hData(pData);
	int neighborCount = args.pop_uint();

	// Find the neighbors
	GKdTree neighborFinder(pData, neighborCount, NULL, true);
	GTEMPBUF(size_t, neighbors, neighborCount);
	GTEMPBUF(double, distances, neighborCount);
	double sumClosest = 0;
	double sumAll = 0;
	for(size_t i = 0; i < pData->rows(); i++)
	{
		neighborFinder.neighbors(neighbors, distances, i);
		neighborFinder.sortNeighbors(neighbors, distances);
		sumClosest += sqrt(distances[0]);
		for(int j = 0; j < neighborCount; j++)
			sumAll += sqrt(distances[j]);
	}
	cout.precision(14);
	cout << "average closest neighbor distance = " << (sumClosest / pData->rows()) << "\n";
	cout << "average neighbor distance = " << (sumAll / (pData->rows() * neighborCount)) << "\n";
}