bcsol.c

Go to the documentation of this file.

/**************************************************************/
/*! \file
  
    \brief This file contains all the routines needed for 
           BCLS optimization. 
*/
/**************************************************************/



#include<slim.h>

/*! \brief Number of iterations */
int bcls_niters = 0; 

/* The following macros are from BCLS */
/*! \brief Compute the maximum of two */
#define CS_MAX(a,b) (((a) > (b)) ? (a) : (b))
/*! \brief Compute the minimum of two */
#define CS_MIN(a,b) (((a) < (b)) ? (a) : (b))
/*! \brief Flip */
#define CS_FLIP(i) (-(i)-2)
/*! \brief Unflip */
#define CS_UNFLIP(i) (((i) < 0) ? CS_FLIP(i) : (i))
/*! \brief Check if marked */
#define CS_MARKED(w,j) (w [j] < 0)
/*! \brief Mark */
#define CS_MARK(w,j) { w [j] = CS_FLIP (w [j]) ; }
/*! \brief CSC */
#define CS_CSC(A) (A && (A->nz == -1))
/*! \brief Triplet*/
#define CS_TRIPLET(A) (A && (A->nz >= 0))


// Type of projected search.                                                                  
static int proj_search = BCLS_PROJ_SEARCH_EXACT;
/* static int proj_search = BCLS_PROJ_SEARCH_APPROX; */
// Method for Newton step computation.                                                        
static int newton_step = BCLS_NEWTON_STEP_LSQR;



/**************************************************************/
/*! \brief call_back function, periodically called by BCLS to 
   test if the user wants to exit. This is from BCLS.
*/
/**************************************************************/                             
int call_back( BCLS *ls, void *UsrWrk ){
  int err;
  err = 0;  /* No error. */
  return err;
}

/**************************************************************/
/*! \brief call_back function, immediately terminate BCLS 
           iterations based on how many iterations it runs
*/
/**************************************************************/
int call_back_it( BCLS *ls, void *UsrWrk ){

  bcls_niters ++; 
  if (bcls_niters == ((worksp*)UsrWrk)->max_bcls_niters)
    return 1; 
  else 
    return 0; 
}



/**************************************************************/
/*! \brief Pretty_printer, this is the print-routine that will 
           be used by BCLS for its output. This is from BCLS.
*/
/**************************************************************/
int pretty_printer( void *io_file, char *msg ) {
  fprintf( io_file, "%s", msg );
  return 0;
}


/**************************************************************/
/*! \brief Wrapper for free
 */
/**************************************************************/
void *cs_free (void *p){
  if (p) free (p) ;     /* free p if it is not already NULL */
  return (NULL) ;       /* return NULL to simplify the use of cs_free */
}


/**************************************************************/
/*! \brief Wrapper for malloc 
 */
/**************************************************************/
void *cs_malloc (int n, size_t size){
  return (malloc (CS_MAX (n,1) * size)) ;
}

/**************************************************************/
/*! \brief Wrapper for calloc 
*/
/**************************************************************/
void *cs_calloc (int n, size_t size){
  return (calloc (CS_MAX (n,1), size)) ;
}



/**************************************************************/
/*! \brief Free a sparse matrix 
 */
/**************************************************************/
cs *cs_spfree (cs *A){
  if (!A) return (NULL) ;     /* do nothing if A already NULL */
  cs_free (A->p) ;
  cs_free (A->i) ;
  cs_free (A->x) ;
  return (cs_free (A)) ;      /* free the cs struct and return NULL */
}

/**************************************************************/
/*! \brief Allocate a sparse matrix (triplet form or 
  compressed-column form) 
*/
/**************************************************************/
cs *cs_spalloc (int m, int n, int nzmax, int values, int triplet){
  cs *A = cs_calloc (1, sizeof (cs)) ;    /* allocate the cs struct */
  if (!A) return (NULL) ;                 /* out of memory */
  A->m = m ;                              /* define dimensions and nzmax */
  A->n = n ;
  A->nzmax = nzmax = CS_MAX (nzmax, 1) ;
  A->nz = triplet ? 0 : -1 ;              /* allocate triplet or comp.col */
  A->p = cs_malloc (triplet ? nzmax : n+1, sizeof (int)) ;
  A->i = cs_malloc (nzmax, sizeof (int)) ;
  A->x = values ? cs_malloc (nzmax, sizeof (double)) : NULL ;
  return ((!A->p || !A->i || (values && !A->x)) ? cs_spfree (A) : A) ;
}


/**************************************************************/
/*! \brief Load a constant into a vector
 */
