C++ SGMatrix类代码示例

void CECOCIHDDecoder::update_delta_cache(const SGMatrix<int32_t> codebook)
{
    if (codebook.matrix == m_codebook.matrix)
        return; // memory address the same

    if (codebook.num_cols == m_codebook.num_cols && codebook.num_rows == m_codebook.num_rows)
    {
        bool the_same = true;
        for (int32_t i=0; i < codebook.num_rows && the_same; ++i)
            for (int32_t j=0; j < codebook.num_cols && the_same; ++j)
                if (codebook(i,j) != m_codebook(i,j))
                    the_same = false;
        if (the_same)
            return; // no need to update delta
    }

    m_codebook = codebook; // operator=
    m_delta = SGMatrix<float64_t>(codebook.num_cols, codebook.num_cols);
    m_delta.zero();
    for (int32_t i=0; i < codebook.num_cols; ++i)
    {
        for (int32_t j=i+1; j < codebook.num_cols; ++j)
        {
            m_delta(i, j) = m_delta(j, i) = 
                CECOCUtil::hamming_distance(codebook.get_column_vector(i), codebook.get_column_vector(j), codebook.num_rows);
        }
    }

    // compute inverse of delta
    SGVector<int32_t> IPIV(m_delta.num_cols);
    clapack_dgetrf(CblasColMajor, m_delta.num_cols, m_delta.num_cols, m_delta.matrix, m_delta.num_cols, IPIV.vector);
    clapack_dgetri(CblasColMajor, m_delta.num_cols, m_delta.matrix, m_delta.num_cols, IPIV.vector);
}

SGVector<T>::SGVector(SGMatrix<T> matrix)
	: SGReferencedData(matrix), vlen(matrix.num_cols * matrix.num_rows),
	  gpu_ptr(NULL)
{
	ASSERT(!matrix.on_gpu())
	vector = matrix.data();
	m_on_gpu.store(false, std::memory_order_release);
}

void build_factor_graph(MultilabelParameter param, SGMatrix<float64_t> feats, SGMatrix<int32_t> labels,
		CFactorGraphFeatures * fg_feats, CFactorGraphLabels * fg_labels,
		const DynArray<CTableFactorType *>& v_ftp_u,
		const DynArray<CTableFactorType *>& v_ftp_t)
{
	int32_t num_sample        = labels.num_cols;
	int32_t num_classes       = labels.num_rows;
	int32_t dim               = feats.num_rows;

	SGMatrix< int32_t > mat_edges = get_edge_list(param.graph_type, num_classes);
	int32_t num_edges = mat_edges.num_rows;

	// prepare features and labels in factor graph
	for (int32_t n = 0; n < num_sample; n++)
	{
		SGVector<int32_t> vc(num_classes);
		SGVector<int32_t>::fill_vector(vc.vector, vc.vlen, NUM_STATUS);

		CFactorGraph * fg = new CFactorGraph(vc);

		float64_t * pfeat = feats.get_column_vector(n);
		SGVector<float64_t> feat_i(dim);
		memcpy(feat_i.vector, pfeat, dim * sizeof(float64_t));

		// add unary factors
		for (int32_t u = 0; u < num_classes; u++)
		{
			SGVector<int32_t> var_index_u(1);
			var_index_u[0] = u;
			CFactor * fac_u = new CFactor(v_ftp_u[u], var_index_u, feat_i);
			fg->add_factor(fac_u);
		}

		// add pairwise factors
		for (int32_t t = 0; t < num_edges; t++)
		{
			SGVector<float64_t> data_t(1);
			data_t[0] = 1.0;
			SGVector<int32_t> var_index_t = mat_edges.get_row_vector(t);
			CFactor * fac_t = new CFactor(v_ftp_t[t], var_index_t, data_t);
			fg->add_factor(fac_t);
		}
		
		// add factor graph instance
		fg_feats->add_sample(fg);

		// add label
		int32_t * plabs = labels.get_column_vector(n);
		SGVector<int32_t> states_gt(num_classes);
		memcpy(states_gt.vector, plabs, num_classes * sizeof(int32_t));
		SGVector<float64_t> loss_weights(num_classes);
		SGVector<float64_t>::fill_vector(loss_weights.vector, loss_weights.vlen, 1.0/num_classes);
		CFactorGraphObservation * fg_obs = new CFactorGraphObservation(states_gt, loss_weights);
		fg_labels->add_label(fg_obs);
	}
}

