updated to latest version of psdlag, compiles well provided libalglib-dev is installed, looks to still be reading headers from the src/inc directory, this could cause problems

src/inc/alglib/alglibmisc.h

@@ -0,0 +1,769 @@
+/*************************************************************************
+Copyright (c) Sergey Bochkanov (ALGLIB project).
+
+>>> SOURCE LICENSE >>>
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation (www.fsf.org); either version 2 of the
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+A copy of the GNU General Public License is available at
+http://www.fsf.org/licensing/licenses
+>>> END OF LICENSE >>>
+*************************************************************************/
+#ifndef _alglibmisc_pkg_h
+#define _alglibmisc_pkg_h
+#include "ap.h"
+#include "alglibinternal.h"
+
+/////////////////////////////////////////////////////////////////////////
+//
+// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES)
+//
+/////////////////////////////////////////////////////////////////////////
+namespace alglib_impl
+{
+typedef struct
+{
+    ae_int_t s1;
+    ae_int_t s2;
+    ae_int_t magicv;
+} hqrndstate;
+typedef struct
+{
+    ae_int_t n;
+    ae_int_t nx;
+    ae_int_t ny;
+    ae_int_t normtype;
+    ae_matrix xy;
+    ae_vector tags;
+    ae_vector boxmin;
+    ae_vector boxmax;
+    ae_vector nodes;
+    ae_vector splits;
+    ae_vector x;
+    ae_int_t kneeded;
+    double rneeded;
+    ae_bool selfmatch;
+    double approxf;
+    ae_int_t kcur;
+    ae_vector idx;
+    ae_vector r;
+    ae_vector buf;
+    ae_vector curboxmin;
+    ae_vector curboxmax;
+    double curdist;
+    ae_int_t debugcounter;
+} kdtree;
+
+}
+
+/////////////////////////////////////////////////////////////////////////
+//
+// THIS SECTION CONTAINS C++ INTERFACE
+//
+/////////////////////////////////////////////////////////////////////////
+namespace alglib
+{
+
+/*************************************************************************
+Portable high quality random number generator state.
+Initialized with HQRNDRandomize() or HQRNDSeed().
+
+Fields:
+    S1, S2      -   seed values
+    V           -   precomputed value
+    MagicV      -   'magic' value used to determine whether State structure
+                    was correctly initialized.
+*************************************************************************/
+class _hqrndstate_owner
+{
+public:
+    _hqrndstate_owner();
+    _hqrndstate_owner(const _hqrndstate_owner &rhs);
+    _hqrndstate_owner& operator=(const _hqrndstate_owner &rhs);
+    virtual ~_hqrndstate_owner();
+    alglib_impl::hqrndstate* c_ptr();
+    alglib_impl::hqrndstate* c_ptr() const;
+protected:
+    alglib_impl::hqrndstate *p_struct;
+};
+class hqrndstate : public _hqrndstate_owner
+{
+public:
+    hqrndstate();
+    hqrndstate(const hqrndstate &rhs);
+    hqrndstate& operator=(const hqrndstate &rhs);
+    virtual ~hqrndstate();
+
+};
+
+/*************************************************************************
+
+*************************************************************************/
+class _kdtree_owner
+{
+public:
+    _kdtree_owner();
+    _kdtree_owner(const _kdtree_owner &rhs);
+    _kdtree_owner& operator=(const _kdtree_owner &rhs);
+    virtual ~_kdtree_owner();
+    alglib_impl::kdtree* c_ptr();
+    alglib_impl::kdtree* c_ptr() const;
+protected:
+    alglib_impl::kdtree *p_struct;
+};
+class kdtree : public _kdtree_owner
+{
+public:
+    kdtree();
+    kdtree(const kdtree &rhs);
+    kdtree& operator=(const kdtree &rhs);
+    virtual ~kdtree();
+
+};
+
+/*************************************************************************
+HQRNDState  initialization  with  random  values  which come from standard
+RNG.
+
+  -- ALGLIB --
+     Copyright 02.12.2009 by Bochkanov Sergey
+*************************************************************************/
+void hqrndrandomize(hqrndstate &state);
+
+
+/*************************************************************************
+HQRNDState initialization with seed values
+
+  -- ALGLIB --
+     Copyright 02.12.2009 by Bochkanov Sergey
+*************************************************************************/
+void hqrndseed(const ae_int_t s1, const ae_int_t s2, hqrndstate &state);
+
+
+/*************************************************************************
+This function generates random real number in (0,1),
+not including interval boundaries
+
+State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
+
+  -- ALGLIB --
+     Copyright 02.12.2009 by Bochkanov Sergey
+*************************************************************************/
+double hqrnduniformr(const hqrndstate &state);
+
+
+/*************************************************************************
+This function generates random integer number in [0, N)
+
+1. State structure must be initialized with HQRNDRandomize() or HQRNDSeed()
+2. N can be any positive number except for very large numbers:
+   * close to 2^31 on 32-bit systems
+   * close to 2^62 on 64-bit systems
+   An exception will be generated if N is too large.
+
+  -- ALGLIB --
+     Copyright 02.12.2009 by Bochkanov Sergey
+*************************************************************************/
+ae_int_t hqrnduniformi(const hqrndstate &state, const ae_int_t n);
+
+
+/*************************************************************************
+Random number generator: normal numbers
+
+This function generates one random number from normal distribution.
+Its performance is equal to that of HQRNDNormal2()
+
+State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
+
+  -- ALGLIB --
+     Copyright 02.12.2009 by Bochkanov Sergey
+*************************************************************************/
+double hqrndnormal(const hqrndstate &state);
+
+
+/*************************************************************************
+Random number generator: random X and Y such that X^2+Y^2=1
+
+State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
+
+  -- ALGLIB --
+     Copyright 02.12.2009 by Bochkanov Sergey
+*************************************************************************/
+void hqrndunit2(const hqrndstate &state, double &x, double &y);
+
+
+/*************************************************************************
+Random number generator: normal numbers
+
+This function generates two independent random numbers from normal
+distribution. Its performance is equal to that of HQRNDNormal()
+
+State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
+
+  -- ALGLIB --
+     Copyright 02.12.2009 by Bochkanov Sergey
+*************************************************************************/
+void hqrndnormal2(const hqrndstate &state, double &x1, double &x2);
+
+
+/*************************************************************************
+Random number generator: exponential distribution
+
+State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
+
+  -- ALGLIB --
+     Copyright 11.08.2007 by Bochkanov Sergey
+*************************************************************************/
+double hqrndexponential(const hqrndstate &state, const double lambdav);
+
+
+/*************************************************************************
+This function generates  random number from discrete distribution given by
+finite sample X.
+
+INPUT PARAMETERS
+    State   -   high quality random number generator, must be
+                initialized with HQRNDRandomize() or HQRNDSeed().
+        X   -   finite sample
+        N   -   number of elements to use, N>=1
+
+RESULT
+    this function returns one of the X[i] for random i=0..N-1
+
+  -- ALGLIB --
+     Copyright 08.11.2011 by Bochkanov Sergey
+*************************************************************************/
+double hqrnddiscrete(const hqrndstate &state, const real_1d_array &x, const ae_int_t n);
+
+
+/*************************************************************************
+This function generates random number from continuous  distribution  given
+by finite sample X.
+
+INPUT PARAMETERS
+    State   -   high quality random number generator, must be
+                initialized with HQRNDRandomize() or HQRNDSeed().
+        X   -   finite sample, array[N] (can be larger, in this  case only
+                leading N elements are used). THIS ARRAY MUST BE SORTED BY
+                ASCENDING.
+        N   -   number of elements to use, N>=1
+
+RESULT
+    this function returns random number from continuous distribution which
+    tries to approximate X as mush as possible. min(X)<=Result<=max(X).
+
+  -- ALGLIB --
+     Copyright 08.11.2011 by Bochkanov Sergey
+*************************************************************************/
+double hqrndcontinuous(const hqrndstate &state, const real_1d_array &x, const ae_int_t n);
+
+/*************************************************************************
+This function serializes data structure to string.
+
+Important properties of s_out:
+* it contains alphanumeric characters, dots, underscores, minus signs
+* these symbols are grouped into words, which are separated by spaces
+  and Windows-style (CR+LF) newlines
+* although  serializer  uses  spaces and CR+LF as separators, you can 
+  replace any separator character by arbitrary combination of spaces,
+  tabs, Windows or Unix newlines. It allows flexible reformatting  of
+  the  string  in  case you want to include it into text or XML file. 
+  But you should not insert separators into the middle of the "words"
+  nor you should change case of letters.
+* s_out can be freely moved between 32-bit and 64-bit systems, little
+  and big endian machines, and so on. You can serialize structure  on
+  32-bit machine and unserialize it on 64-bit one (or vice versa), or
+  serialize  it  on  SPARC  and  unserialize  on  x86.  You  can also 
+  serialize  it  in  C++ version of ALGLIB and unserialize in C# one, 
+  and vice versa.
+*************************************************************************/
+void kdtreeserialize(kdtree &obj, std::string &s_out);
+
+
+/*************************************************************************
+This function unserializes data structure from string.
+*************************************************************************/
+void kdtreeunserialize(std::string &s_in, kdtree &obj);
+
+
+/*************************************************************************
+KD-tree creation
+
+This subroutine creates KD-tree from set of X-values and optional Y-values
+
+INPUT PARAMETERS
+    XY      -   dataset, array[0..N-1,0..NX+NY-1].
+                one row corresponds to one point.
+                first NX columns contain X-values, next NY (NY may be zero)
+                columns may contain associated Y-values
+    N       -   number of points, N>=0.
+    NX      -   space dimension, NX>=1.
+    NY      -   number of optional Y-values, NY>=0.
+    NormType-   norm type:
+                * 0 denotes infinity-norm
+                * 1 denotes 1-norm
+                * 2 denotes 2-norm (Euclidean norm)
+
+OUTPUT PARAMETERS
+    KDT     -   KD-tree
+
+
+NOTES
+
+1. KD-tree  creation  have O(N*logN) complexity and O(N*(2*NX+NY))  memory
+   requirements.
+2. Although KD-trees may be used with any combination of N  and  NX,  they
+   are more efficient than brute-force search only when N >> 4^NX. So they
+   are most useful in low-dimensional tasks (NX=2, NX=3). NX=1  is another
+   inefficient case, because  simple  binary  search  (without  additional
+   structures) is much more efficient in such tasks than KD-trees.
+
+  -- ALGLIB --
+     Copyright 28.02.2010 by Bochkanov Sergey
+*************************************************************************/
+void kdtreebuild(const real_2d_array &xy, const ae_int_t n, const ae_int_t nx, const ae_int_t ny, const ae_int_t normtype, kdtree &kdt);
+void kdtreebuild(const real_2d_array &xy, const ae_int_t nx, const ae_int_t ny, const ae_int_t normtype, kdtree &kdt);
+
+
+/*************************************************************************
+KD-tree creation
+
+This  subroutine  creates  KD-tree  from set of X-values, integer tags and
+optional Y-values
+
+INPUT PARAMETERS
+    XY      -   dataset, array[0..N-1,0..NX+NY-1].
+                one row corresponds to one point.
+                first NX columns contain X-values, next NY (NY may be zero)
+                columns may contain associated Y-values
+    Tags    -   tags, array[0..N-1], contains integer tags associated
+                with points.
+    N       -   number of points, N>=0
+    NX      -   space dimension, NX>=1.
+    NY      -   number of optional Y-values, NY>=0.
+    NormType-   norm type:
+                * 0 denotes infinity-norm
+                * 1 denotes 1-norm
+                * 2 denotes 2-norm (Euclidean norm)
+
+OUTPUT PARAMETERS
+    KDT     -   KD-tree
+
+NOTES
+
+1. KD-tree  creation  have O(N*logN) complexity and O(N*(2*NX+NY))  memory
+   requirements.
+2. Although KD-trees may be used with any combination of N  and  NX,  they
+   are more efficient than brute-force search only when N >> 4^NX. So they
+   are most useful in low-dimensional tasks (NX=2, NX=3). NX=1  is another
+   inefficient case, because  simple  binary  search  (without  additional
+   structures) is much more efficient in such tasks than KD-trees.
+
+  -- ALGLIB --
+     Copyright 28.02.2010 by Bochkanov Sergey
+*************************************************************************/
+void kdtreebuildtagged(const real_2d_array &xy, const integer_1d_array &tags, const ae_int_t n, const ae_int_t nx, const ae_int_t ny, const ae_int_t normtype, kdtree &kdt);
+void kdtreebuildtagged(const real_2d_array &xy, const integer_1d_array &tags, const ae_int_t nx, const ae_int_t ny, const ae_int_t normtype, kdtree &kdt);
+
+
+/*************************************************************************
+K-NN query: K nearest neighbors
+
+INPUT PARAMETERS
+    KDT         -   KD-tree
+    X           -   point, array[0..NX-1].
+    K           -   number of neighbors to return, K>=1
+    SelfMatch   -   whether self-matches are allowed:
+                    * if True, nearest neighbor may be the point itself
+                      (if it exists in original dataset)
+                    * if False, then only points with non-zero distance
+                      are returned
+                    * if not given, considered True
+
+RESULT
+    number of actual neighbors found (either K or N, if K>N).
+
+This  subroutine  performs  query  and  stores  its result in the internal
+structures of the KD-tree. You can use  following  subroutines  to  obtain
+these results:
+* KDTreeQueryResultsX() to get X-values
+* KDTreeQueryResultsXY() to get X- and Y-values
+* KDTreeQueryResultsTags() to get tag values
+* KDTreeQueryResultsDistances() to get distances
+
+  -- ALGLIB --
+     Copyright 28.02.2010 by Bochkanov Sergey
+*************************************************************************/
+ae_int_t kdtreequeryknn(const kdtree &kdt, const real_1d_array &x, const ae_int_t k, const bool selfmatch);
+ae_int_t kdtreequeryknn(const kdtree &kdt, const real_1d_array &x, const ae_int_t k);
+
+
+/*************************************************************************
+R-NN query: all points within R-sphere centered at X
+
+INPUT PARAMETERS
+    KDT         -   KD-tree
+    X           -   point, array[0..NX-1].
+    R           -   radius of sphere (in corresponding norm), R>0
+    SelfMatch   -   whether self-matches are allowed:
+                    * if True, nearest neighbor may be the point itself
+                      (if it exists in original dataset)
+                    * if False, then only points with non-zero distance
+                      are returned
+                    * if not given, considered True
+
+RESULT
+    number of neighbors found, >=0
+
+This  subroutine  performs  query  and  stores  its result in the internal
+structures of the KD-tree. You can use  following  subroutines  to  obtain
+actual results:
+* KDTreeQueryResultsX() to get X-values
+* KDTreeQueryResultsXY() to get X- and Y-values
+* KDTreeQueryResultsTags() to get tag values
+* KDTreeQueryResultsDistances() to get distances
+
+  -- ALGLIB --
+     Copyright 28.02.2010 by Bochkanov Sergey
+*************************************************************************/
+ae_int_t kdtreequeryrnn(const kdtree &kdt, const real_1d_array &x, const double r, const bool selfmatch);
+ae_int_t kdtreequeryrnn(const kdtree &kdt, const real_1d_array &x, const double r);
+
+
+/*************************************************************************
+K-NN query: approximate K nearest neighbors
+
+INPUT PARAMETERS
+    KDT         -   KD-tree
+    X           -   point, array[0..NX-1].
+    K           -   number of neighbors to return, K>=1
+    SelfMatch   -   whether self-matches are allowed:
+                    * if True, nearest neighbor may be the point itself
+                      (if it exists in original dataset)
+                    * if False, then only points with non-zero distance
+                      are returned
+                    * if not given, considered True
+    Eps         -   approximation factor, Eps>=0. eps-approximate  nearest
+                    neighbor  is  a  neighbor  whose distance from X is at
+                    most (1+eps) times distance of true nearest neighbor.
+
+RESULT
+    number of actual neighbors found (either K or N, if K>N).
+
+NOTES
+    significant performance gain may be achieved only when Eps  is  is  on
+    the order of magnitude of 1 or larger.
+
+This  subroutine  performs  query  and  stores  its result in the internal
+structures of the KD-tree. You can use  following  subroutines  to  obtain
+these results:
+* KDTreeQueryResultsX() to get X-values
+* KDTreeQueryResultsXY() to get X- and Y-values
+* KDTreeQueryResultsTags() to get tag values
+* KDTreeQueryResultsDistances() to get distances
+
+  -- ALGLIB --
+     Copyright 28.02.2010 by Bochkanov Sergey
+*************************************************************************/
+ae_int_t kdtreequeryaknn(const kdtree &kdt, const real_1d_array &x, const ae_int_t k, const bool selfmatch, const double eps);
+ae_int_t kdtreequeryaknn(const kdtree &kdt, const real_1d_array &x, const ae_int_t k, const double eps);
+
+
+/*************************************************************************
+X-values from last query
+
+INPUT PARAMETERS
+    KDT     -   KD-tree
+    X       -   possibly pre-allocated buffer. If X is too small to store
+                result, it is resized. If size(X) is enough to store
+                result, it is left unchanged.
+
+OUTPUT PARAMETERS
+    X       -   rows are filled with X-values
+
+NOTES
+1. points are ordered by distance from the query point (first = closest)
+2. if  XY is larger than required to store result, only leading part  will
+   be overwritten; trailing part will be left unchanged. So  if  on  input
+   XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
+   XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
+   you want function  to  resize  array  according  to  result  size,  use
+   function with same name and suffix 'I'.
+
+SEE ALSO
+* KDTreeQueryResultsXY()            X- and Y-values
+* KDTreeQueryResultsTags()          tag values
+* KDTreeQueryResultsDistances()     distances
+
+  -- ALGLIB --
+     Copyright 28.02.2010 by Bochkanov Sergey
+*************************************************************************/
+void kdtreequeryresultsx(const kdtree &kdt, real_2d_array &x);
+
+
+/*************************************************************************
+X- and Y-values from last query
+
+INPUT PARAMETERS
+    KDT     -   KD-tree
+    XY      -   possibly pre-allocated buffer. If XY is too small to store
+                result, it is resized. If size(XY) is enough to store
+                result, it is left unchanged.
+
+OUTPUT PARAMETERS
+    XY      -   rows are filled with points: first NX columns with
+                X-values, next NY columns - with Y-values.
+
+NOTES
+1. points are ordered by distance from the query point (first = closest)
+2. if  XY is larger than required to store result, only leading part  will
+   be overwritten; trailing part will be left unchanged. So  if  on  input
+   XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
+   XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
+   you want function  to  resize  array  according  to  result  size,  use
+   function with same name and suffix 'I'.
+
+SEE ALSO
+* KDTreeQueryResultsX()             X-values
+* KDTreeQueryResultsTags()          tag values
+* KDTreeQueryResultsDistances()     distances
+
+  -- ALGLIB --
+     Copyright 28.02.2010 by Bochkanov Sergey
+*************************************************************************/
+void kdtreequeryresultsxy(const kdtree &kdt, real_2d_array &xy);
+
+
+/*************************************************************************
+Tags from last query
+
+INPUT PARAMETERS
+    KDT     -   KD-tree
+    Tags    -   possibly pre-allocated buffer. If X is too small to store
+                result, it is resized. If size(X) is enough to store
+                result, it is left unchanged.
+
+OUTPUT PARAMETERS
+    Tags    -   filled with tags associated with points,
+                or, when no tags were supplied, with zeros
+
+NOTES
+1. points are ordered by distance from the query point (first = closest)
+2. if  XY is larger than required to store result, only leading part  will
+   be overwritten; trailing part will be left unchanged. So  if  on  input
+   XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
+   XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
+   you want function  to  resize  array  according  to  result  size,  use
+   function with same name and suffix 'I'.
+
+SEE ALSO
+* KDTreeQueryResultsX()             X-values
+* KDTreeQueryResultsXY()            X- and Y-values
+* KDTreeQueryResultsDistances()     distances
+
+  -- ALGLIB --
+     Copyright 28.02.2010 by Bochkanov Sergey
+*************************************************************************/
+void kdtreequeryresultstags(const kdtree &kdt, integer_1d_array &tags);
+
+
+/*************************************************************************
+Distances from last query
+
+INPUT PARAMETERS
+    KDT     -   KD-tree
+    R       -   possibly pre-allocated buffer. If X is too small to store
+                result, it is resized. If size(X) is enough to store
+                result, it is left unchanged.
+
+OUTPUT PARAMETERS
+    R       -   filled with distances (in corresponding norm)
+
+NOTES
+1. points are ordered by distance from the query point (first = closest)
+2. if  XY is larger than required to store result, only leading part  will
+   be overwritten; trailing part will be left unchanged. So  if  on  input
+   XY = [[A,B],[C,D]], and result is [1,2],  then  on  exit  we  will  get
+   XY = [[1,2],[C,D]]. This is done purposely to increase performance;  if
+   you want function  to  resize  array  according  to  result  size,  use
+   function with same name and suffix 'I'.
+
+SEE ALSO
+* KDTreeQueryResultsX()             X-values
+* KDTreeQueryResultsXY()            X- and Y-values
+* KDTreeQueryResultsTags()          tag values
+
+  -- ALGLIB --
+     Copyright 28.02.2010 by Bochkanov Sergey
+*************************************************************************/
+void kdtreequeryresultsdistances(const kdtree &kdt, real_1d_array &r);
+
+
+/*************************************************************************
+X-values from last query; 'interactive' variant for languages like  Python
+which   support    constructs   like  "X = KDTreeQueryResultsXI(KDT)"  and
+interactive mode of interpreter.
+
+This function allocates new array on each call,  so  it  is  significantly
+slower than its 'non-interactive' counterpart, but it is  more  convenient
+when you call it from command line.
+
+  -- ALGLIB --
+     Copyright 28.02.2010 by Bochkanov Sergey
+*************************************************************************/
+void kdtreequeryresultsxi(const kdtree &kdt, real_2d_array &x);
+
+
+/*************************************************************************
+XY-values from last query; 'interactive' variant for languages like Python
+which   support    constructs   like "XY = KDTreeQueryResultsXYI(KDT)" and
+interactive mode of interpreter.
+
+This function allocates new array on each call,  so  it  is  significantly
+slower than its 'non-interactive' counterpart, but it is  more  convenient
+when you call it from command line.
+
+  -- ALGLIB --
+     Copyright 28.02.2010 by Bochkanov Sergey
+*************************************************************************/
+void kdtreequeryresultsxyi(const kdtree &kdt, real_2d_array &xy);
+
+
+/*************************************************************************
+Tags  from  last  query;  'interactive' variant for languages like  Python
+which  support  constructs  like "Tags = KDTreeQueryResultsTagsI(KDT)" and
+interactive mode of interpreter.
+
+This function allocates new array on each call,  so  it  is  significantly
+slower than its 'non-interactive' counterpart, but it is  more  convenient
+when you call it from command line.
+
+  -- ALGLIB --
+     Copyright 28.02.2010 by Bochkanov Sergey
+*************************************************************************/
+void kdtreequeryresultstagsi(const kdtree &kdt, integer_1d_array &tags);
+
+
+/*************************************************************************
+Distances from last query; 'interactive' variant for languages like Python
+which  support  constructs   like  "R = KDTreeQueryResultsDistancesI(KDT)"
+and interactive mode of interpreter.
+
+This function allocates new array on each call,  so  it  is  significantly
+slower than its 'non-interactive' counterpart, but it is  more  convenient
+when you call it from command line.
+
+  -- ALGLIB --
+     Copyright 28.02.2010 by Bochkanov Sergey
+*************************************************************************/
+void kdtreequeryresultsdistancesi(const kdtree &kdt, real_1d_array &r);
+}
+
+/////////////////////////////////////////////////////////////////////////
+//
+// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS)
+//
+/////////////////////////////////////////////////////////////////////////
+namespace alglib_impl
+{
+void hqrndrandomize(hqrndstate* state, ae_state *_state);
+void hqrndseed(ae_int_t s1,
+     ae_int_t s2,
+     hqrndstate* state,
+     ae_state *_state);
+double hqrnduniformr(hqrndstate* state, ae_state *_state);
+ae_int_t hqrnduniformi(hqrndstate* state, ae_int_t n, ae_state *_state);
+double hqrndnormal(hqrndstate* state, ae_state *_state);
+void hqrndunit2(hqrndstate* state, double* x, double* y, ae_state *_state);
+void hqrndnormal2(hqrndstate* state,
+     double* x1,
+     double* x2,
+     ae_state *_state);
+double hqrndexponential(hqrndstate* state,
+     double lambdav,
+     ae_state *_state);
+double hqrnddiscrete(hqrndstate* state,
+     /* Real    */ ae_vector* x,
+     ae_int_t n,
+     ae_state *_state);
+double hqrndcontinuous(hqrndstate* state,
+     /* Real    */ ae_vector* x,
+     ae_int_t n,
+     ae_state *_state);
+ae_bool _hqrndstate_init(void* _p, ae_state *_state, ae_bool make_automatic);
+ae_bool _hqrndstate_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
+void _hqrndstate_clear(void* _p);
+void _hqrndstate_destroy(void* _p);
+void kdtreebuild(/* Real    */ ae_matrix* xy,
+     ae_int_t n,
+     ae_int_t nx,
+     ae_int_t ny,
+     ae_int_t normtype,
+     kdtree* kdt,
+     ae_state *_state);
+void kdtreebuildtagged(/* Real    */ ae_matrix* xy,
+     /* Integer */ ae_vector* tags,
+     ae_int_t n,
+     ae_int_t nx,
+     ae_int_t ny,
+     ae_int_t normtype,
+     kdtree* kdt,
+     ae_state *_state);
+ae_int_t kdtreequeryknn(kdtree* kdt,
+     /* Real    */ ae_vector* x,
+     ae_int_t k,
+     ae_bool selfmatch,
+     ae_state *_state);
+ae_int_t kdtreequeryrnn(kdtree* kdt,
+     /* Real    */ ae_vector* x,
+     double r,
+     ae_bool selfmatch,
+     ae_state *_state);
+ae_int_t kdtreequeryaknn(kdtree* kdt,
+     /* Real    */ ae_vector* x,
+     ae_int_t k,
+     ae_bool selfmatch,
+     double eps,
+     ae_state *_state);
+void kdtreequeryresultsx(kdtree* kdt,
+     /* Real    */ ae_matrix* x,
+     ae_state *_state);
+void kdtreequeryresultsxy(kdtree* kdt,
+     /* Real    */ ae_matrix* xy,
+     ae_state *_state);
+void kdtreequeryresultstags(kdtree* kdt,
+     /* Integer */ ae_vector* tags,
+     ae_state *_state);
+void kdtreequeryresultsdistances(kdtree* kdt,
+     /* Real    */ ae_vector* r,
+     ae_state *_state);
+void kdtreequeryresultsxi(kdtree* kdt,
+     /* Real    */ ae_matrix* x,
+     ae_state *_state);
+void kdtreequeryresultsxyi(kdtree* kdt,
+     /* Real    */ ae_matrix* xy,
+     ae_state *_state);
+void kdtreequeryresultstagsi(kdtree* kdt,
+     /* Integer */ ae_vector* tags,
+     ae_state *_state);
+void kdtreequeryresultsdistancesi(kdtree* kdt,
+     /* Real    */ ae_vector* r,
+     ae_state *_state);
+void kdtreealloc(ae_serializer* s, kdtree* tree, ae_state *_state);
+void kdtreeserialize(ae_serializer* s, kdtree* tree, ae_state *_state);
+void kdtreeunserialize(ae_serializer* s, kdtree* tree, ae_state *_state);
+ae_bool _kdtree_init(void* _p, ae_state *_state, ae_bool make_automatic);
+ae_bool _kdtree_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
+void _kdtree_clear(void* _p);
+void _kdtree_destroy(void* _p);
+
+}
+#endif
+

