C++ real_1d_array::setlength方法代码示例

/*************************************************************************
This function generates  1-dimensional equidistant interpolation task with
moderate Lipshitz constant (close to 1.0)

If N=1 then suborutine generates only one point at the middle of [A,B]

  -- ALGLIB --
     Copyright 02.12.2009 by Bochkanov Sergey
*************************************************************************/
void taskgenint1dequidist(double a,
     double b,
     int n,
     ap::real_1d_array& x,
     ap::real_1d_array& y)
{
    int i;
    double h;

    ap::ap_error::make_assertion(n>=1, "TaskGenInterpolationEqdist1D: N<1!");
    x.setlength(n);
    y.setlength(n);
    if( n>1 )
    {
        x(0) = a;
        y(0) = 2*ap::randomreal()-1;
        h = (b-a)/(n-1);
        for(i = 1; i <= n-1; i++)
        {
            x(i) = a+i*h;
            y(i) = y(i-1)+(2*ap::randomreal()-1)*h;
        }
    }
    else
    {
        x(0) = 0.5*(a+b);
        y(0) = 2*ap::randomreal()-1;
    }
}

/*************************************************************************
Computation of nodes and weights for a Gauss quadrature formula

The algorithm generates the N-point Gauss quadrature formula  with  weight
function given by coefficients alpha and beta  of  a  recurrence  relation
which generates a system of orthogonal polynomials:

P-1(x)   =  0
P0(x)    =  1
Pn+1(x)  =  (x-alpha(n))*Pn(x)  -  beta(n)*Pn-1(x)

and zeroth moment Mu0

Mu0 = integral(W(x)dx,a,b)

INPUT PARAMETERS:
    Alpha   –   array[0..N-1], alpha coefficients
    Beta    –   array[0..N-1], beta coefficients
                Zero-indexed element is not used and may be arbitrary.
                Beta[I]>0.
    Mu0     –   zeroth moment of the weight function.
    N       –   number of nodes of the quadrature formula, N>=1

OUTPUT PARAMETERS:
    Info    -   error code:
                * -3    internal eigenproblem solver hasn't converged
                * -2    Beta[i]<=0
                * -1    incorrect N was passed
                *  1    OK
    X       -   array[0..N-1] - array of quadrature nodes,
                in ascending order.
    W       -   array[0..N-1] - array of quadrature weights.

  -- ALGLIB --
     Copyright 2005-2009 by Bochkanov Sergey
*************************************************************************/
void gqgeneraterec(const ap::real_1d_array& alpha,
     const ap::real_1d_array& beta,
     double mu0,
     int n,
     int& info,
     ap::real_1d_array& x,
     ap::real_1d_array& w)
{
    int i;
    ap::real_1d_array d;
    ap::real_1d_array e;
    ap::real_2d_array z;

    if( n<1 )
    {
        info = -1;
        return;
    }
    info = 1;
    
    //
    // Initialize
    //
    d.setlength(n);
    e.setlength(n);
    for(i = 1; i <= n-1; i++)
    {
        d(i-1) = alpha(i-1);
        if( ap::fp_less_eq(beta(i),0) )
        {
            info = -2;
            return;
        }
        e(i-1) = sqrt(beta(i));
    }
    d(n-1) = alpha(n-1);
    
    //
    // EVD
    //
    if( !smatrixtdevd(d, e, n, 3, z) )
    {
        info = -3;
        return;
    }
    
    //
    // Generate
    //
    x.setlength(n);
    w.setlength(n);
    for(i = 1; i <= n; i++)
    {
        x(i-1) = d(i-1);
        w(i-1) = mu0*ap::sqr(z(0,i-1));
    }
}

