load("matrix.dll");

/********************************************************************
 * Get rid of missing data
 *
 * Input:
 *      V -- data array
 *      vmin, vmax -- data range
 *
 * Output:
 *     Null for no missing or an array with
 *     [0] -- new data array without missing data
 *     [1] -- index of non-missing data of V (int8 array)
 *
 ********************************************************************/

function ridmissing(V, vmin=-1.e30, vmax=1.e30)
{
    assert(isuser(V) & vmax > vmin);

    dim = V:size();
    I = (V >= vmin) && (V <= vmax);
    m = I:sum();
    
    if (m == dim[2]) return null;
    
    U = V[I];
    
    return [U, I];
}

/********************************************************************
 * Least square fitting of
 * 
 *     y = c0 + c1*x1 + c2*x2...
 * 
 * The problem can be formulated to
 * 
 *      C = invert(X~*X)*X~*Y.
 * 
 * Input:
 *      A -- N by M (N > M) design matrix of double type.
 *      Y -- N by 1 matrix of double type.
 * 
 * Output:
 *      Coefficients in a M by 1 double type matrix.
 *
 ********************************************************************/ 

function mlinear(A, Y)
{
    C = A + 0;
    C:trans();
    C:prod(A);
    
    assert(C:invert("ps"));

    A:trans();
    C:prod(A);
    C:prod(Y);
    
    return C;
}

/********************************************************************
 * Trigonometric polynomial curve fitting of
 * 
 *      y = a + b*x + c1*cos(x) + d1*sin(x) + c2*cos(2x) + d2*sin(2x)...
 * 
 * Input:
 *      X -- the independent variable (double array)
 *      Y -- the depend  variable (double array)
 *      period -- period of X
 *      ns -- number of harmonic terms
 * 
 * Output:
 *      coefficients in a double array
 *
 ********************************************************************/

function trigono(X, Y, period, ns)
{
    assert(isuser(X) & isuser(Y) & ns > 0 & period > 0);
    
    dim = X:size();
    n = dim[2];
    xmax = X:max();
    xmin = X:min();
    assert(xmax != 0);
    x = (X - xmin) / xmax;                  // use normalized x for the linear term

    p = 6.28318530717959 / period;
    
    A = matrix("double", n, 2*ns+2);        // the design matrix
    A[null,0] = 1;                          // for constant coefficient
    A[null,1] = x;                            // for trend coefficient

    for (i = 1; i <= ns; i++) {
        P = p * i * X;
        A[null,2*i] = cos(P);               // for cos coefficient
        A[null,2*i+1] = sin(P);             // for sin coefficient
    }

    C = mlinear(A, Y);
    C[1] = C[1] / xmax;
    C[0] = C[0] - xmin * C[1];

    return C;
}

/********************************************************************
 * Evaluate trigonometric fitting
 *
 * input:
 *    X -- the independent variable (double array)
 *    c -- coefficients from trigono() (double array)
 *    period -- period of X
 *    ns -- number of harmonic terms
 *
 * output:
 *    CTK array with
 *    [0] -- the constant
 *    [1] -- the trend (double array)
 *    [2] -- the harmonic term
 * 
 ********************************************************************/

function trigono_eval(X, c, period, ns)
{
    assert(isuser(X) & isarray(c) & ns > 0 & period > 0);
    
    p = __PI / period;

    Y1 = X * 0.0;
    Y1 += c[1] * X;
    
    Y2 = X * 0.0;

    for (i = 1; i <= ns; i++) {
        a = p * i * X;                      // period
        Y2 += c[2*i  ] * cos(a);
        Y2 += c[2*i+1] * sin(a);
    }
    
    return [c[0], Y1, Y2];
}