C++ ArrayXd::resize方法代码示例

void CMATLAB::matGet(string filename,ArrayXd &X,string Xname)
{
    MATFile *pmat=matOpen(filename,string("r"));
    if (!pmat) throw "CMATLAB::importFrom error: matlab.matOpen failed";
    mxArray *pa=matGetVariable(pmat,Xname);
    if (!pa) throw "CMATLAB::importFrom error: matlab.matGetVariable failed";
    int N=mxGetNumberOfElements(pa);
    if (N<=0) throw "CMATLAB::importFrom error: matlab.mxGetNumberOfElements failed";
    X.resize(N);
    memcpy((void *)(mxGetPr(pa)), (void *)X.data(), sizeof((double *)X.data()));
    mxDestroyArray(pa);
    if (!matClose(pmat)) throw "CMATLAB::importFrom error: matlab.matClose failed";
}

double MeanNormalLikelihood::logValue(RefArrayXd modelParameters)
{
    unsigned long n = observations.size();
    double lambda0;
    double lambda;
    ArrayXd argument;
    ArrayXd predictions;

    predictions.resize(n);
    predictions.setZero();
    model.predict(predictions, modelParameters);
    argument = (observations - predictions);
    argument = argument.square()*weights;

    lambda0 = lgammal(n/2.) - log(2) - (n/2.)*log(Functions::PI) + 0.5*weights.log().sum();
    lambda = lambda0 - (n/2.)*log(argument.sum());

    return lambda;
}

//.........这里部分代码省略.........
    secondDerivatives[0] = 0.0;
    partialSolution[0] = 0.0;
    
    
    // Do tridiagonal decomposition for computing second derivatives of observed ordinate
    
    ArrayXd sigma = differenceAbscissa.segment(0,size-2)/differenceAbscissa2.segment(0,size-2);
    double beta;

    
    // Forward computation of partial solutions in tridiagonal system

    for (int i = 1; i < size-1; ++i)
    {
        beta = sigma(i-1) * secondDerivatives[i-1] + 2.0;
        secondDerivatives[i] = (sigma(i-1) - 1.0)/beta;
        partialSolution[i] = differenceOrdinate(i)/differenceAbscissa(i) - differenceOrdinate(i-1)/differenceAbscissa(i-1);
        partialSolution[i] = (6.0*partialSolution[i]/differenceAbscissa2(i-1)-sigma(i-1)*partialSolution[i-1])/beta;
    }

    
    // Upper bound condition for natural spline
    
    secondDerivatives[size-1] = 0.0;


    // Backward substitution

    for (int k = (size-2); k >= 0; --k)
    {
        secondDerivatives[k] = secondDerivatives[k]*secondDerivatives[k+1]+partialSolution[k];
    }


    // Initialize arrays of differences in both ordinate and abscissa
    
    ArrayXd interpolatedOrdinate = ArrayXd::Zero(interpolatedSize);
    ArrayXd remainingInterpolatedAbscissa = interpolatedAbscissa;       // The remaining part of the array of interpolated abscissa
    int cumulatedBinSize = 0;                                           // The cumulated number of interpolated points from the beginning
    int i = 0;                                                          // Bin counter

    while ((i < size-1) && (cumulatedBinSize < interpolatedSize))
    {
        // Find which values of interpolatedAbscissa are containined within the selected bin of observedAbscissa.
        // Since elements in interpolatedAbscissa are monotonically increasing, we cut the input array each time
        // we identify the elements of the current bin. This allows to speed up the process.

        double lowerAbscissa = observedAbscissa(i);
        double upperAbscissa = observedAbscissa(i+1);


        // Find total number of interpolated points falling in the current bin

        int binSize = Functions::countArrayIndicesWithinBoundaries(remainingInterpolatedAbscissa, lowerAbscissa, upperAbscissa);
      

        // Do interpolation only if interpolated points are found within the bin

        if (binSize > 0)
        {
            double lowerOrdinate = observedOrdinate(i);
            double upperOrdinate = observedOrdinate(i+1);
            double denominator = differenceAbscissa(i);
            double upperSecondDerivative = secondDerivatives[i+1];
            double lowerSecondDerivative = secondDerivatives[i];
            ArrayXd interpolatedAbscissaInCurrentBin = remainingInterpolatedAbscissa.segment(0, binSize);
            ArrayXd interpolatedOrdinateInCurrentBin = ArrayXd::Zero(binSize);


            // Compute coefficients for cubic spline interpolation function

            ArrayXd a = (upperAbscissa - interpolatedAbscissaInCurrentBin) / denominator;
            ArrayXd b = 1.0 - a;
            ArrayXd c = (1.0/6.0) * (a.cube() - a)*denominator*denominator;
            ArrayXd d = (1.0/6.0) * (b.cube() - b)*denominator*denominator;
            interpolatedOrdinateInCurrentBin = a*lowerOrdinate + b*upperOrdinate + c*lowerSecondDerivative + d*upperSecondDerivative;

                
            // Merge bin ordinate into total array of ordinate
        
            interpolatedOrdinate.segment(cumulatedBinSize, binSize) = interpolatedOrdinateInCurrentBin;


            // Reduce size of array remainingInterpolatedAbscissa by binSize elements and initialize the array
            // with remaining part of interpolatedAbscissa
        
            int currentRemainingSize = interpolatedSize - cumulatedBinSize;
            remainingInterpolatedAbscissa.resize(currentRemainingSize - binSize);
            cumulatedBinSize += binSize;
            remainingInterpolatedAbscissa = interpolatedAbscissa.segment(cumulatedBinSize, interpolatedSize-cumulatedBinSize);
        }  


        // Move to next bin

        ++i;
    }

    return interpolatedOrdinate;
}