/*************************************************************************
1-dimensional circular real cross-correlation.

For given Pattern/Signal returns corr(Pattern,Signal) (circular).
Algorithm has linearithmic complexity for any M/N.

IMPORTANT:
    for  historical reasons subroutine accepts its parameters in  reversed
    order:   CorrR1DCircular(Signal, Pattern) = Pattern x Signal    (using
    traditional definition of cross-correlation, denoting cross-correlation
    as "x").

INPUT PARAMETERS
    Signal  -   array[0..N-1] - real function to be transformed,
                periodic signal containing pattern
    N       -   problem size
    Pattern -   array[0..M-1] - real function to be transformed,
                non-periodic pattern to search withing signal
    M       -   problem size

OUTPUT PARAMETERS
    R   -   convolution: A*B. array[0..M-1].


  -- ALGLIB --
     Copyright 21.07.2009 by Bochkanov Sergey
*************************************************************************/
void corrr1dcircular(const ap::real_1d_array& signal,
     int m,
     const ap::real_1d_array& pattern,
     int n,
     ap::real_1d_array& c)
{
    ap::real_1d_array p;
    ap::real_1d_array b;
    int i1;
    int i2;
    int i;
    int j2;

    ap::ap_error::make_assertion(n>0&&m>0, "ConvC1DCircular: incorrect N or M!");
    
    //
    // normalize task: make M>=N,
    // so A will be longer (at least - not shorter) that B.
    //
    if( m<n )
    {
        b.setlength(m);
        for(i1 = 0; i1 <= m-1; i1++)
        {
            b(i1) = 0;
        }
        i1 = 0;
        while(i1<n)
        {
            i2 = ap::minint(i1+m-1, n-1);
            j2 = i2-i1;
            ap::vadd(&b(0), &pattern(i1), ap::vlen(0,j2));
            i1 = i1+m;
        }
        corrr1dcircular(signal, m, b, m, c);
        return;
    }
    
    //
    // Task is normalized
    //
    p.setlength(n);
    for(i = 0; i <= n-1; i++)
    {
        p(n-1-i) = pattern(i);
    }
    convr1dcircular(signal, m, p, n, b);
    c.setlength(m);
    ap::vmove(&c(0), &b(n-1), ap::vlen(0,m-n));
    if( m-n+1<=m-1 )
    {
        ap::vmove(&c(m-n+1), &b(0), ap::vlen(m-n+1,m-1));
    }
}

/*************************************************************************
Nonlinear least squares fitting results.

Called after LSFitNonlinearIteration() returned False.

INPUT PARAMETERS:
    State   -   algorithm state (used by LSFitNonlinearIteration).

OUTPUT PARAMETERS:
    Info    -   completetion code:
                    * -1    incorrect parameters were specified
                    *  1    relative function improvement is no more than
                            EpsF.
                    *  2    relative step is no more than EpsX.
                    *  4    gradient norm is no more than EpsG
                    *  5    MaxIts steps was taken
    C       -   array[0..K-1], solution
    Rep     -   optimization report. Following fields are set:
                * Rep.TerminationType completetion code:
                * RMSError          rms error on the (X,Y).
                * AvgError          average error on the (X,Y).
                * AvgRelError       average relative error on the non-zero Y
                * MaxError          maximum error
                                    NON-WEIGHTED ERRORS ARE CALCULATED


  -- ALGLIB --
     Copyright 17.08.2009 by Bochkanov Sergey
*************************************************************************/
void lsfitnonlinearresults(const lsfitstate& state,
     int& info,
     ap::real_1d_array& c,
     lsfitreport& rep)
{

    info = state.repterminationtype;
    if( info>0 )
    {
        c.setlength(state.k);
        ap::vmove(&c(0), 1, &state.c(0), 1, ap::vlen(0,state.k-1));
        rep.rmserror = state.reprmserror;
        rep.avgerror = state.repavgerror;
        rep.avgrelerror = state.repavgrelerror;
        rep.maxerror = state.repmaxerror;
    }
}

