added psdlag source code

psdlag/src/.smhist

@@ -0,0 +1,47 @@
+data check.dat
+read f 1
+read pow 2
+limits f pow
+box
+ptype 10 3
+points f pow
+quit
+data check.dat
+read f 1
+read pow 2
+limits f pow
+ptype 10 3
+points f pow
+box
+quit
+data check.dat
+read f 1
+read pow 2
+limits f pow
+box
+ptype 10 3
+points f pow
+set lgf = log10(f)
+set lgf = alog10(f)
+set lgf = lg(f)
+erase
+box
+erase
+limits lgf pow
+points lgf pow
+box
+erase
+limits lgf 0. 3.
+box
+points lgf pow
+quit
+data check.dat
+read 1 f
+read f 1
+read pow 2
+set lgf = lg(f)
+limits lgf pow
+box
+ptype 10 3
+points lgf pow
+quit

psdlag/src/README

@@ -0,0 +1,32 @@
+
+
+
+To compile, modify the file makefile to so that libdir and incdir point to where
+alglib library and header files are installed respectively, then run: make
+
+
+FILES:
+- README: this file
+- makefile: a make file for compiling
+- main.cpp: the main file, that process the input file, print help message if
+  needed and reads the light curves files. check readLC to see the light curve
+  format.
+- mod.cpp: contains the class mod, which is a base class that uses the likelihood method. This on its own does nothing, but it is meant to be inherited by others, e.g. psd , lag , lag10
+  etc ... It is here where the likelihood is calculated, and maximized.
+
+- psd.cpp: a child class of mod to calculate psds. psd10 is similar by it
+  calculate the log10 of psd. This is better than calculating the psd directly.
+
+- lag.cpp: a child class of mod to calculate cross spectra and phase lags.
+  similar to psd. there is also lag10 that calcualte the log10 of the cross
+  spectrum plus phase lag.
+
+- psdlag.cpp: a child class of mod to calculate psd and cross spec and phase
+  lag at the same time.
+
+- inc: contains the header files.
+
+- inc/mod.hpp: the header for mod.cpp. This file has the Mod class also, which
+  is a class that holds multiple mod's used for fitting multiple light curve
+  simultaneusly, each mod fits for one light curve, and the total log-likelihood
+  is the sum of individual likelihood (i.e. product of probabilities)

psdlag/src/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
+

psdlag/src/check.dat

@@ -0,0 +1,6 @@
+2.5e-5	4.681
+5.5e-5	3.358
+8.5e-5	2.911
+1.5e-4	2.650
+2.5e-4	1.864
+3.5e-4	0.605

psdlag/src/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
+

psdlag/src/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
+

psdlag/src/fittedlags.dat

@@ -0,0 +1,6 @@
+1.022e-01 9.694e-01
+3.260e-01 9.192e-01
+5.118e-01 9.776e-01
+3.254e-01 9.269e-01
+1.000e+00 9.659e-01
+1.000e+00 9.891e-01

psdlag/src/inc/def.hpp

@@ -0,0 +1,36 @@
+/*
+ * def.hpp
+ *
+ *  Created on: May 24, 2013
+ *      Author: abduz
+ */
+
+#ifndef DEF_HPP_
+#define DEF_HPP_
+
+#include <ap.h>
+#include <iostream>
+#include <iomanip>
+
+using namespace alglib;
+using namespace std;
+
+typedef real_1d_array vec;
+typedef real_2d_array vec2;
+typedef integer_1d_array ivec;
+
+
+class lcurve {
+public:
+	int len; double dt; vec lc,lce,t;
+	lcurve(int n=1){t.setlength(n); lc.setlength(n); lce.setlength(n) ;len=n;dt=1.0;}
+	lcurve(vec& t_,vec& lc_,vec& lce_,double dt_){t=t_;lc=lc_;lce=lce_;len=t.length();dt=dt_;}
+	void demean(){double m=0;int i;for(i=0;i<len;i++){m+=lc[i];} m/=len;for(i=0;i<len;i++){lc[i]-=m;}}
+	friend ostream& operator<<(ostream&o,lcurve&l){ o << setprecision(12);
+			for(int i=0;i<l.len;i++){o << l.t[i] << " " << l.lc[i] << " " << l.lce[i] << endl;}
+			return(o);}
+};
+
+
+
+#endif /* DEF_HPP_ */

psdlag/src/inc/lag.hpp

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

psdlag/src/inc/main.hpp

@@ -0,0 +1,41 @@
+/*
+ * main.hpp
+ *
+ *  Created on: May 31, 2013
+ *      Author: azoghbi
+ */
+
+#ifndef MAIN_HPP_
+#define MAIN_HPP_
+
+
+#include <fstream>
+#include <cstring>
+#include <sstream>
+#include <vector>
+
+#include "def.hpp"
+#include "psd.hpp"
+#include "lag.hpp"
+#include "psdlag.hpp"
+#include "mcmc.hpp"
+
+
+void usage();
+void do_work( char* );
+void readLC( vector<vector<lcurve> >&, string , int , int , int , bool );
+
+
+double mcmc_lag10( vec& x ,void*ptr ){
+	int i,np = x.length(); for(i=0;i<np/2;i++){if(x[i]>7){x[i]=7;}if(x[i]<-7){x[i]=-7;}}
+	Mod<lag10> *mod = (Mod<lag10>*) ptr; double logl=mod->loglikelihood(x);
+	return logl;
+}
+
+double mcmc_psdlag10( vec& x ,void*ptr ){
+	int i,np = x.length(); for(i=0;i<np/2;i++){if(x[i]>7){x[i]=7;}if(x[i]<-7){x[i]=-7;}}
+	Mod<psdlag10> *mod = (Mod<psdlag10>*) ptr; double logl=mod->loglikelihood(x);
+	return logl;
+}
+
+#endif /* MAIN_HPP_ */

psdlag/src/inc/mcmc.hpp

@@ -0,0 +1,35 @@
+/*
+ * mcmc.hpp
+ *
+ *  Created on: May 14, 2013
+ *      Author: abduz
+ */
+
+#ifndef MCMC_HPP_
+#define MCMC_HPP_
+
+#include <ap.h>
+#include <alglibmisc.h>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+using namespace alglib;
+using namespace std;
+
+typedef real_1d_array vec;
+typedef real_2d_array vec2;
+typedef integer_1d_array ivec;
+typedef integer_2d_array ivec2;
+
+class mcmc {
+	int npar; double(*loglikelihood)(vec&,void*); void*ptr;
+	hqrndstate rnd;
+public:
+	int nrun,nburn,nwk,ncheck; double avalue;
+	mcmc(int np,double (*f)(vec&,void*),void *p);
+	void setseed(int s1,int s2 ){hqrndseed(s1,s2,rnd);}
+	void run( vec& , vec& , const char*fname="mcmc.dat");
+	virtual ~mcmc();
+};
+
+#endif /* MCMC_HPP_ */

psdlag/src/inc/mod.hpp

