def test_lasso():
    np.random.seed(0)
    print("test lasso")
##############################################
# Decomposition of a large number of signals
##############################################
# data generation
    X = np.asfortranarray(np.random.normal(size=(100,100000)))
    #* X=X./repmat(sqrt(sum(X.^2)),[size(X,1) 1]);
    X = np.asfortranarray(X / np.tile(np.sqrt((X*X).sum(axis=0)),(X.shape[0],1)),dtype= myfloat)
    D = np.asfortranarray(np.random.normal(size=(100,200)))
    D = np.asfortranarray(D / np.tile(np.sqrt((D*D).sum(axis=0)),(D.shape[0],1)),dtype= myfloat)
    # parameter of the optimization procedure are chosen
#param.L=20; # not more than 20 non-zeros coefficients (default: min(size(D,1),size(D,2)))
    param = {
        'lambda1' : 0.15, # not more than 20 non-zeros coefficients
        'numThreads' : -1, # number of processors/cores to use; the default choice is -1
        # and uses all the cores of the machine
        'mode' : spams.PENALTY}        # penalized formulation

    tic = time.time()
    alpha = spams.lasso(X,D = D,return_reg_path = False,**param)
    tac = time.time()
    t = tac - tic
    print("%f signals processed per second\n" %(float(X.shape[1]) / t))
########################################
# Regularization path of a single signal
########################################
    X = np.asfortranarray(np.random.normal(size=(64,1)),dtype= myfloat)
    D = np.asfortranarray(np.random.normal(size=(64,10)))
    D = np.asfortranarray(D / np.tile(np.sqrt((D*D).sum(axis=0)),(D.shape[0],1)),dtype= myfloat)
    (alpha,path) = spams.lasso(X,D = D,return_reg_path = True,**param)
    return None

def test_trainDL_Memory():
    img_file = 'lena.png'
    try:
        img = Image.open(img_file)
    except:
        print("Cannot load image %s : skipping test" %img_file)
        return None
    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
    else:
        A = np.asfortranarray(I)
        rgb = False

    m = 8;n = 8;
    X = spams.im2col_sliding(A,m,n,rgb)

    X = X - np.tile(np.mean(X,0),(X.shape[0],1))
    X = np.asfortranarray(X / np.tile(np.sqrt((X * X).sum(axis=0)),(X.shape[0],1)))
    X = np.asfortranarray(X[:,np.arange(0,X.shape[1],10)],dtype = myfloat)

    param = { 'K' : 200, # learns a dictionary with 100 elements
          'lambda1' : 0.15, 'numThreads' : 4,
          'iter' : 100}

    ############# FIRST EXPERIMENT  ##################
    tic = time.time()
    D = spams.trainDL_Memory(X,**param)
    tac = time.time()
    t = tac - tic
    print('time of computation for Dictionary Learning: %f' %t)

    print('Evaluating cost function...')
    lparam = _extract_lasso_param(param)
    alpha = spams.lasso(X,D = D,**lparam)
    xd = X - D * alpha
    R = np.mean(0.5 * (xd * xd).sum(axis=0) + param['lambda1'] * np.abs(alpha).sum(axis=0))
    print("objective function: %f" %R)
    #* ? DISPLAY

    ############# SECOND EXPERIMENT  ##################
    tic = time.time()
    D = spams.trainDL(X,**param)
    tac = time.time()
    t = tac - tic
    print('time of computation for Dictionary Learning: %f' %t)
    print('Evaluating cost function...')
    alpha = spams.lasso(X,D = D,**lparam)
    xd = X - D * alpha
    R = np.mean(0.5 * (xd * xd).sum(axis=0) + param['lambda1'] * np.abs(alpha).sum(axis=0))
    print("objective function: %f" %R)

    #* ? DISPLAY

    return None

def test_cd():
    np.random.seed(0)
    X = np.asfortranarray(np.random.normal(size = (64,100)))
    X = np.asfortranarray(X / np.tile(np.sqrt((X*X).sum(axis=0)),(X.shape[0],1)),dtype=myfloat)
    D = np.asfortranarray(np.random.normal(size = (64,100)))
    D = np.asfortranarray(D / np.tile(np.sqrt((D*D).sum(axis=0)),(D.shape[0],1)),dtype=myfloat)
    # parameter of the optimization procedure are chosen
    lambda1 = 0.015
    mode = spams.PENALTY
    tic = time.time()
    alpha = spams.lasso(X,D,lambda1 = lambda1,mode = mode,numThreads = 4)
    tac = time.time()
    t = tac - tic
    xd = X - D * alpha
    E = np.mean(0.5 * (xd * xd).sum(axis=0) + lambda1 * np.abs(alpha).sum(axis=0))
    print("%f signals processed per second for LARS" %(X.shape[1] / t))
    print('Objective function for LARS: %g' %E)
    tol = 0.001
    itermax = 1000
    tic = time.time()
