/*
Jason Lawrence
Princeton University

acls.h

Several variants of the Alternating Constrained Least Squares (ACLS)
algorithm described in:

Lawrence, J., Ben-Artzi, A., DeCoro, C., Matusik, W., Pfister, H.,
Ramamoorthi, R., Rusinkiewicz, S. Inverse Shade Trees for
Non-Parametric Material Representation and Editing. ACM Transactions
on Graphics (Proceedings of ACM SIGGRAPH 2006).

This code calls several functions provided by the Numerical Algorithms
Group (NAG) library set.

History:
   4/2006 - Created (jason)
*/

#ifndef __ACLS_H__
#define __ACLS_H__

#include <nag.h>
#include <nag_types.h>
#include <nage04.h>

#define ONE_MB 1048576

typedef struct _acls_workspace_t 
{
  int max_k;

  double *Wt;
  double *Ht;
  double *b;
  double *x;
  double  objf;
  long   *kx;
  double *bl;
  double *bu;
  double *cvec;
  double *sums;
  double *objs;
  int    *inc;
  float  *WH;
  double *A;

} acls_ws_t;


// Initialize workspace variables for future calls to 'acls'.  M and N
// must not change between calls to 'acls', and max_k is the largest
// value of K that 'acls' will be called with.
acls_ws_t * create_acls_ws (int M,
			    int N,
			    int max_k);

// Reclaim memory allocated for acls workspace.
void free_acls_ws(acls_ws_t *ws);

// Main ACLS routine.
float acls (float *V,                                            // Data matrix
	    float *B,                                            // Confidence matrix
	    float *R,                                            // Residual matrix
	    int M, int N,                                        // Dimensions
	    float *W, float *H,                                  // Factorization (must be initialized to lie within feasible region)
	    int K,                                               // Desired term count
	    acls_ws_t *ws,                                       // ACLS Workspace Object
	    float lambda,                                        // Sparsity parameter
	    float mu,                                            // Unity norm parameter
	    void (*cb)(int, float *, float *, int, int, int),    // Visualization callback
	    int min_iters,                                       // Iteration bounds
	    int max_iters,
	    std::vector<double> *t,                              // Vector of time stamps 
	    std::vector<double> *ssds,                           // Vector or error stamps
	    std::vector<double> *lbound,                         // Value bounds for W
	    std::vector<double> *ubound);

#endif