@@ -0,0 +1,283 @@
+/*
+ * mod.hpp
+ *
+ *  Created on: May 31, 2013
+ *      Author: azoghbi
+ */
+
+#ifndef MOD_HPP_
+#define MOD_HPP_
+
+#include "def.hpp"
+#include <specialfunctions.h>
+#include <solvers.h>
+
+class mod {
+
+	void	_CSfq( vec , vector<vec2>& , vector<vec2>& );
+
+public:
+
+	int				n,nfq,npar;
+	double			dt,f1,f2;
+	vec				lc,lce,t,FqL;
+	vector<vec2>	Cfq,Sfq,Cv,Cfq2,Sfq2;
+
+	vec2			C,Ci,I,yyT,yyTmC,dC,C2,C3;
+	vec				Cilc;
+
+	ae_int_t		info; densesolverreport rep;
+
+	mod();
+	virtual ~mod();
+
+	virtual void _C( vec );
+	virtual void _dC( vec , int );
+	virtual void step_pars( int , vec& , vec& );
+	virtual void print_pars( vec& , vec& );
+	virtual void print_pars( vec& , vec& , vec& );
+	virtual void what_pars( int&, int& );
+
+	void	setlc(){ lc.setlength(n); lce.setlength(n); t.setlength(n);}
+	void	init( vec , int );
+	double	loglikelihood( vec );
+	void	dlikelihood( vec , double& , vec& , vec2& );
+	void	optimize( vec& , vec& );
+};
+
+// ----------------------------------------- //
+// ---------- mod container ---------------- //
+// ----------------------------------------- //
+
+template <class M>
+class Mod {
+
+	int			nmod,nfq;
+	vector<M>	mods;
+
+public:
+	int			npar;
+	Mod( vector<lcurve> inlc , vec fqL ){
+		int		i;
+		nmod	=	inlc.size();
+		for( i=0 ; i<nmod ; i++ ){ mods.push_back( M( inlc[i] , fqL ) ); }
+		npar	=	mods[0].npar;
+		nfq		=	mods[0].nfq;
+	}
+	Mod( vector<lcurve> lc1, vector<lcurve> lc2 , vec fqL , vec pars ){
+		int		i;
+		nmod	=	lc1.size();
+		for( i=0 ; i<nmod ; i++ ){ mods.push_back( M( lc1[i] , lc2[i] , fqL , pars ) ); }
+		npar	=	mods[0].npar;
+		nfq		=	mods[0].nfq;
+	}
+	Mod( vector<lcurve> lc1, vector<lcurve> lc2 , vec fqL ){
+		int		i;
+		nmod	=	lc1.size();
+		for( i=0 ; i<nmod ; i++ ){ mods.push_back( M( lc1[i] , lc2[i] , fqL ) ); }
+		npar	=	mods[0].npar;
+		nfq		=	mods[0].nfq;
+	}
+	virtual ~Mod(){}
+
+	double	loglikelihood( vec p ){
+		double l=0; for( int i=0 ; i<nmod ; i++ ){ l += mods[i].loglikelihood(p); }
+		return l;
+	}
+
+	void	dlikelihood( vec p , double& loglike , vec& grad , vec2& hess ){
+		int i,j,nm;vec g; vec2 h; double l;g.setlength(npar);h.setlength(npar,npar);
+		loglike = 0; for( i=0 ; i<npar ; i++ ){ grad[i]=0; for(j=0;j<npar;j++){ hess(i,j) = 0;}}
+
+		for( nm=0 ; nm<nmod ; nm++ ){
+			mods[nm].dlikelihood( p , l , g , h );
+			loglike	+=	l;
+			for( i=0 ; i<npar ; i++ ){ grad[i] += g[i]; for(j=0;j<npar;j++){ hess(i,j) += h(i,j);}}
+		}
+	}
+
+	void	step_pars( int n, vec& dpar , vec& pars ){ mods[0].step_pars( n , dpar , pars );}
+	void	print_pars( vec& pars , vec& errs ){ mods[0].print_pars( pars , errs );}
+	void	print_pars( vec& pars , vec& errs1 , vec& errs2 ){ mods[0].print_pars( pars , errs1 , errs2 );}
+	void	what_pars( int& ip1 , int& ip2 ){ mods[0].what_pars( ip1 , ip2 );}
+
+	double	optimize( vec& pars , vec& errs ){
+		int		nmax = 400;
+		double	tol = 1e-3;
+
+		ae_int_t info; densesolverreport rep;
+		int		i,j,n,sing=0,np;
+		double	loglike,dpmax,dl,prevl=-1e20;
+		vec		tmppars,dpar,grad,prevp; vec2 hess,hessi,ii;
+		tmppars.setlength(npar); dpar.setlength(npar); grad.setlength(npar); prevp.setlength(npar);
+		hess.setlength( npar , npar );ii.setlength(npar,npar);hessi.setlength(npar,npar);
+		for( i=0 ; i<npar ; i++ ){
+			tmppars[i] = pars[i]; dpar[i] = 0.01*pars[i]; ii(i,i)=1;for(j=0;j<i;j++){ii(i,j)=0;ii(j,i)=0;}
+			prevp[i] = pars[i];
+		}
+
+		for( n=1 ; n<= nmax ; n++ ){
+
+			dlikelihood( tmppars , loglike , grad , hess );
+			// ---------- //
+			if( loglike != loglike ){ for( i=0 ; i<npar ; i++ ){ dpar[i]/=10;tmppars[i] = prevp[i] + dpar[i]; } continue;}
+			dl = loglike - prevl;
+			// ---------- //
+			vector<int> ic,iv; sing = 0;
+			for( i=0; i <npar ; i++ ){ if( fabs(grad[i]) < 1e-5 ){ sing=1;ic.push_back(i); }else { iv.push_back(i);}}
+			if( sing == 1 and iv.size()!=0 ){
+				np = iv.size();
+				vec	g; vec2 h,hi,iii; g.setlength(np); h.setlength(np,np); hi.setlength(np,np);iii.setlength(np,np);
+				for( i=0 ; i<np ; i++ ){
+					iii(i,i) = 1; h(i,i) = hess(iv[i],iv[i]);
+					g[i] = grad[iv[i]]; for(j=0;j<i;j++){ h(i,j) = hess(iv[i],iv[j]);h(j,i)=h(i,j);iii(i,j)=0;iii(j,i)=0;}
+				}
+				spdmatrixcholesky( h , np , false );
+				spdmatrixcholeskysolvem( h , np , false , iii , np , info , rep , hi );
+
+				dpmax = -1e20;
+				for( i=0 ; i<np ; i++ ){
+					dpar[iv[i]] = 0; for( j=0 ; j<np ; j++ ){ dpar[iv[i]] += g[j] * hi(i,j);}
+					dpmax	=	max( dpmax , fabs(dpar[iv[i]]) );
+				}
+				for( i=0 ; i<int(ic.size()); i++ ){ dpar[ic[i]] = 0;}
+			}else {
+			// --------- //
+
+				spdmatrixcholesky( hess , npar , false );
+				spdmatrixcholeskysolvem( hess , npar , false , ii , npar , info , rep , hessi );
+
+				dpmax = -1e20;
+				for( i=0 ; i<npar ; i++ ){
+					dpar[i] = 0; for( j=0 ; j<npar ; j++ ){ dpar[i] += grad[j] * hessi(i,j);}
+					dpmax	=	max( dpmax , fabs(dpar[i]) );
+				}
+			}
+			for(i=0;i<npar;i++){prevp[i] = tmppars[i];}
+			prevl = loglike;
+			if ( dl < 2 and n>10 ){for(i=0;i<npar;i++){dpar[i]/=2;}}
+			step_pars( n , dpar , tmppars );
+
+			cout << n << " " << dpmax << " " << loglike << " "; for(j=0;j<npar;j++){cout << tmppars[j] << " ";} cout << endl;
+			if ( dpmax<tol ) break;
+		}
+		for( i=0 ; i<npar ; i++ ){
+			pars[i] = tmppars[i];
+			errs[i]	= sqrt(hessi(i,i));
+		}
+		cout << setfill('-') << setw(40) << "\n";
+		for(i=0;i<npar;i++){cout << pars[i] << " ";} cout << endl;
+		for(i=0;i<npar;i++){cout << errs[i] << " ";} cout << endl;
+		for(i=0;i<npar;i++){cout << grad[i] << " ";} cout << endl;
+		cout << setfill('-') << setw(40) << "\n" << setfill(' ');
+		print_pars( pars, errs );
+		return loglike;
+	}
+
+	double	opti( vec pars , vec errs ,int k ){
+		int		nmax = 400;
+		double	tol = 1e-3;
+
+		ae_int_t info; densesolverreport rep;
+		int		i,j,n,np;
+		double	loglike,dpmax,dl,prevl=-1e20;
+		vec		tmppars,dpar,grad,prevp; vec2 hess,hessi,ii;
+		tmppars.setlength(npar); dpar.setlength(npar); grad.setlength(npar); prevp.setlength(npar);
+		hess.setlength( npar , npar );ii.setlength(npar,npar);hessi.setlength(npar,npar);
+		vector<int> ic,iv; np = npar - 1;
+		vec	g; vec2 h,hi,iii; g.setlength(np); h.setlength(np,np); hi.setlength(np,np);iii.setlength(np,np);
+		ic.push_back(k);
+
+		for( i=0 ; i<npar ; i++ ){
+			tmppars[i] = pars[i]; dpar[i] = 0.01*pars[i]; ii(i,i)=1;for(j=0;j<i;j++){ii(i,j)=0;ii(j,i)=0;}
+			prevp[i] = pars[i];
+			if(i==k){continue;} iv.push_back(i);
+		}
+
+
+		for( n=1 ; n<= nmax ; n++ ){
+
+			dlikelihood( tmppars , loglike , grad , hess );
+			// ---------- //
+			if( loglike != loglike ){ for( i=0 ; i<npar ; i++ ){ tmppars[i] = prevp[i] + dpar[i]/10; } continue;}
+			dl = loglike - prevl;
+			// ---------- //
+
+
+			for( i=0 ; i<np ; i++ ){
+				iii(i,i) = 1; h(i,i) = hess(iv[i],iv[i]);
+				g[i] = grad[iv[i]]; for(j=0;j<i;j++){ h(i,j) = hess(iv[i],iv[j]);h(j,i)=h(i,j);iii(i,j)=0;iii(j,i)=0;}
+			}
+			spdmatrixcholesky( h , np , false );
+			spdmatrixcholeskysolvem( h , np , false , iii , np , info , rep , hi );
+
+			dpmax = -1e20;
+			for( i=0 ; i<np ; i++ ){
+				dpar[iv[i]] = 0; for( j=0 ; j<np ; j++ ){ dpar[iv[i]] += g[j] * hi(i,j);}
+				dpmax	=	max( dpmax , fabs(dpar[iv[i]]) );
+			}
+			for( i=0 ; i<int(ic.size()); i++ ){ dpar[ic[i]] = 0;}
+
+			for(i=0;i<npar;i++){prevp[i] = tmppars[i];}
+			prevl = loglike;
+			if ( dl < 2 and n>10 ){for(i=0;i<npar;i++){dpar[i]/=2;}}
+			step_pars( n , dpar , tmppars );
+
+			cout << n << " " << dpmax << " " << loglike << " "; for(j=0;j<npar;j++){cout << tmppars[j] << " ";} cout << endl;
+			if ( dpmax<tol ) break;
+		}
+		return loglike;
+	}
+
+	void	errors_avg( vec& pars , vec& errs ){
+		vec errs1,errs2;errs1.setlength(npar);errs2.setlength(npar);
+		errors( pars , errs1 , 1 );
+		errors( pars , errs2 , -1 );
+		for(int i=0;i<npar;i++){errs[i] = (errs1[i]+errs2[i])/2;}
+	}
+
+	void	errors( vec& pars , vec& errs1 , vec& errs2 ){
+		errors( pars , errs1 , 1 );
+		errors( pars , errs2 , -1 );
+	}
+
+	void	errors( vec& pars , vec& errs , int sign=-1){
+		double		tol = 1e-3 , delchi = 1;
+		int			i,ip,ip1,ip2;
+		double		pu,pd,phalf=0,bestlike,tmplike,dl;
+		vec			pars0,errs0,tmpp; pars0.setlength(npar);errs0.setlength(npar);tmpp.setlength(npar);
+		bestlike	=	optimize( pars , errs );
+		for(i=0;i<npar;i++){pars0[i] = pars[i];errs0[i]=errs[i];}
+
+		what_pars( ip1 , ip2 );
+
+		for( ip=ip1 ; ip<ip2 ; ip++ ){
+			cout << endl << "******* par " << ip << " *******" << endl << endl;
+			pd	=	pars[ip];
+			for( i=0 ;i<npar ; i++ ){tmpp[i] = pars0[i]; }
+
+			i=2; dl = 0;
+			while( dl < delchi ){
+				pu			=	pars[ip] + sign*i*errs0[ip];
+				tmpp[ip]	=	pu;
+				tmplike		=	opti( tmpp , errs , ip );
+				dl			=	2*(bestlike - tmplike);
+				i++;
+			}
+			while( fabs(dl-delchi)>tol ){
+				phalf		=	(pu+pd)/2.0;
+				tmpp[ip]	=	phalf;
+				tmplike		=	opti( tmpp , errs , ip );
+				dl			=	2*(bestlike - tmplike);
+				cout << "+++ " << bestlike << " " << tmplike << " " << pars[ip] << " "<< phalf << " " << dl << endl;
+				if( dl<delchi ){ pd = phalf;}else{ pu = phalf;}
+			}
+			errs0[ip]		=	sign* ( phalf - pars0[ip] );
+		}
+		for( i=0 ; i<npar ; i++ ){ pars[i] = pars0[i]; errs[i] = errs0[i];}
+		print_pars( pars0 , errs0 );
+		//optimize( pars0 , errs );
+	}
+};
+
+#endif /* MOD_HPP_ */

