static void fillsparsede(ap::real_1d_array& d,
     ap::real_1d_array& e,
     int n,
     double sparcity)
{
    int i;
    int j;

    d.setbounds(0, n-1);
    e.setbounds(0, ap::maxint(0, n-2));
    for(i = 0; i <= n-1; i++)
    {
        if( ap::fp_greater_eq(ap::randomreal(),sparcity) )
        {
            d(i) = 2*ap::randomreal()-1;
        }
        else
        {
            d(i) = 0;
        }
    }
    for(i = 0; i <= n-2; i++)
    {
        if( ap::fp_greater_eq(ap::randomreal(),sparcity) )
        {
            e(i) = 2*ap::randomreal()-1;
        }
        else
        {
            e(i) = 0;
        }
    }
}

/*************************************************************************
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;
    }
}

bool in_out_variable_1D(const ap::boolean_1d_array& in, const ap::real_1d_array& X, ap::real_1d_array& x, ap::real_1d_array& vector,  bool io)
{
// Routine to know the number of variables in/out
	int rows = in.gethighbound(0) + 1;
	int n_invar=0;
	bool flag;
	
	//ap::real_1d_array vector;
	vector.setbounds(0,rows-1);
	unsigned int k=0;
	for (int i=0; i<rows; i++)
	{
		if (in(i)==io) //to know how many variables are in/out
		{
			vector(k) = i;
			k++;
			n_invar++;
		}
	}
	if (n_invar>0)
	{
		// Routine to extract the in/out variables
		x.setbounds(0,n_invar-1);
		for (int i=0; i<n_invar; i++)
		x(i) = X(static_cast<int>(vector(i)));
	
        flag=TRUE;
	}
	else
		flag=FALSE;

	return flag;
}

//===============================
//=====================================================================
void newiteration_residOpticalFlow_mt(int iter, const ap::real_1d_array& x, double f,const ap::real_1d_array& g, void *params)
{
	globs_LBFGS_ *glob_param=(globs_LBFGS_*)params;
	double normG=0.0;
	for(int ii=g.getlowbound();ii<=g.gethighbound();ii++) normG+=g(ii)*g(ii);
	normG=sqrt(normG);
	cout<<"Iter="<<iter<<";fData="<<glob_param->fData<<";fSmooth="<<glob_param->fSmooth<<";fData+lambda*fSmooth="<<f<<";RMS(g)="<<normG/((double)(g.gethighbound()-g.getlowbound()+1))<<endl;
}

/*************************************************************************
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));
    }
}

void lbfgslincomb(const int& n,
     const double& da,
     const ap::real_1d_array& dx,
     int sx,
     ap::real_1d_array& dy,
     int sy)
{
    int fx;
    int fy;

    fx = sx+n-1;
    fy = sy+n-1;
    ap::vadd(dy.getvector(sy, fy), dx.getvector(sx, fx), da);
}

double lbfgsdotproduct(const int& n,
     const ap::real_1d_array& dx,
     int sx,
     const ap::real_1d_array& dy,
     int sy)
{
    double result;
    double v;
    int fx;
    int fy;

    fx = sx+n-1;
    fy = sy+n-1;
    v = ap::vdotproduct(dx.getvector(sx, fx), dy.getvector(sy, fy));
    result = v;
    return result;
}

/*************************************************************************
Serialization of LinearModel strucure

INPUT PARAMETERS:
    LM      -   original

OUTPUT PARAMETERS:
    RA      -   array of real numbers which stores model,
                array[0..RLen-1]
    RLen    -   RA lenght

  -- ALGLIB --
     Copyright 15.03.2009 by Bochkanov Sergey
*************************************************************************/
void lrserialize(const linearmodel& lm, ap::real_1d_array& ra, int& rlen)
{

    rlen = ap::round(lm.w(0))+1;
    ra.setbounds(0, rlen-1);
    ra(0) = lrvnum;
    ap::vmove(&ra(1), &lm.w(0), ap::vlen(1,rlen-1));
}

