C++ RowVectorXd类代码示例

/*******M step, maximazize log-likelihood*/
void Plsa::M_step(MatrixXd &data)
{
    MatrixXd X;
    for (int i=0;i<K;i++)
    {
      Pw_z.col(i)=(Pz_wd[i].cwiseProduct(data)).rowwise().sum();//suma de filas[vector(1XN)]
      Pd_z.col(i)=(Pz_wd[i].cwiseProduct(data)).colwise().sum();//suma de columnas[vector(1XN)]
    }

    //normalize
    RowVectorXd Temp;
    RowVectorXd C;
    VectorXd E;
    VectorXd T;
    Temp=RowVectorXd::Ones(K);
    T=VectorXd::Ones(K); // vector of K with ones
    P_z=Pd_z.colwise().sum();
    cout<<P_z;

    C=Pd_z.colwise().sum(); //suma de columnas[vector(1XN)]
    Temp=Temp.cwiseQuotient(C);
    Pd_z=Pd_z*(Temp.asDiagonal());

    C=Pw_z.colwise().sum(); //
    Temp=Temp.cwiseQuotient(C);
    Pw_z=Pw_z*(Temp.asDiagonal());

    E=P_z.rowwise().sum();
    P_z=P_z.cwiseProduct(T.cwiseQuotient(E));

}

/* compares evaluations on corner cases */
void eval_spline3d_onbrks()
{
  Spline3d spline = spline3d();

  RowVectorXd u = spline.knots();

  MatrixXd pts(11,3);
  pts <<    0.959743958516081,   0.340385726666133,   0.585267750979777,
    0.959743958516081,   0.340385726666133,   0.585267750979777,
    0.959743958516081,   0.340385726666133,   0.585267750979777,
    0.430282980289940,   0.713074680056118,   0.720373307943349,
    0.558074875553060,   0.681617921034459,   0.804417124839942,
    0.407076008291750,   0.349707710518163,   0.617275937419545,
    0.240037008286602,   0.738739390398014,   0.324554153129411,
    0.302434111480572,   0.781162443963899,   0.240177089094644,
    0.251083857976031,   0.616044676146639,   0.473288848902729,
    0.251083857976031,   0.616044676146639,   0.473288848902729,
    0.251083857976031,   0.616044676146639,   0.473288848902729;
  pts.transposeInPlace();

  for (int i=0; i<u.size(); ++i)
  {
    Vector3d pt = spline(u(i));
    VERIFY( (pt - pts.col(i)).norm() < 1e-14 );
  }
}

void Plsa::Normalize(MatrixXd &Mat)
{
  RowVectorXd Temp;
  Temp=RowVectorXd::Ones(K);
  Temp=Temp.cwiseQuotient(Mat.colwise().sum());
  Mat=Mat*Temp.asDiagonal();
}

void RelativeEulerAnglePcaDecoder::Encode(array_view<const DirectX::Quaternion> rots, VectorType & x)
{
	RelativeEulerAngleDecoder::Encode(rots, x);
	RowVectorXd dx = x.transpose().cast<double>();
	dx -= meanY;
	dx *= pcaY;
	x = dx.transpose().cast<float>();
}

QMap<int,QList<QPair<int,double> > > DetectTrigger::detectTriggerFlanksGrad(const MatrixXd& data, const QList<int>& lTriggerChannels, int iOffsetIndex, double dThreshold, bool bRemoveOffset, const QString& type, int iBurstLengthSamp)
{
    QMap<int,QList<QPair<int,double> > > qMapDetectedTrigger;
    RowVectorXd tGradient = RowVectorXd::Zero(data.cols());

    //Find all triggers above threshold in the data block
    for(int i = 0; i < lTriggerChannels.size(); ++i)
    {
//        QTime time;
//        time.start();

        int iChIdx = lTriggerChannels.at(i);

        //Add empty list to map
        QList<QPair<int,double> > temp;
        qMapDetectedTrigger.insert(iChIdx, temp);

        //detect the actual triggers in the current data matrix
        if(iChIdx > data.rows() || iChIdx < 0)
        {
            return qMapDetectedTrigger;
        }

        //Compute gradient
        for(int t = 1; t<tGradient.cols(); t++)
        {
            tGradient(t) = data(iChIdx,t)-data(iChIdx,t-1);
        }

        // If falling flanks are to be detected flip the gradient's sign
        if(type == "Falling")
        {
            tGradient = tGradient * -1;
        }

        //Find positive maximum in gradient vector. This position is equal to the rising trigger flank.
        for(int j = 0; j < tGradient.cols(); ++j)
        {
            double dMatVal = bRemoveOffset ? tGradient(j) - data(iChIdx,0) : tGradient(j);

            if(dMatVal >= dThreshold)
            {
                QPair<int,double> pair;
                pair.first = iOffsetIndex+j;
                pair.second = tGradient(j);

                qMapDetectedTrigger[iChIdx].append(pair);

                j += iBurstLengthSamp;
            }
        }

//        int timeElapsed = time.elapsed();
//        std::cout<<"timeElapsed: "<<timeElapsed<<std::endl;
    }

    return qMapDetectedTrigger;
}