psdlag/src/inc/psd.hpp

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

psdlag/src/inc/psdlag.hpp

@@ -0,0 +1,54 @@
+/*
+ * psdlag.hpp
+ *
+ *  Created on: Jun 1, 2013
+ *      Author: azoghbi
+ */
+
+#include "mod.hpp"
+
+#ifndef PSDLAG_HPP_
+#define PSDLAG_HPP_
+
+class psdlag : public mod {
+
+	int		n1;
+	ivec	ip1,ip2,icx,iphi;
+
+	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( lcurve , lcurve , vec );
+	virtual ~psdlag();
+
+	void	step_pars( int , vec& , vec& );
+	void	print_pars( vec& , vec& );
+};
+
+class psdlag10 : public mod {
+
+	int		n1;
+	ivec	ip1,ip2,icx,iphi;
+
+	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( lcurve , lcurve , vec );
+	virtual ~psdlag10();
+
+	void	step_pars( int , vec& , vec& );
+	void	print_pars( vec& , vec& );
+	void	what_pars( int&, int& );
+};
+
+#endif /* PSDLAG_HPP_ */

psdlag/src/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
+

psdlag/src/lag.cpp

@@ -0,0 +1,227 @@
+/*
+ * lag.cpp
+ *
+ *  Created on: Jun 1, 2013
+ *      Author: azoghbi
+ */
+
+#include "inc/lag.hpp"
+
+lag::lag( lcurve lc1 , lcurve lc2 , vec fqL , vec pars ) {
+
+	// ----------- initial parameters ------------ //
+	n1		=	lc1.len;
+	n		=	n1 + lc2.len;
+	dt		=	lc1.dt;
+	// ------------------------------------------ //
+
+
+	// ----------- light curve setup ------------ //
+	setlc();
+	lc1.demean(); lc2.demean();
+	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];
+		p2[i]	=	pars[i+nfq];
+		icx[i]	=	i;
+		iphi[i]	=	i+nfq;
+	}
+	// ------------------------------------------ //
+
+}
+
+lag::~lag() {}
+
+void lag::_C( vec p ){
+	int	i,j,k; double d;
+	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::_dC( vec p , int k ){
+	if( k<nfq ){	_dCx( p , k );
+	}else{			_dPhi( p , k-nfq );}
+}
+
+void lag::_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::_dPhi( vec p , int k ){
+	int	i,j; double cx=p[k],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::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::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::lag10( lcurve lc1 , lcurve lc2 , vec fqL , vec pars ) {
+
+	// ----------- initial parameters ------------ //
+	n1		=	lc1.len;
+	n		=	n1 + lc2.len;
+	dt		=	lc1.dt;
+	// ------------------------------------------ //
+
+
+	// ----------- light curve setup ------------ //
+	setlc();
+	lc1.demean(); lc2.demean();
+	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]);
+		p2[i]	=	pow(10,pars[i+nfq]);
+		icx[i]	=	i;
+		iphi[i]	=	i+nfq;
+	}
+	// ------------------------------------------ //
+
+}
+
+lag10::~lag10() {}
+
+void lag10::_C( vec p ){
+	int	i,j,k; double d; for(i=0;i<nfq;i++){p[i] = pow(10,p[i]);}
+	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::_dC( vec p , int k ){
+	if( k<nfq ){	_dCx( p , k );
+	}else{			_dPhi( p , k-nfq );}
+}
+
+void lag10::_dCx( vec p , int k ){
+	int	i,j; double phi=p[iphi[k]],cx = log(10)*pow(10,p[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)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
+	}
+}
+
+void lag10::_dPhi( vec p , int k ){
+	int	i,j; double cx=pow(10,p[k]),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::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::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::what_pars( int& ip1 , int& ip2 ){
+	ip1 = nfq; ip2 = npar;
+}

psdlag/src/main.cpp

@@ -0,0 +1,203 @@
+/*
+ * main.cpp
+ *
+ *  Created on: May 31, 2013
+ *      Author: azoghbi
+ */
+
+#include "inc/main.hpp"
+
+/**************************************
+ * Main function
+ * Checks if there is an input and
+ * raises an error.
+ * otherwise, call do_work(fname)
+ **************************************/
+int main( int argc , char* argv[] ){
+	if ( argc < 2 ){
+		cerr << "** Error ** No input file given." << endl;
+		usage();
+	}else {
+		do_work( argv[1] );
+	}
+	return 0;
+}
+
+
+
+
+/**************************************
+ * Main function that does the work
+ * INPUT:	char*	fname: input file name
+ * RETURN:	void
+ **************************************/
+
+void do_work( char* fname ){
+
+	/************   Reading the input file  ***********/
+	ifstream fp(fname);
+	string line,mcmcfile; stringstream ss;
+	int i,nfiles,nfq,mode,npar,bin1,bin2,fit_type,nrun,nburn,nwk; bool strict;
+	getline(fp,line);ss.str(line); ss >> nfiles;								ss.clear();
+	vector<string> files(nfiles);ivec secL;secL.setlength(nfiles);for(i=0;i<nfiles;i++){
+	getline(fp,line); ss.str(line); ss >> files[i] >> secL[i];					ss.clear();}
+	getline(fp,line); ss.str(line); ss >> nfq;vec fqL;fqL.setlength(nfq+1);
+	for(i=0;i<=nfq;i++){ss >> fqL[i];}											ss.clear();
+	getline(fp,line); ss.str(line); ss >> mode;									ss.clear();
+	npar = (mode==0)?4*nfq:nfq; vec pars; pars.setlength(npar);
+	getline(fp,line); ss.str(line); for(i=0;i<npar;i++){ss >> pars[i];}			ss.clear();
+	getline(fp,line); ss.str(line); ss >> bin1 >> bin2;							ss.clear();
+	if(mode!=0){bin1 = (mode==-1)?-1:(mode-1);bin2=-1;}
+	getline(fp,line); ss.str(line); ss >> fit_type;								ss.clear();
+	getline(fp,line); ss.str(line); ss >> i; strict	= i!=0;						ss.clear();
+	getline(fp,line); ss.str(line); ss >> nrun >> nburn >> nwk >> mcmcfile;		ss.clear();
+	fp.close();
+	/*********   END Reading the input file  **********/
+
+	/**********  Read the light curves  **************/
+	vector<vector<lcurve> > LC;
+	for(i=0;i<nfiles;i++) readLC( LC , files[i] , secL[i] , bin1 , bin2 , strict );
+	/********  END Read the light curves  ************/
+
+	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]); }
+		Mod<psd10>	p1( lc1 , fqL );
+
+		p1.optimize( pars , errs );
+		if( fit_type==1 ) p1.errors( pars , errs );
+	}else if( mode == 0 or mode==-1 ){
+		vec ec,e1,e2,errs; ec.setlength(2*nfq);e1.setlength(nfq);e2.setlength(nfq);errs.setlength(4*nfq);
+		vector<lcurve> lc1,lc2; for( i=0 ; i<int(LC.size()) ; i++ ){ lc1.push_back( LC[i][0]); lc2.push_back( LC[i][1]); }
+		Mod<psd10>	p1( lc1 , fqL ), p2( lc2 , fqL );
+		vec		pars1,pars2;pars1.setlength(nfq);pars2.setlength(nfq);for(i=0;i<nfq;i++){pars1[i]=pars[i];pars2[i]=pars[i+nfq];}
+
+		p1.optimize( pars1 , e1 );
+		p2.optimize( pars2 , e2 );
+		if( fit_type==1 ) {
+			p1.errors_avg( pars1 , e1 );
+			p2.errors_avg( pars2 , e2 );
+		}
+
+		vec		pc; pc.setlength(2*nfq);
+		for(i=0;i<nfq;i++){pc[i]=pars[i+2*nfq];pc[i+nfq]=pars[i+3*nfq];pars[i]=pars1[i];pars[i+nfq]=pars2[i];pc[i]=(pars1[i]+pars2[i])/2;}
+		Mod<lag10>		l( lc1 , lc2 , fqL , pars );
+		l.optimize( pc , ec );
+		if( fit_type==1 ) l.errors_avg( pc , ec );
+		if( fit_type==2 ) {
+			mcmc mc( 2*nfq , mcmc_lag10 , (void*)&l );
+			mc.nrun=nrun; mc.nburn=nburn; mc.nwk=nwk;
+			mc.run( pc , errs , mcmcfile.c_str() );
+		}
+
+		for(i=0;i<nfq;i++){ pars[i+2*nfq]=pc[i]; pars[i+3*nfq]=pc[i+nfq];}
+		for(i=0;i<nfq;i++){ errs[i] = e1[i]; errs[i+nfq] = e2[i]; errs[i+2*nfq]=ec[i]; errs[i+3*nfq]=ec[i+nfq];}
+
+		Mod<psdlag10>	pl( lc1 , lc2 , fqL );
+		if( fit_type==3 ) pl.errors_avg( pars , errs );
+		if( fit_type==4 ) {
+			mcmc mc( 4*nfq , mcmc_psdlag10 , (void*)&pl );
+			mc.nrun=nrun; mc.nburn=nburn; mc.nwk=nwk;
+			mc.run( pars , errs , mcmcfile.c_str() );
+		}
+
+		pl.print_pars( pars , errs );
+	}
+}
+
+
+
+
+
+
+
+/**************************************
+ * Usage function
+ * prints a sample input text
+ * INPUT:	none
+ * RETURN:	void
+ **************************************/
+void usage(){
+	cerr << setw(40) << setfill('*') << "\n";
+	cout << setw(30) << left << "1"							<< "# of lc"			<< endl;
+	cout << setw(30) << left << "lc.dat 0"					<< "name1 secL1"		<< endl;
+	cout << setw(30) << left << "3 1e-5 1e-4 5e-4 1e-3"		<< "nfq fqL1 fqL2.." 	<< endl;
+	cout << setw(30) << left << "1"							<< "mode? 0:lag, n:psd" << endl;
+	cout << setw(30) << left << "1 1 1 1 1 1 1 1 1 1 1 1"	<< "init. pars" 		<< endl;
+	cout << setw(30) << left << "0 1"						<< "bin1, bin2,-1: all" << endl;
+	cout << setw(30) << left << "0"							<< "fit_type" 			<< endl;
+	cout << setw(30) << left << "0"							<< "strict" 			<< endl;
+	cout << setw(30) << left << "100 50 50 mcmc.dat"		<< "nrun,nburn,nwk,file"<< endl;
+	cerr << setw(40) << setfill('*') << "\n";
+}
+
+
+
+
+
+/***************************************
+ * READ LIGHTCURVES
+ * Read a standard light curve file that has a description at the top
+ * the first line is something like: # dT nlc b0 b1 b2 ... b_lc (b0 .. are bins boundaries)
+ * Then the content with: time lc1 err1 lc2 err2 ... etc
+ *
+ * INPUT:
+ * 		LC		: a vector of vectors of light curves, where the first dimension is the number of
+ * 				light curves (i.e. segments), this is incremented in this function when reading a new file.
+ * 				The other dimension is either 1 for psd and 2 for lags, which contains either one light curve
+ * 				or two for the lag. Passed by reference
+ * 		fname	: name of the light curve file
+ * 		secL	: the section length to split the light curves to.
+ * 		b1		: the bin1 to read for the first light curve. 0-based, -1 for all
+ * 		b2		: the bin2 to read, 0-based, -1 for none (.i.e just psd)
+ * 		strict	: a boolean on whether to be strict or not in segmenting the light curve
+ * 				e.g. for secL=220: 0 gives 100,100 and 1 gives 100,120
+ *
+ **************************************/
+void readLC( vector<vector<lcurve> >&LC , string fname , int secL , int b1 , int b2 , bool strict ){
+
+	/*		Initialize file stream, define some variables		*/
+	ifstream fp(fname.c_str()); string line,sdum; stringstream ss;
+	int i,j,nlc,n,nsec,ns,sl; double dt;
+
+	/*		read dt from the description line 					*/
+	getline(fp,line); ss.str(line); ss >> sdum >> dt >> nlc; ss.clear();
+
+	/*		get the number of lines 							*/
+	n=0;while(getline(fp,line)){if(line!="")n++;} fp.clear(); fp.seekg(0,ios::beg);getline(fp,line);
+
+	/*		calculate section lengths							*/
+	if(secL==0){secL=n;}
+	nsec = n/secL; vector<int> secl(nsec,secL); if(not strict) secl[nsec-1] += n-nsec*secL;
+
+	/*		Start looping through the number of sections		*/
+	for(ns=0;ns<nsec;ns++){
+		sl	=	secl[ns]; vector<lcurve> Lc;
+		vec t; vector<vec> lc(nlc),lce(nlc); t.setlength( sl ); for(i=0;i<nlc;i++){lc[i].setlength(sl);lce[i].setlength(sl);}
+
+		// Read the first chunk
+		for( i=0 ; i<sl ; i++ ){
+			getline(fp,line);ss.str(line);ss>>t[i];for(j=0;j<nlc;j++){ss>>lc[j][i]>>lce[j][i];} ss.clear();
+		}
+
+		// if we want the whole light curves in bin1
+		if( b1==-1 ){
+			vec lc_,lce_; lc_.setlength( sl );lce_.setlength( sl );for(i=0;i<sl;i++){lc_[i]=0;lce_[i]=0;}
+			for( i=0 ; i<sl ; i++ ){
+				for( j=0 ; j<nlc ; j++ ){ if(j==b2){continue;} lc_[i] += lc[j][i]; lce_[i] += lce[j][i]*lce[j][i]; }
+				lce_[i] = sqrt(lce_[i]);
+			}
+			Lc.push_back(lcurve(t,lc_,lce_,dt));
+
+		// else: i.e. light curve 1 is bin1
+		}else {
+			Lc.push_back(lcurve(t,lc[b1],lce[b1],dt));
+		}
+		// Do we need a second light curve?
+		if(b2!=-1){Lc.push_back(lcurve(t,lc[b2],lce[b2],dt));}
+
+		// push the resulting vector<lcurve> to LC
+		LC.push_back(Lc);
+	}
+	fp.close();
+}