src/inc/alglib/diffequations.h

@@ -0,0 +1,267 @@
+/*************************************************************************
+Copyright (c) Sergey Bochkanov (ALGLIB project).
+
+>>> SOURCE LICENSE >>>
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation (www.fsf.org); either version 2 of the
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+A copy of the GNU General Public License is available at
+http://www.fsf.org/licensing/licenses
+>>> END OF LICENSE >>>
+*************************************************************************/
+#ifndef _diffequations_pkg_h
+#define _diffequations_pkg_h
+#include "ap.h"
+#include "alglibinternal.h"
+
+/////////////////////////////////////////////////////////////////////////
+//
+// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES)
+//
+/////////////////////////////////////////////////////////////////////////
+namespace alglib_impl
+{
+typedef struct
+{
+    ae_int_t n;
+    ae_int_t m;
+    double xscale;
+    double h;
+    double eps;
+    ae_bool fraceps;
+    ae_vector yc;
+    ae_vector escale;
+    ae_vector xg;
+    ae_int_t solvertype;
+    ae_bool needdy;
+    double x;
+    ae_vector y;
+    ae_vector dy;
+    ae_matrix ytbl;
+    ae_int_t repterminationtype;
+    ae_int_t repnfev;
+    ae_vector yn;
+    ae_vector yns;
+    ae_vector rka;
+    ae_vector rkc;
+    ae_vector rkcs;
+    ae_matrix rkb;
+    ae_matrix rkk;
+    rcommstate rstate;
+} odesolverstate;
+typedef struct
+{
+    ae_int_t nfev;
+    ae_int_t terminationtype;
+} odesolverreport;
+
+}
+
+/////////////////////////////////////////////////////////////////////////
+//
+// THIS SECTION CONTAINS C++ INTERFACE
+//
+/////////////////////////////////////////////////////////////////////////
+namespace alglib
+{
+
+/*************************************************************************
+
+*************************************************************************/
+class _odesolverstate_owner
+{
+public:
+    _odesolverstate_owner();
+    _odesolverstate_owner(const _odesolverstate_owner &rhs);
+    _odesolverstate_owner& operator=(const _odesolverstate_owner &rhs);
+    virtual ~_odesolverstate_owner();
+    alglib_impl::odesolverstate* c_ptr();
+    alglib_impl::odesolverstate* c_ptr() const;
+protected:
+    alglib_impl::odesolverstate *p_struct;
+};
+class odesolverstate : public _odesolverstate_owner
+{
+public:
+    odesolverstate();
+    odesolverstate(const odesolverstate &rhs);
+    odesolverstate& operator=(const odesolverstate &rhs);
+    virtual ~odesolverstate();
+    ae_bool &needdy;
+    real_1d_array y;
+    real_1d_array dy;
+    double &x;
+
+};
+
+
+/*************************************************************************
+
+*************************************************************************/
+class _odesolverreport_owner
+{
+public:
+    _odesolverreport_owner();
+    _odesolverreport_owner(const _odesolverreport_owner &rhs);
+    _odesolverreport_owner& operator=(const _odesolverreport_owner &rhs);
+    virtual ~_odesolverreport_owner();
+    alglib_impl::odesolverreport* c_ptr();
+    alglib_impl::odesolverreport* c_ptr() const;
+protected:
+    alglib_impl::odesolverreport *p_struct;
+};
+class odesolverreport : public _odesolverreport_owner
+{
+public:
+    odesolverreport();
+    odesolverreport(const odesolverreport &rhs);
+    odesolverreport& operator=(const odesolverreport &rhs);
+    virtual ~odesolverreport();
+    ae_int_t &nfev;
+    ae_int_t &terminationtype;
+
+};
+
+/*************************************************************************
+Cash-Karp adaptive ODE solver.
+
+This subroutine solves ODE  Y'=f(Y,x)  with  initial  conditions  Y(xs)=Ys
+(here Y may be single variable or vector of N variables).
+
+INPUT PARAMETERS:
+    Y       -   initial conditions, array[0..N-1].
+                contains values of Y[] at X[0]
+    N       -   system size
+    X       -   points at which Y should be tabulated, array[0..M-1]
+                integrations starts at X[0], ends at X[M-1],  intermediate
+                values at X[i] are returned too.
+                SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!!!!
+    M       -   number of intermediate points + first point + last point:
+                * M>2 means that you need both Y(X[M-1]) and M-2 values at
+                  intermediate points
+                * M=2 means that you want just to integrate from  X[0]  to
+                  X[1] and don't interested in intermediate values.
+                * M=1 means that you don't want to integrate :)
+                  it is degenerate case, but it will be handled correctly.
+                * M<1 means error
+    Eps     -   tolerance (absolute/relative error on each  step  will  be
+                less than Eps). When passing:
+                * Eps>0, it means desired ABSOLUTE error
+                * Eps<0, it means desired RELATIVE error.  Relative errors
+                  are calculated with respect to maximum values of  Y seen
+                  so far. Be careful to use this criterion  when  starting
+                  from Y[] that are close to zero.
+    H       -   initial  step  lenth,  it  will  be adjusted automatically
+                after the first  step.  If  H=0,  step  will  be  selected
+                automatically  (usualy  it  will  be  equal  to  0.001  of
+                min(x[i]-x[j])).
+
+OUTPUT PARAMETERS
+    State   -   structure which stores algorithm state between  subsequent
+                calls of OdeSolverIteration. Used for reverse communication.
+                This structure should be passed  to the OdeSolverIteration
+                subroutine.
+
+SEE ALSO
+    AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults.
+
+
+  -- ALGLIB --
+     Copyright 01.09.2009 by Bochkanov Sergey
+*************************************************************************/
+void odesolverrkck(const real_1d_array &y, const ae_int_t n, const real_1d_array &x, const ae_int_t m, const double eps, const double h, odesolverstate &state);
+void odesolverrkck(const real_1d_array &y, const real_1d_array &x, const double eps, const double h, odesolverstate &state);
+
+
+/*************************************************************************
+This function provides reverse communication interface
+Reverse communication interface is not documented or recommended to use.
+See below for functions which provide better documented API
+*************************************************************************/
+bool odesolveriteration(const odesolverstate &state);
+
+
+/*************************************************************************
+This function is used to launcn iterations of ODE solver
+
+It accepts following parameters:
+    diff    -   callback which calculates dy/dx for given y and x
+    ptr     -   optional pointer which is passed to diff; can be NULL
+
+
+  -- ALGLIB --
+     Copyright 01.09.2009 by Bochkanov Sergey
+
+*************************************************************************/
+void odesolversolve(odesolverstate &state,
+    void (*diff)(const real_1d_array &y, double x, real_1d_array &dy, void *ptr),
+    void *ptr = NULL);
+
+
+/*************************************************************************
+ODE solver results
+
+Called after OdeSolverIteration returned False.
+
+INPUT PARAMETERS:
+    State   -   algorithm state (used by OdeSolverIteration).
+
+OUTPUT PARAMETERS:
+    M       -   number of tabulated values, M>=1
+    XTbl    -   array[0..M-1], values of X
+    YTbl    -   array[0..M-1,0..N-1], values of Y in X[i]
+    Rep     -   solver report:
+                * Rep.TerminationType completetion code:
+                    * -2    X is not ordered  by  ascending/descending  or
+                            there are non-distinct X[],  i.e.  X[i]=X[i+1]
+                    * -1    incorrect parameters were specified
+                    *  1    task has been solved
+                * Rep.NFEV contains number of function calculations
+
+  -- ALGLIB --
+     Copyright 01.09.2009 by Bochkanov Sergey
+*************************************************************************/
+void odesolverresults(const odesolverstate &state, ae_int_t &m, real_1d_array &xtbl, real_2d_array &ytbl, odesolverreport &rep);
+}
+
+/////////////////////////////////////////////////////////////////////////
+//
+// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS)
+//
+/////////////////////////////////////////////////////////////////////////
+namespace alglib_impl
+{
+void odesolverrkck(/* Real    */ ae_vector* y,
+     ae_int_t n,
+     /* Real    */ ae_vector* x,
+     ae_int_t m,
+     double eps,
+     double h,
+     odesolverstate* state,
+     ae_state *_state);
+ae_bool odesolveriteration(odesolverstate* state, ae_state *_state);
+void odesolverresults(odesolverstate* state,
+     ae_int_t* m,
+     /* Real    */ ae_vector* xtbl,
+     /* Real    */ ae_matrix* ytbl,
+     odesolverreport* rep,
+     ae_state *_state);
+ae_bool _odesolverstate_init(void* _p, ae_state *_state, ae_bool make_automatic);
+ae_bool _odesolverstate_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
+void _odesolverstate_clear(void* _p);
+void _odesolverstate_destroy(void* _p);
+ae_bool _odesolverreport_init(void* _p, ae_state *_state, ae_bool make_automatic);
+ae_bool _odesolverreport_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
+void _odesolverreport_clear(void* _p);
+void _odesolverreport_destroy(void* _p);
+
+}
+#endif
+

