"maimed"<-
function(data, effmat=NULL, formula = c("name", "type", "size"), niter = 31,
weighted=F, calculateEffort=F, minEff=1/60, jk=F)
{
### wrapper for calling jkfit.  
### data: argument is a data frame with the following records:  
###   name:	name of developer
###     optionally, other variables.  
### effmat: the matrix where rows corrspond to distinct names and
###   columns to months. If matrix is not provided or calcEff=T the
###   matrix is calculated from the data. 
### formula: list of variable names (in dataframe data) to use in
###   the model.
### niter: number of iterations in jkfit
### weighted: if FALSE - equally distribute weights among MRs in a
###   single month done by a single developer.
###   if TRUE alocate proportionally to the time open in that
###   particular month.
### calcEff: do not use unit effort, but calculate it based on the
###   total time all MRs are open during a month.
### minEff: some MRs may be open zero time, minEff specifies minimal
###   time an MR is open.
### jk: do jackknife (i.e., remove individual developers to estimate
###   the variances.

	m <- dim(data)[1];
	naymes <- names(table(data$name));
	o <- length(naymes);
	months <- floor(max(data$close))+1;
	if (!is.matrix(effmat)) { 
	   Effmat <- matrix(1, o, months);
	}else{ 
	   Effmat <- effmat; 
	}
	if (dim(Effmat)[1] != o || dim(Effmat)[2] != months){
		print ("Wrong dimensions of effort matrix");
		return(0);
	}
	MR <- data[,formula];
	if (length(formula)==1){
		print ("Please use at least one factor in the model");
		return(1);
	}


	if (calculateEffort){
		initial1 <- maimed.initial3(data, Effmat, naymes,minEff);
	}else{
		if(weighted){
			initial1 <- maimed.initial2(data, Effmat, naymes, minEff);
		}else{ initial1 <- maimed.initial(data,Effmat,naymes); }
	}	

        tmp1MR <- data.frame(MR[,-match("name",names(MR))]);
	names(tmp1MR) <- names(MR)[-match("name",names(MR))];
	for (i in naymes) tmp1MR[,i] <- as.integer(MR$name==i);
	if (!jk){
		return(maimed.fit(initial1=initial1$initial, 
			Effmat=initial1$Effmat, 
			MR=tmp1MR,data=data,naymes=naymes, 
			niter=niter));
	}else{
		res0 <- maimed.fit(initial1=initial1$initial, 
			Effmat=initial1$Effmat,MR=tmp1MR,data=data,
			naymes=naymes, niter=niter);
		res <- res0;
		for (i in naymes){
			ind <- !(MR$name == i);
			tmp <- maimed.fit(initial1=initial1$initial[ind,],
				Effmat=initial1$Effmat[-match(i, naymes),],
				MR=tmp1MR[ind,],data=data[ind,],
				naymes=naymes[-match(i, naymes)], 
				niter=niter);
			res <-rbind(res,tmp[match(names(res0),names(tmp))]);
		}
		dimnames(res)[[2]] <- names(res0);
		return (res);
	}
}

"maimed.initial"<-
function(data, Effmat, naymes)
{
	newMR <- data[, c("open", "close")];
	m <- dim(data)[1];
	n <- dim(Effmat)[2];
	MAT <- matrix(0, m, n);
	for(i in 1:n) {
            ind <- floor(newMR$open) <=i-1 & floor(newMR$close) >= i-1;
            for (j in 1:length(naymes)){
		ind1 <- data$name[ind]==naymes[j];
		MAT[ind, i][ind1] <- Effmat[j,i]/sum(ind1);
	    }
        }
	return(list(initial=MAT, Effmat=Effmat));
}

"maimed.initial2"<-
function(data,Effmat,naymes, minEff=1/60)
{
	newMR <- data[, c("open", "close")];
	m <- dim(data)[1];
	n <- dim(Effmat)[2];
	MAT <- matrix(0, m, n);
	for(i in 1:n) {
            ind <- floor(newMR$open) <=i-1 & floor(newMR$close) >= i-1;
	    left <- rep(i-1, sum(ind));
	    right <- rep(i, sum(ind));
	    left [left < newMR$open[ind]] <- newMR$open[ind][left < newMR$open[ind]];
	    right [right > newMR$close[ind]] <- newMR$close[ind][right > newMR$close[ind]];
	    ll <- right-left;
	    ll[ll < minEff] <- minEff;
            for (j in 1:length(naymes)){
		ind1 <- data$name[ind]==naymes[j];
		MAT[ind, i][ind1]  <-  ll[ind1]*Effmat[j,i]/sum(ll[ind1]);
	    }
        }
	return(list(initial=MAT, Effmat=Effmat));
}