psdlag/src/makefile

@@ -0,0 +1,7 @@
+
+
+incdir='/home/caes/science/fourier-agn/psdlag/src/'
+
+
+psdlag:
+	g++ *cpp -o psdlag -O3 -Wall  -I${incdir}

psdlag/src/mcmc.cpp

@@ -0,0 +1,73 @@
+/*
+ * mcmc.cpp
+ *
+ *  Created on: May 14, 2013
+ *      Author: abduz
+ */
+
+#include "inc/mcmc.hpp"
+
+mcmc::mcmc(int np,double (*f)(vec&,void*),void *p) {
+	npar	=	np;
+	loglikelihood	=	f;
+	ptr		=	p;
+	hqrndrandomize(rnd);
+	nrun = 100; nburn = 100; nwk = 100; ncheck = 10;avalue = 2.0;
+}
+
+mcmc::~mcmc() {}
+
+void mcmc::run( vec& p, vec& pe, const char*fname){
+	int		i,j,k,nacc,nr,n,nsum;
+	vec2	pars; pars.setlength(nwk,npar);
+	vec		prob,currp,bestp,sum,sum2;
+	ivec2	compj;
+	double	bestlike,tmplike=-1e20,z,q,a1=1./sqrt(avalue),a2=sqrt(avalue);;
+	prob.setlength(nwk);currp.setlength(npar);bestp.setlength(npar);sum.setlength(npar);sum2.setlength(npar);
+	compj.setlength(nwk,nwk-1);
+	nsum=0;for(i=0;i<npar;i++){sum[i]=0;sum2[i]=0;}
+	ofstream fout(fname);
+
+	for(k=0;k<nwk;k++){j=0;for(i=0;i<k;i++){compj[k][j]=i;j++;}for(i=k+1;i<nwk;i++){compj[k][j]=i;j++;}}
+
+	bestlike	=	-1e20;
+	for(k=0;k<nwk;k++){
+		//pars[k][0] = p[0];currp[0] = p[0];
+		for(i=0;i<npar;i++){pars[k][i] = p[i]*(1.0+hqrndnormal(rnd)*0.05); currp[i]=pars[k][i];}
+		prob[k]	=	loglikelihood( currp , ptr );
+		if(bestlike<prob[k]){bestlike=prob[k];for(i=0;i<npar;i++){bestp[i]=pars[k][i];}}
+	}
+
+	nr = 0;nacc = 0;
+	for(n=1;n<=(nrun+nburn);n++){
+		for(k=0;k<nwk;k++){
+			nr++;j = compj[k][hqrnduniformi(rnd,nwk-1)];
+			z	=	hqrnduniformr(rnd)*(a2-a1) + a1; z = z*z;
+			for(i=0;i<npar;i++){currp[i]=pars[j][i] + z*(pars[k][i]-pars[j][i]);}
+			tmplike	=	loglikelihood(currp,ptr);
+			q = pow(z,npar-1)*exp(tmplike-prob[k]);
+			if( hqrnduniformr(rnd)<=q and tmplike>-1e10 ){
+				nacc++;for(i=0;i<npar;i++){pars[k][i]=currp[i];} prob[k] = tmplike;
+				if(bestlike<tmplike){bestlike=tmplike;for(i=0;i<npar;i++){bestp[i]=currp[i];}}
+				if(n>nburn) {
+					for(i=0;i<npar;i++){
+						fout << currp[i] << "  ";
+						sum[i] += currp[i];sum2[i] += currp[i]*currp[i]; nsum++;
+					} fout << tmplike << endl;
+				}
+			}
+		}
+		if(n%ncheck==0){
+			q = 1.0*nacc/nr; nr=0;nacc=0;
+			cout << setw(8) << n << setw(12) << q << setw(12) << tmplike << setw(12) << bestlike << endl;
+		}
+	}
+	for(i=0;i<npar;i++){
+		p[i]=bestp[i];
+		pe[i] = sqrt((sum2[i] - sum[i]*sum[i]/nsum)/(nsum+1));
+	}cout << endl;
+
+
+	fout.close();
+
+}

