OGA <- function(x,y,beta_init){ 
	n <- length(y)
	x <- as.matrix(x)	
	p <- dim(x)[2]
		
	e=y-cbind(rep(1,n),x)%*%beta_init
		
	y_new=y
	correction=0	
## corrections from (E[e])^2 to E[e^2] in computing BIC
	
	m=5*floor(sqrt(n/log(p)))
	oga=OGA2(x,y_new,m)		
	for (k in 1:m){
		U=oga$residuals[,k]
		oga$BIC[k]=log(mean(U^2)+correction/n)+k*log(n)/n
		oga$HDBIC[k]=n*log(mean(U^2)+correction/n)+k*log(n)*log(p)
		oga$HDAIC[k]=n*log(mean(U^2)+correction/n)+k*2.01*log(p)
		oga$HDHQ[k]=n*log(mean(U^2)+correction/n)+2.01*k*log(log(n))*log(p)
	}
	return (oga)
}



OGA2 <- function(x,y,m) {
x <- as.matrix(x)
n <- length(y)	
p <- dim(x)[2]

# centering
mx=colMeans(x)
x=x-matrix(rep(mx,each=n),n,p)
my=mean(y)
y=y-mean(y)

beta_all<-matrix(rep(0,p*m),p,m)
beta_current<-rep(0,p)
alpha_all<-rep(0,m)	# intercept
U_all<-matrix(0,n,m)	# residuals
XO<-matrix(0,n,m)
U<-y
j<-rep(0,m)
a <- b <- matrix(0,m,m)
beta_hat <- RegularBIC <- HiDimBIC <- rep(0,m)
HiDimAIC <- HiDimHQ <- rep(0,m)


for (k in 1:m) {
	corrs<-rep(0,p);
	if (k==1)
		corrs<-((t(U)%*%x)^2)/colSums(x^2)
	else
		corrs[-j]<-((t(U)%*%x[,-j])^2)/colSums(matrix(x[,-j],n,p-k+1)^2) # make a matrix to ensure colSums work for k=2
	j[k]<-which.max(corrs);
	X_o<-x[,j[k]];
	if (k>1) {
		a[k,1:(k-1)]<-(t(X_o)%*%XO[,1:(k-1)])/colSums(matrix(XO[,1:(k-1)],n,(k-1))^2); # needs revision
		X_o<-X_o-XO[,1:(k-1)]%*%matrix(a[k,1:(k-1)],(k-1),1);  # needs revision
	}
	if (sum(X_o^2)>0.00001) { # threshold to be determined
		beta_hat[k]<-(t(U)%*%X_o)/sum(X_o^2);
		XO[,k]<-X_o;
		U=U-beta_hat[k]*X_o;
		U_all[,k]=U;
	}
	RegularBIC[k]=log(mean(U^2))+k*log(n)/n;
	HiDimBIC[k]=n*log(mean(U^2))+k*log(n)*log(p);
	HiDimAIC[k]=n*log(mean(U^2))+k*2.01*log(p);
	HiDimHQ[k]=n*log(mean(U^2))+2.01*k*log(log(n))*log(p);
}
for (k in 2:m)
	b[k,1:(k-1)]=-(a[k,1:(k-1)]+matrix(a[k,1:(k-1)],1,(k-1))%*%matrix(b[1:(k-1),1:(k-1)],(k-1),(k-1))) # needs revision
beta_current[j[1]]=beta_hat[1];
beta_all[j[1],1]=beta_hat[1];
alpha_all[1]=my-beta_current[j[1]]*mx[j[1]]
for (k in 2:m){
	beta_current[j[1:(k-1)]]=beta_current[j[1:(k-1)]]+beta_hat[k]*b[k,1:(k-1)];
	beta_current[j[k]]=beta_hat[k];
	beta_all[,k]=beta_current;
	alpha_all[k]=my-sum(mx[j[1:k]]*beta_current[j[1:k]])
}
beta_all=rbind(alpha_all,beta_all)
return (list(coeff_final=beta_current,coeff_all=beta_all,residuals=U_all,
		BIC=RegularBIC,HDBIC=HiDimBIC,HDAIC=HiDimAIC,HDHQ=HiDimHQ));
}