/*************************************************************************
Dense solver.

Similar to RMatrixSolveM() but solves task with one right part  (where b/x
are vectors, not matrices).

See RMatrixSolveM()  description  for  more  information  about subroutine
parameters.

  -- ALGLIB --
     Copyright 24.08.2009 by Bochkanov Sergey
*************************************************************************/
void rmatrixsolve(const ap::real_2d_array& a,
     int n,
     const ap::real_1d_array& b,
     int& info,
     densesolverreport& rep,
     ap::real_1d_array& x)
{
    ap::real_2d_array bm;
    ap::real_2d_array xm;

    if( n<=0 )
    {
        info = -1;
        return;
    }
    bm.setlength(n, 1);
    ap::vmove(bm.getcolumn(0, 0, n-1), b.getvector(0, n-1));
    rmatrixsolvem(a, n, bm, 1, info, rep, xm);
    x.setlength(n);
    ap::vmove(x.getvector(0, n-1), xm.getcolumn(0, 0, n-1));
}

/*************************************************************************
1-dimensional real cross-correlation.

For given Pattern/Signal returns corr(Pattern,Signal) (non-circular).

Correlation is calculated using reduction to  convolution.  Algorithm with
max(N,N)*log(max(N,N)) complexity is used (see  ConvC1D()  for  more  info
about performance).

IMPORTANT:
    for  historical reasons subroutine accepts its parameters in  reversed
    order: CorrR1D(Signal, Pattern) = Pattern x Signal (using  traditional
    definition of cross-correlation, denoting cross-correlation as "x").

INPUT PARAMETERS
    Signal  -   array[0..N-1] - real function to be transformed,
                signal containing pattern
    N       -   problem size
    Pattern -   array[0..M-1] - real function to be transformed,
                pattern to search withing signal
    M       -   problem size

OUTPUT PARAMETERS
    R       -   cross-correlation, array[0..N+M-2]:
                * positive lags are stored in R[0..N-1],
                  R[i] = sum(pattern[j]*signal[i+j]
                * negative lags are stored in R[N..N+M-2],
                  R[N+M-1-i] = sum(pattern[j]*signal[-i+j]

NOTE:
    It is assumed that pattern domain is [0..M-1].  If Pattern is non-zero
on [-K..M-1],  you can still use this subroutine, just shift result by K.

  -- ALGLIB --
     Copyright 21.07.2009 by Bochkanov Sergey
*************************************************************************/
void corrr1d(const ap::real_1d_array& signal,
     int n,
     const ap::real_1d_array& pattern,
     int m,
     ap::real_1d_array& r)
{
    ap::real_1d_array p;
    ap::real_1d_array b;
    int i;

    ap::ap_error::make_assertion(n>0&&m>0, "CorrR1D: incorrect N or M!");
    p.setlength(m);
    for(i = 0; i <= m-1; i++)
    {
        p(m-1-i) = pattern(i);
    }
    convr1d(p, m, signal, n, b);
    r.setlength(m+n-1);
    ap::vmove(&r(0), &b(m-1), ap::vlen(0,n-1));
    if( m+n-2>=n )
    {
        ap::vmove(&r(n), &b(0), ap::vlen(n,m+n-2));
    }
}