SGMatrix<float64_t> CLeastAngleRegression::cholesky_insert(
		SGMatrix<float64_t> X, SGMatrix<float64_t> R, int32_t i_max_corr)
{
	// diag_k = X[:,i_max_corr]' * X[:,i_max_corr]
	float64_t diag_k = cblas_ddot(X.num_rows, X.get_column_vector(i_max_corr), 1,
		X.get_column_vector(i_max_corr), 1);

	if (m_num_active == 0)
	{ // R isn't allocated yet
		SGMatrix<float64_t> nR(1,1);
		nR(0,0) = CMath::sqrt(diag_k);
		return nR;
	}
	else
	{

		// col_k is the k-th column of (X'X)
		vector<float64_t> col_k(m_num_active);
		for (int32_t i=0; i < m_num_active; ++i)
		{
			// col_k[i] = X[:,i_max_corr]' * X[:,m_active_set[i]]
			col_k[i] = cblas_ddot(X.num_rows, X.get_column_vector(i_max_corr), 1,
				X.get_column_vector(m_active_set[i]), 1);
		}

		// R' * R_k = (X' * X)_k = col_k, solving to get R_k
		vector<float64_t> R_k(col_k);
		cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, m_num_active, 1,
			1, R.matrix, m_num_active, &R_k[0], m_num_active);

		float64_t R_kk = CMath::sqrt(diag_k -
			cblas_ddot(m_num_active, &R_k[0], 1, &R_k[0], 1));

		// new_R = [R R_k; zeros(...) R_kk]
		SGMatrix<float64_t> nR(m_num_active+1, m_num_active+1);
		for (int32_t i=0; i < m_num_active; ++i)
			for (int32_t j=0; j < m_num_active; ++j)
				nR(i,j) = R(i,j);
		for (int32_t i=0; i < m_num_active; ++i)
			nR(i, m_num_active) = R_k[i];
		for (int32_t i=0; i < m_num_active; ++i)
			nR(m_num_active, i) = 0;
		nR(m_num_active, m_num_active) = R_kk;

		return nR;
	}

}

static void plane_rot(float64_t x0, float64_t x1,
	float64_t &y0, float64_t &y1, SGMatrix<float64_t> &G)
{
	G.zero();
	if (x1 == 0)
	{
		G(0, 0) = G(1, 1) = 1;
		y0 = x0;
		y1 = x1;
	}
	else
	{
		float64_t r = CMath::sqrt(x0*x0+x1*x1);
		float64_t sx0 = x0 / r;
		float64_t sx1 = x1 / r;

		G(0,0) = sx0;
		G(1,0) = -sx1;
		G(0,1) = sx1;
		G(1,1) = sx0;

		y0 = r;
		y1 = 0;
	}
}

void gen_rand_data(SGVector<float64_t> lab, SGMatrix<float64_t> feat,
		float64_t dist)
{
	index_t dims=feat.num_rows;
	index_t num=lab.vlen;

	for (int32_t i=0; i<num; i++)
	{
		if (i<num/2)
		{
			lab[i]=-1.0;

			for (int32_t j=0; j<dims; j++)
				feat(j, i)=CMath::random(0.0, 1.0)+dist;
		}
		else
		{
			lab[i]=1.0;

			for (int32_t j=0; j<dims; j++)
				feat(j, i)=CMath::random(0.0, 1.0)-dist;
		}
	}
	lab.display_vector("lab");
	feat.display_matrix("feat");
}

template<class ST> void CDenseFeatures<ST>::copy_feature_matrix(SGMatrix<ST> src)
{
	if (m_subset_stack->has_subsets())
		SG_ERROR("A subset is set, cannot call copy_feature_matrix\n")

	free_feature_matrix();
	feature_matrix = src.clone();
	num_features = src.num_rows;
	num_vectors = src.num_cols;
	initialize_cache();
}