#    A0 = ssp.csc_matrix(np.empty((alpha.shape[0],alpha.shape[1])))
    A0 = ssp.csc_matrix((alpha.shape[0],alpha.shape[1]),dtype=myfloat)
    alpha2 = spams.cd(X,D,A0,lambda1 = lambda1,mode = mode,tol = tol, itermax = itermax,numThreads = 4)
    tac = time.time()
    t = tac - tic
    print("%f signals processed per second for CD" %(X.shape[1] / t))
    xd = X - D * alpha2
    E = np.mean(0.5 * (xd * xd).sum(axis=0) + lambda1 * np.abs(alpha).sum(axis=0))
    print('Objective function for CD: %g' %E)
    print('With Random Design, CD can be much faster than LARS')

    return None

def fit_strf_lasso(input, output, lags, lambda1=1.0, lambda2=1.0, num_threads=10):

    #convert the input into a toeplitz-like matrix
    stime = time.time()
    A = make_toeplitz(input, lags, include_bias=True, fortran_style=True)
    etime = time.time() - stime
    print '[fit_strf_lasso] Time to make Toeplitz matrix: %d seconds' % etime

    fy = np.asfortranarray(output.reshape(len(output), 1))
    #print 'fy.shape=',fy.shape
    #print 'fA.shape=',fA.shape

    #fit the STRF
    stime = time.time()
    fit_params = spams.lasso(fy, A, mode=2, lambda1=lambda1, lambda2=lambda2, numThreads=num_threads)
    etime = time.time() - stime
    print '[fit_strf_lasso] Time to fit STRF: %d seconds' % etime

    #reshape the STRF so that it makes sense
    nt = input.shape[0]
    nf = input.shape[1]
    d = len(lags)
    strf = np.array(fit_params[:-1].todense()).reshape([nf, d])
    bias = fit_params[-1].todense()[0, 0]

    return strf,bias

    def fit( self, y, dirs, KERNELS, params ) :
        nD = dirs.shape[0]
        n1 = len(self.Rs)
        n2 = len(self.ICVFs)
        n3 = len(self.d_ISOs)

        # prepare DICTIONARY from dirs and lookup tables
        A = np.zeros( (len(y), nD*(n1+n2)+n3 ), dtype=np.float64, order='F' )
        o = 0
        for i in xrange(nD) :
            i1, i2 = amico.lut.dir_TO_lut_idx( dirs[i] )
            A[:,o:(o+n1)] = KERNELS['wmr'][:,i1,i2,:].T
            o += n1
        for i in xrange(nD) :
            i1, i2 = amico.lut.dir_TO_lut_idx( dirs[i] )
            A[:,o:(o+n2)] = KERNELS['wmh'][:,i1,i2,:].T
            o += n2
        A[:,o:] = KERNELS['iso'].T

        # empty dictionary
        if A.shape[1] == 0 :
            return [0, 0, 0], None, None, None

        # fit
        x = spams.lasso( np.asfortranarray( y.reshape(-1,1) ), D=A, **params ).todense().A1

        # return estimates
        f1 = x[ :(nD*n1) ].sum()
        f2 = x[ (nD*n1):(nD*(n1+n2)) ].sum()
        v = f1 / ( f1 + f2 + 1e-16 )
        xIC = x[:nD*n1].reshape(-1,n1).sum(axis=0)
        a = 1E6 * 2.0 * np.dot(self.Rs,xIC) / ( f1 + 1e-16 )
        d = (4.0*v) / ( np.pi*a**2 + 1e-16 )
        return [v, a, d], dirs, x, A

    def fit( self, y, dirs, KERNELS, params ) :
        nD = dirs.shape[0]
        if nD > 1 : # model works only with one direction
            raise RuntimeError( '"%s" model requires exactly 1 orientation' % self.name )

        n1 = len(self.d_perps)
        n2 = len(self.d_isos)
        nATOMS = n1+n2
        if nATOMS == 0 : # empty dictionary
            return [0, 0], None, None, None

        # prepare DICTIONARY from dir and lookup tables
        i1, i2 = amico.lut.dir_TO_lut_idx( dirs[0] )
        A = np.zeros( (len(y), nATOMS), dtype=np.float64, order='F' )
        A[:,:(nD*n1)] = KERNELS['D'][:,i1,i2,:].T
        A[:,(nD*n1):] = KERNELS['CSF'].T

        # fit
        x = spams.lasso( np.asfortranarray( y.reshape(-1,1) ), D=A, **params ).todense().A1

        # return estimates
        v = x[ :n1 ].sum() / ( x.sum() + 1e-16 )

        # checking that there is more than 1 isotropic compartment
        if self.type == 'Mouse' :
            v_blood = x[ n1 ] / ( x.sum() + 1e-16 )
            v_csf = x[ n1+1 ] / ( x.sum() + 1e-16 )

            return [ v, 1-v, v_blood, v_csf ], dirs, x, A

        else :
            return [ v, 1-v ], dirs, x, A