//.........这里部分代码省略.........
    double verr;
    bool svdfailed;
    bool zeroa;
    int rfs;
    int nrfs;
    bool terminatenexttime;
    bool smallerr;

    if( nrows<=0||ncols<=0||ap::fp_less(threshold,0) )
    {
        info = -1;
        return;
    }
    if( ap::fp_eq(threshold,0) )
    {
        threshold = 1000*ap::machineepsilon;
    }
    
    //
    // Factorize A first
    //
    svdfailed = !rmatrixsvd(a, nrows, ncols, 1, 2, 2, sv, u, vt);
    zeroa = ap::fp_eq(sv(0),0);
    if( svdfailed||zeroa )
    {
        if( svdfailed )
        {
            info = -4;
        }
        else
        {
            info = 1;
        }
        x.setlength(ncols);
        for(i = 0; i <= ncols-1; i++)
        {
            x(i) = 0;
        }
        rep.n = ncols;
        rep.k = ncols;
        rep.cx.setlength(ncols, ncols);
        for(i = 0; i <= ncols-1; i++)
        {
            for(j = 0; j <= ncols-1; j++)
            {
                if( i==j )
                {
                    rep.cx(i,j) = 1;
                }
                else
                {
                    rep.cx(i,j) = 0;
                }
            }
        }
        rep.r2 = 0;
        return;
    }
    nsv = ap::minint(ncols, nrows);
    if( nsv==ncols )
    {
        rep.r2 = sv(nsv-1)/sv(0);
    }
    else
    {
        rep.r2 = 0;

//.........这里部分代码省略.........
    ap::real_1d_array tau;
    ap::real_2d_array q;
    ap::real_2d_array f2;
    ap::real_1d_array tmp;
    ap::real_1d_array c0;
    double v;

    if( n<1||m<1||k<0 )
    {
        info = -1;
        return;
    }
    if( k>=m )
    {
        info = -3;
        return;
    }
    
    //
    // Solve
    //
    if( k==0 )
    {
        
        //
        // no constraints
        //
        lsfitlinearinternal(y, w, fmatrix, n, m, info, c, rep);
    }
    else
    {
        
        //
        // First, find general form solution of constraints system:
        // * factorize C = L*Q
        // * unpack Q
        // * fill upper part of C with zeros (for RCond)
        //
        // We got C=C0+Q2'*y where Q2 is lower M-K rows of Q.
        //
        rmatrixlq(cmatrix, k, m, tau);
        rmatrixlqunpackq(cmatrix, k, m, tau, m, q);
        for(i = 0; i <= k-1; i++)
        {
            for(j = i+1; j <= m-1; j++)
            {
                cmatrix(i,j) = 0.0;
            }
        }
        if( ap::fp_less(rmatrixlurcondinf(cmatrix, k),1000*ap::machineepsilon) )
        {
            info = -3;
            return;
        }
        tmp.setlength(k);
        for(i = 0; i <= k-1; i++)
        {
            if( i>0 )
            {
                v = ap::vdotproduct(&cmatrix(i, 0), 1, &tmp(0), 1, ap::vlen(0,i-1));
            }
            else
            {
                v = 0;
            }
            tmp(i) = (cmatrix(i,m)-v)/cmatrix(i,i);
        }
        c0.setlength(m);
        for(i = 0; i <= m-1; i++)
        {
            c0(i) = 0;
        }
        for(i = 0; i <= k-1; i++)
        {
            v = tmp(i);
            ap::vadd(&c0(0), 1, &q(i, 0), 1, ap::vlen(0,m-1), v);
        }
        
        //
        // Second, prepare modified matrix F2 = F*Q2' and solve modified task
        //
        tmp.setlength(ap::maxint(n, m)+1);
        f2.setlength(n, m-k);
        matrixvectormultiply(fmatrix, 0, n-1, 0, m-1, false, c0, 0, m-1, -1.0, y, 0, n-1, 1.0);
        matrixmatrixmultiply(fmatrix, 0, n-1, 0, m-1, false, q, k, m-1, 0, m-1, true, 1.0, f2, 0, n-1, 0, m-k-1, 0.0, tmp);
        lsfitlinearinternal(y, w, f2, n, m-k, info, tmp, rep);
        rep.taskrcond = -1;
        if( info<=0 )
        {
            return;
        }
        
        //
        // then, convert back to original answer: C = C0 + Q2'*Y0
        //
        c.setlength(m);
        ap::vmove(&c(0), 1, &c0(0), 1, ap::vlen(0,m-1));
        matrixvectormultiply(q, k, m-1, 0, m-1, true, tmp, 0, m-k-1, 1.0, c, 0, m-1, 1.0);
    }
}

