# CALCULO INFORMAÇÃO MUTUA ENTRE DUAS SERIES TEMPORAIS - A e B
# COM AJUSTE LINEAR DA SÉRIE B

# OBS: Os dados devem estar dispostos em duas colunas com headers A e B


setwd("C:/Users/Sérgio/Documents/Nossos Artigos/Ambiental/Artigo Rotação Espaço")



data <- read.csv('Rn rot 45 45 45 Rebio Jun2008.txt',header = TRUE, sep="\t")

ndat <- nrow(data) - 1

print(ndat)

print(data)

#VALOR MÉDIO DAS SÉRIES A e B:

somaA <- 0.0
somaB <- 0.0

for(i in 1:ndat){

#print(data$AB[i]
#print(data$B[i]

	somaA <- somaA + data$A[i]
	somaB <- somaB + data$B[i]
	
	}

MediaA <- somaA/ndat
MediaB <- somaB/ndat

print(MediaA)
print(MediaB)

if(MediaA*MediaB < 0){

	for(i in 1:ndat){
		data$B[i] <- -data$B[i]
		}
	MediaB <- -MediaB
	}


#DESVIO PADRÃO DAS SÉRIES A e B

somaA <- 0.0
somaB <- 0.0

for(i in 1:ndat){

	somaA <- somaA + (data$A[i] - MediaA)^2
	somaB <- somaB + (data$B[i] - MediaB)^2
	}

DPA <- sqrt(somaA/(ndat-1))
DPB <- sqrt(somaB/(ndat-1))


# REESCALONAMENTO DE B:

CoefLinear <- (MediaA - MediaB)
CoefAng <- DPA/DPB

print("CoefLinear =")
print(CoefLinear)
print("CoefAng =")
print(CoefAng)


for(i in 1:ndat){

	data$B[i] <- data$B[i] + CoefLinear
	}

for(i in 1:ndat){

	data$B[i] <- MediaB + (data$B[i] - MediaB) * CoefAng
	}


# NUMERO DE INTERVALOS e TAO MAXIMO:

ninterv <- 20
taoM <- 20

#as.integer(ndat/1000)

print(nrow(data))


#print(data$A)

Amin <- data$A[1]
Bmin <- data$B[1]
Amax <- data$A[1]
Bmax <- data$A[1]

print(data$A[ndat])
print(data$B[ndat])


# Valores máximos e mínimos de A e B:

for (j in 1:ndat) {
print(j)
	if(data$A[j] < Amin){ Amin <- data$A[j]}
	if(data$B[j] < Bmin){ Bmin <- data$B[j]}
	if(data$A[j] > Amax){ Amax <- data$A[j]}
	if(data$B[j] > Bmax){ Bmax <- data$B[j]}
	}


# Largura dos intervalos:

deltaA <- (Amax - Amin)/ninterv
deltaB <- (Bmax - Bmin)/ninterv


MinfAB <- rep(0.0,times = 2*taoM+1)


for (tao in -taoM:taoM){


ProbA <- rep(0.0,times=ninterv)
ProbB <- rep(0.0,times=ninterv)

ProbAB <- matrix(0.0,nrow=ninterv,ncol=ninterv)

#print(ProbA)

# CALCULO DAS PROBABILIDADES

if(tao >= 0){

for(j in 1:(ndat-tao-1)){

	for(k in 1:ninterv){
		if(data$A[j] >= ((k-1)*deltaA + Amin) && data$A[j] < (k*deltaA + Amin)) ProbA[k] = ProbA[k] + 1
		if(data$B[j+tao] >= ((k-1)*deltaB + Bmin) && data$B[j+tao] < (k*deltaB + Bmin)) ProbB[k] = ProbB[k] + 1
	}
	for(kA in 1:ninterv){
	for(kB in 1:ninterv){
		if(data$A[j] >= ((kA-1)*deltaA + Amin) && data$A[j] < (kA*deltaA + Amin)){
		if(data$B[j+tao] >= ((kB-1)*deltaB + Bmin) && data$B[j+tao] < (kB*deltaB + Bmin)) ProbAB[kA,kB] = ProbAB[kA,kB] + 1
	} # if
	} # for kB
	} # for kA

}	
print(tao)

} #if tao >0

if(tao < 0){

for(j in 1:(ndat+tao-1)){

	for(k in 1:ninterv){
		if(data$A[j-tao] >= ((k-1)*deltaA + Amin) && data$A[j-tao] < (k*deltaA + Amin)) ProbA[k] = ProbA[k] + 1
		if(data$B[j] >= ((k-1)*deltaB + Bmin) && data$B[j] < (k*deltaB + Bmin)) ProbB[k] = ProbB[k] + 1
	}
	for(kA in 1:ninterv){
	for(kB in 1:ninterv){
		if(data$A[j-tao] >= ((kA-1)*deltaA + Amin) && data$A[j-tao] < (kA*deltaA + Amin)){
		if(data$B[j] >= ((kB-1)*deltaB + Bmin) && data$B[j] < (kB*deltaB + Bmin)) ProbAB[kA,kB] = ProbAB[kA,kB] + 1
	} # if
	} # for kB
	} # for kA

}	
print(tao)

} #if tao <0

# NORMALIZAÇÃO

	for(k in 1:ninterv){

		ProbA[k] <- ProbA[k]/(ndat-tao)
		ProbB[k] <- ProbB[k]/(ndat-tao)

	}

	for(kA in 1:ninterv){
	for(kB in 1:ninterv){

	ProbAB[kA,kB] <- ProbAB[kA,kB]/(ndat)
}
}

# CALCULO DA INFORMAÇÃO MÚTUA

soma <- 0.0

	for(kA in 1:ninterv){
	for(kB in 1:ninterv){
	
	if (ProbAB[kA,kB] != 0.0) soma <- soma + ProbAB[kA,kB] * log(ProbAB[kA,kB]/(ProbA[kA]*ProbB[kB]))
}
}
indice <- tao + taoM + 1

MinfAB[indice] <- soma


} #final do for tao

#print(ProbA)
#print(ProbB)
#print(ProbAB)

# Maior Minf:

MaiorMinf <- MinfAB[1]

for(j in 1:(2*taoM+1)){

	if(MinfAB[j] > MaiorMinf){
		MaiorMinf <- MinfAB[j]
		Maiorj <- j
	}
}

MaiorTao <- Maiorj - taoM - 1

print(MaiorMinf)
print(MaiorTao)


# RENORMALIZAÇÃO

if(MaiorMinf > 1.0){
for (j in 1:(2*taoM+1)){

	MinfAB[j] <- MinfAB[j]/MaiorMinf
}
}


# Grafico do resultado:

xs <- seq(from= -taoM,to= taoM,by=1)

plot(xs,MinfAB,xlab="tao",ylab="Mutual Information", type="b")

write.table(MinfAB,file = "Inf Mutua", sep=" ", eol="\n")

write.table(data,file = "Dados Ajustados", sep=" ", eol="\n")