SGVector<float64_t> CWeightedMajorityVote::combine(const SGMatrix<float64_t>& ensemble_result) const
{
	REQUIRE(m_weights.vlen == ensemble_result.num_cols, "The number of results and weights does not match!");
	SGVector<float64_t> mv(ensemble_result.num_rows);
	for (index_t i = 0; i < ensemble_result.num_rows; ++i)
	{
		SGVector<float64_t> rv = ensemble_result.get_row_vector(i);
		mv[i] = combine(rv);
	}

	return mv;
}

void read_data(const char * fname, SGMatrix<int32_t>& labels, SGMatrix<float64_t>& feats)
{
	// sparse data from matrix
	CLibSVMFile * svmfile = new CLibSVMFile(fname);

	SGSparseVector<float64_t>* spv_feats;
	SGVector<float64_t>* pv_labels;
	int32_t dim_feat;
	int32_t num_samples;
	int32_t num_classes;

	svmfile->get_sparse_matrix(spv_feats, dim_feat, num_samples, pv_labels, num_classes);

	SG_SPRINT("Number of the samples: %d\n", num_samples);
	SG_SPRINT("Dimention of the feature: %d\n", dim_feat+1);
	SG_SPRINT("Number of classes: %d\n", num_classes);

	feats  = SGMatrix<float64_t>(dim_feat+1, num_samples);
	labels = SGMatrix<int32_t>(num_classes, num_samples);
	feats.zero();
	labels.zero();

	for (int32_t i = 0; i < num_samples; i++)
	{
		SGVector<float64_t> v_feat = spv_feats[i].get_dense();
		SGVector<float64_t> v_labels = pv_labels[i];

		for (int32_t f = 0; f < v_feat.size(); f++)
			feats(f, i) = v_feat[f];
		
		feats(dim_feat, i) = 1.0; // bias

		for (int32_t l = 0; l < v_labels.size(); l++)
			labels((int32_t)v_labels[l], i) = 1;
	}

	SG_UNREF(svmfile);
	SG_FREE(spv_feats);
	SG_FREE(pv_labels);
}

int32_t CECOCIHDDecoder::decide_label(const SGVector<float64_t> outputs, const SGMatrix<int32_t> codebook)
{
    update_delta_cache(codebook);

    SGVector<float64_t> query = binarize(outputs);
    SGVector<float64_t> L(codebook.num_cols);
    for (int32_t i=0; i < codebook.num_cols; ++i)
        L[i] = CECOCUtil::hamming_distance(query.vector, codebook.get_column_vector(i), query.vlen);

    SGVector<float64_t> res(codebook.num_cols);
    res.zero();
    // res = m_delta * L
    cblas_dgemv(CblasColMajor, CblasNoTrans, m_delta.num_cols, m_delta.num_cols,
            1, m_delta.matrix, m_delta.num_cols, L.vector, 1, 1, res.vector, 1);
    return SGVector<float64_t>::arg_max(res.vector, 1, res.vlen);
}

SGMatrix<float64_t> CJADiagOrth::diagonalize(SGNDArray<float64_t> C, SGMatrix<float64_t> V0,
						double eps, int itermax)
{
	int m = C.dims[0];
	int L = C.dims[2];

	SGMatrix<float64_t> V;
	if (V0.num_rows == m && V0.num_cols == m)
		V = V0.clone();
	else
		V = SGMatrix<float64_t>::create_identity_matrix(m,1);

	bool more = true;
	int rots = 0;

	while (more)
	{
		more = false;

		for (int p = 0; p < m; p++)
		{
			for (int q = p+1; q < m; q++)
			{
				// computation of Givens angle
				float64_t theta = givens_stack(C.array, m, L, p, q);

				// Givens update
				if (fabs(theta) > eps)
				{
					float64_t c = cos(theta);
					float64_t s = sin(theta);
					left_rot_stack (C.array, m, m, L, p, q, c, s);
					right_rot_stack(C.array, m, m, L, p, q, c, s);
					left_rot_simple(V.matrix, m, m, p, q, c, s);
					rots++;
					more = true;
				}
			}
		}
	}

    return V;
}

