/* Last edited on 2024-12-05 11:03:04 by stolfi */

/* ---------------------------------------------------------------------- *
/* jspca.h,c */

  
  /* The array {vin} must contain values for {nd} data vectors and {nv}
  variables. More precisely, the procedure assumes that {vin[id][iv]},
  is the value of variable {iv} in data vector {id}, for {id} in
  {0..nd-1} and {iv} in {0..nv-1}. */
  
  double *Q,  /* Inverse of moment matrix. */
    
  /* The matrix {Q} should be an inverse of {D'*D} as computed by
    {jspca_compute_moment_matrix}. The coefficients
    {ci[0..nd-1]} are computed by the product {Q*D*vi'}.

    The array {Q}, when given, must be a matrix with {nv} rows and {nv}
    columns, linearized by rows, which is the inverse of {D'*D} where {D'}
    is the transpose of {D}.*/

    rmxn_mul_tr(nd, nd, nv, D, D, R);
    rmxn_inv(nd,R,Q);
    /* Check the inverse: */
    double *S = rmxn_alloc(nd,nd); /* Should be the identity. */
    rmxn_mul(nd,nd,nd,R,Q,S);
    for (int32_t i = 0; i < nd; i++)
      { for (int32_t j = 0; j < nd; j++)
          { double Sij = S[i*nd + j];
            double Iij = (i == j ? 1.0 : 0.0);
            if (fabs(Sij - Iij) > 1.0e-6)
              { fprintf(stderr, "  S[%d,%d] = %24.16e\n", i, j, Sij);
                demand(FALSE, "inversion failed, aborted");
              }
          }
      }
    free(S);

void jspca_combine_data_vectors
  ( int32_t nd, 
    int32_t nv, 
    int32_t nd, 
    double *D, 
    double **coef, 
    int32_t mp, 
    int32_t id[], 
    double **vout
  );
  /* Creates a dataset from a selected subset of data vectors from {D}
    and their coefficients from {coef}.

    Assumes that {id[0..mp-1]} are integers in the range {0..nd-1}.
    
    The procedure stores into channels {0..nv-1} of each data vector
    {vout[id]} the EEG variable values corresponding to the
    linear combination of rows {id[0..mp-1]} of {D}, with their
    respective coefficients from {coef[id]}.  That is
    
      {vout[id][iv] = SUM{ coef[id][id[k]]*D[id[k]*nv + iv] : k \in 0..mp-1}}
      
    for {id} in {0..nd-1} and {iv} in {0..nv-1}.  Any excess channels in {vout}
    are left unchanged. */


        
        rmxn_map_col(nd, nv, D, Di, bi);
        rmxn_map_col(nd, nd, Q, bi, ci);

void jspca_combine_data vectors
  ( int32_t nt, 
    int32_t nv, 
    int32_t nd, 
    double *D, 
    double **coef, 
    int32_t mp, 
    int32_t ip[], 
    double **vout
  )
  {
    demand((0 <= nd) && (nd <= nv), "invalid {nd,nv}");
    int32_t it;
    for (it = 0; it < nt; it++)
      { double *cf = coef[it]; /* Fitted data vector coefficients. */
        double *vo = vout[it];        /* Reconstructed variables. */
        /* Combine the requested components: */
        rn_zero(nv, vo);
        int32_t k;
        for (k = 0; k < mp; k++)
          { int32_t ipk = ip[k];
            demand((0 <= ipk) && (ipk < nd), "invalid data vector index");
            double *Pi = &(D[ipk*nv]); /* Selected data vector. */
            double cfi = cf[ipk]; /* Its coefficient. */
            rn_mix_in (nv, cfi, Pi, vo);
          }
      }
  }
  

/* ---------------------------------------------------------------------- *
/* lsq_robust.c */

    auto void gen_case(int it, int nx, double xi[], int nf, double fi[], double *wiP);
      /* Returns row {it} of {X} in {xi}, row {it} of {F} in {fi}, and {W[it]} in {*wiP}. */
      
    void gen_case(int it, int nx, double xi[], int nf, double fi[], double *wiP)
      { double *Xi = &(X[it*nx]);
        double *Fi = &(F[it*nf]);
        int j;
        for (j = 0; j < nx; j++) { xi[j] = Xi[j]; }
        for (j = 0; j < nf; j++) { fi[j] = Fi[j]; }
        (*wiP) = W[it];
      }

/* ---------------------------------------------------------------------- *
/* test_lsq_robust.c */

    int nbez = test_lsq_robust_bez_count(MAX_DEG); /* Number of Bernstein polys with degree at most {MAX_DEG}. */

    int g = (trial < nbez ? test_lsq_robust_bez_deg(trial) : test_lsq_robust_randeg()); /* Poly degree. */

int test_lsq_robust_randeg(void);
  /* Returns a random degree in {0..MAX_DEG}. */
  
int test_lsq_robust_bez_deg(int ibez);
  /* Assumes {ibez} is a global index of a Bernstein polynomial, that
    is, an index in the list of all Bernstein polynomials sorted by
    degree then by Bernstein index within those with same degree.
    Returns the degree of that polynomial. */
  
int test_lsq_robust_bez_count(int g);
  /* Returns the number of Bernstein polynomials with degree {g} or less. */
  
int test_lsq_robust_bez_index(int g, int ibez);
  /* Assumes {ibez} is a global index of a Bernstein polynomial.
    Returns the Bernstein index of that polynomial within those
    of degree {g}. */

int test_lsq_robust_randeg(void)
  { return abrandom(0, MAX_DEG); }
  
int test_lsq_robust_bez_deg(int ibez)
  { int g = 0;
    while (ibez >= test_lsq_robust_bez_count(g)) { g++; }
    return g;
  }
  
int test_lsq_robust_bez_count(int g)
  { return (g+1)*(g+2)/2; }
  
int test_lsq_robust_bez_index(int g, int ibez)
  { int iblo = (g == 0 ? 0 : test_lsq_robust_bez_count(g-1));
    int ibhi = test_lsq_robust_bez_count(g) - 1;
    demand((iblo <= ibez) && (ibez <= ibhi), "invalid degree or global index");
    return ibez - iblo;
  }

/* ---------------------------------------------------------------------- *
/* jspca.c */

    /* Convert {A} to symmetric tridiagnonal form {T}: */
    if (verbose) { fprintf(stderr, "  tridiagonalizing {A}...\n"); }
    double *sT = rn_alloc(nv); /* Subdiagonal of {T}. */
    double *dT = rn_alloc(nv); /* Eigenvalues. */
    double *R = rmxn_alloc(nv,nv); /* Eigenvectors. */
    syei_tridiagonalize(nv, A, dT, sT, R);

    /* Compute the eigensystem for {T}: */
    if (verbose) { fprintf(stderr, "  computing the eigenvalues of the tridiagonal matrix...\n"); }
    int32_t ne; /* Number of eugenvalues actually computed: */
    int32_t absrt = 0; /* Sort eigenvalues by *signed* value. */
    syei_trid_eigen(nv, dT, sT, R, &ne, absrt);