#
#	tri_class.R - A set of R functions for Classical
#		      Item Analysis
#
#Copyright (C) 2005 Fernando Henrique F. P. Rosa <feferraz at ime.usp.br>
#		    Robson Lunardi		 <rlunardi at linux.ime.usp.br>
#	
# 
#This program is free software; you can redistribute it and/or
#modify it under the terms of version 2 of the GNU General Public License
#as published by the Free Software Foundation.  A copy of this license should
#be included in the file COPYING.
#
#This program is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#GNU General Public License for more details.
#
#You should have received a copy of the GNU General Public License
#along with this program; if not, write to the Free Software
#Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
#
#

## Esse conjunto de funções em R implementa análise de 
## teoria clássica para itens -> reproduzindo o comportamento
## do IteMan.
## Para mais detalhes veja:
## http://www.linux.ime.usp.br/~feferraz/files/lista/mae5778-l1.pdf



le.dados <- function(nome.arquivo) {
	arq <- file(nome.arquivo)
	pmlinha <- scan(arq,n=2)
	n.questoes <- pmlinha[1]
	tam.id <- pmlinha[2]
	dados.brutos <- readLines(arq)
	gabarito <- strsplit(dados.brutos[2],split='')[[1]]
	if (length(gabarito) != n.questoes) { 
		stop("Pau! Arquivo inconsisntente.")
	}
        usa.questao <- strsplit(dados.brutos[4],split='')[[1]]

	parse.linha <- function(linha) {
		q.linha <- strsplit(linha,split='')[[1]]
		id <- paste(q.linha[1:tam.id],collapse='')
		respostas <- q.linha[(tam.id+1):(n.questoes+tam.id)]
		res <- c(id,respostas)
		ifelse(res == ' ' | is.na(res),'NULO',res)
	}

	respostas <- (t(sapply(dados.brutos[5:length(dados.brutos)],parse.linha,USE.NAMES=FALSE)))

	usa.questao <- ifelse(usa.questao == "Y",TRUE,FALSE)
	
	obj.res <- list(respostas=respostas,gabarito=gabarito,usa.questao=usa.questao,n.questoes=n.questoes)
	class(obj.res) <- 'tri'
	close(arq)	
	obj.res
}

analise <- function(dados,usa.questao=dados$usa.questao,gabarito=dados$gabarito) {
	score <- apply(dados$respostas[,c(FALSE,usa.questao)],1,function(ind) { sum(ind == dados$gabarito[usa.questao]) })
	N <- length(dados$respostas[,2])	

	# 27% of the highest
        ind.high <- N-ceiling(N*0.27)+1
        highest <- score >= sort(score)[ind.high]
	# 27% of the lowest
	ind.low <- ceiling(N*0.27)
	lowest <- score <= sort(score)[ind.low]


	analise.questao <- function(coluna) {
  		ind.dif <- sapply(sort(unique(coluna)),function(resp) { sum(resp == coluna) })/N
		end.low <- sapply(sort(unique(coluna)),function(resp) { sum(resp == coluna[lowest]) })/sum(lowest)
		end.high <- sapply(sort(unique(coluna)),function(resp) { sum(resp == coluna[highest]) })/sum(highest)
		pt.biss <- sapply(sort(unique(coluna)),function(resp) { cor(score,coluna == resp) })
		biss <- sapply(names(pt.biss),function(resp) { pt.biss[resp] * sqrt(ind.dif[resp]*(1-ind.dif[resp])) / dnorm(qnorm(ind.dif[resp])) },USE.NAMES=FALSE )
		res <- t(rbind(ind.dif,end.low,end.high,pt.biss,biss))
		res
	}

	an.questoes <- lapply((2:(dados$n.questoes+1))[usa.questao],function(questao) analise.questao(dados$respostas[,questao]))
	names(an.questoes) <- (1:dados$n.questoes)[usa.questao]
	an.questoes <-lapply((1:dados$n.questoes)[usa.questao],function(q) { res <- an.questoes[[as.character(q)]]; attr(res,'gabarito') <- gabarito[q]; attr(res,'questao') <- q ;res  })
        names(an.questoes) <- (1:dados$n.questoes)[usa.questao]


	#calculo do alpha

	prop.acerto <- sapply(1:sum(usa.questao),function(i) an.questoes[[i]][gabarito[usa.questao][i],1])
	names(prop.acerto) <- names(an.questoes)

	K <- sum(usa.questao)
	alpha <- K/(K-1) * (1 - sum(prop.acerto * (1-prop.acerto))/var(score))
	
	SEM = sd(score)*sqrt(1-alpha)

	an.global <- list(questoes.efetivas=K,N=N,media=mean(score),variancia=var(score),desvio=sd(score),mediana=median(score),alpha=alpha,SEM=SEM,max.low=sort(score)[ind.low],N.low=sum(lowest),min.high=sort(score)[ind.high],N.high=sum(highest))
	
	names(score) <- dados$respostas[,1]
	obj.res <- list(an.questoes=an.questoes,an.global=an.global,dados=dados,outros=list(score=score,prop.acerto=prop.acerto,usa.questao=usa.questao))
	class(obj.res) <- 'trianal'
	obj.res


}