void CConvolutionalFeatureMap::compute_weight_gradients(
	SGMatrix< float64_t > inputs, 
	SGMatrix< float64_t > local_gradients, 
	SGMatrix< float64_t > weight_gradients,
	int32_t inputs_row_offset,
 	int32_t local_gradients_row_offset)
{
	weight_gradients.zero();
	for (int32_t i=0; i<local_gradients.num_cols; i++)
	{
		SGMatrix<float64_t> image(
			inputs.matrix+i*inputs.num_rows + inputs_row_offset, 
			m_input_height, m_input_width, false);
		
		SGMatrix<float64_t> LG_image(
			local_gradients.matrix+i*local_gradients.num_rows 
			+ local_gradients_row_offset, m_output_height, m_output_width, false);
		
		for (int32_t x=0; x<m_input_width; x+=m_stride_x)
		{
			for (int32_t y=0; y<m_input_height; y+=m_stride_y)
			{
				for (int32_t x1=x-m_radius_x; x1<=x+m_radius_x; x1++)
				{
					for (int32_t y1=y-m_radius_y; y1<=y+m_radius_y; y1++)
					{
						if (x1>=0 && y1>=0 && x1<image.num_cols && y1<image.num_rows)
						{
							if (m_autoencoder_position == NLAP_NONE)
								weight_gradients(m_radius_y-y1+y,m_radius_x-x1+x) +=
									LG_image(y/m_stride_y,x/m_stride_x)*image(y1,x1);
							else
								weight_gradients(m_radius_y-y1+y,m_radius_x-x1+x) +=
									LG_image(y,x)*image(y1,x1);
						}
					}
				}
			}
		}
	}
}

void CTwoStateModel::reshape_transmission_params(
		SGMatrix< float64_t >& transmission_weights, SGVector< float64_t > w)
{
	transmission_weights.set_const(-CMath::INFTY);

	// Legend for state indices:
	// 0 -> start state
	// 1 -> stop state
	// 2 -> negative state (label == 0)
	// 3 -> positive state (label == 1)

	// From start
	transmission_weights(0,2) = 0;    // to negative
	transmission_weights(0,3) = 0;    // to positive
	// From negative
	transmission_weights(2,1) = 0;    // to stop
	transmission_weights(2,2) = w[0]; // to negative
	transmission_weights(2,3) = w[1]; // to positive
	// From positive
	transmission_weights(3,1) = 0;    // to stop
	transmission_weights(3,2) = w[3]; // to positive
	transmission_weights(3,3) = w[2]; // to negative
}