src/inc/alglib/fasttransforms.h

@@ -0,0 +1,691 @@
+/*************************************************************************
+Copyright (c) Sergey Bochkanov (ALGLIB project).
+
+>>> SOURCE LICENSE >>>
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation (www.fsf.org); either version 2 of the
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+A copy of the GNU General Public License is available at
+http://www.fsf.org/licensing/licenses
+>>> END OF LICENSE >>>
+*************************************************************************/
+#ifndef _fasttransforms_pkg_h
+#define _fasttransforms_pkg_h
+#include "ap.h"
+#include "alglibinternal.h"
+
+/////////////////////////////////////////////////////////////////////////
+//
+// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES)
+//
+/////////////////////////////////////////////////////////////////////////
+namespace alglib_impl
+{
+
+}
+
+/////////////////////////////////////////////////////////////////////////
+//
+// THIS SECTION CONTAINS C++ INTERFACE
+//
+/////////////////////////////////////////////////////////////////////////
+namespace alglib
+{
+
+
+/*************************************************************************
+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
+*************************************************************************/
+void fftc1d(complex_1d_array &a, const ae_int_t n);
+void fftc1d(complex_1d_array &a);
+
+
+/*************************************************************************
+1-dimensional complex inverse FFT.
+
+Array size N may be arbitrary number (composite or prime).  Algorithm  has
+O(N*logN) complexity for any N (composite or prime).
+
+See FFTC1D() description for more information about algorithm performance.
+
+INPUT PARAMETERS
+    A   -   array[0..N-1] - complex array to be transformed
+    N   -   problem size
+
+OUTPUT PARAMETERS
+    A   -   inverse DFT of a input array, array[0..N-1]
+            A_out[j] = SUM(A_in[k]/N*exp(+2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
+
+
+  -- ALGLIB --
+     Copyright 29.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void fftc1dinv(complex_1d_array &a, const ae_int_t n);
+void fftc1dinv(complex_1d_array &a);
+
+
+/*************************************************************************
+1-dimensional real FFT.
+
+Algorithm has O(N*logN) complexity for any N (composite or prime).
+
+INPUT PARAMETERS
+    A   -   array[0..N-1] - real function to be transformed
+    N   -   problem size
+
+OUTPUT PARAMETERS
+    F   -   DFT of a input array, array[0..N-1]
+            F[j] = SUM(A[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
+
+NOTE:
+    F[] satisfies symmetry property F[k] = conj(F[N-k]),  so just one half
+of  array  is  usually needed. But for convinience subroutine returns full
+complex array (with frequencies above N/2), so its result may be  used  by
+other FFT-related subroutines.
+
+
+  -- ALGLIB --
+     Copyright 01.06.2009 by Bochkanov Sergey
+*************************************************************************/
+void fftr1d(const real_1d_array &a, const ae_int_t n, complex_1d_array &f);
+void fftr1d(const real_1d_array &a, complex_1d_array &f);
+
+
+/*************************************************************************
+1-dimensional real inverse FFT.
+
+Algorithm has O(N*logN) complexity for any N (composite or prime).
+
+INPUT PARAMETERS
+    F   -   array[0..floor(N/2)] - frequencies from forward real FFT
+    N   -   problem size
+
+OUTPUT PARAMETERS
+    A   -   inverse DFT of a input array, array[0..N-1]
+
+NOTE:
+    F[] should satisfy symmetry property F[k] = conj(F[N-k]), so just  one
+half of frequencies array is needed - elements from 0 to floor(N/2).  F[0]
+is ALWAYS real. If N is even F[floor(N/2)] is real too. If N is odd,  then
+F[floor(N/2)] has no special properties.
+
+Relying on properties noted above, FFTR1DInv subroutine uses only elements
+from 0th to floor(N/2)-th. It ignores imaginary part of F[0],  and in case
+N is even it ignores imaginary part of F[floor(N/2)] too.
+
+When you call this function using full arguments list - "FFTR1DInv(F,N,A)"
+- you can pass either either frequencies array with N elements or  reduced
+array with roughly N/2 elements - subroutine will  successfully  transform
+both.
+
+If you call this function using reduced arguments list -  "FFTR1DInv(F,A)"
+- you must pass FULL array with N elements (although higher  N/2 are still
+not used) because array size is used to automatically determine FFT length
+
+
+  -- ALGLIB --
+     Copyright 01.06.2009 by Bochkanov Sergey
+*************************************************************************/
+void fftr1dinv(const complex_1d_array &f, const ae_int_t n, real_1d_array &a);
+void fftr1dinv(const complex_1d_array &f, real_1d_array &a);
+
+/*************************************************************************
+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
+*************************************************************************/
+void convc1d(const complex_1d_array &a, const ae_int_t m, const complex_1d_array &b, const ae_int_t n, complex_1d_array &r);
+
+
+/*************************************************************************
+1-dimensional complex non-circular deconvolution (inverse of ConvC1D()).
+
+Algorithm has M*log(M)) complexity for any M (composite or prime).
+
+INPUT PARAMETERS
+    A   -   array[0..M-1] - convolved signal, A = conv(R, B)
+    M   -   convolved signal length
+    B   -   array[0..N-1] - response
+    N   -   response length, N<=M
+
+OUTPUT PARAMETERS
+    R   -   deconvolved signal. array[0..M-N].
+
+NOTE:
+    deconvolution is unstable process and may result in division  by  zero
+(if your response function is degenerate, i.e. has zero Fourier coefficient).
+
+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
+*************************************************************************/
+void convc1dinv(const complex_1d_array &a, const ae_int_t m, const complex_1d_array &b, const ae_int_t n, complex_1d_array &r);
+
+
+/*************************************************************************
+1-dimensional circular complex convolution.
+
+For given S/R returns conv(S,R) (circular). Algorithm has linearithmic
+complexity for any M/N.
+
+IMPORTANT:  normal convolution is commutative,  i.e.   it  is symmetric  -
+conv(A,B)=conv(B,A).  Cyclic convolution IS NOT.  One function - S - is  a
+signal,  periodic function, and another - R - is a response,  non-periodic
+function with limited length.
+
+INPUT PARAMETERS
+    S   -   array[0..M-1] - complex periodic signal
+    M   -   problem size
+    B   -   array[0..N-1] - complex non-periodic response
+    N   -   problem size
+
+OUTPUT PARAMETERS
+    R   -   convolution: A*B. array[0..M-1].
+
+NOTE:
+    It is assumed that B is zero at T<0. If  it  has  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
+*************************************************************************/
+void convc1dcircular(const complex_1d_array &s, const ae_int_t m, const complex_1d_array &r, const ae_int_t n, complex_1d_array &c);
+
+
+/*************************************************************************
+1-dimensional circular complex deconvolution (inverse of ConvC1DCircular()).
+
+Algorithm has M*log(M)) complexity for any M (composite or prime).
+
+INPUT PARAMETERS
+    A   -   array[0..M-1] - convolved periodic signal, A = conv(R, B)
+    M   -   convolved signal length
+    B   -   array[0..N-1] - non-periodic response
+    N   -   response length
+
+OUTPUT PARAMETERS
+    R   -   deconvolved signal. array[0..M-1].
+
+NOTE:
+    deconvolution is unstable process and may result in division  by  zero
+(if your response function is degenerate, i.e. has zero Fourier coefficient).
+
+NOTE:
+    It is assumed that B is zero at T<0. If  it  has  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
+*************************************************************************/
+void convc1dcircularinv(const complex_1d_array &a, const ae_int_t m, const complex_1d_array &b, const ae_int_t n, complex_1d_array &r);
+
+
+/*************************************************************************
+1-dimensional real convolution.
+
+Analogous to ConvC1D(), see ConvC1D() comments for more details.
+
+INPUT PARAMETERS
+    A   -   array[0..M-1] - real function to be transformed
+    M   -   problem size
+    B   -   array[0..N-1] - real 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
+*************************************************************************/
+void convr1d(const real_1d_array &a, const ae_int_t m, const real_1d_array &b, const ae_int_t n, real_1d_array &r);
+
+
+/*************************************************************************
+1-dimensional real deconvolution (inverse of ConvC1D()).
+
+Algorithm has M*log(M)) complexity for any M (composite or prime).
+
+INPUT PARAMETERS
+    A   -   array[0..M-1] - convolved signal, A = conv(R, B)
+    M   -   convolved signal length
+    B   -   array[0..N-1] - response
+    N   -   response length, N<=M
+
+OUTPUT PARAMETERS
+    R   -   deconvolved signal. array[0..M-N].
+
+NOTE:
+    deconvolution is unstable process and may result in division  by  zero
+(if your response function is degenerate, i.e. has zero Fourier coefficient).
+
+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
+*************************************************************************/
+void convr1dinv(const real_1d_array &a, const ae_int_t m, const real_1d_array &b, const ae_int_t n, real_1d_array &r);
+
+
+/*************************************************************************
+1-dimensional circular real convolution.
+
+Analogous to ConvC1DCircular(), see ConvC1DCircular() comments for more details.
+
+INPUT PARAMETERS
+    S   -   array[0..M-1] - real signal
+    M   -   problem size
+    B   -   array[0..N-1] - real response
+    N   -   problem size
+
+OUTPUT PARAMETERS
+    R   -   convolution: A*B. array[0..M-1].
+
+NOTE:
+    It is assumed that B is zero at T<0. If  it  has  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
+*************************************************************************/
+void convr1dcircular(const real_1d_array &s, const ae_int_t m, const real_1d_array &r, const ae_int_t n, real_1d_array &c);
+
+
+/*************************************************************************
+1-dimensional complex deconvolution (inverse of ConvC1D()).
+
+Algorithm has M*log(M)) complexity for any M (composite or prime).
+
+INPUT PARAMETERS
+    A   -   array[0..M-1] - convolved signal, A = conv(R, B)
+    M   -   convolved signal length
+    B   -   array[0..N-1] - response
+    N   -   response length
+
+OUTPUT PARAMETERS
+    R   -   deconvolved signal. array[0..M-N].
+
+NOTE:
+    deconvolution is unstable process and may result in division  by  zero
+(if your response function is degenerate, i.e. has zero Fourier coefficient).
+
+NOTE:
+    It is assumed that B is zero at T<0. If  it  has  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
+*************************************************************************/
+void convr1dcircularinv(const real_1d_array &a, const ae_int_t m, const real_1d_array &b, const ae_int_t n, real_1d_array &r);
+
+/*************************************************************************
+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
+*************************************************************************/
+void corrc1d(const complex_1d_array &signal, const ae_int_t n, const complex_1d_array &pattern, const ae_int_t m, complex_1d_array &r);
+
+
+/*************************************************************************
+1-dimensional circular complex 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:   CorrC1DCircular(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,
+                periodic signal containing pattern
+    N       -   problem size
+    Pattern -   array[0..M-1] - complex 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 corrc1dcircular(const complex_1d_array &signal, const ae_int_t m, const complex_1d_array &pattern, const ae_int_t n, complex_1d_array &c);
+
+
+/*************************************************************************
+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 real_1d_array &signal, const ae_int_t n, const real_1d_array &pattern, const ae_int_t m, real_1d_array &r);
+
+
+/*************************************************************************
+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 real_1d_array &signal, const ae_int_t m, const real_1d_array &pattern, const ae_int_t n, real_1d_array &c);
+
+/*************************************************************************
+1-dimensional Fast Hartley Transform.
+
+Algorithm has O(N*logN) complexity for any N (composite or prime).
+
+INPUT PARAMETERS
+    A   -   array[0..N-1] - real function to be transformed
+    N   -   problem size
+
+OUTPUT PARAMETERS
+    A   -   FHT of a input array, array[0..N-1],
+            A_out[k] = sum(A_in[j]*(cos(2*pi*j*k/N)+sin(2*pi*j*k/N)), j=0..N-1)
+
+
+  -- ALGLIB --
+     Copyright 04.06.2009 by Bochkanov Sergey
+*************************************************************************/
+void fhtr1d(real_1d_array &a, const ae_int_t n);
+
+
+/*************************************************************************
+1-dimensional inverse FHT.
+
+Algorithm has O(N*logN) complexity for any N (composite or prime).
+
+INPUT PARAMETERS
+    A   -   array[0..N-1] - complex array to be transformed
+    N   -   problem size
+
+OUTPUT PARAMETERS
+    A   -   inverse FHT of a input array, array[0..N-1]
+
+
+  -- ALGLIB --
+     Copyright 29.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void fhtr1dinv(real_1d_array &a, const ae_int_t n);
+}
+
+/////////////////////////////////////////////////////////////////////////
+//
+// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS)
+//
+/////////////////////////////////////////////////////////////////////////
+namespace alglib_impl
+{
+void fftc1d(/* Complex */ ae_vector* a, ae_int_t n, ae_state *_state);
+void fftc1dinv(/* Complex */ ae_vector* a, ae_int_t n, ae_state *_state);
+void fftr1d(/* Real    */ ae_vector* a,
+     ae_int_t n,
+     /* Complex */ ae_vector* f,
+     ae_state *_state);
+void fftr1dinv(/* Complex */ ae_vector* f,
+     ae_int_t n,
+     /* Real    */ ae_vector* a,
+     ae_state *_state);
+void fftr1dinternaleven(/* Real    */ ae_vector* a,
+     ae_int_t n,
+     /* Real    */ ae_vector* buf,
+     fasttransformplan* plan,
+     ae_state *_state);
+void fftr1dinvinternaleven(/* Real    */ ae_vector* a,
+     ae_int_t n,
+     /* Real    */ ae_vector* buf,
+     fasttransformplan* plan,
+     ae_state *_state);
+void convc1d(/* Complex */ ae_vector* a,
+     ae_int_t m,
+     /* Complex */ ae_vector* b,
+     ae_int_t n,
+     /* Complex */ ae_vector* r,
+     ae_state *_state);
+void convc1dinv(/* Complex */ ae_vector* a,
+     ae_int_t m,
+     /* Complex */ ae_vector* b,
+     ae_int_t n,
+     /* Complex */ ae_vector* r,
+     ae_state *_state);
+void convc1dcircular(/* Complex */ ae_vector* s,
+     ae_int_t m,
+     /* Complex */ ae_vector* r,
+     ae_int_t n,
+     /* Complex */ ae_vector* c,
+     ae_state *_state);
+void convc1dcircularinv(/* Complex */ ae_vector* a,
+     ae_int_t m,
+     /* Complex */ ae_vector* b,
+     ae_int_t n,
+     /* Complex */ ae_vector* r,
+     ae_state *_state);
+void convr1d(/* Real    */ ae_vector* a,
+     ae_int_t m,
+     /* Real    */ ae_vector* b,
+     ae_int_t n,
+     /* Real    */ ae_vector* r,
+     ae_state *_state);
+void convr1dinv(/* Real    */ ae_vector* a,
+     ae_int_t m,
+     /* Real    */ ae_vector* b,
+     ae_int_t n,
+     /* Real    */ ae_vector* r,
+     ae_state *_state);
+void convr1dcircular(/* Real    */ ae_vector* s,
+     ae_int_t m,
+     /* Real    */ ae_vector* r,
+     ae_int_t n,
+     /* Real    */ ae_vector* c,
+     ae_state *_state);
+void convr1dcircularinv(/* Real    */ ae_vector* a,
+     ae_int_t m,
+     /* Real    */ ae_vector* b,
+     ae_int_t n,
+     /* Real    */ ae_vector* r,
+     ae_state *_state);
+void convc1dx(/* Complex */ ae_vector* a,
+     ae_int_t m,
+     /* Complex */ ae_vector* b,
+     ae_int_t n,
+     ae_bool circular,
+     ae_int_t alg,
+     ae_int_t q,
+     /* Complex */ ae_vector* r,
+     ae_state *_state);
+void convr1dx(/* Real    */ ae_vector* a,
+     ae_int_t m,
+     /* Real    */ ae_vector* b,
+     ae_int_t n,
+     ae_bool circular,
+     ae_int_t alg,
+     ae_int_t q,
+     /* Real    */ ae_vector* r,
+     ae_state *_state);
+void corrc1d(/* Complex */ ae_vector* signal,
+     ae_int_t n,
+     /* Complex */ ae_vector* pattern,
+     ae_int_t m,
+     /* Complex */ ae_vector* r,
+     ae_state *_state);
+void corrc1dcircular(/* Complex */ ae_vector* signal,
+     ae_int_t m,
+     /* Complex */ ae_vector* pattern,
+     ae_int_t n,
+     /* Complex */ ae_vector* c,
+     ae_state *_state);
+void corrr1d(/* Real    */ ae_vector* signal,
+     ae_int_t n,
+     /* Real    */ ae_vector* pattern,
+     ae_int_t m,
+     /* Real    */ ae_vector* r,
+     ae_state *_state);
+void corrr1dcircular(/* Real    */ ae_vector* signal,
+     ae_int_t m,
+     /* Real    */ ae_vector* pattern,
+     ae_int_t n,
+     /* Real    */ ae_vector* c,
+     ae_state *_state);
+void fhtr1d(/* Real    */ ae_vector* a, ae_int_t n, ae_state *_state);
+void fhtr1dinv(/* Real    */ ae_vector* a, ae_int_t n, ae_state *_state);
+
+}
+#endif
+