RowVectorXd FilterData::applyFFTFilter(const RowVectorXd& data, bool keepOverhead, CompensateEdgeEffects compensateEdgeEffects) const
{
    #ifdef EIGEN_FFTW_DEFAULT
        fftw_make_planner_thread_safe();
    #endif

    if(data.cols()<m_dCoeffA.cols() && compensateEdgeEffects==MirrorData) {
        qDebug()<<QString("Error in FilterData: Number of filter taps(%1) bigger then data size(%2). Not enough data to perform mirroring!").arg(m_dCoeffA.cols()).arg(data.cols());
        return data;
    }

//    std::cout<<"m_iFFTlength: "<<m_iFFTlength<<std::endl;
//    std::cout<<"2*m_dCoeffA.cols() + data.cols(): "<<2*m_dCoeffA.cols() + data.cols()<<std::endl;

    if(2*m_dCoeffA.cols() + data.cols()>m_iFFTlength) {
        qDebug()<<"Error in FilterData: Number of mirroring/zeropadding size plus data size is bigger then fft length!";
        return data;
    }

    //Do zero padding or mirroring depending on user input
    RowVectorXd t_dataZeroPad = RowVectorXd::Zero(m_iFFTlength);

    switch(compensateEdgeEffects) {
        case MirrorData:
            t_dataZeroPad.head(m_dCoeffA.cols()) = data.head(m_dCoeffA.cols()).reverse();   //front
            t_dataZeroPad.segment(m_dCoeffA.cols(), data.cols()) = data;                    //middle
            t_dataZeroPad.tail(m_dCoeffA.cols()) = data.tail(m_dCoeffA.cols()).reverse();   //back
            break;

        case ZeroPad:
            t_dataZeroPad.head(data.cols()) = data;
            break;

        default:
            t_dataZeroPad.head(data.cols()) = data;
            break;
    }

    //generate fft object
    Eigen::FFT<double> fft;
    fft.SetFlag(fft.HalfSpectrum);

    //fft-transform data sequence
    RowVectorXcd t_freqData;
    fft.fwd(t_freqData,t_dataZeroPad);

    //perform frequency-domain filtering
    RowVectorXcd t_filteredFreq = m_dFFTCoeffA.array()*t_freqData.array();

    //inverse-FFT
    RowVectorXd t_filteredTime;
    fft.inv(t_filteredTime,t_filteredFreq);

    //Return filtered data
    if(!keepOverhead)
        return t_filteredTime.segment(m_dCoeffA.cols()/2, data.cols());

    return t_filteredTime.head(data.cols()+m_dCoeffA.cols());
}

QPair< int,QList<double> > BCI::applyFeatureCalcConcurrentlyOnSensorLevel(const QPair<int, RowVectorXd> &chdata)
{
    RowVectorXd data = chdata.second;
    QList<double> features;

    // TODO: Divide into subsignals
    //features << data.squaredNorm(); // Compute variance
    features << abs(log10(data.squaredNorm())); // Compute log of variance

    return QPair< int,QList<double> >(chdata.first, features);
}

RowVectorXd DataPackage::cutData(const RowVectorXd &originalData, int cutFront, int cutBack)
{
    if(originalData.cols()-cutFront-cutBack < 0 || cutFront>originalData.cols()) {
        qDebug()<<"DataPackage::cutData - cutFront or cutBack do not fit. Aborting mapping and returning original data.";
        RowVectorXd returnVec = originalData;
        return returnVec;
    }

    //Cut original data using segment
    return (RowVectorXd)originalData.segment(cutFront, originalData.cols()-cutFront-cutBack);
}

