swh:1:snp:3aec91c51d538d62f3f51f6a0af59fe452f330ab
Tip revision: 71c7015946b000eca777d0192d64b689d8c1732d authored by Zhu Wang on 19 February 2020, 20:10:02 UTC
version 0.3-22
version 0.3-22
Tip revision: 71c7015
nclreg_ad.f
C used in nclreg.R for active=TRUE and decreasing=TRUE
C output: nlambdacal: number of lambda actually computed including
C repeated ones
C
subroutine nclreg_ad(x, y, weights, n,m,start, etastart,
+ mustart, offset, iter, nlambda, lambda, alpha, gam,
+ standardize, intercept, penaltyfactor, maxit, eps, epscycle,
+ penalty, trace, del,rfamily, B, s, thresh, cost,
+ decreasing, beta, b0, yhat, los, pll, nlambdacal)
implicit none
integer n,m,i,ii,j,jj,penalty,nlambda, standardize, maxit,
+ trace, iter, rfamily, jk, activeset(m), intercept,
+ activeset_old(m), m_act, nlambdacal, uturn, decreasing,
+ cutpoint, nact, conv, jc, fakejk, AllocateStatus
double precision x(n, m), y(n), weights(n),start(m+1),etastart(n),
+ mustart(n), offset(n), lambda(nlambda), alpha, gam, eps,
+ penaltyfactor(m), thresh, beta(m, nlambda), epscycle,
+ betaall(m), b0all, b0(nlambda), b0_1, yhat(n), del,
+ lambda_i, s, B, fk(n), los(nlambda),
+ pll(nlambda), cost, penval
double precision, dimension(:, :), allocatable :: x_act
double precision, dimension(:), allocatable :: start_act, beta_1,
+ penaltyfactor_act
b0all=start(1)
do ii=1, m
activeset(ii)=ii
activeset_old(ii)=ii
betaall(ii)=start(ii+1)
enddo
call find_activeset(m, start(2:(m+1)), eps, activeset, jk)
fakejk=0
if(jk==0)then
jk=1
activeset(1)=1
fakejk=1
endif
do ii=1, jk
activeset_old(ii)=activeset(ii)
enddo
m_act = jk
allocate(start_act(m_act+1), stat=AllocateStatus)
if(AllocateStatus .NE. 0)then
return
endif
allocate(penaltyfactor_act(m_act), stat=AllocateStatus)
allocate(x_act(n, m_act), stat=AllocateStatus)
if(AllocateStatus .NE. 0)then
return
endif
do jj=1, n
do ii=1, m_act
x_act(jj, ii)=x(jj, activeset(ii))
enddo
enddo
allocate(beta_1(m_act), stat=AllocateStatus)
if(AllocateStatus .NE. 0)then
return
endif
start_act(1)=start(1)
do j=1, m_act
beta_1(j)=0
start_act(j+1)=start(1+activeset(j))
penaltyfactor_act(j)=penaltyfactor(activeset(j))
enddo
i=1
nlambdacal=0
uturn=0
cutpoint=1
10 if(i .LE. nlambda)then
if(trace .EQ. 1)then
call intpr("nclreg_ad lambda iteration i=", -1, i, 1)
endif
lambda_i=lambda(i)/B
j=1
nact=2
conv=0
13000 if(j <= nact .AND. conv==0)then
if(trace .EQ. 1)then
call intpr("begin activeset nclreg_onelambda", -1, 1, 1)
endif
call nclreg_onelambda(x_act, y,weights, n,m_act,start_act,
+ etastart, mustart, yhat, offset, lambda_i, alpha, gam,
+ intercept, penaltyfactor_act, maxit, eps, penalty,
+ trace,iter,del,rfamily, B, s, thresh, beta_1, b0_1, fk)
start(1)=b0_1
do ii=1, m_act
start(activeset(ii)+1)=beta_1(ii)
enddo
if(j .NE. nact)then
if(trace .EQ. 1)then
C call dblepr("b0_1=", -1, b0_1, 1)
call dblepr("beta_1=", -1, beta_1, m_act)
call intpr("begin fullset nclreg_onelambda", -1, 1, 1)
call dblepr("start=", -1, start, m+1)
endif
call nclreg_onelambda(x, y,weights, n,m,start, etastart,
+ mustart, yhat, offset, lambda_i, alpha, gam,
+ intercept,penaltyfactor, maxit, eps,penalty,trace,1,
+ del, rfamily, B, s, thresh, betaall, b0all, fk)
call find_activeset(m, betaall, eps, activeset, jk)
if(jk==0)then
jk=1
activeset(1)=1
fakejk=1
endif
endif
if(jk .NE. m_act)then
conv=0
else
jc=0
do ii=1, max(jk, m_act)
if(activeset(ii)==activeset_old(ii))then
jc=jc+1
endif
enddo
if(jk==jc)then
conv=1
endif
endif
if(jk .NE. m_act .OR. jk .NE. jc)then
deallocate(beta_1, start_act, penaltyfactor_act, x_act)
allocate(beta_1(jk), stat=AllocateStatus)
allocate(start_act(jk+1), stat=AllocateStatus)
allocate(penaltyfactor_act(jk), stat=AllocateStatus)
allocate(x_act(n, jk), stat=AllocateStatus)
start_act(1)=b0all
do 11135 ii=1, jk
beta_1(ii)=0
start_act(ii+1)=betaall(activeset(ii))
penaltyfactor_act(ii)=penaltyfactor(activeset(ii))
activeset_old(ii)=activeset(ii)
11135 continue
do 11155 jj=1, n
do 11145 ii=1, jk
x_act(jj, ii)=x(jj, activeset(ii))
11145 continue
11155 continue
m_act = jk
else
if(fakejk==1)then
start_act(1)=b0all
start_act(2)=betaall(1)
endif
endif
j=j+1
goto 13000
endif
nlambdacal=nlambdacal+1
b0(i)=b0_1
if(jk .GT. 0)then
do ii=1, jk
beta(activeset(ii), i)=beta_1(ii)
enddo
endif
call loss(n, y, fk, cost, rfamily, s, los(i))
call penGLM(beta_1, m_act, lambda_i*penaltyfactor_act,
+ alpha, gam, penalty, penval)
if(standardize .EQ. 1)then
pll(i)=los(i) + n*penval
else
pll(i)=los(i) + penval
endif
C redo (i-1)-lambda estimates with the current start_act for the
C i-th lambda. Note, i=i-2 not i-1 will do this since i=i+1 is
C computed before the end of the loop.
if(decreasing==1)then
if(i > 1)then
if(abs(pll(i)-pll(i-1))/pll(i) > epscycle)then
if(cutpoint==1)then
cutpoint = i
endif
i = i - 2
uturn=1
else
uturn=0
endif
endif
endif
if(cutpoint > 1 .AND. uturn==0)then
i = cutpoint
cutpoint = 1
endif
i = i + 1
goto 10
endif
deallocate(beta_1, start_act, x_act, penaltyfactor_act)
return
end