src/inc/alglib/integration.h

@@ -0,0 +1,837 @@
+/*************************************************************************
+Copyright (c) Sergey Bochkanov (ALGLIB project).
+
+>>> SOURCE LICENSE >>>
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation (www.fsf.org); either version 2 of the
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+A copy of the GNU General Public License is available at
+http://www.fsf.org/licensing/licenses
+>>> END OF LICENSE >>>
+*************************************************************************/
+#ifndef _integration_pkg_h
+#define _integration_pkg_h
+#include "ap.h"
+#include "alglibinternal.h"
+#include "linalg.h"
+#include "specialfunctions.h"
+
+/////////////////////////////////////////////////////////////////////////
+//
+// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES)
+//
+/////////////////////////////////////////////////////////////////////////
+namespace alglib_impl
+{
+typedef struct
+{
+    ae_int_t terminationtype;
+    ae_int_t nfev;
+    ae_int_t nintervals;
+} autogkreport;
+typedef struct
+{
+    double a;
+    double b;
+    double eps;
+    double xwidth;
+    double x;
+    double f;
+    ae_int_t info;
+    double r;
+    ae_matrix heap;
+    ae_int_t heapsize;
+    ae_int_t heapwidth;
+    ae_int_t heapused;
+    double sumerr;
+    double sumabs;
+    ae_vector qn;
+    ae_vector wg;
+    ae_vector wk;
+    ae_vector wr;
+    ae_int_t n;
+    rcommstate rstate;
+} autogkinternalstate;
+typedef struct
+{
+    double a;
+    double b;
+    double alpha;
+    double beta;
+    double xwidth;
+    double x;
+    double xminusa;
+    double bminusx;
+    ae_bool needf;
+    double f;
+    ae_int_t wrappermode;
+    autogkinternalstate internalstate;
+    rcommstate rstate;
+    double v;
+    ae_int_t terminationtype;
+    ae_int_t nfev;
+    ae_int_t nintervals;
+} autogkstate;
+
+}
+
+/////////////////////////////////////////////////////////////////////////
+//
+// THIS SECTION CONTAINS C++ INTERFACE
+//
+/////////////////////////////////////////////////////////////////////////
+namespace alglib
+{
+
+
+
+
+
+/*************************************************************************
+Integration report:
+* TerminationType = completetion code:
+    * -5    non-convergence of Gauss-Kronrod nodes
+            calculation subroutine.
+    * -1    incorrect parameters were specified
+    *  1    OK
+* Rep.NFEV countains number of function calculations
+* Rep.NIntervals contains number of intervals [a,b]
+  was partitioned into.
+*************************************************************************/
+class _autogkreport_owner
+{
+public:
+    _autogkreport_owner();
+    _autogkreport_owner(const _autogkreport_owner &rhs);
+    _autogkreport_owner& operator=(const _autogkreport_owner &rhs);
+    virtual ~_autogkreport_owner();
+    alglib_impl::autogkreport* c_ptr();
+    alglib_impl::autogkreport* c_ptr() const;
+protected:
+    alglib_impl::autogkreport *p_struct;
+};
+class autogkreport : public _autogkreport_owner
+{
+public:
+    autogkreport();
+    autogkreport(const autogkreport &rhs);
+    autogkreport& operator=(const autogkreport &rhs);
+    virtual ~autogkreport();
+    ae_int_t &terminationtype;
+    ae_int_t &nfev;
+    ae_int_t &nintervals;
+
+};
+
+
+/*************************************************************************
+This structure stores state of the integration algorithm.
+
+Although this class has public fields,  they are not intended for external
+use. You should use ALGLIB functions to work with this class:
+* autogksmooth()/AutoGKSmoothW()/... to create objects
+* autogkintegrate() to begin integration
+* autogkresults() to get results
+*************************************************************************/
+class _autogkstate_owner
+{
+public:
+    _autogkstate_owner();
+    _autogkstate_owner(const _autogkstate_owner &rhs);
+    _autogkstate_owner& operator=(const _autogkstate_owner &rhs);
+    virtual ~_autogkstate_owner();
+    alglib_impl::autogkstate* c_ptr();
+    alglib_impl::autogkstate* c_ptr() const;
+protected:
+    alglib_impl::autogkstate *p_struct;
+};
+class autogkstate : public _autogkstate_owner
+{
+public:
+    autogkstate();
+    autogkstate(const autogkstate &rhs);
+    autogkstate& operator=(const autogkstate &rhs);
+    virtual ~autogkstate();
+    ae_bool &needf;
+    double &x;
+    double &xminusa;
+    double &bminusx;
+    double &f;
+
+};
+
+/*************************************************************************
+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 real_1d_array &alpha, const real_1d_array &beta, const double mu0, const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &w);
+
+
+/*************************************************************************
+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(const real_1d_array &alpha, const real_1d_array &beta, const double mu0, const double a, const double b, const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &w);
+
+
+/*************************************************************************
+Computation of nodes and weights for a Gauss-Radau quadrature formula
+
+The algorithm generates the N-point Gauss-Radau  quadrature  formula  with
+weight function given by the 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-1], beta coefficients
+                Zero-indexed element is not used.
+                Beta[I]>0
+    Mu0     –   zeroth moment of the weighting function.
+    A       –   left boundary of the integration interval.
+    N       –   number of nodes of the quadrature formula, N>=2
+                (including the left boundary node).
+
+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 gqgenerategaussradaurec(const real_1d_array &alpha, const real_1d_array &beta, const double mu0, const double a, const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &w);
+
+
+/*************************************************************************
+Returns nodes/weights for Gauss-Legendre quadrature on [-1,1] with N
+nodes.
+
+INPUT PARAMETERS:
+    N           -   number of nodes, >=1
+
+OUTPUT PARAMETERS:
+    Info        -   error code:
+                    * -4    an  error   was   detected   when  calculating
+                            weights/nodes.  N  is  too  large   to  obtain
+                            weights/nodes  with  high   enough   accuracy.
+                            Try  to   use   multiple   precision  version.
+                    * -3    internal eigenproblem solver hasn't  converged
+                    * -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 12.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void gqgenerategausslegendre(const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &w);
+
+
+/*************************************************************************
+Returns  nodes/weights  for  Gauss-Jacobi quadrature on [-1,1] with weight
+function W(x)=Power(1-x,Alpha)*Power(1+x,Beta).
+
+INPUT PARAMETERS:
+    N           -   number of nodes, >=1
+    Alpha       -   power-law coefficient, Alpha>-1
+    Beta        -   power-law coefficient, Beta>-1
+
+OUTPUT PARAMETERS:
+    Info        -   error code:
+                    * -4    an  error  was   detected   when   calculating
+                            weights/nodes. Alpha or  Beta  are  too  close
+                            to -1 to obtain weights/nodes with high enough
+                            accuracy, or, may be, N is too large.  Try  to
+                            use multiple precision version.
+                    * -3    internal eigenproblem solver hasn't converged
+                    * -1    incorrect N/Alpha/Beta 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 12.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void gqgenerategaussjacobi(const ae_int_t n, const double alpha, const double beta, ae_int_t &info, real_1d_array &x, real_1d_array &w);
+
+
+/*************************************************************************
+Returns  nodes/weights  for  Gauss-Laguerre  quadrature  on  [0,+inf) with
+weight function W(x)=Power(x,Alpha)*Exp(-x)
+
+INPUT PARAMETERS:
+    N           -   number of nodes, >=1
+    Alpha       -   power-law coefficient, Alpha>-1
+
+OUTPUT PARAMETERS:
+    Info        -   error code:
+                    * -4    an  error  was   detected   when   calculating
+                            weights/nodes. Alpha is too  close  to  -1  to
+                            obtain weights/nodes with high enough accuracy
+                            or, may  be,  N  is  too  large.  Try  to  use
+                            multiple precision version.
+                    * -3    internal eigenproblem solver hasn't converged
+                    * -1    incorrect N/Alpha 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 12.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void gqgenerategausslaguerre(const ae_int_t n, const double alpha, ae_int_t &info, real_1d_array &x, real_1d_array &w);
+
+
+/*************************************************************************
+Returns  nodes/weights  for  Gauss-Hermite  quadrature on (-inf,+inf) with
+weight function W(x)=Exp(-x*x)
+
+INPUT PARAMETERS:
+    N           -   number of nodes, >=1
+
+OUTPUT PARAMETERS:
+    Info        -   error code:
+                    * -4    an  error  was   detected   when   calculating
+                            weights/nodes.  May be, N is too large. Try to
+                            use multiple precision version.
+                    * -3    internal eigenproblem solver hasn't converged
+                    * -1    incorrect N/Alpha 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 12.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void gqgenerategausshermite(const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &w);
+
+/*************************************************************************
+Computation of nodes and weights of a Gauss-Kronrod quadrature formula
+
+The algorithm generates the N-point Gauss-Kronrod 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 zero moment Mu0
+
+    Mu0 = integral(W(x)dx,a,b)
+
+
+INPUT PARAMETERS:
+    Alpha       –   alpha coefficients, array[0..floor(3*K/2)].
+    Beta        –   beta coefficients,  array[0..ceil(3*K/2)].
+                    Beta[0] is not used and may be arbitrary.
+                    Beta[I]>0.
+    Mu0         –   zeroth moment of the weight function.
+    N           –   number of nodes of the Gauss-Kronrod quadrature formula,
+                    N >= 3,
+                    N =  2*K+1.
+
+OUTPUT PARAMETERS:
+    Info        -   error code:
+                    * -5    no real and positive Gauss-Kronrod formula can
+                            be created for such a weight function  with  a
+                            given number of nodes.
+                    * -4    N is too large, task may be ill  conditioned -
+                            x[i]=x[i+1] found.
+                    * -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.
+    WKronrod    -   array[0..N-1] - Kronrod weights
+    WGauss      -   array[0..N-1] - Gauss weights (interleaved with zeros
+                    corresponding to extended Kronrod nodes).
+
+  -- ALGLIB --
+     Copyright 08.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void gkqgeneraterec(const real_1d_array &alpha, const real_1d_array &beta, const double mu0, const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &wkronrod, real_1d_array &wgauss);
+
+
+/*************************************************************************
+Returns   Gauss   and   Gauss-Kronrod   nodes/weights  for  Gauss-Legendre
+quadrature with N points.
+
+GKQLegendreCalc (calculation) or  GKQLegendreTbl  (precomputed  table)  is
+used depending on machine precision and number of nodes.
+
+INPUT PARAMETERS:
+    N           -   number of Kronrod nodes, must be odd number, >=3.
+
+OUTPUT PARAMETERS:
+    Info        -   error code:
+                    * -4    an  error   was   detected   when  calculating
+                            weights/nodes.  N  is  too  large   to  obtain
+                            weights/nodes  with  high   enough   accuracy.
+                            Try  to   use   multiple   precision  version.
+                    * -3    internal eigenproblem solver hasn't converged
+                    * -1    incorrect N was passed
+                    * +1    OK
+    X           -   array[0..N-1] - array of quadrature nodes, ordered in
+                    ascending order.
+    WKronrod    -   array[0..N-1] - Kronrod weights
+    WGauss      -   array[0..N-1] - Gauss weights (interleaved with zeros
+                    corresponding to extended Kronrod nodes).
+
+
+  -- ALGLIB --
+     Copyright 12.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void gkqgenerategausslegendre(const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &wkronrod, real_1d_array &wgauss);
+
+
+/*************************************************************************
+Returns   Gauss   and   Gauss-Kronrod   nodes/weights   for   Gauss-Jacobi
+quadrature on [-1,1] with weight function
+
+    W(x)=Power(1-x,Alpha)*Power(1+x,Beta).
+
+INPUT PARAMETERS:
+    N           -   number of Kronrod nodes, must be odd number, >=3.
+    Alpha       -   power-law coefficient, Alpha>-1
+    Beta        -   power-law coefficient, Beta>-1
+
+OUTPUT PARAMETERS:
+    Info        -   error code:
+                    * -5    no real and positive Gauss-Kronrod formula can
+                            be created for such a weight function  with  a
+                            given number of nodes.
+                    * -4    an  error  was   detected   when   calculating
+                            weights/nodes. Alpha or  Beta  are  too  close
+                            to -1 to obtain weights/nodes with high enough
+                            accuracy, or, may be, N is too large.  Try  to
+                            use multiple precision version.
+                    * -3    internal eigenproblem solver hasn't converged
+                    * -1    incorrect N was passed
+                    * +1    OK
+                    * +2    OK, but quadrature rule have exterior  nodes,
+                            x[0]<-1 or x[n-1]>+1
+    X           -   array[0..N-1] - array of quadrature nodes, ordered in
+                    ascending order.
+    WKronrod    -   array[0..N-1] - Kronrod weights
+    WGauss      -   array[0..N-1] - Gauss weights (interleaved with zeros
+                    corresponding to extended Kronrod nodes).
+
+
+  -- ALGLIB --
+     Copyright 12.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void gkqgenerategaussjacobi(const ae_int_t n, const double alpha, const double beta, ae_int_t &info, real_1d_array &x, real_1d_array &wkronrod, real_1d_array &wgauss);
+
+
+/*************************************************************************
+Returns Gauss and Gauss-Kronrod nodes for quadrature with N points.
+
+Reduction to tridiagonal eigenproblem is used.
+
+INPUT PARAMETERS:
+    N           -   number of Kronrod nodes, must be odd number, >=3.
+
+OUTPUT PARAMETERS:
+    Info        -   error code:
+                    * -4    an  error   was   detected   when  calculating
+                            weights/nodes.  N  is  too  large   to  obtain
+                            weights/nodes  with  high   enough   accuracy.
+                            Try  to   use   multiple   precision  version.
+                    * -3    internal eigenproblem solver hasn't converged
+                    * -1    incorrect N was passed
+                    * +1    OK
+    X           -   array[0..N-1] - array of quadrature nodes, ordered in
+                    ascending order.
+    WKronrod    -   array[0..N-1] - Kronrod weights
+    WGauss      -   array[0..N-1] - Gauss weights (interleaved with zeros
+                    corresponding to extended Kronrod nodes).
+
+  -- ALGLIB --
+     Copyright 12.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void gkqlegendrecalc(const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &wkronrod, real_1d_array &wgauss);
+
+
+/*************************************************************************
+Returns Gauss and Gauss-Kronrod nodes for quadrature with N  points  using
+pre-calculated table. Nodes/weights were  computed  with  accuracy  up  to
+1.0E-32 (if MPFR version of ALGLIB is used). In standard double  precision
+accuracy reduces to something about 2.0E-16 (depending  on your compiler's
+handling of long floating point constants).
+
+INPUT PARAMETERS:
+    N           -   number of Kronrod nodes.
+                    N can be 15, 21, 31, 41, 51, 61.
+
+OUTPUT PARAMETERS:
+    X           -   array[0..N-1] - array of quadrature nodes, ordered in
+                    ascending order.
+    WKronrod    -   array[0..N-1] - Kronrod weights
+    WGauss      -   array[0..N-1] - Gauss weights (interleaved with zeros
+                    corresponding to extended Kronrod nodes).
+
+
+  -- ALGLIB --
+     Copyright 12.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void gkqlegendretbl(const ae_int_t n, real_1d_array &x, real_1d_array &wkronrod, real_1d_array &wgauss, double &eps);
+
+/*************************************************************************
+Integration of a smooth function F(x) on a finite interval [a,b].
+
+Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result
+is calculated with accuracy close to the machine precision.
+
+Algorithm works well only with smooth integrands.  It  may  be  used  with
+continuous non-smooth integrands, but with  less  performance.
+
+It should never be used with integrands which have integrable singularities
+at lower or upper limits - algorithm may crash. Use AutoGKSingular in such
+cases.
+
+INPUT PARAMETERS:
+    A, B    -   interval boundaries (A<B, A=B or A>B)
+
+OUTPUT PARAMETERS
+    State   -   structure which stores algorithm state
+
+SEE ALSO
+    AutoGKSmoothW, AutoGKSingular, AutoGKResults.
+
+
+  -- ALGLIB --
+     Copyright 06.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void autogksmooth(const double a, const double b, autogkstate &state);
+
+
+/*************************************************************************
+Integration of a smooth function F(x) on a finite interval [a,b].
+
+This subroutine is same as AutoGKSmooth(), but it guarantees that interval
+[a,b] is partitioned into subintervals which have width at most XWidth.
+
+Subroutine  can  be  used  when  integrating nearly-constant function with
+narrow "bumps" (about XWidth wide). If "bumps" are too narrow, AutoGKSmooth
+subroutine can overlook them.
+
+INPUT PARAMETERS:
+    A, B    -   interval boundaries (A<B, A=B or A>B)
+
+OUTPUT PARAMETERS
+    State   -   structure which stores algorithm state
+
+SEE ALSO
+    AutoGKSmooth, AutoGKSingular, AutoGKResults.
+
+
+  -- ALGLIB --
+     Copyright 06.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void autogksmoothw(const double a, const double b, const double xwidth, autogkstate &state);
+
+
+/*************************************************************************
+Integration on a finite interval [A,B].
+Integrand have integrable singularities at A/B.
+
+F(X) must diverge as "(x-A)^alpha" at A, as "(B-x)^beta" at B,  with known
+alpha/beta (alpha>-1, beta>-1).  If alpha/beta  are  not known,  estimates
+from below can be used (but these estimates should be greater than -1 too).
+
+One  of  alpha/beta variables (or even both alpha/beta) may be equal to 0,
+which means than function F(x) is non-singular at A/B. Anyway (singular at
+bounds or not), function F(x) is supposed to be continuous on (A,B).
+
+Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result
+is calculated with accuracy close to the machine precision.
+
+INPUT PARAMETERS:
+    A, B    -   interval boundaries (A<B, A=B or A>B)
+    Alpha   -   power-law coefficient of the F(x) at A,
+                Alpha>-1
+    Beta    -   power-law coefficient of the F(x) at B,
+                Beta>-1
+
+OUTPUT PARAMETERS
+    State   -   structure which stores algorithm state
+
+SEE ALSO
+    AutoGKSmooth, AutoGKSmoothW, AutoGKResults.
+
+
+  -- ALGLIB --
+     Copyright 06.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void autogksingular(const double a, const double b, const double alpha, const double beta, autogkstate &state);
+
+
+/*************************************************************************
+This function provides reverse communication interface
+Reverse communication interface is not documented or recommended to use.
+See below for functions which provide better documented API
+*************************************************************************/
+bool autogkiteration(const autogkstate &state);
+
+
+/*************************************************************************
+This function is used to launcn iterations of the 1-dimensional integrator
+
+It accepts following parameters:
+    func    -   callback which calculates f(x) for given x
+    ptr     -   optional pointer which is passed to func; can be NULL
+
+
+  -- ALGLIB --
+     Copyright 07.05.2009 by Bochkanov Sergey
+
+*************************************************************************/
+void autogkintegrate(autogkstate &state,
+    void (*func)(double x, double xminusa, double bminusx, double &y, void *ptr),
+    void *ptr = NULL);
+
+
+/*************************************************************************
+Adaptive integration results
+
+Called after AutoGKIteration returned False.
+
+Input parameters:
+    State   -   algorithm state (used by AutoGKIteration).
+
+Output parameters:
+    V       -   integral(f(x)dx,a,b)
+    Rep     -   optimization report (see AutoGKReport description)
+
+  -- ALGLIB --
+     Copyright 14.11.2007 by Bochkanov Sergey
+*************************************************************************/
+void autogkresults(const autogkstate &state, double &v, autogkreport &rep);
+}
+
+/////////////////////////////////////////////////////////////////////////
+//
+// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS)
+//
+/////////////////////////////////////////////////////////////////////////
+namespace alglib_impl
+{
+void gqgeneraterec(/* Real    */ ae_vector* alpha,
+     /* Real    */ ae_vector* beta,
+     double mu0,
+     ae_int_t n,
+     ae_int_t* info,
+     /* Real    */ ae_vector* x,
+     /* Real    */ ae_vector* w,
+     ae_state *_state);
+void gqgenerategausslobattorec(/* Real    */ ae_vector* alpha,
+     /* Real    */ ae_vector* beta,
+     double mu0,
+     double a,
+     double b,
+     ae_int_t n,
+     ae_int_t* info,
+     /* Real    */ ae_vector* x,
+     /* Real    */ ae_vector* w,
+     ae_state *_state);
+void gqgenerategaussradaurec(/* Real    */ ae_vector* alpha,
+     /* Real    */ ae_vector* beta,
+     double mu0,
+     double a,
+     ae_int_t n,
+     ae_int_t* info,
+     /* Real    */ ae_vector* x,
+     /* Real    */ ae_vector* w,
+     ae_state *_state);
+void gqgenerategausslegendre(ae_int_t n,
+     ae_int_t* info,
+     /* Real    */ ae_vector* x,
+     /* Real    */ ae_vector* w,
+     ae_state *_state);
+void gqgenerategaussjacobi(ae_int_t n,
+     double alpha,
+     double beta,
+     ae_int_t* info,
+     /* Real    */ ae_vector* x,
+     /* Real    */ ae_vector* w,
+     ae_state *_state);
+void gqgenerategausslaguerre(ae_int_t n,
+     double alpha,
+     ae_int_t* info,
+     /* Real    */ ae_vector* x,
+     /* Real    */ ae_vector* w,
+     ae_state *_state);
+void gqgenerategausshermite(ae_int_t n,
+     ae_int_t* info,
+     /* Real    */ ae_vector* x,
+     /* Real    */ ae_vector* w,
+     ae_state *_state);
+void gkqgeneraterec(/* Real    */ ae_vector* alpha,
+     /* Real    */ ae_vector* beta,
+     double mu0,
+     ae_int_t n,
+     ae_int_t* info,
+     /* Real    */ ae_vector* x,
+     /* Real    */ ae_vector* wkronrod,
+     /* Real    */ ae_vector* wgauss,
+     ae_state *_state);
+void gkqgenerategausslegendre(ae_int_t n,
+     ae_int_t* info,
+     /* Real    */ ae_vector* x,
+     /* Real    */ ae_vector* wkronrod,
+     /* Real    */ ae_vector* wgauss,
+     ae_state *_state);
+void gkqgenerategaussjacobi(ae_int_t n,
+     double alpha,
+     double beta,
+     ae_int_t* info,
+     /* Real    */ ae_vector* x,
+     /* Real    */ ae_vector* wkronrod,
+     /* Real    */ ae_vector* wgauss,
+     ae_state *_state);
+void gkqlegendrecalc(ae_int_t n,
+     ae_int_t* info,
+     /* Real    */ ae_vector* x,
+     /* Real    */ ae_vector* wkronrod,
+     /* Real    */ ae_vector* wgauss,
+     ae_state *_state);
+void gkqlegendretbl(ae_int_t n,
+     /* Real    */ ae_vector* x,
+     /* Real    */ ae_vector* wkronrod,
+     /* Real    */ ae_vector* wgauss,
+     double* eps,
+     ae_state *_state);
+void autogksmooth(double a,
+     double b,
+     autogkstate* state,
+     ae_state *_state);
+void autogksmoothw(double a,
+     double b,
+     double xwidth,
+     autogkstate* state,
+     ae_state *_state);
+void autogksingular(double a,
+     double b,
+     double alpha,
+     double beta,
+     autogkstate* state,
+     ae_state *_state);
+ae_bool autogkiteration(autogkstate* state, ae_state *_state);
+void autogkresults(autogkstate* state,
+     double* v,
+     autogkreport* rep,
+     ae_state *_state);
+ae_bool _autogkreport_init(void* _p, ae_state *_state, ae_bool make_automatic);
+ae_bool _autogkreport_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
+void _autogkreport_clear(void* _p);
+void _autogkreport_destroy(void* _p);
+ae_bool _autogkinternalstate_init(void* _p, ae_state *_state, ae_bool make_automatic);
+ae_bool _autogkinternalstate_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
+void _autogkinternalstate_clear(void* _p);
+void _autogkinternalstate_destroy(void* _p);
+ae_bool _autogkstate_init(void* _p, ae_state *_state, ae_bool make_automatic);
+ae_bool _autogkstate_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
+void _autogkstate_clear(void* _p);
+void _autogkstate_destroy(void* _p);
+
+}
+#endif
+

