Demo entry 6341048

怡然代码

   

Submitted by anonymous on Dec 30, 2016 at 06:35
Language: Python. Code size: 5.9 kB.

import numpy as np
import numpy.ma as ma
import math
import csv
from scipy import linalg



    ##定义函数

class RESULT(object):
    def __init__(self,xurl,yurl,lamb2=0):
        self.x=np.loadtxt(xurl,delimiter=',',skiprows=0)
        self.y=np.loadtxt(yurl,delimiter=',',skiprows=0)
        self.n=self.x.shape[0]
        self.p=self.x.shape[1]
        self.lamb2=lamb2
        self.xy=np.dot(self.x.T,self.y)
        self.yy=np.inner(self.y,self.y)
        self.xx=np.dot(self.x.T,self.x)
        self.re=[]
        self.index=[]
    
    ##定义变量
    def lar(self):
        beta=np.zeros(self.p,dtype=float)
        df=0
        RSS=self.yy
        A=np.array([False]*self.p)
        cc=self.xy
        cc_abs=np.absolute(cc)
        CCI=np.argmax(cc_abs)
        CC=cc_abs[CCI]
        lamb1=CC
        A[CCI]=True
        self.index.append(CCI)
        self.re.append([np.copy(A),np.copy(beta),df,RSS,lamb1,self.index])
        XA=np.copy(self.x[:,A])
        LA=np.array([[np.sqrt(np.inner(XA[:,0],XA[:,0])+self.lamb2)]],dtype=float)
       
    ##循环内
        while(True):
            print [np.copy(A),np.copy(beta),RSS,lamb1,self.index]
            S=np.sign(cc)
            LLA=np.sum(A)
            SA=[] 
            for i in range(LLA):
                SA.append(S[self.index[i]])
            SA=np.array(SA)
            SA1=SA.reshape(LLA,1)
            SAM=np.repeat(SA1,LLA,axis=1)
       
            LA=LA*SAM
            one=np.ones(np.sum(A),dtype=float)
            lb = linalg.solve_triangular(LA,one,lower=True)
            RA=linalg.solve_triangular(LA.T,lb)
            AA=1/np.sqrt((np.dot(one,RA)))
            WA=AA*RA
            UA=SA*WA
            XA=[]
            index2=np.array(np.copy(self.index)).reshape((LLA,1))
            for i in range(LLA):
                XA.append(self.x[:,index2[i][0]])
            XA=np.array(XA).reshape((LLA,self.n))
            a=np.dot(self.x.T,np.dot(XA.T,UA))
            gam=np.zeros(self.p,dtype=float)
            LA=LA*SAM
            
    ##接下来是赋值的阶段,针对两种集合,即活跃的集合和不活跃的集合      
    ##首先,处理不活跃的集合
    ##将更新的步长进行赋值
    
            if(LLA<self.p):
                gamN=(CC+cc[~A])/(AA+a[~A])
                gamP=(CC-cc[~A])/(AA-a[~A])
                gamt=np.where(((gamP>gamN)&(gamN>0))|((gamN>0)&(gamP<0)),gamN,gamP)
                gam[~A]=gamt
                
    ##然后,处理活的跃集合
    ##将更新的步长进行赋值
    
            if(LLA>1):
                beta_pos=[]
                for i in range(LLA):
                    beta_pos.append(beta[self.index[i]])
                beta_pos=np.array(beta_pos)
                gamt=-beta_pos/(WA*SA)
                for i in range(LLA):
                    gam[self.index[i]]=gamt[i]
            if(LLA==self.p):
                gamCA=CC/AA
                gam=np.where(gam<gamCA,gam,0)
            if(np.all(gam<=0)):
                gammin=CC/AA
            else:
                gam1=ma.masked_array(gam, mask=(gam<=0))
                gammina=np.argmin(gam1)
                gammin=gam1[gammina]
                
    ##完成变量选择的过程,从而找出下一步可以进行更新的变量      
    ##对于参数的进一步update
    
            updater=gammin*WA*SA
            for i in range(len(updater)):
                beta[self.index[i]]=beta[self.index[i]]+updater[i]
            df=LLA
            xxb=np.dot(self.xx,beta)
            RSS=self.yy-2*np.inner(self.xy,beta)+np.inner(xxb,beta)
            
    ##正式进入引入变量的过程
    ##该过程使用增量算法的新方法
            
            if(gammina not in self.index):    
                x1=self.x[:,gammina]
                xx1_pre=[]
                for i in range(LLA):
                    xx1_pre.append(self.x[:,self.index[i]])
                xx1_pre=np.array(xx1_pre).T
                xx1=np.dot(xx1_pre.T,x1)
                l1=linalg.solve_triangular(LA,xx1,lower=True)
                l11=np.sqrt(np.dot(x1.T,x1)+self.lamb2-np.dot(l1.T,l1))
                l1=np.hstack((l1,l11)).reshape(1,LLA+1)
                LA=np.column_stack([LA,np.zeros(LLA)])  
                LA=np.row_stack([LA,l1])
                self.index.append(gammina)
                
    ##进入变量退出的过程
    ##该过程我们使用Givens算法,简便易行
            
            else:
                index_minus=range(LLA)
                index_minus.remove(self.index.index(gammina))
                LA=LA[index_minus]
                IA=index_minus
                LLA=LLA-1
                for i in range(LLA):
                    ia=IA[i]
                    if(i!=ia):
                        v=np.copy(LA[i,i:(ia+1)])
                        lv=np.linalg.norm(v)
                        LA[i,i]=lv
                        LA[i,(i+1):(ia+1)]=0                    
                        if(i<(LLA-1)):
                            TLA=LA[(i+1):,i:(ia+1)]
                            Epsilon_1=v[0]
                            for j in range(ia-i):
                                Epsilon_2=v[j+1]
                                mod=math.sqrt(pow(Epsilon_1,2)+pow(Epsilon_2,2))
                                c=Epsilon_2/mod
                                s=Epsilon_1/mod
                                Epsilon_1=mod
                                TLA[:,[0,j+1]]=np.dot(TLA[:,[0,j+1]],np.array([[s,c],[c,-s]]))                       
                LA=LA[:,:LLA]
                self.index.remove(gammina)

    ##分别对lambda,CC,cc,beta进行更新
    ##循环中进行判断
               
            if(np.any(gam>0)):
                A[gammina]=~A[gammina]
                cc=self.xy-xxb-self.lamb2*beta
                CC=abs(cc[gammina]) 
                lamb1=CC
            else:
                lamb1=0
            self.re.append([np.copy(A),np.copy(beta),RSS,lamb1,self.index])
            
            if(np.all(gam<=0)):break
                
        return self.re

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).