print.trianal <- function(obj,mostrar.questoes=names(obj$an.questoes),global=TRUE) {
	if (length(mostrar.questoes) > 0) {
	mostrar.questoes <- as.character(mostrar.questoes)
	if (sum(is.na(match(mostrar.questoes,names(obj$an.questoes)))) > 0) {
		stop("Trying to show questions not considered for analysis.")
	}
	imprime.quest <- function(quest) {
		gabarito <- attr(quest,'gabarito')
		cat(paste('\n\nQuestion',attr(quest,'questao')),'\n\n') 
		disc.index <- quest[gabarito,3] - quest[gabarito,2]
		print(data.frame(key=gabarito,prop.correct=quest[gabarito,1],disc.index=disc.index,row.names=' ',pt.biss=quest[gabarito,4],biss=quest[gabarito,5]))
		key <- rep(' ',length(quest[,1]))
		key[which(gabarito == row.names(quest))] <- "*"
		cat("\n")
		print(data.frame(ind.dif=quest[,1],end.low=quest[,2],end.high=quest[,3],pt.biss=quest[,4],biss=quest[,5],key))
		cat("\n")

	}

	sapply(mostrar.questoes,function(quest) imprime.quest(obj$an.questoes[[quest]]) )
	}
	if (global) {
		global <- t(data.frame(obj$an.global))
		colnames(global) <- "Global Statistics"
		cat('\n')
		print(global)
		cat('\n')
	}

	return(invisible(obj))

}

triplot <- function(quest,obj,classes=7) {
	K <- obj$an.global$questoes.efetivas 
	tam.class <- (K+1)/classes
	tam.efet <- round(tam.class)
	if (tam.efet == floor(tam.class)) {
		niveis <- sort(rep(1:classes,tam.efet,each=T))
		if (length(niveis) != K+1) {
			falt <- K+1 - length(niveis)
			niveis[(length(niveis)+1):(K+1)] <- niveis[length(niveis)]
		} 
	}
	else {
                niveis <- sort(rep(1:classes,tam.efet,each=T))
                if (length(niveis) != K+1) {
                        excess <- length(niveis) - (K+1)
			niveis <- niveis[1:(K+1)]
                } 
	}

	processa.med <- function(i,quest,obj) {
		nivel.scr <- which(niveis == i)-1
		ind <- obj$outros$score %in% nivel.scr
		coluna <- obj$dados$respostas[,quest+1]	
		gabarito <- obj$dados$gabarito[quest]
		prop.acerto <- sum(coluna[ind]== gabarito)/sum(ind)
		media.score <- mean(obj$outros$score[ind])
		c(media.score,prop.acerto)

	}
	res <- list(res=t(sapply(1:niveis[K+1],processa.med,quest,obj)),niveis=niveis,questao=quest)
	attr(res,'class') <- 'triplot'
	res
}

print.triplot <- function(obj,...,ylim=c(0,1)) {
	plot(obj$res,main=paste('Diagrama de Dispersão do Escore Médio e Proporção de Acerto\nQuestão',obj$questao),xlab=paste('Escore médio\nClasses: ',obj$niveis[length(obj$niveis)]),ylab='Proporção de Acerto',type='l',col='blue',ylim=ylim,...)
}