CDataObject::CDataObject(int total_pts, enumStreetIndices br_ndx, int times_acted, int ndealt, eBetType action_type) : 
	m_npoints(total_pts), 
	m_br(br_ndx), 
	m_nacted(times_acted),
	m_ndealt(ndealt),
	m_action(action_type)
{
	assert(br_ndx >= ePreflopIndex && br_ndx < eRoundIndices);
	assert(times_acted >= 0);

	m_ndims = times_acted + (br_ndx == ePreflopIndex ? num_prior_dims_preflop : num_prior_dims_postflop) + 1;

	// this is only used for ann
	m_data = new double*[m_npoints];
	// allocate hand_ids for corresponding data points
	m_hand_ids = new long[m_npoints];
	for(long i = 0; i < m_npoints; i++)
		m_hand_ids[i] = -1L;

	m_profits	= VectorXd::Zero(m_npoints);
	m_points	= Matrix<double, Dynamic, Dynamic, RowMajor>::Zero(m_npoints, m_ndims);

	CDatabase p_db;
	GetData(&p_db);

	// Set the weights
	diag_factor	= Matrix<double, Dynamic, Dynamic, RowMajor>::Identity(m_ndims, m_ndims);

	#ifdef KASPER_WEIGHTS

		RowVectorXd tmp;
		if(ePreflopIndex == br_ndx)
		{
			tmp = RowVectorXd::Zero(9);
			tmp << 10, 20, 100, 10, 0, 0, 2, 2, 20;
		}
		else
		{
			tmp = RowVectorXd::Zero(11);
			tmp << 10, 80, 10, 0, 0, 0, 2, 2, 50, 1, 1;
		}
		diag_factor.bottomRightCorner(tmp.cols(), tmp.cols()) = tmp.asDiagonal();

	#else

		FeatureNormalize();
		CRegressionObject regress(m_points, m_profits);
		diag_factor = regress.get_theta().asDiagonal();

	#endif

	gLog.WriteLog(eSeverityInfo, eCatPerformance, "br%d_%d: %8s - %8d points\n", br_ndx+1, times_acted, bets_str[action_type], m_npoints);
}

void IKoptions::setAdditionaltSamples(const RowVectorXd &t_samples) {
  if (t_samples.size() > 0) {
    set<double> unique_sort_t(t_samples.data(),
                              t_samples.data() + t_samples.size());
    this->additional_tSamples.resize(unique_sort_t.size());
    int t_idx = 0;
    for (auto it = unique_sort_t.begin(); it != unique_sort_t.end(); it++) {
      this->additional_tSamples(t_idx) = *it;
      t_idx++;
    }
  } else {
    this->additional_tSamples.resize(0);
  }
}

QPair<int,double> ConnectivityMeasures::calcCrossCorrelation(const RowVectorXd& vecFirst, const RowVectorXd& vecSecond)
{
    Eigen::FFT<double> fft;

    int N = std::max(vecFirst.cols(), vecSecond.cols());

    //Compute the FFT size as the "next power of 2" of the input vector's length (max)
    int b = ceil(log2(2.0 * N - 1));
    int fftsize = pow(2,b);
//    int end = fftsize - 1;
//    int maxlag = N - 1;

    //Zero Padd
    RowVectorXd xCorrInputVecFirst = RowVectorXd::Zero(fftsize);
    xCorrInputVecFirst.head(vecFirst.cols()) = vecFirst;

    RowVectorXd xCorrInputVecSecond = RowVectorXd::Zero(fftsize);
    xCorrInputVecSecond.head(vecSecond.cols()) = vecSecond;

    //FFT for freq domain to both vectors
    RowVectorXcd freqvec;
    RowVectorXcd freqvec2;

    fft.fwd(freqvec, xCorrInputVecFirst);
    fft.fwd(freqvec2, xCorrInputVecSecond);

    //Create conjugate complex
    freqvec2.conjugate();

    //Main step of cross corr
    for (int i = 0; i < fftsize; i++) {
        freqvec[i] = freqvec[i] * freqvec2[i];
    }

    RowVectorXd result;
    fft.inv(result, freqvec);

    //Will get rid of extra zero padding
    RowVectorXd result2 = result;//.segment(maxlag, N);

    QPair<int,int> minMaxRange;
    int idx = 0;
    result2.minCoeff(&idx);
    minMaxRange.first = idx;
    result2.maxCoeff(&idx);
    minMaxRange.second = idx;

//    std::cout<<"result2(minMaxRange.first)"<<result2(minMaxRange.first)<<std::endl;
//    std::cout<<"result2(minMaxRange.second)"<<result2(minMaxRange.second)<<std::endl;
//    std::cout<<"b"<<b<<std::endl;
//    std::cout<<"fftsize"<<fftsize<<std::endl;
//    std::cout<<"end"<<end<<std::endl;
//    std::cout<<"maxlag"<<maxlag<<std::endl;

    //Return val
    int resultIndex = minMaxRange.second;
    double maxValue = result2(resultIndex);

    return QPair<int,double>(resultIndex, maxValue);
}