/*************************************************************************
Solving a system of linear equations with a system matrix given by its
LU decomposition.

The algorithm solves a system of linear equations whose matrix is given by
its LU decomposition. In case of a singular matrix, the algorithm  returns
False.

The algorithm solves systems with a square matrix only.

Input parameters:
    A       -   LU decomposition of a system matrix in compact  form  (the
                result of the RMatrixLU subroutine).
    Pivots  -   row permutation table (the result of a
                RMatrixLU subroutine).
    B       -   right side of a system.
                Array whose index ranges within [0..N-1].
    N       -   size of matrix A.

Output parameters:
    X       -   solution of a system.
                Array whose index ranges within [0..N-1].

Result:
    True, if the matrix is not singular.
    False, if the matrux is singular. In this case, X doesn't contain a
solution.

  -- ALGLIB --
     Copyright 2005-2008 by Bochkanov Sergey
*************************************************************************/
bool rmatrixlusolve(const ap::real_2d_array& a,
     const ap::integer_1d_array& pivots,
     ap::real_1d_array b,
     int n,
     ap::real_1d_array& x)
{
    bool result;
    ap::real_1d_array y;
    int i;
    int j;
    double v;

    y.setbounds(0, n-1);
    x.setbounds(0, n-1);
    result = true;
    for(i = 0; i <= n-1; i++)
    {
        if( a(i,i)==0 )
        {
            result = false;
            return result;
        }
    }
    
    //
    // pivots
    //
    for(i = 0; i <= n-1; i++)
    {
        if( pivots(i)!=i )
        {
            v = b(i);
            b(i) = b(pivots(i));
            b(pivots(i)) = v;
        }
    }
    
    //
    // Ly = b
    //
    y(0) = b(0);
    for(i = 1; i <= n-1; i++)
    {
        v = ap::vdotproduct(&a(i, 0), &y(0), ap::vlen(0,i-1));
        y(i) = b(i)-v;
    }
    
    //
    // Ux = y
    //
    x(n-1) = y(n-1)/a(n-1,n-1);
    for(i = n-2; i >= 0; i--)
    {
        v = ap::vdotproduct(&a(i, i+1), &x(i+1), ap::vlen(i+1,n-1));
        x(i) = (y(i)-v)/a(i,i);
    }
    return result;
}

/*************************************************************************
Unpacks coefficients of linear model.

INPUT PARAMETERS:
    LM          -   linear model in ALGLIB format

OUTPUT PARAMETERS:
    V           -   coefficients, array[0..NVars]
    NVars       -   number of independent variables (one less than number
                    of coefficients)

  -- ALGLIB --
     Copyright 30.08.2008 by Bochkanov Sergey
*************************************************************************/
void lrunpack(const linearmodel& lm, ap::real_1d_array& v, int& nvars)
{
    int offs;

    ap::ap_error::make_assertion(ap::round(lm.w(1))==lrvnum, "LINREG: Incorrect LINREG version!");
    nvars = ap::round(lm.w(2));
    offs = ap::round(lm.w(3));
    v.setbounds(0, nvars);
    ap::vmove(&v(0), &lm.w(offs), ap::vlen(0,nvars));
}

/*************************************************************************
Conversion of a series of Chebyshev polynomials to a power series.

Represents A[0]*T0(x) + A[1]*T1(x) + ... + A[N]*Tn(x) as
B[0] + B[1]*X + ... + B[N]*X^N.

Input parameters:
    A   -   Chebyshev series coefficients
    N   -   degree, N>=0
    
Output parameters
    B   -   power series coefficients
*************************************************************************/
void fromchebyshev(const ap::real_1d_array& a,
     const int& n,
     ap::real_1d_array& b)
{
    int i;
    int k;
    double e;
    double d;

    b.setbounds(0, n);
    for(i = 0; i <= n; i++)
    {
        b(i) = 0;
    }
    d = 0;
    i = 0;
    do
    {
        k = i;
        do
        {
            e = b(k);
            b(k) = 0;
            if( i<=1&&k==i )
            {
                b(k) = 1;
            }
            else
            {
                if( i!=0 )
                {
                    b(k) = 2*d;
                }
                if( k>i+1 )
                {
                    b(k) = b(k)-b(k-2);
                }
            }
            d = e;
            k = k+1;
        }
        while(k<=n);
        d = b(i);
        e = 0;
        k = i;
        while(k<=n)
        {
            e = e+b(k)*a(k);
            k = k+2;
        }
        b(i) = e;
        i = i+1;
    }
    while(i<=n);
}