src/inc/alglib/stdafx.h

@@ -0,0 +1,2 @@
+
+

src/inc/lag_rms.hpp

@@ -0,0 +1,58 @@
+/*
+ * lag.hpp
+ *
+ *  Created on: Jun 1, 2013
+ *      Author: azoghbi
+ */
+
+#ifndef LAG_RMS_HPP_
+#define LAG_RMS_HPP_
+
+#include "mod.hpp"
+
+class lag_rms : public mod {
+
+	int		n1;
+	vec		p1,p2;
+	ivec	icx,iphi;
+	double	mean1,mean2,mean1sq,mean2sq,mean12;
+
+	void	_C( vec );
+	void	_dC( vec , int );
+	void	_dCx( vec , int );
+	void	_dPhi( vec , int );
+
+
+public:
+	lag_rms( lcurve , lcurve , vec , vec );
+	virtual ~lag_rms();
+
+	void	step_pars( int , vec& , vec& );
+	void	print_pars( vec& , vec& );
+};
+
+// ----------------------------------- //
+
+class lag10_rms : public mod {
+
+	int		n1;
+	vec		p1,p2;
+	ivec	icx,iphi;
+	double	mean1,mean2,mean1sq,mean2sq,mean12;
+
+	void	_C( vec );
+	void	_dC( vec , int );
+	void	_dCx( vec , int );
+	void	_dPhi( vec , int );
+
+
+public:
+	lag10_rms( lcurve , lcurve , vec , vec );
+	virtual ~lag10_rms();
+
+	void	step_pars( int , vec& , vec& );
+	void	print_pars( vec& , vec& );
+	void	what_pars( int& , int& );
+};
+
+#endif /* LAG_HPP_ */