psdlag/src/mod.cpp

@@ -0,0 +1,235 @@
+/*
+ * mod.cpp
+ *
+ *  Created on: May 31, 2013
+ *      Author: azoghbi
+ */
+
+#include "inc/mod.hpp"
+
+mod::mod() {
+	n = 0 ; nfq = 0 ; npar = 0; dt = 0; info = 0;f1=0;f2=0;
+
+}
+
+mod::~mod() {}
+
+void mod::init( vec fqL , int fac ){
+
+	// --------------- Freq stuff----------------- //
+	int		i,j;
+	nfq		=	fqL.length() - 1;
+	npar	=	nfq * fac ;
+	FqL.setlength(nfq+1); for( i=0 ; i<=nfq ; i++ ){ FqL[i] = fqL[i]; }
+	// ------------------------------------------ //
+
+
+	// ---------------- Integrals ---------------- //
+	Cfq.resize( nfq );Sfq.resize( nfq );
+	for( i=0 ; i<nfq ; i++ ){ Cfq[i].setlength( n , n ); Sfq[i].setlength( n , n );}
+	_CSfq( fqL , Cfq , Sfq );
+
+	Cfq2.resize( 3 );Sfq2.resize( 3 );
+	for( i=0 ; i<3 ; i++ ){ Cfq2[i].setlength( n , n ); Sfq2[i].setlength( n , n );}
+	vec		fqL2;fqL2.setlength(4); fqL2[0] = 1e-1*fqL[0]; fqL2[1] = fqL[0]; fqL2[2] = fqL[nfq]; fqL2[3] = 10*fqL[nfq];
+	f1 = 2; f2 = 0.5;
+	_CSfq( fqL2 , Cfq2 , Sfq2 );
+	// ------------------------------------------ //
+
+
+	// -------------- Containers ---------------- //
+	C.setlength( n , n );
+	Ci.setlength( n , n );
+	I.setlength( n , n );
+	for( i=0 ; i<n ; i++ ){ I(i,i) = 1.0; for( j=0 ; j<i ; j++ ){ I(i,j)=0; I(j,i)=0;}}
+	Cilc.setlength( n );
+	yyT.setlength( n , n );
+	for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ yyT(i,j) = lc[i]*lc[j];}}
+	yyTmC.setlength( n , n );
+	dC.setlength( n , n );
+	Cv.resize( npar ); for( i=0 ; i<npar ; i++ ){ Cv[i].setlength( n , n ); }
+	C2.setlength( n , n );
+	C3.setlength( n , n );
+
+	// ------------------------------------------ //
+
+}
+
+void mod::_CSfq( vec fqL , vector<vec2>&cfq , vector<vec2>& sfq ){
+	//fqL[0] *= 1e-6; fqL[nfq]*=2;
+	int		i,j,k, nfq=fqL.length()-1;
+	double	tt,norm,t1,t2,pi=M_PI,pidt=pi*dt,sdt2,sdt1,cdt2,cdt1,s2dt2,s2dt1,c2dt2,c2dt1,dtmt,dtpt;
+	double	st2,st1,ct2,ct1, sm2,sm1,cm2,cm1,  sp2,sp1,cp2,cp1;
+	vec		w; w.setlength(nfq+1); for( k=0 ; k<=nfq ; k++ ){ w[k] = 2*pi*fqL[k]; }
+	norm	=	1./( pidt*pidt );
+
+	for( i=0 ; i<n ; i++ ){ for( j=0 ; j<n ; j++ ){
+		tt	=	t[j] - t[i];
+		dtmt = dt - tt; dtpt = dt+tt;
+		for( k=0 ; k<nfq ; k++ ){
+			if( tt==0 ){
+				sinecosineintegrals( w[k+1]*dt , sdt2 , cdt2 );
+				sinecosineintegrals( w[k]  *dt , sdt1 , cdt1 );
+				t1	=	( cos( w[k+1]*dt ) - 1 ) / ( 2 * fqL[k+1] ) + pidt * sdt2;
+				t1	-=	( cos( w[k]  *dt ) - 1 ) / ( 2 * fqL[k]   ) + pidt * sdt1;
+				cfq[k](i,j)		=	norm * t1;
+				sfq[k](i,j)		=	0;
+			}else if( fabs(tt)==dt ){
+				sinecosineintegrals( w[k+1]*dt , sdt2 , cdt2 );
+				sinecosineintegrals( w[k]  *dt , sdt1 , cdt1 );
+				sinecosineintegrals( 2*w[k+1]*dt , s2dt2 , c2dt2 );
+				sinecosineintegrals( 2*w[k]  *dt , s2dt1 , c2dt1 );
+				t1	=	-( cos( w[k+1]*dt )*pow( sin( 0.5*w[k+1]*dt ) , 2 ) ) / fqL[k+1];
+				t1	-=	-( cos( w[k]  *dt )*pow( sin( 0.5*w[k]  *dt ) , 2 ) ) / fqL[k]  ;
+				t2	=	pidt * ( s2dt2 - sdt2 );
+				t2	-=	pidt * ( s2dt1 - sdt1 );
+				cfq[k](i,j)		=	norm * ( t1 + t2 );
+
+				t1	=	-( 2*cos( 0.5*w[k+1]*dt )*pow( sin( 0.5*w[k+1]*dt ) , 3 ) ) / fqL[k+1];
+				t1	-=	-( 2*cos( 0.5*w[k]  *dt )*pow( sin( 0.5*w[k]  *dt ) , 3 ) ) / fqL[k]  ;
+				t2	=	pidt * ( cdt2 - c2dt2 );
+				t2	-=	pidt * ( cdt1 - c2dt1 );
+				sfq[k](i,j)		=	norm * fabs(tt)/tt * ( t1 + t2 );
+			}else{
+				sinecosineintegrals( w[k+1]*tt , st2 , ct2 );
+				sinecosineintegrals( w[k]  *tt , st1 , ct1 );
+				sinecosineintegrals( w[k+1]*dtmt , sm2 , cm2 );
+				sinecosineintegrals( w[k]  *dtmt , sm1 , cm1 );
+				sinecosineintegrals( w[k+1]*dtpt , sp2 , cp2 );
+				sinecosineintegrals( w[k]  *dtpt , sp1 , cp1 );
+				t1	=	( cos( w[k+1]*dtmt ) - 2*cos( w[k+1]*tt ) + cos( w[k+1]*dtpt ) )/(4*fqL[k+1]);
+				t1	-=	( cos( w[k]  *dtmt ) - 2*cos( w[k]  *tt ) + cos( w[k]  *dtpt ) )/(4*fqL[k]  );
+				t2	=	( dtmt*sm2 - 2*tt*st2 + dtpt*sp2 ) * (pi/2);
+				t2	-=	( dtmt*sm1 - 2*tt*st1 + dtpt*sp1 ) * (pi/2);
+				cfq[k](i,j)		=	norm * ( t1 + t2 );
+
+				t1	=	( -sin( w[k+1]*dtmt ) - 2*sin( w[k+1]*tt ) + sin( w[k+1]*dtpt ) )/(4*fqL[k+1]);
+				t1	-=	( -sin( w[k]  *dtmt ) - 2*sin( w[k]  *tt ) + sin( w[k]  *dtpt ) )/(4*fqL[k]  );
+				t2	=	( dtmt*cm2 + 2*tt*ct2 - dtpt*cp2 ) * (pi/2);
+				t2	-=	( dtmt*cm1 + 2*tt*ct1 - dtpt*cp1 ) * (pi/2);
+				sfq[k](i,j)		=	norm * ( t1 + t2 );
+			}
+		}
+
+	}}
+}
+
+
+void mod::_C( vec ){
+	cerr << endl << "mod::_C is an empty function. It should not be called directly" << endl << endl;;
+	exit(1);
+}
+
+void mod::_dC( vec , int ){
+	cerr << endl << "mod::_C is an empty function. It should not be called directly" << endl << endl;;
+	exit(1);
+}
+
+double mod::loglikelihood( vec p ){
+	double		logdet , chi2 , loglike ; int i;
+	_C( p );
+	if (not spdmatrixcholesky( C , n , false )) {return -1e20;};
+	logdet = 0 ; for( i=0 ; i<n ; i++ ){ logdet += log( C(i,i) ); } logdet *= 2;
+	spdmatrixcholeskysolvem( C , n , false , I , n , info , rep , Ci );
+	rmatrixmv( n , n , Ci , 0 , 0 , 0 , lc , 0 , Cilc , 0 );
+	chi2 = 0 ; for( i=0 ; i<n ; i++ ){ chi2 += Cilc[i]*lc[i] ; }
+	loglike	=	-0.5 * ( chi2 + logdet + n*log(2*M_PI) );
+	return loglike;
+}
+
+void mod::dlikelihood( vec p , double& loglike , vec& grad , vec2& hess ){
+	double		logdet , chi2 , res; int i,j,k;
+	_C( p );
+	for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ yyTmC(i,j) = yyT(i,j) - C(i,j); yyTmC(j,i) = yyTmC(i,j);}}
+	spdmatrixcholesky( C , n , false );
+	logdet = 0 ; for( i=0 ; i<n ; i++ ){ logdet += log( C(i,i) ); } logdet *= 2;
+	spdmatrixcholeskysolvem( C , n , false , I , n , info , rep , Ci );
+	rmatrixmv( n , n , Ci , 0 , 0 , 0 , lc , 0 , Cilc , 0 );
+	chi2 = 0 ; for( i=0 ; i<n ; i++ ){ chi2 += Cilc[i]*lc[i] ; }
+	loglike	=	-0.5 * ( chi2 + logdet + n*log(2*M_PI) );
+
+
+	for( i=0 ; i<npar ; i++ ){
+		_dC( p , i );
+		rmatrixgemm( n , n , n , 1.0 , Ci , 0 , 0 , 0 , dC , 0 , 0 , 0 , 0.0 , Cv[i] , 0 , 0 );
+	}
+
+	for( i=0 ; i<npar ; i++ ){
+		rmatrixgemm( n , n , n , 1.0 , Cv[i] , 0 , 0 , 0 , Ci , 0 , 0 , 0 , 0.0 , C2 , 0 , 0 );
+		rmatrixgemm( n , n , n , 1.0 , yyTmC , 0 , 0 , 0 , C2 , 0 , 0 , 0 , 0.0 , C3 , 0 , 0 );
+		res = 0; for( k=0 ; k<n ; k++ ){ res += C3(k,k); }
+		grad[i]		=	0.5 * res;
+		for( j=0 ; j<=i ; j++ ){
+			rmatrixgemm( n , n , n , 1.0 , Cv[i] , 0 , 0 , 0 , Cv[j] , 0 , 0 , 0 , 0.0 , C3 , 0 , 0 );
+			res = 0; for( k=0 ; k<n ; k++ ){ res += C3(k,k); }
+			hess(i,j)	=	0.5 * res; hess(j,i) = 0.5*res;
+		}
+	}
+}
+
+
+void mod::optimize( vec& pars , vec& errs ){
+	int		nmax = 200;
+	double	tol = 1e-3;
+
+	int		i,j,n;
+	double	loglike,dpmax;
+	vec		tmppars,dpar,grad; vec2 hess,hessi,ii;
+	tmppars.setlength(npar); dpar.setlength(npar); grad.setlength(npar);
+	hess.setlength( npar , npar );ii.setlength(npar,npar);hessi.setlength(npar,npar);
+	for( i=0 ; i<npar ; i++ ){
+		tmppars[i] = pars[i]; dpar[i] = 0.01*pars[i]; ii(i,i)=1;for(j=0;j<i;j++){ii(i,j)=0;ii(j,i)=0;}
+	}
+
+	for( n=1 ; n<= nmax ; n++ ){
+
+		dlikelihood( tmppars , loglike , grad , hess );
+		spdmatrixcholesky( hess , npar , false );
+		spdmatrixcholeskysolvem( hess , npar , false , ii , npar , info , rep , hessi );
+
+		dpmax = -1e20;
+		for( i=0 ; i<npar ; i++ ){
+			dpar[i] = 0; for( j=0 ; j<npar ; j++ ){ dpar[i] += grad[j] * hessi(i,j);}
+			dpmax	=	max( dpmax , fabs(dpar[i]) );
+		}
+		step_pars( n , dpar , tmppars );
+
+		cout << n << " " << dpmax << " " << loglike << " "; for(j=0;j<npar;j++){cout << tmppars[j] << " ";} cout << endl;
+		if ( dpmax<tol ) break;
+	}
+	for( i=0 ; i<npar ; i++ ){
+		pars[i] = tmppars[i];
+		errs[i]	= sqrt(hessi(i,i));
+	}
+	cout << setfill('-') << setw(40) << "\n";
+	for(i=0;i<npar;i++){cout << pars[i] << " ";} cout << endl;
+	for(i=0;i<npar;i++){cout << errs[i] << " ";} cout << endl;
+	for(i=0;i<npar;i++){cout << grad[i] << " ";} cout << endl;
+	cout << setfill('-') << setw(40) << "\n" << setfill(' ');
+	print_pars( pars, errs );
+}
+
+
+void mod::step_pars( int n , vec& dpar , vec& pars ){
+	for( int i=0 ; i<npar ; i++ ){
+		pars[i] += dpar[i];
+	}
+}
+
+void mod::print_pars( vec& pars, vec& errs ){
+	cout << setfill('*') << setw(40) << "\n" << setfill(' ');
+	cout << endl << "# "; for(int i=0;i<=nfq;i++){cout << FqL[i] << " ";} cout << endl;
+	for(int i=0;i<npar;i++){ printf( "%3.3e %3.3e\n" , pars[i] , errs[i] );}
+	cout << setfill('*') << setw(40) << "\n" << setfill(' ');
+}
+
+void mod::print_pars( vec& pars, vec& errs1 , vec& errs2 ){
+	cout << setfill('*') << setw(40) << "\n" << setfill(' ');
+	cout << endl << "# "; for(int i=0;i<=nfq;i++){cout << FqL[i] << " ";} cout << endl;
+	for(int i=0;i<npar;i++){ printf( "%3.3e %3.3e %3.3e\n" , pars[i] , errs1[i] , errs2[i] );}
+	cout << setfill('*') << setw(40) << "\n" << setfill(' ');
+}
+
+void mod::what_pars( int& ip1 , int& ip2 ){
+	ip1 = 0; ip2 = npar;
+}