def dictEval( X, D, param, lam=None, dsfactor=None, patchSize=None, patchFnGrp=None, kind='avg'):
    if dsfactor is not None:
        X_useme,dsz  = downsamplePatchList( X, patchSize, dsfactor, kind=kind )
        D_useme,Ddsz = downsamplePatchList( D, patchSize, dsfactor, kind=kind )

        if patchFnGrp:
            patchFnGrp.create_dataset('patchesDown', data=X_useme)
    else:
        X_useme = X
        D_useme = D

    if lam is None:
        lam = param['lambda1']

    alpha = spams.lasso( np.asfortranarray(X_useme), D = np.asfortranarray(D_useme), **param )
    Xre = ( D * alpha )

    if patchFnGrp:
        patchFnGrp.create_dataset('patchesRecon', data=Xre)

    xd = X - Xre 

    R = np.mean( (xd * xd).sum(axis=0))

    if lam > 0:
        print "   dictEval - lambda: ", lam
        R = R + lam * np.mean( np.abs(alpha).sum(axis=0))

    return R

    def predict(self, imgs, neuron_idx=None, penalty_lambda=None, algorithm=None):
        """ get neuron response to images

        Parameters
        ----------
        imgs

        Returns
        -------

        """
        imgs_array = make_2d_input_matrix(imgs)
        if neuron_idx is None:
            dict_to_use = self.w
        else:
            dict_to_use = self.w[neuron_idx:(neuron_idx + 1), :]

        if penalty_lambda is None:
            _lambda = self._lambda
        else:
            _lambda = penalty_lambda
        assert np.isscalar(_lambda)

        if algorithm is None:
            _algorithm = self.algorithm
        else:
            _algorithm = algorithm


        # let's call sparse encoder to do it!
        # no scaling at all!
        # having /nsample in objective function is exactly the same as sovling each problem separately.
        # the underlying function called is elastic net, and that function fits each column of y separately.
        # each column of y is each stimulus. This is because when passing imgs_array and dict_to_use to Elastic Net,
        # they are transposed. That is, y = imgs_array.T
        #
        # in the code there's also a subtle detail, where alpha is divided by number of pixels in each stimulus.
        # I haven't figured that out, but seems that's simply a detail for using ElasticNet to do this.
        if _algorithm in ['lasso_lars', 'lasso_cd']:
            response = sparse_encode(imgs_array, dict_to_use, alpha=_lambda, algorithm=_algorithm, max_iter=10000)
        else:
            assert _algorithm == 'spams'
            #print(imgs_array.dtype, dict_to_use.dtype, _lambda.shape)
            response = lasso(np.asfortranarray(imgs_array.T), D=np.asfortranarray(dict_to_use.T), lambda1=_lambda,
                             mode=2)
            response = response.T.toarray()  # because lasso returns sparse matrix...
        # this can be used for debugging, for comparison with SPAMS.
        # notice here I give per sample cost.
        self.last_cost_recon = 0.5 * np.sum((imgs_array - np.dot(response, dict_to_use)) ** 2, axis=1)
        self.last_cost_sparsity = _lambda * np.abs(response).sum(axis=1)
        assert self.last_cost_sparsity.shape == (imgs_array.shape[0], )
        assert self.last_cost_recon.shape == (imgs_array.shape[0],)
        self.last_cost = np.mean(self.last_cost_recon + self.last_cost_sparsity)

        return response

    def __call__(self, X, D):
        import spams
        lasso_params = {
            'lambda1': self._lambda,
            'lambda2': 0,
            'numThreads': self.n_jobs,
            'mode': 2
        }

        return np.array(spams.lasso(np.asfortranarray(X, np.float64), D=np.asfortranarray(D, np.float64),
                                    return_reg_path=False, **lasso_params).todense())

