From 2a4780f4106071da6cfc0e9503bfd655c99eeff0 Mon Sep 17 00:00:00 2001 From: Samuel St-Jean <3030760+samuelstjean@users.noreply.github.com> Date: Tue, 28 Jan 2025 11:20:40 -0500 Subject: [PATCH 1/9] add pos threshold to corelars2w --- spams_wrap/spams/decomp/decomp.h | 292 ++++++++++++++++--------------- 1 file changed, 148 insertions(+), 144 deletions(-) diff --git a/spams_wrap/spams/decomp/decomp.h b/spams_wrap/spams/decomp/decomp.h index 25b4e0b5..a4e31559 100644 --- a/spams_wrap/spams/decomp/decomp.h +++ b/spams_wrap/spams/decomp/decomp.h @@ -1,5 +1,5 @@ /*! - * Software SPAMS v2.2 - Copyright 2009-2011 Julien Mairal + * Software SPAMS v2.2 - Copyright 2009-2011 Julien Mairal * * This file is part of SPAMS. * @@ -24,34 +24,34 @@ * julien.mairal@inria.fr * * File decomp.h - * \brief Contains sparse decomposition algorithms + * \brief Contains sparse decomposition algorithms * It requires the toolbox linalg */ #ifndef DECOMP_H #define DECOMP_H #include "utils.h" - + static char low='l'; static char nonUnit='n'; /* ************************** - * Greedy Forward Selection + * Greedy Forward Selection * **************************/ -/// Forward Selection (or Orthogonal matching pursuit) +/// Forward Selection (or Orthogonal matching pursuit) /// Address the problem of: -/// \forall i, \min_{\alpha_i} ||X_i-D\alpha_i||_2^2 +/// \forall i, \min_{\alpha_i} ||X_i-D\alpha_i||_2^2 /// s.t. ||\alphai||_0 <= L or -/// \forall i, \min_{\alpha_i} ||\alpha_i||_0 +/// \forall i, \min_{\alpha_i} ||\alpha_i||_0 /// s.t. ||\X_i-D\alpha_i||_2^2 <= epsilon -/// This function is +/// This function is /// * based on Cholesky decompositions /// * parallel /// * optimized for a large number of signals (precompute the Gramm matrix template -void omp(const Matrix& X, const Matrix& D, SpMatrix& spalpha, +void omp(const Matrix& X, const Matrix& D, SpMatrix& spalpha, const int *L, const T* eps, const T* lambda, const bool vecL = false, const bool vecEps = false, const bool Lambda=false, const int numThreads=-1, Matrix* path = NULL); @@ -64,16 +64,16 @@ void omp_mask(const Matrix& X, const Matrix& D, SpMatrix& spalpha, cons /// Auxiliary function of omp template -void coreORMP(Vector& scores, Vector& norm, Vector& tmp, - Matrix& Un, Matrix& Undn, Matrix& Unds, Matrix& Gs, - Vector& Rdn, const AbstractMatrix& G, Vector& ind, +void coreORMP(Vector& scores, Vector& norm, Vector& tmp, + Matrix& Un, Matrix& Undn, Matrix& Unds, Matrix& Gs, + Vector& Rdn, const AbstractMatrix& G, Vector& ind, Vector& RUn, T& normX, const T* eps, const int* L, const T* lambda, T* path = NULL); /// Auxiliary function of omp template -void coreORMPB(Vector& RtD, const AbstractMatrix& G, Vector& ind, +void coreORMPB(Vector& RtD, const AbstractMatrix& G, Vector& ind, Vector& coeffs, T& normX, const int L, const T eps, const T lambda = 0); @@ -85,39 +85,39 @@ void coreORMPWeighted(Vector& scores, Vector& weights, Vector& norm, ind, Vector& RUn, T& normX, const T eps, const int L, const T lambda);*/ /* ************** - * LARS - Lasso + * LARS - Lasso * **************/ /// Defines different types of problem, /// - constraint on the l1 norm of the coefficients /// - constraint on the reconstruction error -/// - l1-sparsity penalty +/// - l1-sparsity penalty enum constraint_type { L1COEFFS, L2ERROR, PENALTY, SPARSITY, L2ERROR2, PENALTY2,FISTAMODE}; /// Implementation of LARS-Lasso for solving -/// \forall i, \min_{\alpha_i} ||X_i-D\alpha_i||_2^2 +/// \forall i, \min_{\alpha_i} ||X_i-D\alpha_i||_2^2 /// s.t. ||\alphai||_1 <= constraint or -/// \forall i, \min_{\alpha_i} ||\alpha_i||_1 +/// \forall i, \min_{\alpha_i} ||\alpha_i||_1 /// s.t. ||\X_i-D\alpha_i||_2^2 <= constraint or /// \forall i, \min_{\alpha_i} constraint*||\alpha_i||_1 + ... /// ... ||\X_i-D\alpha_i||_2^2 <= T -/// Optionally, the solution might be positive (boolean pos), and a +/// Optionally, the solution might be positive (boolean pos), and a /// Least-Square can be solved as a post-processing step. /// L is a maximum number of coefficients. -/// This function is +/// This function is /// * efficient (Cholesky-based) /// * parallel /// * optimized for a big number of signals (precompute the Gramm matrix template -void lasso(const Matrix& X, const Matrix& D, - SpMatrix& spalpha, +void lasso(const Matrix& X, const Matrix& D, + SpMatrix& spalpha, int L, const T constraint, const T lambda2 = 0, constraint_type mode = PENALTY, const bool pos = false, const bool ols = false, const int numThreads=-1, Matrix* path = NULL, const int length_path=-1); template void lasso(const Data& X, const AbstractMatrix& G, const AbstractMatrix& DtX, - SpMatrix& spalpha, + SpMatrix& spalpha, int L, const T constraint, constraint_type mode = PENALTY, const bool pos = false, const bool ols = false, const int numThreads=-1, Matrix* path = NULL, const int length_path=-1); @@ -149,15 +149,15 @@ void lassoReweighted(const Matrix& X, const Matrix& D, SpMatrix& spalph /// Auxiliary function for lasso template -void coreLARS(Vector& Rdn, Vector& Xdn, Vector& A, +void coreLARS(Vector& Rdn, Vector& Xdn, Vector& A, Vector& u, Vector& sig, - Vector& av, Vector& RUn, Matrix& Un, + Vector& av, Vector& RUn, Matrix& Un, Matrix& Unds, Matrix& Gs, Matrix& Gsa, Matrix& workT, Matrix& R, - const AbstractMatrix& G,T& normX, + const AbstractMatrix& G,T& normX, Vector& ind,Vector& coeffs,const T constraint, const bool ols = false, - const bool pos =false, + const bool pos =false, constraint_type mode = L1COEFFS, T* path = NULL, int length_path=-1); @@ -215,31 +215,31 @@ void downDateLasso(int& j,int& minBasis,T& normX,const bool ols, * ************************/ /// Implementation of IST for solving -/// \forall i, \min_{\alpha_i} ||\alpha_i||_1 +/// \forall i, \min_{\alpha_i} ||\alpha_i||_1 /// s.t. ||\X_i-D\alpha_i||_2^2 <= constraint or /// \forall i, \min_{\alpha_i} constraint*||\alpha_i||_1 + ... /// ... ||\X_i-D\alpha_i||_2^2 <= T template -void ist(const Matrix& X, const Matrix& D, +void ist(const Matrix& X, const Matrix& D, SpMatrix& spalpha, T lambda, constraint_type mode, - const int itermax=500, + const int itermax=500, const T tol = 0.5, const int numThreads = -1); template -void ist(const Matrix& X, const Matrix& D, +void ist(const Matrix& X, const Matrix& D, Matrix& spalpha, T lambda, constraint_type mode, - const int itermax=500, + const int itermax=500, const T tol = 0.5, const int numThreads=-1); /// coreIST template void coreIST(const AbstractMatrix& G, Vector& DtR, Vector& coeffs, - const T thrs, const int itermax = 500, + const T thrs, const int itermax = 500, const T tol = 0.5); template void coreISTW(const AbstractMatrix& G, Vector& DtR, Vector& coeffs, const Vector& weights, - const T thrs, const int itermax = 500, + const T thrs, const int itermax = 500, const T tol = 0.5); @@ -247,13 +247,13 @@ void coreISTW(const AbstractMatrix& G, Vector& DtR, Vector& coeffs, con template void coreISTconstrained(const AbstractMatrix& G, Vector& DtR, Vector& coeffs, const T normX2, - const T thrs, const int itermax = 500, + const T thrs, const int itermax = 500, const T tol = 0.5); /// ist for group Lasso template void ist_groupLasso(const Matrix* XT, const Matrix& D, - Matrix* alphaT, const int Ngroups, + Matrix* alphaT, const int Ngroups, const T lambda, const constraint_type mode, const int itermax = 500, const T tol = 0.5, const int numThreads = -1); @@ -287,15 +287,15 @@ T computeError(const T normX2, SpVector& coeffs_tmp); /* ****************** - * Simultaneous OMP + * Simultaneous OMP * *****************/ template -void somp(const Matrix* X, const Matrix& D, SpMatrix* spalpha, +void somp(const Matrix* X, const Matrix& D, SpMatrix* spalpha, const int Ngroups, const int L, const T* pr_eps, const bool adapt=false, const int numThreads=-1); template -void somp(const Matrix* X, const Matrix& D, SpMatrix* spalpha, +void somp(const Matrix* X, const Matrix& D, SpMatrix* spalpha, const int Ngroups, const int L, const T eps, const int numThreads=-1); @@ -308,20 +308,20 @@ void coreSOMP(const Matrix& X, const Matrix& D, const Matrix& G, * Implementation of OMP * *********************/ -/// Forward Selection (or Orthogonal matching pursuit) +/// Forward Selection (or Orthogonal matching pursuit) /// Address the problem of: -/// \forall i, \min_{\alpha_i} ||X_i-D\alpha_i||_2^2 +/// \forall i, \min_{\alpha_i} ||X_i-D\alpha_i||_2^2 /// s.t. ||\alphai||_0 <= L or -/// \forall i, \min_{\alpha_i} ||\alpha_i||_0 +/// \forall i, \min_{\alpha_i} ||\alpha_i||_0 /// s.t. ||\X_i-D\alpha_i||_2^2 <= epsilon -/// This function is +/// This function is /// * efficient (Cholesky-based) /// * parallel /// * optimized for a big number of signals (precompute the Gramm matrix template -void omp(const Matrix& X, const Matrix& D, SpMatrix& spalpha, - const int* pL, const T* peps, const T* pLambda, +void omp(const Matrix& X, const Matrix& D, SpMatrix& spalpha, + const int* pL, const T* peps, const T* pLambda, const bool vecL, const bool vecEps, const bool vecLambda, const int numThreads, Matrix* path) { int L; @@ -364,7 +364,7 @@ void omp(const Matrix& X, const Matrix& D, SpMatrix& spalpha, } int i; -#pragma omp parallel for private(i) +#pragma omp parallel for private(i) for (i = 0; i< M; ++i) { #ifdef _OPENMP int numT=omp_get_thread_num(); @@ -386,7 +386,7 @@ void omp(const Matrix& X, const Matrix& D, SpMatrix& spalpha, D.multTrans(Xi,Rdn); coreORMP(scoresT[numT],normT[numT],tmpT[numT],UnT[numT],UndnT[numT],UndsT[numT], GsT[numT],Rdn,G,ind,RUn, normX, vecEps ? peps+i : peps, - vecL ? pL+i : pL, vecLambda ? pLambda+i : pLambda, + vecL ? pL+i : pL, vecLambda ? pLambda+i : pLambda, path && i==0 ? path->rawX() : NULL); } @@ -453,7 +453,7 @@ void omp_mask(const Matrix& X, const Matrix& D, SpMatrix& spalpha, cons } int i; -#pragma omp parallel for private(i) +#pragma omp parallel for private(i) for (i = 0; i< M; ++i) { #ifdef _OPENMP int numT=omp_get_thread_num(); @@ -479,7 +479,7 @@ void omp_mask(const Matrix& X, const Matrix& D, SpMatrix& spalpha, cons T normX = Xi.nrm2sq(); coreORMP(scoresT[numT],normT[numT],tmpT[numT],UnT[numT],UndnT[numT],UndsT[numT], GsT[numT],Rdn,G,ind,RUn, normX, vecEps ? peps+i : peps, - vecL ? pL+i : pL, vecLambda ? pLambda+i : pLambda, + vecL ? pL+i : pL, vecLambda ? pLambda+i : pLambda, path && i==0 ? path->rawX() : NULL); } else { D.copyMask(DmaskT[numT],maski); @@ -492,8 +492,8 @@ void omp_mask(const Matrix& X, const Matrix& D, SpMatrix& spalpha, cons coreORMP(scoresT[numT],normT[numT],tmpT[numT], UnT[numT],UndnT[numT],UndsT[numT], GsT[numT],Rdn,GT[numT],ind,RUn, - normX, &eps_mask, vecL ? pL+i : pL, - vecLambda ? pLambda+i : pLambda, + normX, &eps_mask, vecL ? pL+i : pL, + vecLambda ? pLambda+i : pLambda, path && i==0 ? path->rawX() : NULL); DmaskT[numT].setm(D.m()); @@ -520,7 +520,7 @@ void omp_mask(const Matrix& X, const Matrix& D, SpMatrix& spalpha, cons /// Auxiliary function of omp template -void coreORMPB(Vector& RtD, const AbstractMatrix& G, Vector& ind, +void coreORMPB(Vector& RtD, const AbstractMatrix& G, Vector& ind, Vector& coeffs, T& normX, const int L, const T eps, const T lambda) { const int K = G.n(); Vector scores(K); @@ -539,7 +539,7 @@ template void coreORMP(Vector& scores, Vector& norm, Vector& tmp, Matrix& Un, Matrix& Undn, Matrix& Unds, Matrix& Gs, Vector& Rdn, const AbstractMatrix& G, - Vector& ind, Vector& RUn, + Vector& ind, Vector& RUn, T& normX, const T* peps, const int* pL, const T* plambda, T* path) { const T eps = abs(*peps); @@ -579,9 +579,9 @@ void coreORMP(Vector& scores, Vector& norm, Vector& tmp, Matrix& Un, ind[j]=currentInd; //for (int k = 0; k& scores, Vector& norm, Vector& tmp, Matrix& Un, cblas_copy(j,prUndn+currentInd,K,prUn+j*L,1); cblas_trmv(CblasColMajor,CblasUpper,CblasNoTrans,CblasNonUnit,j,prUn,L,prUn+j*L,1); cblas_scal(j+1,-invNorm,prUn+j*L,1); - + if (j == L-1 || (normX <= eps)) { ++j; break; @@ -606,7 +606,7 @@ void coreORMP(Vector& scores, Vector& norm, Vector& tmp, Matrix& Un, T* last_path=path+(L-1)*K; cblas_copy(j+1,prRUn,1,last_path,1); cblas_trmv(CblasColMajor,CblasUpper,CblasNoTrans,CblasNonUnit, - j+1,prUn,L,last_path,1); + j+1,prUn,L,last_path,1); for (int k = 0; k<=j; ++k) { path[j*K+ind[k]]=last_path[k]; } @@ -628,9 +628,9 @@ void coreORMP(Vector& scores, Vector& norm, Vector& tmp, Matrix& Un, scores.div(norm); for (int k = 0; k<=j; ++k) scores[ind[k]]=T(); } - // compute the final coefficients + // compute the final coefficients cblas_trmv(CblasColMajor,CblasUpper,CblasNoTrans,CblasNonUnit, - j,prUn,L,prRUn,1); + j,prUn,L,prRUn,1); if (path) { memset(path+(L-1)*K,0,L*sizeof(T)); for (int k = 0; k& scores, Vector& norm, Vector& tmp, Matrix& Un, void coreORMPWeighted(Vector& scores, Vector& weights, Vector& norm, Vector& tmp, Matrix& Un, Matrix& Undn, Matrix& Unds, Matrix& Gs, Vector& Rdn, const AbstractMatrix& G, - Vector& ind, Vector& RUn, + Vector& ind, Vector& RUn, T& normX, const T peps, const int pL, const T plambda) { const T eps = abs(*peps); const int L = MIN(*pL,Gs.n()); @@ -684,7 +684,7 @@ void coreORMPWeighted(Vector& scores, Vector& weights, Vector& norm, Ve cblas_copy(j,prUndn+currentInd,K,prUn+j*L,1); cblas_trmv(CblasColMajor,CblasUpper,CblasNoTrans,CblasNonUnit,j,prUn,L,prUn+j*L,1); cblas_scal(j+1,-invNorm,prUn+j*L,1); - + if (j == L-1 || (normX <= eps)) { ++j; break; @@ -706,35 +706,35 @@ void coreORMPWeighted(Vector& scores, Vector& weights, Vector& norm, Ve scores.div(weights); for (int k = 0; k<=j; ++k) scores[ind[k]]=T(); } - // compute the final coefficients + // compute the final coefficients cblas_trmv(CblasColMajor,CblasUpper,CblasNoTrans,CblasNonUnit, - j,prUn,L,prRUn,1); + j,prUn,L,prRUn,1); };*/ /* ************** - * LARS - Lasso + * LARS - Lasso * **************/ /// Implementation of LARS-Lasso for solving -/// \forall i, \min_{\alpha_i} ||X_i-D\alpha_i||_2^2 +/// \forall i, \min_{\alpha_i} ||X_i-D\alpha_i||_2^2 /// s.t. ||\alphai||_1 <= constraint or -/// \forall i, \min_{\alpha_i} ||\alpha_i||_1 +/// \forall i, \min_{\alpha_i} ||\alpha_i||_1 /// s.t. ||\X_i-D\alpha_i||_2^2 <= constraint or /// \forall i, \min_{\alpha_i} constraint*||\alpha_i||_1 + ... /// ... ||\X_i-D\alpha_i||_2^2 <= T -/// Optionally, the solution might be positive (boolean pos), and a +/// Optionally, the solution might be positive (boolean pos), and a /// Least-Square can be solved as a post-processing step. /// L is a maximum number of coefficients. -/// This function is +/// This function is /// * efficient (Cholesky-based) /// * parallel /// * optimized for a big number of signals (precompute the Gramm matrix template -void lasso(const Matrix& X, const Matrix& D, SpMatrix& spalpha, - int L, const T lambda, const T lambda2, constraint_type mode, +void lasso(const Matrix& X, const Matrix& D, SpMatrix& spalpha, + int L, const T lambda, const T lambda2, constraint_type mode, const bool pos, const bool ols, const int numThreads, Matrix* path, const int length_path) { ProdMatrix G(D, X.n() > 10 && D.n() < 50000); @@ -744,9 +744,9 @@ void lasso(const Matrix& X, const Matrix& D, SpMatrix& spalpha, } template -void lasso(const Data& X, const AbstractMatrix& G, - const AbstractMatrix& DtX, SpMatrix& spalpha, - int L, const T lambda, constraint_type mode, +void lasso(const Data& X, const AbstractMatrix& G, + const AbstractMatrix& DtX, SpMatrix& spalpha, + int L, const T lambda, constraint_type mode, const bool pos, const bool ols, const int numThreads, Matrix* path, const int length_path) { @@ -760,7 +760,7 @@ void lasso(const Data& X, const AbstractMatrix& G, if (L <= 0) return; if (path) path->setZeros(); - + int NUM_THREADS=init_omp(numThreads); //ProdMatrix G(D, K < 25000 && M > 10); @@ -799,14 +799,14 @@ void lasso(const Data& X, const AbstractMatrix& G, Vector norms; X.norm_2sq_cols(norms); int i; -#pragma omp parallel for private(i) +#pragma omp parallel for private(i) for (i = 0; i< M; ++i) { #ifdef _OPENMP int numT=omp_get_thread_num(); #else int numT=0; #endif - T normX = norms[i]; + T normX = norms[i]; Vector ind; rM.refCol(i,ind); @@ -817,7 +817,7 @@ void lasso(const Data& X, const AbstractMatrix& G, Vector& Rdn=RdnT[numT]; DtX.copyCol(i,Rdn); coreLARS(Rdn,XdnT[numT], AT[numT], uT[numT], sigT[numT], avT[numT], - RUnT[numT], UnT[numT], UndsT[numT], GsT[numT], GsaT[numT], + RUnT[numT], UnT[numT], UndsT[numT], GsT[numT], GsaT[numT], workT[numT],RT[numT],G,normX, ind,coeffs,lambda,ols,pos, mode,path && i==0 ? path->rawX() : NULL, length_path); } @@ -845,8 +845,8 @@ template void coreLARS(Vector& Rdnv, Vector& Xdnv, Vector& Av, Vector& uv, Vector& sigv, Vector& avv, Vector& RUnv, Matrix& Unm, Matrix& Undsm, Matrix& Gsm, - Matrix& Gsam, Matrix& workm, Matrix& Rm, - const AbstractMatrix& Gm,T& normX, + Matrix& Gsam, Matrix& workm, Matrix& Rm, + const AbstractMatrix& Gm,T& normX, Vector& indv,Vector& coeffsv,const T constraint, const bool ols,const bool pos, constraint_type mode, T* path, int length_path) { @@ -925,7 +925,7 @@ void coreLARS(Vector& Rdnv, Vector& Xdnv, Vector& Av, // cerr << "bad exit" << endl; break; } - + // int iter2 = norm2 < 0.5 ? 2 : 1; // for(int k = 0; k& Rdnv, Vector& Xdnv, Vector& Av, work[k]= diff <= 0 ? INFINITY : (Cmax-Rdn[k])/diff; } for (int k = 0; k<=j; ++k) { - work[ind[k]]=INFINITY; + work[ind[k]]=INFINITY; } - for (int k = 0; k(K,work,1); } else { memset(work,0,2*K*sizeof(T)); for (int k = 0; k<=j; ++k) { const int index=2*ind[k]; - work[index]=INFINITY; - work[index+1]=INFINITY; + work[index]=INFINITY; + work[index+1]=INFINITY; } for (int k = 0; k& Rdnv, Vector& Xdnv, Vector& Av, vDiv(j+1,coeffs,u,work); cblas_scal(j+1,-T(1.0),work,1); /// voir pour petites valeurs - for (int k=0; k<=j; ++k) + for (int k=0; k<=j; ++k) if (coeffs[k]==0 || work[k] <=0) work[k]=INFINITY; minBasis=cblas_iamin(j+1,work,1); gammaMin=work[minBasis]; @@ -1014,7 +1014,7 @@ void coreLARS(Vector& Rdnv, Vector& Xdnv, Vector& Av, T Tu = 0.0; for (int k = 0; k<=j; ++k) Tu += u[k]; - if (Tu > EPSILON) + if (Tu > EPSILON) gamma= MIN(gamma,(constraint-thrs)/Tu); thrs+=gamma*Tu; } @@ -1051,7 +1051,7 @@ void coreLARS(Vector& Rdnv, Vector& Xdnv, Vector& Av, cblas_axpy(K,-gamma,A,1,Rdn,1); if (!pos) currentInd/= 2; if (path) { - for (int k = 0; k<=j; ++k) + for (int k = 0; k<=j; ++k) path[iter*K+ind[k]]=coeffs[k]*sig[k]; } @@ -1070,9 +1070,9 @@ void coreLARS(Vector& Rdnv, Vector& Xdnv, Vector& Av, thrs=abs(Rdn[ind[0]]); } - if ((j == L-1) || + if ((j == L-1) || (mode == PENALTY && (thrs - constraint < 1e-15)) || - (mode == L1COEFFS && (thrs - constraint > -1e-15)) || + (mode == L1COEFFS && (thrs - constraint > -1e-15)) || (newAtom && mode == L2ERROR && (normX - constraint < 1e-15)) || (normX < 1e-15) || (iter >= length_path)) { @@ -1151,13 +1151,13 @@ inline void downDateLasso(int& j,int& minBasis,T& normX,const bool ols, } // Update Unds - for (int k = minBasis+1; k<=j; ++k) + for (int k = minBasis+1; k<=j; ++k) cblas_axpy(j-minBasis,av[k-minBasis-1],Unds+minBasis*L+minBasis+1,1, Unds+k*L+minBasis+1,1); - for (int k = 0; k(j-minBasis,Unds+k*L+minBasis+1,1,Unds+(k-1)*L+minBasis,1); if (num > 0) cblas_trmm(CblasColMajor,CblasRight,CblasLower,CblasTrans,CblasNonUnit, @@ -1218,7 +1218,7 @@ void lassoReweighted(const Matrix& X, const Matrix& D, SpMatrix& spalph vM.resize(L,M); rM.resize(L,M); const int iterR = 30; - + if (L <= 0) return; int NUM_THREADS=init_omp(numThreads); @@ -1254,7 +1254,7 @@ void lassoReweighted(const Matrix& X, const Matrix& D, SpMatrix& spalph } int i; -#pragma omp parallel for private(i) +#pragma omp parallel for private(i) for (i = 0; i< M; ++i) { #ifdef _OPENMP int numT=omp_get_thread_num(); @@ -1314,7 +1314,7 @@ void lassoReweighted(const Matrix& X, const Matrix& D, SpMatrix& spalph template void lassoWeight(const Matrix& X, const Matrix& D, const Matrix& weights, - SpMatrix& spalpha, + SpMatrix& spalpha, int L, const T constraint, constraint_type mode, const bool pos, const int numThreads) { @@ -1325,7 +1325,7 @@ void lassoWeight(const Matrix& X, const Matrix& D, const Matrix& weight Matrix rM; vM.resize(L,M); rM.resize(L,M); - + if (L <= 0) return; int NUM_THREADS=init_omp(numThreads); @@ -1354,7 +1354,7 @@ void lassoWeight(const Matrix& X, const Matrix& D, const Matrix& weight } int i; -#pragma omp parallel for private(i) +#pragma omp parallel for private(i) for (i = 0; i< M; ++i) { #ifdef _OPENMP int numT=omp_get_thread_num(); @@ -1392,7 +1392,7 @@ void lassoWeight(const Matrix& X, const Matrix& D, const Matrix& weight template void lassoWeightPreComputed(const Matrix& X, const Matrix& G, const Matrix& DtR, const Matrix& weights, - SpMatrix& spalpha, + SpMatrix& spalpha, int L, const T constraint, constraint_type mode, const bool pos, const int numThreads) { @@ -1403,7 +1403,7 @@ void lassoWeightPreComputed(const Matrix& X, const Matrix& G, const Matrix Matrix rM; vM.resize(L,M); rM.resize(L,M); - + if (L <= 0) return; int NUM_THREADS=init_omp(numThreads); @@ -1426,7 +1426,7 @@ void lassoWeightPreComputed(const Matrix& X, const Matrix& G, const Matrix } int i; -#pragma omp parallel for private(i) +#pragma omp parallel for private(i) for (i = 0; i< M; ++i) { #ifdef _OPENMP int numT=omp_get_thread_num(); @@ -1505,7 +1505,7 @@ void lasso_mask(const Matrix& X, const Matrix& D, SpMatrix& spalpha, co } int i; -#pragma omp parallel for private(i) +#pragma omp parallel for private(i) for (i = 0; i< M; ++i) { #ifdef _OPENMP int numT=omp_get_thread_num(); @@ -1561,7 +1561,7 @@ void lasso_mask(const Matrix& X, const Matrix& D, SpMatrix& spalpha, co }; template -void lasso2(const Matrix& X, const Matrix& D, SpMatrix& spalpha, +void lasso2(const Matrix& X, const Matrix& D, SpMatrix& spalpha, int L, const T constraint, const T lambda2, constraint_type mode, const bool pos, const int numThreads, Matrix* path, int length_path) { ProdMatrix G(D,X.n() > 10 && D.n() < 50000); @@ -1573,7 +1573,7 @@ void lasso2(const Matrix& X, const Matrix& D, SpMatrix& spalpha, template void lasso2(const Data& X, const AbstractMatrix& G, const AbstractMatrix& DtX, - SpMatrix& spalpha, + SpMatrix& spalpha, int L, const T constraint, constraint_type mode, const bool pos, const int numThreads, Matrix* path, int length_path) { spalpha.clear(); @@ -1608,7 +1608,7 @@ void lasso2(const Data& X, const AbstractMatrix& G, const AbstractMatrix norms; X.norm_2sq_cols(norms); -#pragma omp parallel for private(i) +#pragma omp parallel for private(i) for (i = 0; i< M; ++i) { #ifdef _OPENMP int numT=omp_get_thread_num(); @@ -1650,7 +1650,7 @@ void coreLARS2W(Vector& DtR, const AbstractMatrix& G, Vector& coeffs, const Vector& weights, T normX, const constraint_type mode, const T constraint, const bool pos) { - const INTM p = G.m(); + const INTM p = G.m(); const INTM L = p; Vector v; v.resize(L); @@ -1679,7 +1679,7 @@ void coreLARS2(Vector& DtR, const AbstractMatrix& G, Vector& coeffs, T normX, const constraint_type mode, const T constraint, const bool pos) { - const INTM p = G.m(); + const INTM p = G.m(); const INTM L = p; Vector v; v.resize(L); @@ -1703,7 +1703,7 @@ void coreLARS2(Vector& DtR, const AbstractMatrix& G, }; }; -/// Auxiliary function for lasso +/// Auxiliary function for lasso template void coreLARS2(Vector& DtR, const AbstractMatrix& G, Matrix& Gs, @@ -1770,7 +1770,7 @@ void coreLARS2(Vector& DtR, const AbstractMatrix& G, } } - // Compute the path direction + // Compute the path direction for (int j = 0; j<=i; ++j) pr_work[j]= pr_DtR[pr_ind[j]] > 0 ? T(1.0) : T(-1.0); cblas_symv(CblasColMajor,CblasUpper,i+1,T(1.0),pr_invGs,LL, @@ -1866,7 +1866,7 @@ void coreLARS2(Vector& DtR, const AbstractMatrix& G, thrs += step*coeff1; if (path) { - for (int k = 0; k<=i; ++k) + for (int k = 0; k<=i; ++k) path[iter*K+ind[k]]=pr_coeffs[k]; } @@ -1914,7 +1914,7 @@ void coreLARS2(Vector& DtR, const AbstractMatrix& G, } } -/// Auxiliary function for lasso +/// Auxiliary function for lasso template void coreLARS2W(Vector& DtR, const AbstractMatrix& G, Matrix& Gs, @@ -1982,7 +1982,7 @@ void coreLARS2W(Vector& DtR, const AbstractMatrix& G, } } - // Compute the path direction + // Compute the path direction for (int j = 0; j<=i; ++j) pr_work[j]= pr_DtR[pr_ind[j]] > 0 ? weights[pr_ind[j]] : -weights[pr_ind[j]]; cblas_symv(CblasColMajor,CblasUpper,i+1,T(1.0),pr_invGs,LL, @@ -2030,7 +2030,7 @@ void coreLARS2W(Vector& DtR, const AbstractMatrix& G, // compute the coefficients of the polynome representing normX^2 T coeff1 = 0; for (int j = 0; j<=i; ++j) - coeff1 += pr_DtR[pr_ind[j]] > 0 ? pr_weights[pr_ind[j]]*pr_u[j] : + coeff1 += pr_DtR[pr_ind[j]] > 0 ? pr_weights[pr_ind[j]]*pr_u[j] : -pr_weights[pr_ind[j]]*pr_u[j]; T coeff2 = 0; for (int j = 0; j<=i; ++j) @@ -2055,6 +2055,11 @@ void coreLARS2W(Vector& DtR, const AbstractMatrix& G, // Update coefficients cblas_axpy(i+1,step,pr_u,1,pr_coeffs,1); + if (pos) { + for (int j = 0; j(K,-step,pr_work+2*K,1,pr_DtR,1); @@ -2113,14 +2118,14 @@ void coreLARS2W(Vector& DtR, const AbstractMatrix& G, * ************************/ /// Implementation of IST for solving -/// \forall i, \min_{\alpha_i} ||\alpha_i||_1 +/// \forall i, \min_{\alpha_i} ||\alpha_i||_1 /// s.t. ||\X_i-D\alpha_i||_2^2 <= constraint or /// \forall i, \min_{\alpha_i} constraint*||\alpha_i||_1 + ... -/// ... ||\X_i-D\alpha_i||_2^2 <= lambda +/// ... ||\X_i-D\alpha_i||_2^2 <= lambda template -void ist(const Matrix& X, const Matrix& D, +void ist(const Matrix& X, const Matrix& D, SpMatrix& spalpha, T lambda, constraint_type mode, - const int itermax, + const int itermax, const T tol, const int numThreads) { Matrix alpha; @@ -2131,9 +2136,9 @@ void ist(const Matrix& X, const Matrix& D, } template -void ist(const Matrix& X, const Matrix& D, +void ist(const Matrix& X, const Matrix& D, Matrix& alpha, T lambda, constraint_type mode, - const int itermax, + const int itermax, const T tol, const int numThreads) { if (mode == L1COEFFS) { @@ -2168,7 +2173,7 @@ void ist(const Matrix& X, const Matrix& D, }; int i; -#pragma omp parallel for private(i) +#pragma omp parallel for private(i) for (i = 0; i< M; ++i) { #ifdef _OPENMP int numT=omp_get_thread_num(); @@ -2184,7 +2189,7 @@ void ist(const Matrix& X, const Matrix& D, DtX.copyCol(i,DtR); Vector Xi; X.refCol(i,Xi); - T normX2 = Xi.nrm2sq(); + T normX2 = Xi.nrm2sq(); if (norm1 > EPSILON) { coeffs.toSparse(spAlpha); @@ -2196,7 +2201,7 @@ void ist(const Matrix& X, const Matrix& D, } else { coreISTconstrained(G,DtR,coeffs,normX2,lambda,itermax,tol); } - } + } delete[](DtRT); delete[](spAlphaT); @@ -2212,7 +2217,7 @@ inline void generalCD(const AbstractMatrix& G, Vector& DtRv, Vector& co const int K = G.n(); T* const coeffs = coeffsv.rawX(); T* const DtR = DtRv.rawX(); - + for (int iter=0; iter < itermax; ++iter) { if (iter % 5 == 0) { T eps1=DtRv.fmaxval()/lambda-1; @@ -2226,7 +2231,7 @@ inline void generalCD(const AbstractMatrix& G, Vector& DtRv, Vector& co } } eps2=-(eps2/lambda-1); - if (eps2 <= tol) + if (eps2 <= tol) break; } } @@ -2253,7 +2258,7 @@ inline void generalCD(const AbstractMatrix& G, Vector& DtRv, Vector& co template inline void coreIST(const AbstractMatrix& G, Vector& DtRv, Vector& coeffsv, - const T thrs, const int itermax, + const T thrs, const int itermax, const T tol) { const int K = G.n(); @@ -2292,7 +2297,7 @@ inline void coreIST(const AbstractMatrix& G, Vector& DtRv, Vector& coef } } if (iter % 5 == 1) { - vSub(K,DtR,coeffs,DtR); + vSub(K,DtR,coeffs,DtR); maxDtR = DtRv.fmaxval(); norm1 =T(); T DtRa = T(); @@ -2302,7 +2307,7 @@ inline void coreIST(const AbstractMatrix& G, Vector& DtRv, Vector& coef DtRa += DtR[j]*coeffs[j]; } } - vAdd(K,DtR,coeffs,DtR); + vAdd(K,DtR,coeffs,DtR); const T kappa = -DtRa+norm1*maxDtR; if (abs(lambda - maxDtR) < tol && kappa <= tol) break; @@ -2312,7 +2317,7 @@ inline void coreIST(const AbstractMatrix& G, Vector& DtRv, Vector& coef template inline void coreISTW(const Matrix& G, Vector& DtRv, Vector& coeffsv,const Vector& weightsv, - const T lambda, const int itermax, + const T lambda, const int itermax, const T tol) { T opt=0; @@ -2363,7 +2368,7 @@ inline void coreISTW(const Matrix& G, Vector& DtRv, Vector& coeffsv,con /*template inline void coreIST_unnormalized(const AbstractMatrix& G, Vector& DtRv, Vector& coeffsv, - const T thrs, const int itermax, + const T thrs, const int itermax, const T tol) { const int K = G.n(); @@ -2384,7 +2389,7 @@ inline void coreIST_unnormalized(const AbstractMatrix& G, Vector& DtRv, Ve T diff=coeffs[j]; coeffs[j]=DtR[j]-lambda; diff-=coeffs[j]; - + DtR[j]-=diff; G.add_rawCol(j,DtR,diff); } else if (DtR[j] < -lambda) { @@ -2401,7 +2406,7 @@ inline void coreIST_unnormalized(const AbstractMatrix& G, Vector& DtRv, Ve } } if (iter % 5 == 1) { - vSub(K,DtR,coeffs,DtR); + vSub(K,DtR,coeffs,DtR); maxDtR = DtRv.fmaxval(); norm1 =T(); T DtRa = T(); @@ -2506,7 +2511,7 @@ void coreISTconstrained(const AbstractMatrix& G, Vector& DtRv, Vector& /// ist for group Lasso template void ist_groupLasso(const Matrix* XT, const Matrix& D, - Matrix* alphaT, const int Ngroups, + Matrix* alphaT, const int Ngroups, const T lambda, const constraint_type mode, const int itermax, const T tol, const int numThreads) { @@ -2534,7 +2539,7 @@ void ist_groupLasso(const Matrix* XT, const Matrix& D, Matrix* alphatT = new Matrix[NUM_THREADS]; int i; -#pragma omp parallel for private(i) +#pragma omp parallel for private(i) for (i = 0; i< Ngroups; ++i) { #ifdef _OPENMP int numT=omp_get_thread_num(); @@ -2564,14 +2569,14 @@ void ist_groupLasso(const Matrix* XT, const Matrix& D, } else { Vector meanVec(n); X.meanCol(meanVec); - normX2=meanVec.nrm2sq(); + normX2=meanVec.nrm2sq(); coreISTconstrained(G,DtR_mean,coeffs_mean,normX2, lambda,itermax,tol); SpVector spalpha(K); normX2-=computeError(normX2,G,DtR_mean,coeffs_mean,spalpha); normX2=X.normFsq()-M*normX2; } - alphat.fillRow(coeffs_mean); + alphat.fillRow(coeffs_mean); } if (M > 1) { @@ -2638,7 +2643,7 @@ void coreGroupIST(const Matrix& G, Matrix& RtDm, cblas_ger(CblasColMajor,M,K,T(1.0),old_coeff,1,prG+j*K,1,RtD,M); activate[j]=5; } else { - memset(coeffs+j*M,0,M*sizeof(T)); + memset(coeffs+j*M,0,M*sizeof(T)); norms[j]=T(); cblas_ger(CblasColMajor,M,K,T(1.0),old_coeff,1,prG+j*K,1,RtD,M); --activate[j]; @@ -2735,7 +2740,7 @@ void coreGroupISTConstrained(const Matrix& G, Matrix& RtDm, cblas_ger(CblasColMajor,M,K,T(1.0),old_coeff,1,prG+j*K,1,RtD,M); activate[j]=3; } else { - memset(coeffs+j*M,0,M*sizeof(T)); + memset(coeffs+j*M,0,M*sizeof(T)); norms[j]=T(); err += cblas_dot(M,old_coeff,1,old_coeff,1) +2*cblas_dot(M,old_coeff,1,RtD+j*M,1); @@ -2813,7 +2818,7 @@ T computeError(const T normX2,const Vector& norms, return err2; } -/// auxiliary function for +/// auxiliary function for template T computeError(const T normX2, const Matrix& G,const Vector& DtR,const Vector& coeffs, @@ -2823,17 +2828,17 @@ T computeError(const T normX2, }; /* ****************** - * Simultaneous OMP + * Simultaneous OMP * *****************/ template -void somp(const Matrix* X, const Matrix& D, SpMatrix* spalpha, +void somp(const Matrix* X, const Matrix& D, SpMatrix* spalpha, const int Ngroups, const int L, const T eps,const int numThreads) { somp(X,D,spalpha,Ngroups,L,&eps,false,numThreads); } template -void somp(const Matrix* XT, const Matrix& D, SpMatrix* spalphaT, +void somp(const Matrix* XT, const Matrix& D, SpMatrix* spalphaT, const int Ngroups, const int LL, const T* eps, const bool adapt, const int numThreads) { if (LL <= 0) return; @@ -2852,7 +2857,7 @@ void somp(const Matrix* XT, const Matrix& D, SpMatrix* spalphaT, init_omp(numThreads); int i; -#pragma omp parallel for private(i) +#pragma omp parallel for private(i) for (i = 0; i< Ngroups; ++i) { const Matrix& X = XT[i]; const INTM M = X.n(); @@ -2862,7 +2867,7 @@ void somp(const Matrix* XT, const Matrix& D, SpMatrix* spalphaT, Matrix vM; T thrs = adapt ? eps[i] : M*(*eps); coreSOMP(X,D,G,vM,rv,L,thrs); - spalpha.convert2(vM,rv,K); + spalpha.convert2(vM,rv,K); } } @@ -3023,7 +3028,7 @@ void coreSOMP(const Matrix& X, const Matrix& D, const Matrix& G, cblas_gemv(CblasColMajor,CblasNoTrans,K,j+1,T(1.0),prFs,K,prS+j*L,1, T(0.0),prB+j,L); for (int k = 0; k(j,prS[j*L+k],prB+r[k]*L,1,pr_c,1); f.add(tmp,f[currentInd]*invNorm*invNorm); if (j > 0) { @@ -3059,4 +3064,3 @@ void coreSOMP(const Matrix& X, const Matrix& D, const Matrix& G, #endif // DECOMP_H - From a6d027d0d2502454d8003f45db4ff5ed019239d0 Mon Sep 17 00:00:00 2001 From: Samuel St-Jean <3030760+samuelstjean@users.noreply.github.com> Date: Mon, 24 Feb 2025 11:43:26 -0500 Subject: [PATCH 2/9] relax convergence for LassoW like in Lasso regular --- spams_wrap/spams/decomp/decomp.h | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/spams_wrap/spams/decomp/decomp.h b/spams_wrap/spams/decomp/decomp.h index a4e31559..2c8ccec5 100644 --- a/spams_wrap/spams/decomp/decomp.h +++ b/spams_wrap/spams/decomp/decomp.h @@ -2044,9 +2044,11 @@ void coreLARS2W(Vector& DtR, const AbstractMatrix& G, /// L2ERROR const T delta = coeff2*coeff2-coeff1*coeff3; step_max2 = delta < 0 ? INFINITY : (coeff2-sqrt(delta))/coeff1; + step_max2 = MIN(current_correlation,step_max2); } else { /// L1COEFFS step_max2 = coeff1 < 0 ? INFINITY : (constraint-thrs)/coeff1; + step_max2 = MIN(current_correlation,step_max2); } step = MIN(MIN(step,step_max2),step_max); @@ -2101,11 +2103,11 @@ void coreLARS2W(Vector& DtR, const AbstractMatrix& G, newAtom=true; } // Choose next action - if (iter > 4*L || abs(step) < 1e-10 || - step == step_max2 || (normX < 1e-10) || + if (iter > 4*L || abs(step) < 1e-15 || + step == step_max2 || (normX < 1e-15) || (i == (L-1)) || - (mode == L2ERROR && normX - constraint < 1e-10) || - (mode == L1COEFFS && (constraint-thrs < 1e-10))) { + (mode == L2ERROR && normX - constraint < 1e-15) || + (mode == L1COEFFS && (constraint-thrs < 1e-15))) { break; } } From 1a4af446aa46e6c834e92361fd403c5af34be07d Mon Sep 17 00:00:00 2001 From: Samuel St-Jean <3030760+samuelstjean@users.noreply.github.com> Date: Mon, 28 Apr 2025 15:53:40 -0400 Subject: [PATCH 3/9] add weigthed lasso test --- test/test_dictLearn.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/test/test_dictLearn.py b/test/test_dictLearn.py index 398eb45d..cf21619c 100644 --- a/test/test_dictLearn.py +++ b/test/test_dictLearn.py @@ -15,6 +15,8 @@ imgpath = os.path.dirname(os.path.realpath(__file__)) import pytest +from scipy.optimize import nnls + def _extract_lasso_param(f_param): lst = [ 'L','lambda1','lambda2','mode','pos','ols','numThreads','length_path','verbose','cholesky'] l_param = {'return_reg_path' : False} @@ -460,3 +462,32 @@ def test_archetypalAnalysis(myfloat): tac = time.time() t = tac - tic print('time of computation for Robust Archetypal Dictionary Learning: %f' %t) + +@pytest.mark.parametrize("myfloat", [np.float64]) +def test_lasso_weigthed_pos(myfloat): + + myfloat = np.float64 + rng = np.random.default_rng(123456) + m, n, p = 50, 100, 200 + + alpha = rng.standard_normal([n, p]) + alpha[alpha < 0] = 0 + D = rng.standard_normal([m, n]) + D = np.abs(D) * 10 + X = D @ alpha + + X = np.asfortranarray(X).astype(myfloat) + D = np.asfortranarray(D).astype(myfloat) + W = np.ones(alpha.shape, dtype=myfloat, order='F') + + alpha1 = spams.lasso(X, D, lambda1=0.0, pos=True).toarray() + alpha2 = spams.lassoWeighted(X, D, W=W, lambda1=0.0, pos=True).toarray() + alpha3 = [nnls(D, X[:, i])[0] for i in range(X.shape[1])] + alpha3 = np.asarray(alpha3).T + + assert np.all(alpha1 >= 0) + assert np.all(alpha2 >= 0) + assert np.all(alpha3 >= 0) + + np.testing.assert_allclose(alpha1, alpha2) + np.testing.assert_allclose(alpha2, alpha3) From 7429f09702e5cb6f255468078a87e14ece134d51 Mon Sep 17 00:00:00 2001 From: Samuel St-Jean <3030760+samuelstjean@users.noreply.github.com> Date: Mon, 28 Apr 2025 16:08:01 -0400 Subject: [PATCH 4/9] add weigthed lasso test --- test/test_dictLearn.py | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_dictLearn.py b/test/test_dictLearn.py index cf21619c..d1ec1c28 100644 --- a/test/test_dictLearn.py +++ b/test/test_dictLearn.py @@ -471,7 +471,6 @@ def test_lasso_weigthed_pos(myfloat): m, n, p = 50, 100, 200 alpha = rng.standard_normal([n, p]) - alpha[alpha < 0] = 0 D = rng.standard_normal([m, n]) D = np.abs(D) * 10 X = D @ alpha From 1d11e54479c80fe0ab0fea289f5dc407ee02058a Mon Sep 17 00:00:00 2001 From: Samuel St-Jean <3030760+samuelstjean@users.noreply.github.com> Date: Wed, 30 Apr 2025 15:48:36 -0400 Subject: [PATCH 5/9] tag version --- meson.build | 2 +- pyproject.toml | 2 +- test/test_dictLearn.py | 12 ++++++------ 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/meson.build b/meson.build index 0d9dff63..c9c2d8a0 100644 --- a/meson.build +++ b/meson.build @@ -1,6 +1,6 @@ project('_spams_wrap', 'cpp', - version: '2.6.11', + version: '2.6.12', meson_version: '>= 1.2.3', default_options : [ 'warning_level=1', diff --git a/pyproject.toml b/pyproject.toml index e22a7a4d..2f69c507 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -7,7 +7,7 @@ build-backend = "mesonpy" [project] name='spams-bin' requires-python='>=3.9.0' -version='2.6.11' +version='2.6.12' dependencies=['numpy>=1.21.3', 'scipy>=1.5'] license = {text = "GPLv2"} diff --git a/test/test_dictLearn.py b/test/test_dictLearn.py index d1ec1c28..5512c795 100644 --- a/test/test_dictLearn.py +++ b/test/test_dictLearn.py @@ -49,7 +49,7 @@ def test_trainDL(myfloat): except: print("Cannot load image %s : skipping test" %img_file) return None - I = np.array(img) / 255. + I = np.array(img) / 255 if I.ndim == 3: A = np.asfortranarray(I.reshape((I.shape[0],I.shape[1] * I.shape[2]))) rgb = True @@ -141,7 +141,7 @@ def test_trainDL_Memory(myfloat): except: print("Cannot load image %s : skipping test" %img_file) return None - I = np.array(img) / 255. + I = np.array(img) / 255 if I.ndim == 3: A = np.asfortranarray(I.reshape((I.shape[0],I.shape[1] * I.shape[2]))) rgb = True @@ -199,7 +199,7 @@ def test_structTrainDL(myfloat): except Exception as e: print("Cannot load image %s (%s) : skipping test" %(img_file,e)) return None - I = np.array(img) / 255. + I = np.array(img) / 255 if I.ndim == 3: A = np.asfortranarray(I.reshape((I.shape[0],I.shape[1] * I.shape[2])),dtype = myfloat) rgb = True @@ -371,7 +371,7 @@ def test_nmf(myfloat): except: print("Cannot load image %s : skipping test" %img_file) return None - I = np.array(img) / 255. + I = np.array(img) / 255 if I.ndim == 3: A = np.asfortranarray(I.reshape((I.shape[0],I.shape[1] * I.shape[2])),dtype = myfloat) rgb = True @@ -406,7 +406,7 @@ def test_archetypalAnalysis(myfloat): except Exception as e: print("Cannot load image %s (%s) : skipping test" %(img_file,e)) return None - I = np.array(img) / 255. + I = np.array(img) / 255 if I.ndim == 3: A = np.asfortranarray(I.reshape((I.shape[0],I.shape[1] * I.shape[2])),dtype = myfloat) rgb = True @@ -463,7 +463,7 @@ def test_archetypalAnalysis(myfloat): t = tac - tic print('time of computation for Robust Archetypal Dictionary Learning: %f' %t) -@pytest.mark.parametrize("myfloat", [np.float64]) +@pytest.mark.parametrize("myfloat", [np.float32, np.float64]) def test_lasso_weigthed_pos(myfloat): myfloat = np.float64 From c70f0785f0e59646072b5dbd39758b07d8bcf217 Mon Sep 17 00:00:00 2001 From: Samuel St-Jean <3030760+samuelstjean@users.noreply.github.com> Date: Wed, 30 Apr 2025 15:50:47 -0400 Subject: [PATCH 6/9] fix typo --- test/test_dictLearn.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_dictLearn.py b/test/test_dictLearn.py index 5512c795..f361483e 100644 --- a/test/test_dictLearn.py +++ b/test/test_dictLearn.py @@ -464,7 +464,7 @@ def test_archetypalAnalysis(myfloat): print('time of computation for Robust Archetypal Dictionary Learning: %f' %t) @pytest.mark.parametrize("myfloat", [np.float32, np.float64]) -def test_lasso_weigthed_pos(myfloat): +def test_lasso_weighted_pos(myfloat): myfloat = np.float64 rng = np.random.default_rng(123456) From abfc0743e431f90cfa4f243f3f6a9fd66b26f187 Mon Sep 17 00:00:00 2001 From: Samuel St-Jean <3030760+samuelstjean@users.noreply.github.com> Date: Thu, 1 May 2025 08:23:45 -0400 Subject: [PATCH 7/9] Update test_dictLearn.py --- test/test_dictLearn.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/test_dictLearn.py b/test/test_dictLearn.py index f361483e..b937a9d4 100644 --- a/test/test_dictLearn.py +++ b/test/test_dictLearn.py @@ -465,8 +465,6 @@ def test_archetypalAnalysis(myfloat): @pytest.mark.parametrize("myfloat", [np.float32, np.float64]) def test_lasso_weighted_pos(myfloat): - - myfloat = np.float64 rng = np.random.default_rng(123456) m, n, p = 50, 100, 200 From 67911b9bbdd9a22fe017bf517981f4066152999d Mon Sep 17 00:00:00 2001 From: Samuel St-Jean <3030760+samuelstjean@users.noreply.github.com> Date: Thu, 1 May 2025 08:33:57 -0400 Subject: [PATCH 8/9] Update test_dictLearn.py --- test/test_dictLearn.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_dictLearn.py b/test/test_dictLearn.py index b937a9d4..52b84bc6 100644 --- a/test/test_dictLearn.py +++ b/test/test_dictLearn.py @@ -480,7 +480,7 @@ def test_lasso_weighted_pos(myfloat): alpha1 = spams.lasso(X, D, lambda1=0.0, pos=True).toarray() alpha2 = spams.lassoWeighted(X, D, W=W, lambda1=0.0, pos=True).toarray() alpha3 = [nnls(D, X[:, i])[0] for i in range(X.shape[1])] - alpha3 = np.asarray(alpha3).T + alpha3 = np.asarray(alpha3, dtype=myfloat).T assert np.all(alpha1 >= 0) assert np.all(alpha2 >= 0) From 3e6331d89d4c0ded715b527f060d5b21b0603d29 Mon Sep 17 00:00:00 2001 From: Samuel St-Jean <3030760+samuelstjean@users.noreply.github.com> Date: Thu, 1 May 2025 09:08:45 -0400 Subject: [PATCH 9/9] fix precision test issue --- test/test_dictLearn.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/test/test_dictLearn.py b/test/test_dictLearn.py index 52b84bc6..ca519e75 100644 --- a/test/test_dictLearn.py +++ b/test/test_dictLearn.py @@ -480,11 +480,15 @@ def test_lasso_weighted_pos(myfloat): alpha1 = spams.lasso(X, D, lambda1=0.0, pos=True).toarray() alpha2 = spams.lassoWeighted(X, D, W=W, lambda1=0.0, pos=True).toarray() alpha3 = [nnls(D, X[:, i])[0] for i in range(X.shape[1])] - alpha3 = np.asarray(alpha3, dtype=myfloat).T + alpha3 = np.asarray(alpha3).T assert np.all(alpha1 >= 0) assert np.all(alpha2 >= 0) assert np.all(alpha3 >= 0) np.testing.assert_allclose(alpha1, alpha2) - np.testing.assert_allclose(alpha2, alpha3) + + # only tested on double precision because nnls casts internally to double + # while spams runs at the given precision + if myfloat is np.float64: + np.testing.assert_allclose(alpha2, alpha3)