int main(int argc, char **argv)
{
	init_shogun_with_defaults();


	/* create some data and labels */
	SGMatrix<float64_t> matrix =
			SGMatrix<float64_t>(dim_vectors, num_vectors);

	SGMatrix<float64_t> matrix2 =
			SGMatrix<float64_t>(dim_vectors, num_vectors);
			
	CRegressionLabels* labels=new CRegressionLabels(num_vectors);

	build_matrices(matrix2, matrix, labels);
	
	/* create training features */
	CDenseFeatures<float64_t>* features=new CDenseFeatures<float64_t> ();
	features->set_feature_matrix(matrix);

	/* create testing features */
	CDenseFeatures<float64_t>* features2=new CDenseFeatures<float64_t> ();
	features2->set_feature_matrix(matrix2);

	SG_REF(features);
	SG_REF(features2);

	SG_REF(labels);
	
	/*Allocate our Kernel*/
	CGaussianKernel* test_kernel = new CGaussianKernel(10, 2);

	test_kernel->init(features, features);

	/*Allocate our mean function*/
	CZeroMean* mean = new CZeroMean();
	
	/*Allocate our likelihood function*/
	CGaussianLikelihood* lik = new CGaussianLikelihood();

	/*Allocate our inference method*/
	CExactInferenceMethod* inf =
			new CExactInferenceMethod(test_kernel, 
						  features, mean, labels, lik);

	SG_REF(inf);

	/*Finally use these to allocate the Gaussian Process Object*/
	CGaussianProcessRegression* gp =
			new CGaussianProcessRegression(inf, features, labels);

	SG_REF(gp);
	
	/*Build the parameter tree for model selection*/
	CModelSelectionParameters* root = build_tree(inf, lik, test_kernel);

	/*Criterion for gradient search*/
	CGradientCriterion* crit = new CGradientCriterion();

	/*This will evaluate our inference method for its derivatives*/
	CGradientEvaluation* grad=new CGradientEvaluation(gp, features, labels,
			crit);

	grad->set_function(inf);

	gp->print_modsel_params();

	root->print_tree();

	/* handles all of the above structures in memory */
	CGradientModelSelection* grad_search=new CGradientModelSelection(
			root, grad);

	/* set autolocking to false to get rid of warnings */
	grad->set_autolock(false);

	/*Search for best parameters*/
	CParameterCombination* best_combination=grad_search->select_model(true);

	/*Output all the results and information*/
	if (best_combination)
	{
		SG_SPRINT("best parameter(s):\n");
		best_combination->print_tree();

		best_combination->apply_to_machine(gp);
	}

	CGradientResult* result=(CGradientResult*)grad->evaluate();

	if(result->get_result_type() != GRADIENTEVALUATION_RESULT)
		SG_SERROR("Evaluation result not a GradientEvaluationResult!");

	result->print_result();

	SGVector<float64_t> alpha = inf->get_alpha();
	SGVector<float64_t> labe = labels->get_labels();
	SGVector<float64_t> diagonal = inf->get_diagonal_vector();
	SGMatrix<float64_t> cholesky = inf->get_cholesky();
	gp->set_return_type(CGaussianProcessRegression::GP_RETURN_COV);
//.........这里部分代码省略.........

int main(int argc, char** argv)
{
	int32_t num_vectors = 0;
	int32_t num_feats   = 0;

	init_shogun_with_defaults();

	const char*fname_train = "../data/7class_example4_train.dense";
	CStreamingAsciiFile *train_file = new CStreamingAsciiFile(fname_train);
	SG_REF(train_file);

	CStreamingDenseFeatures<float64_t> *stream_features = new CStreamingDenseFeatures<float64_t>(train_file, true, 1024);
	SG_REF(stream_features);

	SGMatrix<float64_t> mat;
	SGVector<float64_t> labvec(1000);

	stream_features->start_parser();
	SGVector< float64_t > vec;
	while (stream_features->get_next_example())
	{
		vec = stream_features->get_vector();
		if (num_feats == 0)
		{
			num_feats = vec.vlen;
			mat = SGMatrix<float64_t>(num_feats, 1000);
		}
		std::copy(vec.vector, vec.vector+vec.vlen, mat.get_column_vector(num_vectors));
		labvec[num_vectors] = stream_features->get_label();
		num_vectors++;
		stream_features->release_example();
	}
	stream_features->end_parser();
	mat.num_cols = num_vectors;
	labvec.vlen = num_vectors;
	
	CMulticlassLabels* labels = new CMulticlassLabels(labvec);
	SG_REF(labels);

	// Create features with the useful values from mat
	CDenseFeatures< float64_t >* features = new CDenseFeatures<float64_t>(mat);
	SG_REF(features);

	SG_SPRINT("Performing ShareBoost on a %d-class problem\n", labels->get_num_classes());

	// Create ShareBoost Machine
	CShareBoost *machine = new CShareBoost(features, labels, 10);
	SG_REF(machine);

	machine->train();

	SGVector<int32_t> activeset = machine->get_activeset();
	SG_SPRINT("%d out of %d features are selected:\n", activeset.vlen, mat.num_rows);
	for (int32_t i=0; i < activeset.vlen; ++i)
		SG_SPRINT("activeset[%02d] = %d\n", i, activeset[i]);

	CDenseSubsetFeatures<float64_t> *subset_fea = new CDenseSubsetFeatures<float64_t>(features, machine->get_activeset());
	SG_REF(subset_fea);
	CMulticlassLabels* output = CMulticlassLabels::obtain_from_generic(machine->apply(subset_fea));

	int32_t correct = 0;
	for (int32_t i=0; i < output->get_num_labels(); ++i)
		if (output->get_int_label(i) == labels->get_int_label(i))
			correct++;
	SG_SPRINT("Accuracy = %.4f\n", float64_t(correct)/labels->get_num_labels());

	// Free resources
	SG_UNREF(machine);
	SG_UNREF(output);
	SG_UNREF(subset_fea);
	SG_UNREF(features);
	SG_UNREF(labels);
	SG_UNREF(train_file);
	SG_UNREF(stream_features);
	exit_shogun();

	return 0;
}