C++ Matrix::GetRows方法代码示例

/**
 * returns a diagonal matrix
 * @param  v a vector
 * @return a diagonal matrix with the given vector v on the diagonal
 */
Matrix Diag(const Matrix& v)
{
  Matrix res;
  if (v.GetCols() == 1)
  {
    // the given matrix is a vector n x 1
    int rows = v.GetRows();
    res = Matrix(rows, rows);

    // copy the values of the vector to the matrix
    for (int r=1; r <= rows; r++)
    {
      res(r, r) = v.get(r, 1);
    }
  }
  else if (v.GetRows() == 1)
  {
    // the given matrix is a vector 1 x n
    int cols = v.GetCols();
    res = Matrix(cols, cols);

    // copy the values of the vector to the matrix
    for (int c=1; c <= cols; c++)
    {
      res(c, c) = v.get(1, c);
    }
  }
  else
  {
    throw Exception("Parameter for diag must be a vector");
  }
  return res;
}

norm_tools::dist_stat calc_dist_matrix_impl(const Matrix& data, Matrix& dist, norm_tools::vv_norm_fn_t norm2_fn)
{
	//some constant values used
	const ulong points_num = data.row_num();
	const double mult = 2.0/(points_num * (points_num - 1));
	const double mult1 = 1.0/(points_num - 1);

	//localy used matrices
	Matrix dv, norm; //, dist_row(1, points_num);
	//statistics
	norm_tools::dist_stat stat;
	//meand distance^2 & meand distance^2 between nearest neighbours
	double meand2 = 0, meand2_nn = 0;
	//current distance & minimum distances
	double cur_dist, cur_mind = 0, cur_mind2 = 0;

	//resize distance matrix
	dist.Resize(points_num, points_num);
	//zero distance matrix
	dist = 0;
	//start distances calculation
	for(ulong i = 0; i < points_num - 1; ++i) {
		dv <<= data.GetRows(i);
		//calc dist^2 to all other rows
		for(ulong j = i + 1; j < points_num; ++j) {
			//calc dist^2
			//norm <<= dv - data.GetRows(j);
			//transform(norm, norm, multiplies<double>());
			//cur_dist = distance^2
			//cur_dist = norm.Sum();
			cur_dist = norm2_fn(dv, data.GetRows(j));

			//update mean of distance^2
			meand2 += mult * cur_dist;
			//update current min distance^2
			if(j == i + 1 || cur_dist < cur_mind2) cur_mind2 = cur_dist;
			//cur_dist = pure l2 distance
			cur_dist = sqrt(cur_dist);
			//update current minimum distance (to nearest neighbour)
			if(j == i + 1 || cur_dist < cur_mind) cur_mind = cur_dist;
			//save it into matrix elements (i,j) and (j, i)
			dist(i, j) = cur_dist; dist(j, i) = cur_dist;
			//update global mean distance
			stat.mean_ += cur_dist * mult;
		}
		//update global min distance
		if(i == 0 || cur_mind < stat.min_) stat.min_ = cur_mind;
		//update mean nearest neighbour distance
		stat.mean_nn_ += cur_mind * mult1;
		//update mean nearest neighbour distance^2
		meand2_nn += cur_mind2 * mult1;
	}

	//calc mse
	stat.mse_ = sqrt(stat.mean_*(meand2/stat.mean_ - stat.mean_));
	stat.mse_nn_ = sqrt(stat.mean_nn_*(meand2_nn/stat.mean_nn_ - stat.mean_nn_));
	return stat;
}

// / Calculates the mean value of the elements in a matrix.
double mean(Matrix < double >&m)
{
    double sum = 0.;

    for (int r = 0; r < m.GetRows(); r++)
        for (int c = 0; c < m.GetCols(); c++)
            sum += m[r][c];
    return sum / double (m.GetRows() * m.GetCols());
}

Complex scalprod(Matrix < Complex > &m1, Matrix < double >&m2)
{
    assert((m1.GetRows() == m2.GetRows()) && (m1.GetCols() == m2.GetCols()));
    Complex sum(0, 0);

    for (int r = 0; r < m1.GetRows(); r++)
        for (int c = 0; c < m1.GetCols(); c++)
            sum += m1[r][c] * m2[r][c];
    return sum;
}

