/*** seq_ls.h  Manuel Loth 2010 **********************/

#ifndef SEQ_LS_H
#define SEQ_LS_H

/**
  Notation: M' is the transpose of M

  A seq_ls structure maintains in some form a squared
  matrix A=MM' troughout line appendings/removals in M.
  It provides a function 'solve' which, given a vector b,
  solves Ax = b  for x.

  The typical use of solving (MM')x = b is least-squares 
  regression:
    Given a vector y of responses in R and a matrix X of 
    observations  (columns = observations associated to 
    responses, lines = n predictors/attributes/features), 
    one wishes to approximate y as X'w, by finding w that minimizes

     1/2 |y - X'*w|2² = Sum_i [y_i - Sum_j (wj*Xji)]²

    , which amounts to zeroing the gradient:
      X (y - X'w) = (0)
      XX'w = Xy

     (hence   A = XX',   b = Xy)
  
  The appending / removal of lines in X corresponds to 
  adding / removing an attribute.

***********************************************************/

typedef struct seq_ls {
  
  int n;        /* number of lines in M (observations) */
  double ** M;

  /** 
     Appends line u at bottom of M.
     /!\ CAUTION: u itself (not a copy) will be used in
     subsequent steps, so it should not be changed or
     deallocated until being removed from M.  Returns
     non-zero value in case of success and zero in case of
     failure, which consists in A being no more
     positive-definite after insertion of u, in which case
     it is not inserted (the function does nothing).
  */
  int (*append_line) (struct seq_ls *self, double * u);

  /** 
     Removes line i (i>=0) from M.
  */
  double * (*remove_line) (struct seq_ls *self, int i);

  /**
     Solves MM' * x = b. Stores the result at given address
     if not null, otherwise allocates as needed (and returns
     the new address).
  */
  double * (*solve) (struct seq_ls *self, double * b, double * x);

  /**
     Multiply A(=MM') or M by a vector v (A*v or
     v'*M). Store the result at given address if not null,
     otherwise allocate as needed (and return the new
     address).
  */
  double * (*mul_Av) (struct seq_ls *self, 
		      double * v, double * dst);
  double * (*mul_vM) (struct seq_ls *self, 
		      double * v, double * dst);

  /**
     Frees this structure and its workspace
  */
  int (*delete) (struct seq_ls *self);

  void *private;

} seq_ls_t;



/**
 Constructors return new seq_ls structures, with empty M / A.
 M_max_height is not a strict limit, but an indication for
 intial memory allocation 
*/

/* maintains the Cholesky decomposition of A */
seq_ls_t *new_seq_Cholesky (int M_width, int M_max_height);

/* maintains the QR decomposition of A */
/* seq_ls_t *new_seq_QR (int M_width, int M_max_height); */

#endif