def sparseCoding(X):
    X = np.asfortranarray(X)
    param = { 'K' : NCLUSTER,	# size of the dictionary 
          'lambda1' : 0.15, 
          #'posD' : True,	# dictionary positive constrain
          #'modeD' : 1,	# L1 regulization regularization on D
          'iter' : ITER} # runtime limit 15mins
    
    D = spams.trainDL_Memory(X,**param)
    lparam = _extract_lasso_param(param)
    print 'genrating codes...'
    alpha = spams.lasso(X,D = D,**lparam)
    return D, alpha

def findSparseRep(A, y, epsilon):
    # A: data in columns
    # y: image to match
    # epsilon: allowed error
    yCol = np.asfortranarray(np.reshape(y, (len(y), 1)), dtype=float)
    A = np.asfortranarray(np.array(A), dtype=float)
    param = {'lambda1' : epsilon,
            'numThreads' : 4,
            'pos' : True,
            'mode' : 1}
    x = sp.lasso(yCol, D=A, **param).toarray()
    residuals = [np.linalg.norm(A[:,i]*x[i] - y) for i in range(0,np.shape(A)[1])]
    return (x, residuals)

 def fit(self, x, y):
     self.coef_ = (
         spams.lasso(
             np.asfortranarray(y, dtype=np.float),
             D=np.asfortranarray(x, dtype=np.float),
             pos=self.positive,
             lambda1=self.lambda_,
             **self.params
         )
         .toarray()
         .flatten()
     )
     return self

    def coding_series(self, segment_list, D, a=None, batch=False, iter=-5):
        k = D.shape[1]
        bow_data = numpy.zeros([len(segment_list), k])

        # BoW for data
        for index, item in enumerate(segment_list):
            code = spams.lasso(numpy.asfortranarray(item), D, lambda1=a, pos=True)
            code = numpy.sum(code.todense(), axis=1)
            bow_data[index : index + 1, :] += code.reshape([1, k])
            div = numpy.linalg.norm(bow_data[index, :])
            if div > 0:
                bow_data[index, :] = bow_data[index, :] / div
        return bow_data

    def solve(self, A, y, x0, as_signs):

        #print 'dense_solver: y.shape=',y.shape
        #print 'dense_solver: A.shape=',A.shape
        fy = np.asfortranarray(y.reshape(len(y), 1))
        fA = np.asfortranarray(A)
        #print 'fy.shape=',fy.shape
        #print 'fA.shape=',fA.shape
        xnew = lasso(fy, fA, mode=2, lambda1=self.lambda1, lambda2=self.lambda2)
        xnew = np.array(xnew.todense()).reshape(x0.shape)

        #print 'dense_solver: xnew.shape=',xnew.shape

        return xnew

def _objective(X,D,param,imgname = None):
    print('Evaluating cost function...')
    lparam = _extract_lasso_param(param)
    alpha = spams.lasso(X,D = D,**lparam)
    # NB : as alpha is sparse, D*alpha is the dot product
    xd = X - D * alpha
    R = np.mean(0.5 * (xd * xd).sum(axis=0) + param['lambda1'] * np.abs(alpha).sum(axis=0))
    print("objective function: %f" %R)
    #* display ?
    if imgname != None:
        img = spams.displayPatches(D)
        print("IMG %s" %str(img.shape))
        x = np.uint8(img[:,:,0] * 255.)
        image = Image.fromarray(x,mode = 'L')
        image.save("%s.png" %imgname)

def generateCode(X,D):
    X = np.asfortranarray(X)
    print X.shape
    D = np.asfortranarray(D)
    print D.shape
    param = { 'K' : NCLUSTER,	# size of the dictionary 
          'lambda1' : 0.05, 
          #'posD' : True,	# dictionary positive constrain
          #'modeD' : 1,	# L1 regulization regularization on D
          'iter' : ITER} # runtime limit 15mins

    lparam = _extract_lasso_param(param)
    print 'genrating codes...'
    alpha = spams.lasso(X,D = D,**lparam)
    return alpha

