# --- Jan. 29th, 2007 ---
# --- Meta TGDR, estimation ---

#     March 19th, 2007
#     simulated data for plot

print(Sys.time());

tau<-1.0;
step<-1000;
increase<-0.01;

d1.y<-read.table("out-1.dat");
d2.y<-read.table("out-2.dat");
d3.y<-read.table("out-3.dat");

all.y<-rbind(d1.y, d2.y, d3.y);
all.y<-all.y$V1;
samplesize<-length(all.y);


d1.x<-read.table("exp-1.dat");
d2.x<-read.table("exp-2.dat");
d3.x<-read.table("exp-3.dat");

all.x<-cbind(d1.x, d2.x, d3.x);


size<-c(ncol(d1.x), ncol(d2.x), ncol(d3.x));
group<-rep(3, sum(size));
for(i in 1:size[1]) group[i]<-1;
for(i in (size[1]+1):(size[1]+size[2])) group[i]<-2;


# --- remember to perturbe group also ---
random<-read.table("random.dat")$V1;
random.w<-random[1:length(all.y)];
rank.w<-(1:length(all.y))[order(random.w)];

all.y<-all.y[rank.w];
all.x<-all.x[ ,rank.w];
group.status<-group[rank.w];
###


group.name<-sort(unique(group.status));
num.group<-length(unique(group.status));
group.size<-rep(0, num.group);
for(i in 1:num.group){
   work.group<-group.name[i];
   group.size[i]<-length(group.status[group.status==work.group]);
}

aug.x<-t(cbind(rep(1, samplesize), t(all.x)));
num.cov<-nrow(all.x);



# --- estimation ---
gradient.mat<-matrix(0, (num.cov+1), num.group); # -- intercept
beta.mat<-matrix(0, (num.cov+1), num.group);

aug.z.dat<-aug.x;
indicator<-all.y;


for(s in 1:step){

  # --- for data from each group ---
  for(g in 1:num.group){
     work.z.dat<-aug.z.dat[ ,group.status==group.name[g]];
     work.mat<-matrix(0, nrow(work.z.dat), ncol(work.z.dat));
     for(mm in 1:ncol(work.z.dat))
        work.mat[ ,mm]<-work.z.dat[ ,mm];
     work.beta<-beta.mat[ ,g];
     work.y<-indicator[group.status==group.name[g]];

     beta.z<-t(work.mat)%*%work.beta;
     beta.p<-1/(1+exp(-beta.z));
     term.1<-work.y/beta.p-(1-work.y)/(1-beta.p);
     term.2<-exp(beta.z)/((1+exp(beta.z))^2);
     work.gradient<-(t(term.1*term.2))%*%t(work.mat);

     gradient.mat[ ,g]<-work.gradient[1, ];  ###/group.size[g];
  }

  sum.gradient<-abs(apply(gradient.mat, 1, sum));
  max.gradient<-max(sum.gradient[2:length(sum.gradient)]);

  gradient<-gradient.mat/max.gradient;

  f.v<-ifelse(sum.gradient>=(tau*max(sum.gradient)), 1, 0);
  f.v[1]<-1;

  g.v<-gradient*f.v;

  beta.mat<-beta.mat+increase*g.v;

  beta.out<-cbind(t(beta.mat[2:5,1]), t(beta.mat[2:5 ,2]),
		t(beta.mat[2:5,3]));
  write.table(round(beta.out, digits=3), 
	file="beta-path.txt", quote=F, append=T, row.name=F, col.name=F); 
}

beta.sum<-apply(abs(beta.mat), 1, sum);
size<-length(beta.sum[beta.sum!=0]);

score<-rep(0, samplesize);
for(i in 1:samplesize){
   x<-aug.z.dat[ ,i];
   temp.group<-group.status[i];
   score[i]<-sum(x*beta.mat[ ,temp.group]);
}

cut<-rep(0, num.group);
for(i in 1:num.group){
   score.group<-score[group.status==group.name[i]];
   group.ind<-indicator[group.status==group.name[i]];
   cut[i]<-sort(score.group, decreasing=T)[sum(group.ind)];
}

pred<-rep(0, samplesize);
for(i in 1:samplesize){
   pred[i]<-ifelse(score[i]>=cut[group.status[i]], 1, 0);
}

error<-sum(abs(pred-all.y))/length(all.y);

result<-cbind(tau, s, increase, size, error);

print(Sys.time());