/**************************************************************/                              
void dload( const int n, const double alpha, double x[] ) {

  int j;
  for (j = 0; j < n; j++) x[j] = alpha;
  return;
}




/**************************************************************/ 
/*! \brief Aprod, matrix-vector products. This is from BCLS.

  \details
  If     mode == BCLS_PROD_A  (0), compute y <- A *x, with x untouched;
  and if mode == BCLS_PROD_At (1), compute x <- A'*y, with y untouched.

*/
/**************************************************************/ 
int Aprod( int mode, int m, int n, int nix, int ix[],
       double x[], double y[], void *UsrWrk ) {
  
  int     i, j, k, l;
  double  aij;
  double xj, sum;
  worksp * Wrk = (worksp *)UsrWrk;
  cs *A = (cs *)Wrk->A;
  int * Ai = A->i; 
  int * Ap = A->p; 
  float * Ax = A->x; 
  int * acol = Wrk->acol; 


  if (mode == BCLS_PROD_A) {

    gk_dset(m, 0.0, y); 

    for (l = 0; l < nix; l++) {
      j = ix[l];

      /* skip the inactive column */
      if (!acol[j]) continue; 
      
      xj = x[j];
      if (xj == 0.0)
    ; // Relax.
      else
    for (k = Ap[j]; k < Ap[j+1]; k++) {
      aij   = Ax[k];


      /* this is to handle float-valued A matrix */
      i     = Ai[k];
      y[i] += aij * xj; 
 
 
    }
    }
  }

  else if (mode == BCLS_PROD_At) {
    for (l = 0; l < nix; l++) {
      j = ix[l];
      sum = 0;

      /* skip the inactive column */
      if (!acol[j]){
    x[j] = sum; 
    continue; 
      }

      for (k = Ap[j]; k < Ap[j+1]; k++) {
    aij  = Ax[k];

    /* this is to handle float-valued A matrix */
    i    = Ai[k];
    sum += aij * y[i]; 

    
      }
      x[j] = sum;
    }
  }


  return 0;
}



/**************************************************************/ 
/*! \brief BCLS learning. This is from BCLS

    \details This is to solve the problem
    \f[
    \begin{array}{ll} 
    \displaystyle\mathop{\hbox{minimize}}_x & \frac12\|Ax-a_i\|_2^2 + \frac{1}{2}\beta\|x\|_2^2 + \lambda \|x\|_1 \\ 
    \hbox{subject to} & 0 \le x \\
                      & x_i = 0 \\
    \end{array} 
    \f]
 */
/**************************************************************/ 
void  bcsol(ctrl_t * ctrl, gk_csr_t * AA, double * bb, double * x, worksp * Wrk, 
        double * bl, double * bu, 
        double beta, double * c){

  bcls_niters = 0; 

  /*  Problem dimensions. */
  int m = AA->nrows; 
  int n = AA->ncols; 

  /* init a bcls problem */
  BCLS   *ls = bcls_create_prob( m, n ); 

  bcls_set_problem_data(ls, m, n, Aprod,  Wrk, beta, x, bb, c, bl, bu); 


  /* set up tolerance */
  ls->optTol = ctrl->optTol;  
    
  /* whatever */
  bcls_set_print_hook( ls, stdout, pretty_printer );
  ls->proj_search    = proj_search; 
  ls->newton_step    = newton_step; 
  /* if (ctrl->test_bcls) */
  if (ctrl->max_bcls_niters > 0)
    ls->CallBack       = call_back_it; 
  else
    ls->CallBack       = call_back; 

  /* call the solver */
  int err = bcls_solve_prob( ls ); 

  /* solution */
  if (ctrl->dbglvl > 1){
    int nnzx = 0; 
    printf("\n Solution\n --------\n");
    printf("%4s  %18s %1s %18s %1s %18s  %18s\n",
       "Var","Lower","","Value","","Upper","Gradient");
    for (int j = 0; j < n; j++) {
      if (x[j] > 1e-10){
    nnzx ++; 
    char * blActiv = "";
    char * buActiv = "";
    if (x[j] - bl[j] < ls->epsx) blActiv = "=";
    if (bu[j] - x[j] < ls->epsx) buActiv = "=";
    printf("%4d  %18.11e %1s %18.11e %1s %18.11e  %18.11e\n",
           j+1, bl[j], blActiv, x[j], buActiv, bu[j], (ls->g)[j]);
      }
    }
    printf("%d nnz solution values\n", nnzx); 
  }
    

  /* free the problem */
  err = bcls_free_prob( ls );

}