def bow(trajectory, D, a, w_len, interval=1):
    k = D.shape[1]
    histogram = numpy.zeros([len(trajectory), k])
    for index, item in enumerate(trajectory):
        temp = list()
        stamp_index = range(0, item.shape[0]-w_len+1, interval)
        for i in stamp_index:
            temp.append(item[i:i+w_len, :1])
        temp = numpy.hstack(temp)
        code = spams.lasso(numpy.asfortranarray(temp), D, lambda1=a, pos=True)
        code = numpy.sum(code.todense(), axis=1)
        histogram[index:index+1, :] += code.reshape([1, k])
        div = numpy.linalg.norm(histogram[index, :])
        if div > 0:
            histogram[index, :] = histogram[index, :] / div
    return histogram

 def doFeatureEncoding(self, features):
     """ do feature encoding to original features"""
     encodedFeatures = None
     whitenedFeatures = whiten(features)
     
     if self._featureEncodingMethod == 'vector-quantization':
         # Vector quantization
         # each row is a feature vector
         index, _ = vq(whitenedFeatures, self._codebook)
         row, _ = features.shape
         col = config.codebookSize
         encodedFeatures = np.zeros((row, col))
         
         for i in xrange(len(index)):
             encodedFeatures[i, index[i]] = 1
             
     elif self._featureEncodingMethod == 'sparse-coding':
         # Sparse coding
         # each column is a feature vector
         X = np.asfortranarray(whitenedFeatures.transpose())
         X = np.asfortranarray(X / np.tile(np.sqrt((X*X).sum(axis=0)),
                                           (X.shape[0],1)),
                               dtype= X.dtype)
         D = np.asfortranarray(self._codebook.transpose())
         D = np.asfortranarray(D / np.tile(np.sqrt((D*D).sum(axis=0)),
                                           (D.shape[0],1)), 
                               dtype = D.dtype)
         
         # Parameters of the optimization are chosen
         param = {
             'lambda1': 0.15, 
             'numThreads': -1,
             'mode': 0    
             }
         
         alpha = spams.lasso(X, D, **param)   # alpha is sparse matrix
         
         alphaShape = (D.shape[1], X.shape[1])
         denseMatrix = csc_matrix(alpha, shape = alphaShape).todense()
         encodedFeatures = np.asarray(denseMatrix).transpose()
             
     return encodedFeatures

    def solve(self, A, y, x0, as_signs):

        #print 'dense_solver: y.shape=',y.shape
        #print 'dense_solver: A.shape=',A.shape
        fy = np.asfortranarray(y.reshape(len(y), 1))
        fA = np.asfortranarray(A)
        #print 'fy.shape=',fy.shape
        #print 'fA.shape=',fA.shape
        if not self.positive:
            xnew = spams.lasso(fy, fA, mode=2, lambda1=self.lambda1, lambda2=self.lambda2)
        else:
            W = np.ones(len(x0))
            params = {'lambda1':self.lambda1, 'pos':True}
            xnew = spams.lassoWeighted(fy, fA, W, **params)

        xnew = np.array(xnew.todense()).reshape(x0.shape)

        #print 'dense_solver: xnew.shape=',xnew.shape

        return xnew

    # for each fixed dictionary K, we will repeat dictionary
    # learning for 100 times, each with a different initial value
    test_cases = 10
    R = np.zeros((test_cases, ))
    for i in xrange(0, test_cases):
        if randomStart == 1:
            D0 = util.dictLearnInit(X, K, 'random', 0)
        param['D'] = np.asfortranarray(D0)
        lparam = {'lambda1': Lambda,
                  'pos': True,
                  'mode': 2,
                  'numThreads': -1
        }
        D = spams.trainDL(X, **param)
        alpha = spams.lasso(X, D, **lparam)
        R[i] = np.mean(0.5 * sum((X - D * alpha) ** 2) + param['lambda1'] * sum(abs(alpha)))
        print R[i]

        #(permInd, cost) = util.dissimilarityDict(Dtemplate, D, 'euclidean')
        #D = D[:, permInd] # permute the columns of D to match the template Dtemplate

        if i >= 1:
            if R[i] < Rbest:
                Dbest = D
                Rbest = R[i]
        else:
            Dbest = D
            Rbest = R[0]

    # print path