src/inc/psd_rms.hpp

@@ -0,0 +1,41 @@
+/*
+ * psd.hpp
+ *
+ *  Created on: May 31, 2013
+ *      Author: azoghbi
+ */
+
+#ifndef PSD_RMS_HPP_
+#define PSD_RMS_HPP_
+
+#include "mod.hpp"
+
+class psd_rms: public mod {
+
+	double	mean,mean2;
+
+	void	_C( vec );
+	void	_dC( vec , int );
+
+public:
+	psd_rms( lcurve , vec );
+	virtual ~psd_rms();
+};
+
+// ---------------------- //
+
+class psd10_rms: public mod {
+
+	double	mean,mean2;
+
+	void	_C( vec );
+	void	_dC( vec , int );
+
+public:
+	psd10_rms( lcurve , vec );
+	virtual ~psd10_rms();
+
+	void	step_pars( int , vec& , vec& );
+};
+
+#endif /* PSD_RMS_HPP_ */

src/inc/psdlag_rms.hpp

@@ -0,0 +1,56 @@
+/*
+ * psdlag.hpp
+ *
+ *  Created on: Jun 1, 2013
+ *      Author: azoghbi
+ */
+
+#include "mod.hpp"
+
+#ifndef PSDLAG_RMS_HPP_
+#define PSDLAG_RMS_HPP_
+
+class psdlag_rms : public mod {
+
+	int		n1;
+	ivec	ip1,ip2,icx,iphi;
+	double	mean1,mean2,mean1sq,mean2sq,mean12;
+
+	void	_C( vec );
+	void	_dC( vec , int );
+	void	_dP1( vec , int );
+	void	_dP2( vec , int );
+	void	_dCx( vec , int );
+	void	_dPhi( vec , int );
+
+public:
+	psdlag_rms( lcurve , lcurve , vec );
+	virtual ~psdlag_rms();
+
+	void	step_pars( int , vec& , vec& );
+	void	print_pars( vec& , vec& );
+};
+
+class psdlag10_rms : public mod {
+
+	int		n1;
+	ivec	ip1,ip2,icx,iphi;
+	double	mean1,mean2,mean1sq,mean2sq,mean12;
+
+	void	_C( vec );
+	void	_dC( vec , int );
+	void	_dP1( vec , int );
+	void	_dP2( vec , int );
+	void	_dCx( vec , int );
+	void	_dPhi( vec , int );
+
+public:
+	psdlag10_rms( lcurve , lcurve , vec );
+	virtual ~psdlag10_rms();
+
+	void	step_pars( int , vec& , vec& );
+	void	print_pars( vec& , vec& );
+	void	what_pars( int&, int& );
+};
+
+#endif /* PSDLAG_RMS_HPP_ */