/*************************************************************************
Representation of Ln as C[0] + C[1]*X + ... + C[N]*X^N

Input parameters:
    N   -   polynomial degree, n>=0

Output parameters:
    C   -   coefficients
*************************************************************************/
void laguerrecoefficients(const int& n, ap::real_1d_array& c)
{
    int i;

    c.setbounds(0, n);
    c(0) = 1;
    for(i = 0; i <= n-1; i++)
    {
        c(i+1) = -c(i)*(n-i)/(i+1)/(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));
    }
}

/*************************************************************************
L-BFGS algorithm results

Called after MinLBFGSIteration() returned False.

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

OUTPUT PARAMETERS:
    X       -   array[0..N-1], solution
    Rep     -   optimization report:
                * Rep.TerminationType completetion code:
                    * -2    rounding errors prevent further improvement.
                            X contains best point found.
                    * -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
                    *  7    stopping conditions are too stringent,
                            further improvement is impossible
                * Rep.IterationsCount contains iterations count
                * NFEV countains number of function calculations

  -- ALGLIB --
     Copyright 02.04.2010 by Bochkanov Sergey
*************************************************************************/
void minlbfgsresults(const minlbfgsstate& state,
     ap::real_1d_array& x,
     minlbfgsreport& rep)
{

    x.setbounds(0, state.n-1);
    ap::vmove(&x(0), 1, &state.x(0), 1, ap::vlen(0,state.n-1));
    rep.iterationscount = state.repiterationscount;
    rep.nfev = state.repnfev;
    rep.terminationtype = state.repterminationtype;
}

/*************************************************************************
Serialization of DecisionForest strucure

INPUT PARAMETERS:
    DF      -   original

OUTPUT PARAMETERS:
    RA      -   array of real numbers which stores decision forest,
                array[0..RLen-1]
    RLen    -   RA lenght

  -- ALGLIB --
     Copyright 13.02.2009 by Bochkanov Sergey
*************************************************************************/
void dfserialize(const decisionforest& df, ap::real_1d_array& ra, int& rlen)
{

    ra.setbounds(0, df.bufsize+5-1);
    ra(0) = dfvnum;
    ra(1) = df.nvars;
    ra(2) = df.nclasses;
    ra(3) = df.ntrees;
    ra(4) = df.bufsize;
    ap::vmove(&ra(5), 1, &df.trees(0), 1, ap::vlen(5,5+df.bufsize-1));
    rlen = 5+df.bufsize;
}

/*************************************************************************
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));
}

/*************************************************************************
Obsolete 1-based subroutine.
See RMatrixBDUnpackDiagonals for 0-based replacement.
*************************************************************************/
void unpackdiagonalsfrombidiagonal(const ap::real_2d_array& b,
     int m,
     int n,
     bool& isupper,
     ap::real_1d_array& d,
     ap::real_1d_array& e)
{
    int i;

    isupper = m>=n;
    if( m==0||n==0 )
    {
        return;
    }
    if( isupper )
    {
        d.setbounds(1, n);
        e.setbounds(1, n);
        for(i = 1; i <= n-1; i++)
        {
            d(i) = b(i,i);
            e(i) = b(i,i+1);
        }
        d(n) = b(n,n);
    }
    else
    {
        d.setbounds(1, m);
        e.setbounds(1, m);
        for(i = 1; i <= m-1; i++)
        {
            d(i) = b(i,i);
            e(i) = b(i+1,i);
        }
        d(m) = b(m,m);
    }
}