void
PatchSample::sample_frustum(double hfov, double vfov, int nh, int nv,
                            Vector3d p, Vector3d u,
                            MatrixXd &mx, MatrixXd &my, MatrixXd &mz)
{
  //TBD: check input args

  // pointing, up, left
  p = p/p.norm();
  u = u/u.norm();
  Vector3d l = u.cross(p);

  // sample like a camera would
  double ll = (tan(hfov/2.0)*((double)(nh-1)))/nh;
  double uu = (tan(vfov/2.0)*((double)(nv-1)))/nv;
  RowVectorXd y;
  y.setLinSpaced(nh,-ll,ll);
  MatrixXd yy;
  yy = y.replicate(nv,1);
  VectorXd z;
  z.setLinSpaced(nv,-uu,uu);
  MatrixXd zz;
  zz = z.replicate(1,nh);
  MatrixXd xx = MatrixXd::Ones(nv,nh);
  MatrixXd nn = (xx.array().square() + yy.array().square() + zz.array().square()).cwiseSqrt();
  xx = xx.array() / nn.array();
  yy = yy.array() / nn.array();
  zz = zz.array() / nn.array();

  // rotation matrix
  Matrix3d rr;
  rr << p, l, u;

  // rotate points
  MatrixXd xyz;
  MatrixXd cam (3,xx.rows()*xx.cols());
  cam.row(0) = vectorizeColWise(xx);
  cam.row(1) = vectorizeColWise(yy);
  cam.row(2) = vectorizeColWise(zz);
  xyz = rr*cam;

  // extract coordinates
  xx = xyz.row(0);
  yy = xyz.row(1);
  zz = xyz.row(2);
  mx = Map<MatrixXd>(xx.data(),nv,nh);
  my = Map<MatrixXd>(yy.data(),nv,nh);
  mz = Map<MatrixXd>(zz.data(),nv,nh);
}

VectorXd probutils::mahaldist (
    const MatrixXd& X,
    const RowVectorXd& mu,
    const MatrixXd& A
    )
{
  // Check for same number of dimensions, D
  if((X.cols() != mu.cols()) || (X.cols() != A.cols()))
    throw invalid_argument("Arguments do not have the same dimensionality");

  // Check if A is square
  if (A.rows() != A.cols())
    throw invalid_argument("Matrix A must be square!");

  // Decompose A
  LDLT<MatrixXd> Aldl(A);

  // Check if A is PD
  if ((Aldl.vectorD().array() <= 0).any() == true)
    throw invalid_argument("Matrix A is not positive definite");

  // Do the Mahalanobis distance for each sample (N times)
  MatrixXd X_mu = (X.rowwise() - mu).transpose();
  return ((X_mu.array() * (Aldl.solve(X_mu)).array())
          .colwise().sum()).transpose();
}

void FrequencySpectrumDelegate::createPlotPath(const QModelIndex &index, const QStyleOptionViewItem &option, QPainterPath& path, RowVectorXd& data) const
{
    const FrequencySpectrumModel* t_pModel = static_cast<const FrequencySpectrumModel*>(index.model());

    float fMaxValue = data.maxCoeff();

    float fValue;
    float fScaleY = option.rect.height()/(fMaxValue*0.5);

    float y_base = path.currentPosition().y();
    QPointF qSamplePosition;

    qint32 lowerIdx = t_pModel->getLowerFrqBound();
    qint32 upperIdx = t_pModel->getUpperFrqBound();

    //Move to initial starting point
    if(data.size() > 0)
    {
        float val = 0;
        fValue = val*fScaleY;

        float newY = y_base+fValue;

        qSamplePosition.setY(newY);
        qSamplePosition.setX((double)option.rect.width()*t_pModel->getFreqScaleBound()[lowerIdx]);

        path.moveTo(qSamplePosition);
    }


    //create lines from one to the next sample
    qint32 i;
    for(i = lowerIdx+1; i <= upperIdx; ++i) {
        float val = data[i]-data[0]; //remove first sample data[0] as offset
        fValue = val*fScaleY;

        float newY = y_base+fValue;

        qSamplePosition.setY(newY);
        qSamplePosition.setX((double)option.rect.width()*t_pModel->getFreqScaleBound()[i]);

        path.lineTo(qSamplePosition);
    }
}

//target function przyjmuje macierz, wiec mozemy oba porownywane punkty polaczyc i przeslac.
bool NelderMead::Compare(RowVectorXd first, RowVectorXd second)
{
    RowVector2d val;
    MatrixXd temp(2, first.cols());
    temp.row(0) = first;
    temp.row(1) = second;
    val = this->TargetFunction(temp);

    return val[0] < val[1];
}