src/lag_rms.cpp

@@ -0,0 +1,230 @@
+/*
+ * lag.cpp
+ *
+ *  Created on: Jun 1, 2013
+ *      Author: azoghbi
+ */
+
+#include "inc/lag_rms.hpp"
+
+lag_rms::lag_rms( lcurve lc1 , lcurve lc2 , vec fqL , vec pars ) {
+
+	// ----------- initial parameters ------------ //
+	n1		=	lc1.len;
+	n		=	n1 + lc2.len;
+	dt		=	lc1.dt;
+	// ------------------------------------------ //
+
+
+	// ----------- light curve setup ------------ //
+	setlc();
+	mean1	=	lc1.demean(); mean2	=	lc2.demean();
+	mean1sq	=	mean1*mean1; mean2sq = mean2*mean2; mean12 = mean1*mean2;
+	int		i;
+	for( i=0 ; i<n1 ; i++ ){
+		lc[i]	=	lc1.lc[i];
+		lce[i]	=	lc1.lce[i]*lc1.lce[i];
+		t[i]	=	lc1.t[i];
+	}
+	for( i=n1 ; i<n ; i++ ){
+		lc[i]	=	lc2.lc[i-n1];
+		lce[i]	=	lc2.lce[i-n1]*lc2.lce[i-n1];
+		t[i]	=	lc2.t[i-n1];
+	}
+	// ------------------------------------------ //
+
+	// ----------- parent initilizer ------------ //
+	init( fqL , 2 );
+	// ------------------------------------------ //
+
+	// --------- constants and indices ---------- //
+	p1.setlength(nfq); p2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
+	for( i=0 ; i<nfq ; i++ ){
+		p1[i]	=	pars[i] * mean1sq;
+		p2[i]	=	pars[i+nfq] * mean2sq;
+		icx[i]	=	i;
+		iphi[i]	=	i+nfq;
+	}
+	// ------------------------------------------ //
+
+}
+
+lag_rms::~lag_rms() {}
+
+void lag_rms::_C( vec p ){
+	int	i,j,k; double d;
+	for( i=0 ; i<nfq ; i++ ){ p[icx[i]] *= mean12;}
+	for( i=0 ; i<n1 ; i++ ){
+		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
+		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,j); } C(i,j) = d;}
+	}
+	for( i=n1 ; i<n ; i++ ){
+		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
+		for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,j); } C(i,j) = d;}
+		for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
+			d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
+	}
+}
+
+void lag_rms::_dC( vec p , int k ){
+	if( k<nfq ){	_dCx( p , k );
+	}else{			_dPhi( p , k-nfq );}
+}
+
+void lag_rms::_dCx( vec p , int k ){
+	int	i,j; double phi=p[iphi[k]];
+	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
+	for( i=n1 ; i<n ; i++ ){
+		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
+		for( j=0 ; j<n1 ; j++ ){
+			dC(i,j) = ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
+	}
+}
+
+void lag_rms::_dPhi( vec p , int k ){
+	int	i,j; double cx=p[k]*mean12,phi=p[iphi[k]];
+	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
+	for( i=n1 ; i<n ; i++ ){
+		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
+		for( j=0 ; j<n1 ; j++ ){
+			dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
+	}
+}
+
+
+void lag_rms::step_pars( int n , vec& dpar , vec& pars ){
+	for( int i=0 ; i<npar ; i++ ){
+		if( dpar[i]>3000 ){dpar[i] = 3000;} if( dpar[i]<-3000 ){dpar[i] = -3000;}
+		pars[i] += dpar[i]/((n<10)?10:1);
+	}
+}
+
+
+void lag_rms::print_pars( vec& pars , vec& errs ){
+	for( int i=0 ; i<nfq ; i++ ){
+		if( pars[i]<0 ){ pars[i]*=-1; pars[i+nfq] += M_PI;}
+		while( pars[i+nfq] > M_PI ){ pars[i+nfq] -= 2*M_PI; }
+		while( pars[i+nfq] <-M_PI ){ pars[i+nfq] += 2*M_PI; }
+	}
+	mod::print_pars( pars , errs );
+}
+
+
+// ****************************************** //
+
+
+lag10_rms::lag10_rms( lcurve lc1 , lcurve lc2 , vec fqL , vec pars ) {
+
+	// ----------- initial parameters ------------ //
+	n1		=	lc1.len;
+	n		=	n1 + lc2.len;
+	dt		=	lc1.dt;
+	// ------------------------------------------ //
+
+
+	// ----------- light curve setup ------------ //
+	setlc();
+	mean1	=	lc1.demean(); mean2	=	lc2.demean();
+	mean1sq	=	mean1*mean1; mean2sq = mean2*mean2; mean12 = mean1*mean2;
+	int		i;
+	for( i=0 ; i<n1 ; i++ ){
+		lc[i]	=	lc1.lc[i];
+		lce[i]	=	lc1.lce[i]*lc1.lce[i];
+		t[i]	=	lc1.t[i];
+	}
+	for( i=n1 ; i<n ; i++ ){
+		lc[i]	=	lc2.lc[i-n1];
+		lce[i]	=	lc2.lce[i-n1]*lc2.lce[i-n1];
+		t[i]	=	lc2.t[i-n1];
+	}
+	// ------------------------------------------ //
+
+	// ----------- parent initilizer ------------ //
+	init( fqL , 2 );
+	// ------------------------------------------ //
+
+	// --------- constants and indices ---------- //
+	p1.setlength(nfq); p2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
+	for( i=0 ; i<nfq ; i++ ){
+		p1[i]	=	pow(10,pars[i]) * mean1sq;
+		p2[i]	=	pow(10,pars[i+nfq]) * mean2sq;
+		icx[i]	=	i;
+		iphi[i]	=	i+nfq;
+	}
+	// ------------------------------------------ //
+
+}
+
+lag10_rms::~lag10_rms() {}
+
+void lag10_rms::_C( vec p ){
+	int	i,j,k; double d; for(i=0;i<nfq;i++){p[icx[i]] = pow(10,p[i]) * mean12;}
+	for( i=0 ; i<n1 ; i++ ){
+		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,i); }
+		d += f1*p1[0]*Cfq2[0](i,i) + f2*p1[nfq-1]*Cfq2[2](i,i);
+		d += lce[i]; C(i,i) = d;
+		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,j); }
+		d += f1*p1[0]*Cfq2[0](i,j) + f2*p1[nfq-1]*Cfq2[2](i,j);
+		C(i,j) = d;}
+	}
+	for( i=n1 ; i<n ; i++ ){
+		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,i); }
+		d += f1*p2[0]*Cfq2[0](i,i) + f2*p2[nfq-1]*Cfq2[2](i,i);
+		d += lce[i]; C(i,i) = d;
+		for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,j); }
+		d += f1*p2[0]*Cfq2[0](i,j) + f2*p2[nfq-1]*Cfq2[2](i,j);
+		C(i,j) = d;}
+		for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
+			d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); }
+		d += f1*p[icx[0]] * ( Cfq2[0](i,j)*cos(p[iphi[0]]) - Sfq2[0](i,j)*sin(p[iphi[0]]) );
+		d += f2*p[icx[nfq-1]] * ( Cfq2[2](i,j)*cos(p[iphi[nfq-1]]) - Sfq2[2](i,j)*sin(p[iphi[nfq-1]]) );
+		C(i,j) = d;}
+	}
+}
+
+void lag10_rms::_dC( vec p , int k ){
+	if( k<nfq ){	_dCx( p , k );
+	}else{			_dPhi( p , k-nfq );}
+}
+
+void lag10_rms::_dCx( vec p , int k ){
+	int	i,j; double phi=p[iphi[k]],cx = log(10)*pow(10,p[k])*mean12;
+	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
+	for( i=n1 ; i<n ; i++ ){
+		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
+		for( j=0 ; j<n1 ; j++ ){
+			dC(i,j) = cx * ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
+	}
+}
+
+void lag10_rms::_dPhi( vec p , int k ){
+	int	i,j; double cx=pow(10,p[k])*mean12,phi=p[iphi[k]];
+	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
+	for( i=n1 ; i<n ; i++ ){
+		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
+		for( j=0 ; j<n1 ; j++ ){
+			dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
+	}
+}
+
+
+void lag10_rms::step_pars( int n , vec& dpar , vec& pars ){
+	for( int i=0 ; i<npar ; i++ ){
+		if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
+		//pars[i] += dpar[i];
+		pars[i] += dpar[i]/((n<5)?10:1);
+	}
+}
+
+
+void lag10_rms::print_pars( vec& pars , vec& errs ){
+	for( int i=0 ; i<nfq ; i++ ){
+		while( pars[i+nfq] > M_PI ){ pars[i+nfq] -= 2*M_PI; }
+		while( pars[i+nfq] <-M_PI ){ pars[i+nfq] += 2*M_PI; }
+	}
+	mod::print_pars( pars , errs );
+}
+
+void lag10_rms::what_pars( int& ip1 , int& ip2 ){
+	ip1 = nfq; ip2 = npar;
+}

src/main.cpp

@@ -59,20 +59,11 @@ void do_work( char* fname ){
 	fp.close();
 	/*********   END Reading the input file  **********/
 
-	std::cout << "Completed reading file." << bin1 << bin2 
-			<< files[0]
-			//<< files[1]
-			<< std::endl;
-
 	/**********  Read the light curves  **************/
-	//std::cout << "hello";
-	//std::cout << bin1 << bin2;
 	vector<vector<lcurve> > LC,LC_ref;
 	if ( ref == -10 ){
-		//std::cout << "yo";
 		for(i=0;i<nfiles;i++)  readLC( LC , files[i] , secL[i] , bin1 , bin2 , strict );
 	}else {
-		//std::cout << "yo2";
 		for(i=0;i<nfiles;i++) {
 			if( i!=ref ) { readLC( LC , files[i] , secL[i] , bin2 , -1 , strict ); }
 		}
@@ -80,8 +71,6 @@ void do_work( char* fname ){
 	}
 	/********  END Read the light curves  ************/
 
-	//std::cout << "Read lcs." << std::endl;
-
 	if( mode > 0 or mode==-1 ){
 		vec errs; errs.setlength(nfq);
 		vector<lcurve> lc1; for( i=0 ; i<int(LC.size()) ; i++ ){ lc1.push_back( LC[i][0]); }

src/makefile

@@ -6,5 +6,6 @@ libdir='/home/caes/science/psdlag-agn/src'
 incdir='/home/caes/science/psdlag-agn/src/inc'
 
 
+
 psdlag:
 	g++ *cpp -o psdlag -O3 -Wall -lalglib -I${incdir} -L${libdir}

src/psd_rms.cpp

@@ -0,0 +1,108 @@
+/*
+ * psd.cpp
+ *
+ *  Created on: May 31, 2013
+ *      Author: azoghbi
+ */
+
+#include "inc/psd_rms.hpp"
+
+psd_rms::psd_rms( lcurve inlc , vec fqL ) {
+
+	// ----------- initial parameters ------------ //
+	n		=	inlc.len;
+	dt		=	inlc.dt;
+	// ------------------------------------------ //
+
+
+	// ----------- light curve setup ------------ //
+	setlc();
+	mean	=	inlc.demean();
+	mean2	=	mean*mean;
+	int		i;
+	for( i=0 ; i<n ; i++ ){
+		lc[i]	=	inlc.lc[i];
+		lce[i]	=	inlc.lce[i]*inlc.lce[i];
+		t[i]	=	inlc.t[i];
+	}
+	// ------------------------------------------ //
+
+	// ----------- parent initilizer ------------ //
+	init( fqL , 1 );
+	// ------------------------------------------ //
+}
+
+psd_rms::~psd_rms() {}
+
+void psd_rms::_C( vec p ){
+	int	i,j,k; double d;
+	for( i=0 ; i<npar ; i++ ){ p[i] *= mean2; }
+	for( i=0 ; i<n ; i++ ){
+		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,i); }
+		d += f1*p[0]*Cfq2[0](i,i) + f2*p[nfq-1]*Cfq2[2](i,i);
+		d += lce[i]; C(i,i) = d;
+		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,j); }
+		d += f1*p[0]*Cfq2[0](i,j) + f2*p[nfq-1]*Cfq2[2](i,j);
+		C(i,j) = d;}
+	}
+}
+
+void psd_rms::_dC( vec p , int k ){
+	int	i,j;
+	for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = Cfq[k](i,j); dC(j,i) = dC(i,j);} }
+}
+
+
+// ********************************************** //
+// ********************************************** //
+
+psd10_rms::psd10_rms( lcurve inlc , vec fqL ) {
+
+	// ----------- initial parameters ------------ //
+	n		=	inlc.len;
+	dt		=	inlc.dt;
+	// ------------------------------------------ //
+
+
+	// ----------- light curve setup ------------ //
+	setlc();
+	mean		=	inlc.demean();
+	mean2	=	mean*mean;
+	int		i;
+	for( i=0 ; i<n ; i++ ){
+		lc[i]	=	inlc.lc[i];
+		lce[i]	=	inlc.lce[i]*inlc.lce[i];
+		t[i]	=	inlc.t[i];
+	}
+	// ------------------------------------------ //
+
+	// ----------- parent initilizer ------------ //
+	init( fqL , 1 );
+	// ------------------------------------------ //
+}
+
+psd10_rms::~psd10_rms() {}
+
+void psd10_rms::_C( vec p ){
+	int	i,j,k; double d; for(i=0;i<npar;i++){p[i]=pow(10,p[i])*mean2;}
+	for( i=0 ; i<n ; i++ ){
+		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,i); }
+		d += f1*p[0]*Cfq2[0](i,i) + f2*p[nfq-1]*Cfq2[2](i,i);
+		d += lce[i]; C(i,i) = d;
+		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,j); }
+		d += f1*p[0]*Cfq2[0](i,j) + f2*p[nfq-1]*Cfq2[2](i,j);
+		C(i,j) = d;}
+	}
+}
+
+void psd10_rms::_dC( vec p , int k ){
+	int	i,j; double d = log(10)*pow(10,p[k])*mean2;
+	for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = d*Cfq[k](i,j); dC(j,i) = dC(i,j);} }
+}
+
+void psd10_rms::step_pars( int n , vec& dpar , vec& pars ){
+	for( int i=0 ; i<npar ; i++ ){
+		if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
+		pars[i] += dpar[i];
+	}
+}