/*************************************************************************
Unpacking of the main and secondary diagonals of bidiagonal decomposition
of matrix A.

Input parameters:
    B   -   output of RMatrixBD subroutine.
    M   -   number of rows in matrix B.
    N   -   number of columns in matrix B.

Output parameters:
    IsUpper -   True, if the matrix is upper bidiagonal.
                otherwise IsUpper is False.
    D       -   the main diagonal.
                Array whose index ranges within [0..Min(M,N)-1].
    E       -   the secondary diagonal (upper or lower, depending on
                the value of IsUpper).
                Array index ranges within [0..Min(M,N)-1], the last
                element is not used.

  -- ALGLIB --
     Copyright 2005-2007 by Bochkanov Sergey
*************************************************************************/
void rmatrixbdunpackdiagonals(const ap::real_2d_array& b,
     int m,
     int n,
     bool& isupper,
     ap::real_1d_array& d,
     ap::real_1d_array& e)
{
    int i;

    isupper = m>=n;
    if( m<=0||n<=0 )
    {
        return;
    }
    if( isupper )
    {
        d.setbounds(0, n-1);
        e.setbounds(0, n-1);
        for(i = 0; i <= n-2; i++)
        {
            d(i) = b(i,i);
            e(i) = b(i,i+1);
        }
        d(n-1) = b(n-1,n-1);
    }
    else
    {
        d.setbounds(0, m-1);
        e.setbounds(0, m-1);
        for(i = 0; i <= m-2; i++)
        {
            d(i) = b(i,i);
            e(i) = b(i+1,i);
        }
        d(m-1) = b(m-1,m-1);
    }
}

/*************************************************************************
Representation of Hn as C[0] + C[1]*X + ... + C[N]*X^N

Input parameters:
    N   -   polynomial degree, n>=0

Output parameters:
    C   -   coefficients
*************************************************************************/
void hermitecoefficients(const int& n, ap::real_1d_array& c)
{
    int i;

    c.setbounds(0, n);
    for(i = 0; i <= n; i++)
    {
        c(i) = 0;
    }
    c(n) = exp(n*log(double(2)));
    for(i = 0; i <= n/2-1; i++)
    {
        c(n-2*(i+1)) = -c(n-2*i)*(n-2*i)*(n-2*i-1)/4/(i+1);
    }
}

/*************************************************************************
Multiclass Fisher LDA

Subroutine finds coefficients of linear combination which optimally separates
training set on classes.

INPUT PARAMETERS:
    XY          -   training set, array[0..NPoints-1,0..NVars].
                    First NVars columns store values of independent
                    variables, next column stores number of class (from 0
                    to NClasses-1) which dataset element belongs to. Fractional
                    values are rounded to nearest integer.
    NPoints     -   training set size, NPoints>=0
    NVars       -   number of independent variables, NVars>=1
    NClasses    -   number of classes, NClasses>=2


OUTPUT PARAMETERS:
    Info        -   return code:
                    * -4, if internal EVD subroutine hasn't converged
                    * -2, if there is a point with class number
                          outside of [0..NClasses-1].
                    * -1, if incorrect parameters was passed (NPoints<0,
                          NVars<1, NClasses<2)
                    *  1, if task has been solved
                    *  2, if there was a multicollinearity in training set,
                          but task has been solved.
    W           -   linear combination coefficients, array[0..NVars-1]

  -- ALGLIB --
     Copyright 31.05.2008 by Bochkanov Sergey
*************************************************************************/
void fisherlda(const ap::real_2d_array& xy,
     int npoints,
     int nvars,
     int nclasses,
     int& info,
     ap::real_1d_array& w)
{
    ap::real_2d_array w2;

    fisherldan(xy, npoints, nvars, nclasses, info, w2);
    if( info>0 )
    {
        w.setbounds(0, nvars-1);
        ap::vmove(&w(0), 1, &w2(0, 0), w2.getstride(), ap::vlen(0,nvars-1));
    }
}