// / Calculates the varaince of the elements in a matrix.
double variance(Matrix < double >&m)
{
    double sum = 0.;

    double average = mean(m);

    for (int r = 0; r < m.GetRows(); r++)
        for (int c = 0; c < m.GetCols(); c++)
            sum += (m[r][c] - average) * (m[r][c] - average);
    return sum / double (m.GetRows() * m.GetCols() - 1);
}

Complex scalconj(Matrix < Complex > &m1, Matrix < Complex > &m2)        // total(A.B*)
{
    assert((m1.GetRows() == m2.GetRows()) && (m1.GetCols() == m2.GetCols()));
    Complex sum(0, 0);

    for (int r = 0; r < m1.GetRows(); r++)
        for (int c = 0; c < m1.GetCols(); c++)
            sum += m1[r][c] * conj(m2[r][c]);
    return sum;

}

void GA::nn_addon::_build_surf(ulong net_ind, const Matrix& input, const Matrix& targets)
{
	/*
	Matrix x(1, 100), y(100, 1);
	x = abs(_ga.opt_.initRange(1, 0) - _ga.opt_.initRange(0, 0))/(x.size() - 1);
	x[0] = 0; x <<= !(!x).CumSum(); x += _ga.opt_.initRange(0, 0);
	y = abs(_ga.opt_.initRange(1, 1) - _ga.opt_.initRange(0, 1))/(y.size() - 1);
	y[0] = 0; y <<= y.CumSum(); y += _ga.opt_.initRange(0, 1);

	Matrix yrow(1, 100);
	Matrix surf;
	for(ulong i = 0; i < y.row_num(); ++i) {
		yrow = y[i];
		if(opt_.netType == matrix_nn)
			surf &= _net.Sim(x & yrow);
		else
			surf &= _onet[net_ind]->sim(x & yrow);
	}l
	DumpMatrix(x, "x.txt");
	DumpMatrix(y, "y.txt");
	DumpMatrix(surf, "z.txt");
	*/

	Matrix points(ga_.opt_.initRange.col_num(), 10000);
	generate(points.begin(), points.end(), prg::rand01);
	const ulong pnum = points.col_num();
	//range matrix
	Matrix range = input.minmax(false);
	//lower bound
	Matrix a = range.GetRows(0);
	//difference between bounds
	Matrix scale = range.GetRows(1) - a;
	//start moving points
	Matrix::r_iterator pos = points.begin();
	Matrix::r_iterator p_a = a.begin();
	Matrix::r_iterator p_sc = scale.begin();
	for(ulong i = 0; i < points.row_num(); ++i, ++p_a, ++p_sc)
	{
		//x = x*(b - a)
		transform(pos, pos + pnum, pos, bind2nd(multiplies<double>(), *p_sc));
		pos = transform(pos, pos + pnum, pos, bind2nd(plus<double>(), *p_a));
		//pos += points.col_num();
	}
	DumpMatrix(points, "y.txt");

	Matrix surf;
	if(opt_.netType == matrix_nn)
		surf <<= net_.Sim(points);
	else
		surf <<= _onet[net_ind]->sim(points);
	DumpMatrix(surf, "z.txt");
}

Matrix nn_addon::_best_filter(const Matrix& p, const Matrix& f, Matrix& lp, Matrix& lf)
{
	lf = f;
	indMatrix mInd = lf.RawSort();
	lp.NewMatrix(min<ulong>(opt_.bestCount, p.row_num()), p.col_num());
	for(ulong i = 0; i < lp.row_num(); ++i)
		lp.SetRows(p.GetRows(mInd[i]), i);
	if(lp.row_num() < lf.row_num())
		lf.DelRows(lp.row_num(), lf.row_num() - lp.row_num());

	if(opt_.search_samples_num > 0)
		return lp.GetRows(0, opt_.search_samples_num).minmax();
	else
		return lp.minmax();
}