src/psdlag_rms.cpp

@@ -0,0 +1,255 @@
+/*
+ * psdlag.cpp
+ *
+ *  Created on: Jun 1, 2013
+ *      Author: azoghbi
+ */
+
+#include "inc/psdlag_rms.hpp"
+
+psdlag_rms::psdlag_rms( lcurve lc1, lcurve lc2 , vec fqL ) {
+
+	// ----------- initial parameters ------------ //
+	n1		=	lc1.len;
+	n		=	n1 + lc2.len;
+	dt		=	lc1.dt;
+	// ------------------------------------------ //
+
+
+	// ----------- light curve setup ------------ //
+	setlc();
+	mean1	=	lc1.demean(); mean2	=	lc2.demean();
+	mean1sq	=	mean1*mean1; mean2sq = mean2*mean2; mean12 = mean1*mean2;
+	int		i;
+	for( i=0 ; i<n1 ; i++ ){
+		lc[i]	=	lc1.lc[i];
+		lce[i]	=	lc1.lce[i]*lc1.lce[i];
+		t[i]	=	lc1.t[i];
+	}
+	for( i=n1 ; i<n ; i++ ){
+		lc[i]	=	lc2.lc[i-n1];
+		lce[i]	=	lc2.lce[i-n1]*lc2.lce[i-n1];
+		t[i]	=	lc2.t[i-n1];
+	}
+	// ------------------------------------------ //
+
+	// ----------- parent initilizer ------------ //
+	init( fqL , 4 );
+	// ------------------------------------------ //
+
+	// --------- constants and indices ---------- //
+	ip1.setlength(nfq); ip2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
+	for( i=0 ; i<nfq ; i++ ){
+		ip1[i]	=	i;
+		ip2[i]	=	i+nfq;
+		icx[i]	=	i+2*nfq;
+		iphi[i]	=	i+3*nfq;
+	}
+	// ------------------------------------------ //
+
+}
+
+psdlag_rms::~psdlag_rms() {}
+
+void psdlag_rms::_C( vec p ){
+	int	i,j,k; double d;
+	for( i=0 ; i<nfq ; i++ ){ p[ip1[i]] *= mean1sq; p[ip2[i]] *= mean2sq; p[icx[i]] *= mean12;}
+	for( i=0 ; i<n1 ; i++ ){
+		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
+		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,j); } C(i,j) = d;}
+	}
+	for( i=n1 ; i<n ; i++ ){
+		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
+		for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,j); } C(i,j) = d;}
+		for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
+			d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
+	}
+}
+
+void psdlag_rms::_dC( vec p , int k ){
+	if( k<nfq ){			_dP1( p , k );
+	}else if( k<2*nfq ){	_dP2( p , k%nfq );
+	}else if( k<3*nfq ){	_dCx( p , k%nfq );
+	}else{					_dPhi( p , k%nfq );}
+}
+
+void psdlag_rms::_dP1( vec p , int k ){
+	int	i,j;
+	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = Cfq[k](i,j); dC(j,i) = dC(i,j);} }
+	for( i=n1 ; i<n ; i++ ){
+		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
+		for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
+	}
+}
+
+void psdlag_rms::_dP2( vec p , int k ){
+	int	i,j;
+	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);} }
+	for( i=n1 ; i<n ; i++ ){
+		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = Cfq[k](i,j); dC(j,i) = dC(i,j);}
+		for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
+	}
+}
+
+void psdlag_rms::_dCx( vec p , int k ){
+	int	i,j; double phi=p[iphi[k]];
+	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
+	for( i=n1 ; i<n ; i++ ){
+		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
+		for( j=0 ; j<n1 ; j++ ){
+			dC(i,j) = ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
+	}
+}
+
+void psdlag_rms::_dPhi( vec p , int k ){
+	int	i,j; double cx=p[icx[k]]*mean12,phi=p[iphi[k]];
+	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
+	for( i=n1 ; i<n ; i++ ){
+		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
+		for( j=0 ; j<n1 ; j++ ){
+			dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
+	}
+}
+
+void psdlag_rms::step_pars( int n , vec& dpar , vec& pars ){
+	for( int i=0 ; i<npar ; i++ ){
+		if( dpar[i]>3000 ){dpar[i] = 3000;} if( dpar[i]<-3000 ){dpar[i] = -3000;}
+		pars[i] += dpar[i]/((n<10)?10:1);
+	}
+}
+
+
+void psdlag_rms::print_pars( vec& pars , vec& errs ){
+	for( int i=0 ; i<nfq ; i++ ){
+		if( pars[i+2*nfq]<0 ){ pars[i+2*nfq]*=-1; pars[i+3*nfq] += M_PI;}
+		while( pars[i+3*nfq] > M_PI ){ pars[i+3*nfq] -= 2*M_PI; }
+		while( pars[i+3*nfq] <-M_PI ){ pars[i+3*nfq] += 2*M_PI; }
+		mod::print_pars( pars , errs );
+	}
+}
+
+// ++++++++++++++++++++++++++++++++++++++++++++++++++++++ //
+
+psdlag10_rms::psdlag10_rms( lcurve lc1, lcurve lc2 , vec fqL ) {
+
+	// ----------- initial parameters ------------ //
+	n1		=	lc1.len;
+	n		=	n1 + lc2.len;
+	dt		=	lc1.dt;
+	// ------------------------------------------ //
+
+
+	// ----------- light curve setup ------------ //
+	setlc();
+	mean1	=	lc1.demean(); mean2	=	lc2.demean();
+	mean1sq	=	mean1*mean1; mean2sq = mean2*mean2; mean12 = mean1*mean2;
+	int		i;
+	for( i=0 ; i<n1 ; i++ ){
+		lc[i]	=	lc1.lc[i];
+		lce[i]	=	lc1.lce[i]*lc1.lce[i];
+		t[i]	=	lc1.t[i];
+	}
+	for( i=n1 ; i<n ; i++ ){
+		lc[i]	=	lc2.lc[i-n1];
+		lce[i]	=	lc2.lce[i-n1]*lc2.lce[i-n1];
+		t[i]	=	lc2.t[i-n1];
+	}
+	// ------------------------------------------ //
+
+	// ----------- parent initilizer ------------ //
+	init( fqL , 4 );
+	// ------------------------------------------ //
+
+	// --------- constants and indices ---------- //
+	ip1.setlength(nfq); ip2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
+	for( i=0 ; i<nfq ; i++ ){
+		ip1[i]	=	i;
+		ip2[i]	=	i+nfq;
+		icx[i]	=	i+2*nfq;
+		iphi[i]	=	i+3*nfq;
+	}
+	// ------------------------------------------ //
+
+}
+
+psdlag10_rms::~psdlag10_rms() {}
+
+void psdlag10_rms::_C( vec p ){
+	int	i,j,k; double d; for(i=0;i<nfq;i++){ p[ip1[i]] = pow(10,p[ip1[i]])*mean1sq; p[ip2[i]] = pow(10,p[ip2[i]])*mean2sq; p[icx[i]] = pow(10,p[icx[i]])*mean12;}
+	for( i=0 ; i<n1 ; i++ ){
+		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
+		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,j); } C(i,j) = d;}
+	}
+	for( i=n1 ; i<n ; i++ ){
+		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
+		for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,j); } C(i,j) = d;}
+		for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
+			d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
+	}
+}
+
+void psdlag10_rms::_dC( vec p , int k ){
+	if( k<nfq ){			_dP1( p , k );
+	}else if( k<2*nfq ){	_dP2( p , k%nfq );
+	}else if( k<3*nfq ){	_dCx( p , k%nfq );
+	}else{					_dPhi( p , k%nfq );}
+}
+
+void psdlag10_rms::_dP1( vec p , int k ){
+	int	i,j; double p1 = log(10)*pow(10,p[k])*mean1sq;
+	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = p1 * Cfq[k](i,j); dC(j,i) = dC(i,j);} }
+	for( i=n1 ; i<n ; i++ ){
+		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
+		for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
+	}
+}
+
+void psdlag10_rms::_dP2( vec p , int k ){
+	int	i,j; double p2 = log(10) * pow(10,p[ip2[k]])*mean2sq;
+	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);} }
+	for( i=n1 ; i<n ; i++ ){
+		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = p2 * Cfq[k](i,j); dC(j,i) = dC(i,j);}
+		for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
+	}
+}
+
+void psdlag10_rms::_dCx( vec p , int k ){
+	int	i,j; double phi=p[iphi[k]], cx=log(10)*pow(10,p[icx[k]])*mean12;
+	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
+	for( i=n1 ; i<n ; i++ ){
+		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
+		for( j=0 ; j<n1 ; j++ ){
+			dC(i,j) = cx * ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
+	}
+}
+
+void psdlag10_rms::_dPhi( vec p , int k ){
+	int	i,j; double cx=pow(10,p[icx[k]])*mean12,phi=p[iphi[k]];
+	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
+	for( i=n1 ; i<n ; i++ ){
+		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
+		for( j=0 ; j<n1 ; j++ ){
+			dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
+	}
+}
+
+void psdlag10_rms::step_pars( int n , vec& dpar , vec& pars ){
+	for( int i=0 ; i<npar ; i++ ){
+		if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
+		//pars[i] += dpar[i];
+		pars[i] += dpar[i]/((n<5)?10:1);
+	}
+}
+
+
+void psdlag10_rms::print_pars( vec& pars , vec& errs ){
+	for( int i=0 ; i<nfq ; i++ ){
+		while( pars[i+3*nfq] > M_PI ){ pars[i+3*nfq] -= 2*M_PI; }
+		while( pars[i+3*nfq] <-M_PI ){ pars[i+3*nfq] += 2*M_PI; }
+	}
+	mod::print_pars( pars , errs );
+}
+
+void psdlag10_rms::what_pars( int& ip1 , int& ip2 ){
+	ip1 = 3*nfq; ip2 = 4*nfq;
+}