/*************************************************************************
        1-dimensional complex 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: CorrC1D(Signal, Pattern) = Pattern x Signal (using  traditional
            definition of cross-correlation, denoting cross-correlation as "x").

        INPUT PARAMETERS
            Signal  -   array[0..N-1] - complex function to be transformed,
                        signal containing pattern
            N       -   problem size
            Pattern -   array[0..M-1] - complex 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(conj(pattern[j])*signal[i+j]
                        * negative lags are stored in R[N..N+M-2],
                          R[N+M-1-i] = sum(conj(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
        *************************************************************************/
        public static void corrc1d(ref AP.Complex[] signal,
            int n,
            ref AP.Complex[] pattern,
            int m,
            ref AP.Complex[] r)
        {
            AP.Complex[] p = new AP.Complex[0];
            AP.Complex[] b = new AP.Complex[0];
            int i = 0;
            int i_ = 0;
            int i1_ = 0;

            System.Diagnostics.Debug.Assert(n>0 & m>0, "CorrC1D: incorrect N or M!");
            p = new AP.Complex[m];
            for(i=0; i<=m-1; i++)
            {
                p[m-1-i] = AP.Math.Conj(pattern[i]);
            }
            conv.convc1d(ref p, m, ref signal, n, ref b);
            r = new AP.Complex[m+n-1];
            i1_ = (m-1) - (0);
            for(i_=0; i_<=n-1;i_++)
            {
                r[i_] = b[i_+i1_];
            }
            if( m+n-2>=n )
            {
                i1_ = (0) - (n);
                for(i_=n; i_<=m+n-2;i_++)
                {
                    r[i_] = b[i_+i1_];
                }
            }
        }

        /*************************************************************************
        Block size for complex subroutines.

          -- ALGLIB routine --
             15.12.2009
             Bochkanov Sergey
        *************************************************************************/
        public static int ablascomplexblocksize(ref AP.Complex[,] a)
        {
            int result = 0;

            result = 24;
            return result;
        }

 internal CompletionResult(string name, string completion, string doc, PythonMemberType memberType, AP.CompletionValue[] values) {
     _name = name;
     _memberType = memberType;
     _completion = completion;
     _doc = doc;
     _values = values;
 }

 public PythonParameter(ISignature signature, AP.Parameter param, Span locus, Span ppLocus, AnalysisVariable[] variables) {
     _signature = signature;
     _param = param;
     _locus = locus;
     _ppLocus = ppLocus;
     _documentation = _param.doc.LimitLines(15, stopAtFirstBlankLine: true);
     _variables = variables;
 }

        /*************************************************************************
        Calculation of the determinant of a general matrix

        Input parameters:
            A       -   matrix, array[0..N-1, 0..N-1]
            N       -   size of matrix A.

        Result: determinant of matrix A.

          -- ALGLIB --
             Copyright 2005 by Bochkanov Sergey
        *************************************************************************/
        public static AP.Complex cmatrixdet(AP.Complex[,] a,
            int n)
        {
            AP.Complex result = 0;
            int[] pivots = new int[0];

            a = (AP.Complex[,])a.Clone();

            trfac.cmatrixlu(ref a, n, n, ref pivots);
            result = cmatrixludet(ref a, ref pivots, n);
            return result;
        }

        /*************************************************************************
        Complex ABLASSplitLength

          -- ALGLIB routine --
             15.12.2009
             Bochkanov Sergey
        *************************************************************************/
        public static void ablascomplexsplitlength(ref AP.Complex[,] a,
            int n,
            ref int n1,
            ref int n2)
        {
            if( n>ablascomplexblocksize(ref a) )
            {
                ablasinternalsplitlength(n, ablascomplexblocksize(ref a), ref n1, ref n2);
            }
            else
            {
                ablasinternalsplitlength(n, ablasmicroblocksize(), ref n1, ref n2);
            }
        }

        /*************************************************************************
        Fast kernel

          -- ALGLIB routine --
             19.01.2010
             Bochkanov Sergey
        *************************************************************************/
        public static bool cmatrixrank1f(int m,
            int n,
            ref AP.Complex[,] a,
            int ia,
            int ja,
            ref AP.Complex[] u,
            int iu,
            ref AP.Complex[] v,
            int iv)
        {
            bool result = new bool();

            result = false;
            return result;
        }

        /*************************************************************************
        Application of an elementary reflection to a rectangular matrix of size MxN

        The  algorithm  pre-multiplies  the  matrix  by  an  elementary reflection
        transformation  which  is  given  by  column  V  and  scalar  Tau (see the
        description of the GenerateReflection). Not the whole matrix  but  only  a
        part of it is transformed (rows from M1 to M2, columns from N1 to N2). Only
        the elements of this submatrix are changed.

        Note: the matrix is multiplied by H, not by H'.   If  it  is  required  to
        multiply the matrix by H', it is necessary to pass Conj(Tau) instead of Tau.

        Input parameters:
            C       -   matrix to be transformed.
            Tau     -   scalar defining transformation.
            V       -   column defining transformation.
                        Array whose index ranges within [1..M2-M1+1]
            M1, M2  -   range of rows to be transformed.
            N1, N2  -   range of columns to be transformed.
            WORK    -   working array whose index goes from N1 to N2.

        Output parameters:
            C       -   the result of multiplying the input matrix C by the
                        transformation matrix which is given by Tau and V.
                        If N1>N2 or M1>M2, C is not modified.

          -- LAPACK auxiliary routine (version 3.0) --
             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
             Courant Institute, Argonne National Lab, and Rice University
             September 30, 1994
        *************************************************************************/
        public static void complexapplyreflectionfromtheleft(ref AP.Complex[,] c,
            AP.Complex tau,
            ref AP.Complex[] v,
            int m1,
            int m2,
            int n1,
            int n2,
            ref AP.Complex[] work)
        {
            AP.Complex t = 0;
            int i = 0;
            int vm = 0;
            int i_ = 0;

            if (tau == 0 | n1 > n2 | m1 > m2)
            {
                return;
            }

            //
            // w := C^T * conj(v)
            //
            vm = m2 - m1 + 1;
            for (i = n1; i <= n2; i++)
            {
                work[i] = 0;
            }
            for (i = m1; i <= m2; i++)
            {
                t = AP.Math.Conj(v[i + 1 - m1]);
                for (i_ = n1; i_ <= n2; i_++)
                {
                    work[i_] = work[i_] + t * c[i, i_];
                }
            }

            //
            // C := C - tau * v * w^T
            //
            for (i = m1; i <= m2; i++)
            {
                t = v[i - m1 + 1] * tau;
                for (i_ = n1; i_ <= n2; i_++)
                {
                    c[i, i_] = c[i, i_] - t * work[i_];
                }
            }
        }

        /*************************************************************************
        Fast kernel

          -- ALGLIB routine --
             19.01.2010
             Bochkanov Sergey
        *************************************************************************/
        public static bool cmatrixmvf(int m,
            int n,
            ref AP.Complex[,] a,
            int ia,
            int ja,
            int opa,
            ref AP.Complex[] x,
            int ix,
            ref AP.Complex[] y,
            int iy)
        {
            bool result = new bool();

            result = false;
            return result;
        }

        /*************************************************************************
        1-dimensional complex FFT.

        Array size N may be arbitrary number (composite or prime).  Composite  N's
        are handled with cache-oblivious variation of  a  Cooley-Tukey  algorithm.
        Small prime-factors are transformed using hard coded  codelets (similar to
        FFTW codelets, but without low-level  optimization),  large  prime-factors
        are handled with Bluestein's algorithm.

        Fastests transforms are for smooth N's (prime factors are 2, 3,  5  only),
        most fast for powers of 2. When N have prime factors  larger  than  these,
        but orders of magnitude smaller than N, computations will be about 4 times
        slower than for nearby highly composite N's. When N itself is prime, speed
        will be 6 times lower.

        Algorithm has O(N*logN) complexity for any N (composite or prime).

        INPUT PARAMETERS
            A   -   array[0..N-1] - complex function to be transformed
            N   -   problem size

        OUTPUT PARAMETERS
            A   -   DFT of a input array, array[0..N-1]
                    A_out[j] = SUM(A_in[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)

          -- ALGLIB --
             Copyright 29.05.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void fftc1d(ref AP.Complex[] a,
            int n)
        {
            ftbase.ftplan plan = new ftbase.ftplan();
            int i = 0;
            double[] buf = new double[0];

            System.Diagnostics.Debug.Assert(n>0, "FFTC1D: incorrect N!");

            //
            // Special case: N=1, FFT is just identity transform.
            // After this block we assume that N is strictly greater than 1.
            //
            if( n==1 )
            {
                return;
            }

            //
            // convert input array to the more convinient format
            //
            buf = new double[2*n];
            for(i=0; i<=n-1; i++)
            {
                buf[2*i+0] = a[i].x;
                buf[2*i+1] = a[i].y;
            }

            //
            // Generate plan and execute it.
            //
            // Plan is a combination of a successive factorizations of N and
            // precomputed data. It is much like a FFTW plan, but is not stored
            // between subroutine calls and is much simpler.
            //
            ftbase.ftbasegeneratecomplexfftplan(n, ref plan);
            ftbase.ftbaseexecuteplan(ref buf, 0, n, ref plan);

            //
            // result
            //
            for(i=0; i<=n-1; i++)
            {
                a[i].x = buf[2*i+0];
                a[i].y = buf[2*i+1];
            }
        }

        /*************************************************************************
        Fast kernel

          -- ALGLIB routine --
             19.01.2010
             Bochkanov Sergey
        *************************************************************************/
        public static bool cmatrixlefttrsmf(int m,
            int n,
            ref AP.Complex[,] a,
            int i1,
            int j1,
            bool isupper,
            bool isunit,
            int optype,
            ref AP.Complex[,] x,
            int i2,
            int j2)
        {
            bool result = new bool();

            result = false;
            return result;
        }

        /*************************************************************************
        1-dimensional complex convolution.

        For given A/B returns conv(A,B) (non-circular). Subroutine can automatically
        choose between three implementations: straightforward O(M*N)  formula  for
        very small N (or M), overlap-add algorithm for  cases  where  max(M,N)  is
        significantly larger than min(M,N), but O(M*N) algorithm is too slow,  and
        general FFT-based formula for cases where two previois algorithms are  too
        slow.

        Algorithm has max(M,N)*log(max(M,N)) complexity for any M/N.

        INPUT PARAMETERS
            A   -   array[0..M-1] - complex function to be transformed
            M   -   problem size
            B   -   array[0..N-1] - complex function to be transformed
            N   -   problem size

        OUTPUT PARAMETERS
            R   -   convolution: A*B. array[0..N+M-2].

        NOTE:
            It is assumed that A is zero at T<0, B is zero too.  If  one  or  both
        functions have non-zero values at negative T's, you  can  still  use  this
        subroutine - just shift its result correspondingly.

          -- ALGLIB --
             Copyright 21.07.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void convc1d(ref AP.Complex[] a,
            int m,
            ref AP.Complex[] b,
            int n,
            ref AP.Complex[] r)
        {
            System.Diagnostics.Debug.Assert(n>0 & m>0, "ConvC1D: incorrect N or M!");
            
            //
            // normalize task: make M>=N,
            // so A will be longer that B.
            //
            if( m<n )
            {
                convc1d(ref b, n, ref a, m, ref r);
                return;
            }
            convc1dx(ref a, m, ref b, n, false, -1, 0, ref r);
        }

        /// <summary>
        /// Updates the code applying the changes to the existing text buffer and updates the version.
        /// </summary>
        public static string UpdateCode(this IProjectEntry entry, AP.VersionChanges[] versions, int buffer, int version) {
            lock (_currentCodeKey) {
                CurrentCode curCode = GetCurrentCode(entry, buffer);
                var strBuffer = curCode.Text;

                foreach (var versionChange in versions) {
                    int delta = 0;

                    foreach (var change in versionChange.changes) {
                        strBuffer.Remove(change.start + delta, change.length);
                        strBuffer.Insert(change.start + delta, change.newText);

                        delta += change.newText.Length - change.length;
                    }
                }

                curCode.Version = version;
                return strBuffer.ToString();
            }
        }

        /*************************************************************************
        Determinant calculation of the matrix given by its LU decomposition.

        Input parameters:
            A       -   LU decomposition of the matrix (output of
                        RMatrixLU subroutine).
            Pivots  -   table of permutations which were made during
                        the LU decomposition.
                        Output of RMatrixLU subroutine.
            N       -   size of matrix A.

        Result: matrix determinant.

          -- ALGLIB --
             Copyright 2005 by Bochkanov Sergey
        *************************************************************************/
        public static AP.Complex cmatrixludet(ref AP.Complex[,] a,
            ref int[] pivots,
            int n)
        {
            AP.Complex result = 0;
            int i = 0;
            int s = 0;

            result = 1;
            s = 1;
            for (i = 0; i <= n - 1; i++)
            {
                result = result * a[i, i];
                if (pivots[i] != i)
                {
                    s = -s;
                }
            }
            result = result * s;
            return result;
        }

        /*************************************************************************
        Fast kernel

          -- ALGLIB routine --
             19.01.2010
             Bochkanov Sergey
        *************************************************************************/
        public static bool cmatrixgemmf(int m,
            int n,
            int k,
            AP.Complex alpha,
            ref AP.Complex[,] a,
            int ia,
            int ja,
            int optypea,
            ref AP.Complex[,] b,
            int ib,
            int jb,
            int optypeb,
            AP.Complex beta,
            ref AP.Complex[,] c,
            int ic,
            int jc)
        {
            bool result = new bool();

            result = false;
            return result;
        }

        /*************************************************************************
        Inversion of a Hermitian positive definite matrix.

        Given an upper or lower triangle of a Hermitian positive definite matrix,
        the algorithm generates matrix A^-1 and saves the upper or lower triangle
        depending on the input.

        Input parameters:
            A       -   matrix to be inverted (upper or lower triangle).
                        Array with elements [0..N-1,0..N-1].
            N       -   size of matrix A.
            IsUpper -   storage format.
                        If IsUpper = True, then the upper triangle of matrix A is
                        given, otherwise the lower triangle is given.

        Output parameters:
            Info    -   return code, same as in RMatrixLUInverse
            Rep     -   solver report, same as in RMatrixLUInverse
            A       -   inverse of matrix A, same as in RMatrixLUInverse

          -- ALGLIB routine --
             10.02.2010
             Bochkanov Sergey
        *************************************************************************/
        public static void hpdmatrixinverse(ref AP.Complex[,] a,
            int n,
            bool isupper,
            ref int info,
            ref matinvreport rep)
        {
            if (n < 1)
            {
                info = -1;
                return;
            }
            info = 1;
            if (trfac.hpdmatrixcholesky(ref a, n, isupper))
            {
                hpdmatrixcholeskyinverse(ref a, n, isupper, ref info, ref rep);
            }
            else
            {
                info = -3;
            }
        }

        /*************************************************************************
        Inversion of a Hermitian positive definite matrix which is given
        by Cholesky decomposition.

        Input parameters:
            A       -   Cholesky decomposition of the matrix to be inverted:
                        A=U�*U or A = L*L'.
                        Output of  HPDMatrixCholesky subroutine.
            N       -   size of matrix A.
            IsUpper �   storage format.
                        If IsUpper = True, then matrix A is given as A = U'*U
                        (matrix contains upper triangle).
                        Similarly, if IsUpper = False, then A = L*L'.

        Output parameters:
            Info    -   return code, same as in RMatrixLUInverse
            Rep     -   solver report, same as in RMatrixLUInverse
            A       -   inverse of matrix A, same as in RMatrixLUInverse

          -- ALGLIB routine --
             10.02.2010
             Bochkanov Sergey
        *************************************************************************/
        public static void hpdmatrixcholeskyinverse(ref AP.Complex[,] a,
            int n,
            bool isupper,
            ref int info,
            ref matinvreport rep)
        {
            int i = 0;
            int j = 0;
            int info2 = 0;
            matinvreport rep2 = new matinvreport();
            AP.Complex[] tmp = new AP.Complex[0];
            AP.Complex v = 0;

            if (n < 1)
            {
                info = -1;
                return;
            }
            info = 1;

            //
            // calculate condition numbers
            //
            rep.r1 = rcond.hpdmatrixcholeskyrcond(ref a, n, isupper);
            rep.rinf = rep.r1;
            if ((double)(rep.r1) < (double)(rcond.rcondthreshold()) | (double)(rep.rinf) < (double)(rcond.rcondthreshold()))
            {
                if (isupper)
                {
                    for (i = 0; i <= n - 1; i++)
                    {
                        for (j = i; j <= n - 1; j++)
                        {
                            a[i, j] = 0;
                        }
                    }
                }
                else
                {
                    for (i = 0; i <= n - 1; i++)
                    {
                        for (j = 0; j <= i; j++)
                        {
                            a[i, j] = 0;
                        }
                    }
                }
                rep.r1 = 0;
                rep.rinf = 0;
                info = -3;
                return;
            }

            //
            // Inverse
            //
            tmp = new AP.Complex[n];
            hpdmatrixcholeskyinverserec(ref a, 0, n, isupper, ref tmp);
        }

        /*************************************************************************
        Inversion of a general matrix.

        Input parameters:
            A   -   matrix, array[0..N-1,0..N-1].
            N   -   size of A.

        Output parameters:
            Info    -   return code, same as in RMatrixLUInverse
            Rep     -   solver report, same as in RMatrixLUInverse
            A       -   inverse of matrix A, same as in RMatrixLUInverse

          -- ALGLIB --
             Copyright 2005 by Bochkanov Sergey
        *************************************************************************/
        public static void cmatrixinverse(ref AP.Complex[,] a,
            int n,
            ref int info,
            ref matinvreport rep)
        {
            int[] pivots = new int[0];

            trfac.cmatrixlu(ref a, n, n, ref pivots);
            cmatrixluinverse(ref a, ref pivots, n, ref info, ref rep);
        }

        /*************************************************************************
        Inversion of a matrix given by its LU decomposition.

        INPUT PARAMETERS:
            A       -   LU decomposition of the matrix (output of CMatrixLU subroutine).
            Pivots  -   table of permutations which were made during the LU decomposition
                        (the output of CMatrixLU subroutine).
            N       -   size of matrix A.

        OUTPUT PARAMETERS:
            Info    -   return code, same as in RMatrixLUInverse
            Rep     -   solver report, same as in RMatrixLUInverse
            A       -   inverse of matrix A, same as in RMatrixLUInverse

          -- ALGLIB routine --
             05.02.2010
             Bochkanov Sergey
        *************************************************************************/
        public static void cmatrixluinverse(ref AP.Complex[,] a,
            ref int[] pivots,
            int n,
            ref int info,
            ref matinvreport rep)
        {
            AP.Complex[] work = new AP.Complex[0];
            int i = 0;
            int j = 0;
            int k = 0;
            AP.Complex v = 0;

            info = 1;

            //
            // Quick return if possible
            //
            if (n == 0)
            {
                info = -1;
                return;
            }
            for (i = 0; i <= n - 1; i++)
            {
                if (pivots[i] > n - 1 | pivots[i] < i)
                {
                    info = -1;
                    return;
                }
            }

            //
            // calculate condition numbers
            //
            rep.r1 = rcond.cmatrixlurcond1(ref a, n);
            rep.rinf = rcond.cmatrixlurcondinf(ref a, n);
            if ((double)(rep.r1) < (double)(rcond.rcondthreshold()) | (double)(rep.rinf) < (double)(rcond.rcondthreshold()))
            {
                for (i = 0; i <= n - 1; i++)
                {
                    for (j = 0; j <= n - 1; j++)
                    {
                        a[i, j] = 0;
                    }
                }
                rep.r1 = 0;
                rep.rinf = 0;
                info = -3;
                return;
            }

            //
            // Call cache-oblivious code
            //
            work = new AP.Complex[n];
            cmatrixluinverserec(ref a, 0, n, ref work, ref info, ref rep);

            //
            // apply permutations
            //
            for (i = 0; i <= n - 1; i++)
            {
                for (j = n - 2; j >= 0; j--)
                {
                    k = pivots[j];
                    v = a[i, j];
                    a[i, j] = a[i, k];
                    a[i, k] = v;
                }
            }
        }

        /*************************************************************************
        Recursive subroutine for HPD inversion.

          -- ALGLIB routine --
             10.02.2010
             Bochkanov Sergey
        *************************************************************************/
        private static void hpdmatrixcholeskyinverserec(ref AP.Complex[,] a,
            int offs,
            int n,
            bool isupper,
            ref AP.Complex[] tmp)
        {
            int i = 0;
            int j = 0;
            AP.Complex v = 0;
            int n1 = 0;
            int n2 = 0;
            int info2 = 0;
            matinvreport rep2 = new matinvreport();
            int i_ = 0;
            int i1_ = 0;

            if (n < 1)
            {
                return;
            }

            //
            // Base case
            //
            if (n <= ablas.ablascomplexblocksize(ref a))
            {
                cmatrixtrinverserec(ref a, offs, n, isupper, false, ref tmp, ref info2, ref rep2);
                if (isupper)
                {

                    //
                    // Compute the product U * U'.
                    // NOTE: we never assume that diagonal of U is real
                    //
                    for (i = 0; i <= n - 1; i++)
                    {
                        if (i == 0)
                        {

                            //
                            // 1x1 matrix
                            //
                            a[offs + i, offs + i] = AP.Math.Sqr(a[offs + i, offs + i].x) + AP.Math.Sqr(a[offs + i, offs + i].y);
                        }
                        else
                        {

                            //
                            // (I+1)x(I+1) matrix,
                            //
                            // ( A11  A12 )   ( A11^H        )   ( A11*A11^H+A12*A12^H  A12*A22^H )
                            // (          ) * (              ) = (                                )
                            // (      A22 )   ( A12^H  A22^H )   ( A22*A12^H            A22*A22^H )
                            //
                            // A11 is IxI, A22 is 1x1.
                            //
                            i1_ = (offs) - (0);
                            for (i_ = 0; i_ <= i - 1; i_++)
                            {
                                tmp[i_] = AP.Math.Conj(a[i_ + i1_, offs + i]);
                            }
                            for (j = 0; j <= i - 1; j++)
                            {
                                v = a[offs + j, offs + i];
                                i1_ = (j) - (offs + j);
                                for (i_ = offs + j; i_ <= offs + i - 1; i_++)
                                {
                                    a[offs + j, i_] = a[offs + j, i_] + v * tmp[i_ + i1_];
                                }
                            }
                            v = AP.Math.Conj(a[offs + i, offs + i]);
                            for (i_ = offs; i_ <= offs + i - 1; i_++)
                            {
                                a[i_, offs + i] = v * a[i_, offs + i];
                            }
                            a[offs + i, offs + i] = AP.Math.Sqr(a[offs + i, offs + i].x) + AP.Math.Sqr(a[offs + i, offs + i].y);
                        }
                    }
                }
                else
                {

                    //
                    // Compute the product L' * L
                    // NOTE: we never assume that diagonal of L is real
                    //
                    for (i = 0; i <= n - 1; i++)
                    {
                        if (i == 0)
                        {

                            //
                            // 1x1 matrix
//.........这里部分代码省略.........