Matrix Matrix::operator*(Matrix oper)
{
	float *B = oper.GetMatrix();
	int RB = oper.GetRows();
	int CB = oper.GetColumns();
	Matrix C(oper.GetColumns(), GetColumns());
	int CSize = C.GetColumns()*C.GetRows();
	int Columns = C.GetColumns();
	float *Result = C.GetMatrix();

	Result[0] = A[0]*B[0] + A[1]*B[3] + A[2]*B[6]; 
	Result[1] = A[0]*B[1] + A[1]*B[4] + A[2]*B[7]; 
	Result[2] = A[0]*B[2] + A[1]*B[5] + A[2]*B[8]; 

	Result[3] = A[3]*B[0] + A[4]*B[3] + A[5]*B[6]; 
	Result[4] = A[3]*B[1] + A[4]*B[4] + A[5]*B[7]; 
	Result[5] = A[3]*B[2] + A[4]*B[5] + A[5]*B[8]; 

	Result[6] = A[6]*B[0] + A[7]*B[3] + A[8]*B[6]; 
	Result[7] = A[6]*B[1] + A[7]*B[4] + A[8]*B[7]; 
	Result[8] = A[6]*B[2] + A[7]*B[5] + A[8]*B[8]; 
	

	return C;
}

Matrix fft2(Matrix &m, int isign)
{
    assert(m.bComplex == TRUE);
    int ndim = 2;	
    long row = m.GetRows(); 
    long col = m.GetCols();
	
    Matrix ret(m);
    unsigned long nn[3];
		
    nn[1] = row;
    nn[2] = col / 2;
		
#if 1
    float *pdata = new float[row * col + 1];
    memcpy( pdata + 1, m.p, sizeof(float) * row * col);
    fourn( pdata, nn, ndim, isign );
    memcpy( ret.p, pdata + 1, sizeof(float) * row * col);
    delete[] pdata;
#else
	float* pdata = new float[row * col];
	memcpy(pdata, m.p, sizeof(float) * row * col);
	CUFFT::fft(pdata, col / 2, row, isign == 1);
	CUFFT::fft(pdata, col / 2, row, isign != 1);
	memcpy(ret.p, pdata, sizeof(float) * row * col);
	delete [] pdata;
#endif
    ret.bComplex = TRUE;
		
    return ret;	
}

Matrix WhitenFrame(Matrix &m)
{
    int nRow = m.GetRows();
    int nCol = m.GetCols();
    assert(nRow == nCol);
	
    Matrix fx(nRow), fy(nRow);
    Meshgrid(fx, -nRow / 2, nRow / 2 - 1, fy, -nRow / 2, nRow / 2 - 1);
    Matrix theta = fx;
    Matrix rho = fx;
    
    cart2pol( fx,fy,theta,rho );
	
#if 1
    Matrix complex = Complex(m);
    Matrix fftm = fft2(complex, 1);
    Matrix imF = fftshift(fftm);
    Matrix time = times(rho, imF);
    Matrix imW = fftshift(time);
#else
    Matrix complex = Complex(m);
    Matrix fftm = fft2(complex, 1);
    Matrix imF = fftshift(fftm);
    Matrix time = times(rho, imF);
    Matrix imW = fftshift(time);
#endif

    return imW;
}

Matrix nn_addon::GetChromMM(const Matrix& pop)
{
	Matrix res(pop.row_num(), 2);
	for(ulong i=0; i<pop.row_num(); ++i) {
		res(i, 0) = pop.GetRows(i).Min();
		res(i, 1) = pop.GetRows(i).Max();
	}
	return res;
}

Complex total(Matrix < Complex > &m)
{
    Complex sum(0, 0);

    for (int r = 0; r < m.GetRows(); r++)
        for (int c = 0; c < m.GetCols(); c++)
            sum += m[r][c];
    return sum;
}

double norm2(Matrix < Complex > &m)
{
    double sum = 0.;

    for (int r = 0; r < m.GetRows(); r++)
        for (int c = 0; c < m.GetCols(); c++)
            sum += norm(m[r][c]);
    return sum;

}

template< > _CLASS_DECLSPEC
Matrix norm_tools::vm_norm2< norm_tools::l2 >(const Matrix& v_from, const Matrix& m_to) {
	Matrix norm(1, m_to.row_num());
	Matrix diff;
	for(ulong i = 0; i < norm.size(); ++i) {
		diff <<= v_from - m_to.GetRows(i);
		norm[i] = diff.Mul(diff).Sum();
	}
	return norm;
}