"maimed.initial3"<-
function(data, Effmat, naymes, minEff=1/60)
{
	newMR <- data[, c("open", "close")];
	m <- dim(data)[1];
	n <- dim(Effmat)[2];
	MAT <- matrix(0, m, n);
	tmpEffmat <- Effmat;
	for(i in 1:n) {
            ind<-floor(newMR$open) <=i-1 & floor(newMR$close) >= i-1;
	    left <- rep(i-1, sum(ind));
	    right <- rep(i, sum(ind));
	    left [left < newMR$open[ind]] <- newMR$open[ind][left < newMR$open[ind]];
	    right [right > newMR$close[ind]] <- newMR$close[ind][right > newMR$close[ind]];
	    ll <- right-left;
	    ll[ll < minEff] <- minEff;
            for (j in 1:length(naymes)){
		ind1 <- data$name[ind]==naymes[j];
		MAT[ind, i][ind1]  <-  ll[ind1];
		tmpEffmat[j,i] <- sum(ll[ind1]);
	    }
	}
        return(list(initial=MAT, Effmat=tmpEffmat));
}

maimed.fit<-
function(initial1, Effmat, MR, data, naymes, niter=10){
	m <- dim(MR)[1];
	n <- dim(Effmat)[2];

        tmpMR <- MR;
	for (i in 1:niter){
		tmpMR$mrEff <- initial1%*%rep(1, n);
		assign("tmpMR", tmpMR, inherits=T);
		mod <- glm(mrEff~.-1,family=poisson, data=tmpMR);
		assign("mod", mod, inherits=T);
		summm<-summary(mod);
		coef<-summm$coefficients;
		devi<-summm$deviance;
		print (paste("Iteration:",i," deviance:",devi,sep=""));print (coef[,1]);
		mrEff1<-fitted(mod);
		initial1<-initial1*as.vector(mrEff1/tmpMR$mrEff);
		for (j in 1:length(naymes)){
       	    		ind <- data$name==naymes[j];
			if (sum (ind) == 1) {
 				print (paste ("developer ", naymes[j], "has only one MR, can not proceed."));
 				exit();
			}

	    		monthEff<-rep(1, sum(ind))%*%initial1[ind,];
	    		res<-as.vector(Effmat[j,]/monthEff);
	    		res[monthEff <= 0] <- 1;
	    		initial1[ind,]<-t(res*t(initial1[ind,]))
		}
	}
	coef[,1];
}
maimed.fit1<-
function(initial1, Effmat, MR, niter=10){
	m <- dim(MR)[1];
	n <- length(Effmat);

        tmpMR <- MR;
	for (i in 1:niter){
		tmpMR$mrEff <- initial1%*%rep(1, n);
		assign("tmpMR", tmpMR, inherits = T);
		mod <- glm(mrEff~.-1,family=poisson, data=tmpMR);
		assign("mod", mod, inherits=T);
		summm<-summary(mod);
		coef<-summm$coefficients;
		devi<-summm$deviance;
		print (paste("Iteration:",i," deviance:",devi,sep=""));print (coef[,1]);
		mrEff1<-fitted(mod);
		initial1<-initial1*as.vector(mrEff1/tmpMR$mrEff);
	    	monthEff<-apply(initial1, 2, sum);
	    	res<-as.vector(Effmat/monthEff);
	    	res[monthEff <= 0] <- 1;
	    	initial1<-t(res*t(initial1))
	}
	coef[,1];
}
"maimed.sd"<-function(x){
 x <- x[!is.na(x)];
 n <- length(x);
 (n-1)*sqrt(var(x)/n);
}
"maimed.validate"<-function(model)
{
### jackknife estimates of standard errors 
### mode: model estimated using maimed with jk=T
### outputs estimates, standard deviations, t-statistic, 
### two-sided p-values, and 95% CI
  mm <- model[1,];
  sd <- apply(model, 2, maimed.sd);
  tt <- mm/sd;
  pp <- tt;pp[pp>0]<- -pp[pp>0];
  pp <- pt (pp, sum(!is.na(mm)))*2;
  ci <- cbind(mm + sd*qt(.025, sum(!is.na(mm))), 
        mm + sd*qt(.975, sum(!is.na(mm))));
  res <- cbind(mm, sd, tt, pp,ci);
  dimnames(res) <- list(names(mm), c("estimate","stdev","t-statitic", "p-value", "95%CI", "95%CI"));
  res;
}