/*************************************************************************
Internal fitting subroutine
*************************************************************************/
static void lsfitlinearinternal(const ap::real_1d_array& y,
     const ap::real_1d_array& w,
     const ap::real_2d_array& fmatrix,
     int n,
     int m,
     int& info,
     ap::real_1d_array& c,
     lsfitreport& rep)
{
    double threshold;
    ap::real_2d_array ft;
    ap::real_2d_array q;
    ap::real_2d_array l;
    ap::real_2d_array r;
    ap::real_1d_array b;
    ap::real_1d_array wmod;
    ap::real_1d_array tau;
    int i;
    int j;
    double v;
    ap::real_1d_array sv;
    ap::real_2d_array u;
    ap::real_2d_array vt;
    ap::real_1d_array tmp;
    ap::real_1d_array utb;
    ap::real_1d_array sutb;
    int relcnt;

    if( n<1||m<1 )
    {
        info = -1;
        return;
    }
    info = 1;
    threshold = sqrt(ap::machineepsilon);
    
    //
    // Degenerate case, needs special handling
    //
    if( n<m )
    {
        
        //
        // Create design matrix.
        //
        ft.setlength(n, m);
        b.setlength(n);
        wmod.setlength(n);
        for(j = 0; j <= n-1; j++)
        {
            v = w(j);
            ap::vmove(&ft(j, 0), 1, &fmatrix(j, 0), 1, ap::vlen(0,m-1), v);
            b(j) = w(j)*y(j);
            wmod(j) = 1;
        }
        
        //
        // LQ decomposition and reduction to M=N
        //
        c.setlength(m);
        for(i = 0; i <= m-1; i++)
        {
            c(i) = 0;
        }
        rep.taskrcond = 0;
        rmatrixlq(ft, n, m, tau);
        rmatrixlqunpackq(ft, n, m, tau, n, q);
        rmatrixlqunpackl(ft, n, m, l);
        lsfitlinearinternal(b, wmod, l, n, n, info, tmp, rep);
        if( info<=0 )
        {
            return;
        }
        for(i = 0; i <= n-1; i++)
        {
            v = tmp(i);
            ap::vadd(&c(0), 1, &q(i, 0), 1, ap::vlen(0,m-1), v);
        }
        return;
    }
    
    //
    // N>=M. Generate design matrix and reduce to N=M using
    // QR decomposition.
    //
    ft.setlength(n, m);
    b.setlength(n);
    for(j = 0; j <= n-1; j++)
    {
        v = w(j);
        ap::vmove(&ft(j, 0), 1, &fmatrix(j, 0), 1, ap::vlen(0,m-1), v);
        b(j) = w(j)*y(j);
    }
    rmatrixqr(ft, n, m, tau);
    rmatrixqrunpackq(ft, n, m, tau, m, q);
    rmatrixqrunpackr(ft, n, m, r);
    tmp.setlength(m);
//.........这里部分代码省略.........

void lsfitscalexy(ap::real_1d_array& x,
     ap::real_1d_array& y,
     int n,
     ap::real_1d_array& xc,
     ap::real_1d_array& yc,
     const ap::integer_1d_array& dc,
     int k,
     double& xa,
     double& xb,
     double& sa,
     double& sb,
     ap::real_1d_array& xoriginal,
     ap::real_1d_array& yoriginal)
{
    double xmin;
    double xmax;
    int i;

    ap::ap_error::make_assertion(n>=1, "LSFitScaleXY: incorrect N");
    ap::ap_error::make_assertion(k>=0, "LSFitScaleXY: incorrect K");
    
    //
    // Calculate xmin/xmax.
    // Force xmin<>xmax.
    //
    xmin = x(0);
    xmax = x(0);
    for(i = 1; i <= n-1; i++)
    {
        xmin = ap::minreal(xmin, x(i));
        xmax = ap::maxreal(xmax, x(i));
    }
    for(i = 0; i <= k-1; i++)
    {
        xmin = ap::minreal(xmin, xc(i));
        xmax = ap::maxreal(xmax, xc(i));
    }
    if( ap::fp_eq(xmin,xmax) )
    {
        if( ap::fp_eq(xmin,0) )
        {
            xmin = -1;
            xmax = +1;
        }
        else
        {
            xmin = 0.5*xmin;
        }
    }
    
    //
    // Transform abscissas: map [XA,XB] to [0,1]
    //
    // Store old X[] in XOriginal[] (it will be used
    // to calculate relative error).
    //
    xoriginal.setlength(n);
    ap::vmove(&xoriginal(0), 1, &x(0), 1, ap::vlen(0,n-1));
    xa = xmin;
    xb = xmax;
    for(i = 0; i <= n-1; i++)
    {
        x(i) = 2*(x(i)-0.5*(xa+xb))/(xb-xa);
    }
    for(i = 0; i <= k-1; i++)
    {
        ap::ap_error::make_assertion(dc(i)>=0, "LSFitScaleXY: internal error!");
        xc(i) = 2*(xc(i)-0.5*(xa+xb))/(xb-xa);
        yc(i) = yc(i)*pow(0.5*(xb-xa), double(dc(i)));
    }
    
    //
    // Transform function values: map [SA,SB] to [0,1]
    // SA = mean(Y),
    // SB = SA+stddev(Y).
    //
    // Store old Y[] in YOriginal[] (it will be used
    // to calculate relative error).
    //
    yoriginal.setlength(n);
    ap::vmove(&yoriginal(0), 1, &y(0), 1, ap::vlen(0,n-1));
    sa = 0;
    for(i = 0; i <= n-1; i++)
    {
        sa = sa+y(i);
    }
    sa = sa/n;
    sb = 0;
    for(i = 0; i <= n-1; i++)
    {
        sb = sb+ap::sqr(y(i)-sa);
    }
    sb = sqrt(sb/n)+sa;
    if( ap::fp_eq(sb,sa) )
    {
        sb = 2*sa;
    }
    if( ap::fp_eq(sb,sa) )
    {
        sb = sa+1;
//.........这里部分代码省略.........

/*************************************************************************
Computation of nodes and weights for a Gauss-Lobatto quadrature formula

The algorithm generates the N-point Gauss-Lobatto quadrature formula  with
weight function given by coefficients alpha and beta of a recurrence which
generates a system of orthogonal polynomials.

P-1(x)   =  0
P0(x)    =  1
Pn+1(x)  =  (x-alpha(n))*Pn(x)  -  beta(n)*Pn-1(x)

and zeroth moment Mu0

Mu0 = integral(W(x)dx,a,b)

INPUT PARAMETERS:
    Alpha   –   array[0..N-2], alpha coefficients
    Beta    –   array[0..N-2], beta coefficients.
                Zero-indexed element is not used, may be arbitrary.
                Beta[I]>0
    Mu0     –   zeroth moment of the weighting function.
    A       –   left boundary of the integration interval.
    B       –   right boundary of the integration interval.
    N       –   number of nodes of the quadrature formula, N>=3
                (including the left and right boundary nodes).

OUTPUT PARAMETERS:
    Info    -   error code:
                * -3    internal eigenproblem solver hasn't converged
                * -2    Beta[i]<=0
                * -1    incorrect N was passed
                *  1    OK
    X       -   array[0..N-1] - array of quadrature nodes,
                in ascending order.
    W       -   array[0..N-1] - array of quadrature weights.

  -- ALGLIB --
     Copyright 2005-2009 by Bochkanov Sergey
*************************************************************************/
void gqgenerategausslobattorec(ap::real_1d_array alpha,
     ap::real_1d_array beta,
     double mu0,
     double a,
     double b,
     int n,
     int& info,
     ap::real_1d_array& x,
     ap::real_1d_array& w)
{
    int i;
    ap::real_1d_array d;
    ap::real_1d_array e;
    ap::real_2d_array z;
    double pim1a;
    double pia;
    double pim1b;
    double pib;
    double t;
    double a11;
    double a12;
    double a21;
    double a22;
    double b1;
    double b2;
    double alph;
    double bet;

    if( n<=2 )
    {
        info = -1;
        return;
    }
    info = 1;
    
    //
    // Initialize, D[1:N+1], E[1:N]
    //
    n = n-2;
    d.setlength(n+2);
    e.setlength(n+1);
    for(i = 1; i <= n+1; i++)
    {
        d(i-1) = alpha(i-1);
    }
    for(i = 1; i <= n; i++)
    {
        if( ap::fp_less_eq(beta(i),0) )
        {
            info = -2;
            return;
        }
        e(i-1) = sqrt(beta(i));
    }
    
    //
    // Caclulate Pn(a), Pn+1(a), Pn(b), Pn+1(b)
    //
    beta(0) = 0;
    pim1a = 0;
    pia = 1;
//.........这里部分代码省略.........