psdlag/src/psd.cpp

@@ -0,0 +1,105 @@
+/*
+ * psd.cpp
+ *
+ *  Created on: May 31, 2013
+ *      Author: azoghbi
+ */
+
+#include "inc/psd.hpp"
+
+psd::psd( lcurve inlc , vec fqL ) {
+
+	// ----------- initial parameters ------------ //
+	n		=	inlc.len;
+	dt		=	inlc.dt;
+	// ------------------------------------------ //
+
+
+	// ----------- light curve setup ------------ //
+	setlc();
+	inlc.demean();
+	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::~psd() {}
+
+void psd::_C( vec p ){
+	int	i,j,k; double d;
+	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::_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::psd10( lcurve inlc , vec fqL ) {
+
+	// ----------- initial parameters ------------ //
+	n		=	inlc.len;
+	dt		=	inlc.dt;
+	// ------------------------------------------ //
+
+
+	// ----------- light curve setup ------------ //
+	setlc();
+	inlc.demean();
+	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::~psd10() {}
+
+void psd10::_C( vec p ){
+	int	i,j,k; double d; for(i=0;i<npar;i++){p[i]=pow(10,p[i]);}
+	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::_dC( vec p , int k ){
+	int	i,j; double d = log(10)*pow(10,p[k]);
+	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::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];
+	}
+}

psdlag/src/psdlag.cpp

@@ -0,0 +1,252 @@
+/*
+ * psdlag.cpp
+ *
+ *  Created on: Jun 1, 2013
+ *      Author: azoghbi
+ */
+
+#include "inc/psdlag.hpp"
+
+psdlag::psdlag( lcurve lc1, lcurve lc2 , vec fqL ) {
+
+	// ----------- initial parameters ------------ //
+	n1		=	lc1.len;
+	n		=	n1 + lc2.len;
+	dt		=	lc1.dt;
+	// ------------------------------------------ //
+
+
+	// ----------- light curve setup ------------ //
+	setlc();
+	lc1.demean(); lc2.demean();
+	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::~psdlag() {}
+
+void psdlag::_C( vec p ){
+	int	i,j,k; double d;
+	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::_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::_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::_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::_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::_dPhi( vec p , int k ){
+	int	i,j; double cx=p[icx[k]],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::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::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::psdlag10( lcurve lc1, lcurve lc2 , vec fqL ) {
+
+	// ----------- initial parameters ------------ //
+	n1		=	lc1.len;
+	n		=	n1 + lc2.len;
+	dt		=	lc1.dt;
+	// ------------------------------------------ //
+
+
+	// ----------- light curve setup ------------ //
+	setlc();
+	lc1.demean(); lc2.demean();
+	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::~psdlag10() {}
+
+void psdlag10::_C( vec p ){
+	int	i,j,k; double d; for(i=0;i<3*nfq;i++){ p[i] = pow(10,p[i]);}
+	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::_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::_dP1( vec p , int k ){
+	int	i,j; double p1 = log(10)*pow(10,p[k]);
+	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::_dP2( vec p , int k ){
+	int	i,j; double p2 = log(10) * pow(10,p[ip2[k]]);
+	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::_dCx( vec p , int k ){
+	int	i,j; double phi=p[iphi[k]], cx=log(10)*pow(10,p[icx[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)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
+	}
+}
+
+void psdlag10::_dPhi( vec p , int k ){
+	int	i,j; double cx=pow(10,p[icx[k]]),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::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::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::what_pars( int& ip1 , int& ip2 ){
+	ip1 = 3*nfq; ip2 = 4*nfq;
+}

psdlag/src/stdafx.h

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