#***************************************
# Creado por Mary Pili Vargas 27-09-05
#***************************************

# Este algoritmo realiza la estimación de valores faltantes en una estación particular
# considerando para la estimación cualquiera de los siguientes métodos  
#(REGRESIÓN LINEAL, PROMEDIO SIMPLE Y NORMAL DE  PRECIPITACIÓN), la selección del método 
# se deja a criterio del usuario. 

# La base de datos relacionada con cada estación debe existir en el directorio de trabajo como
# un objeto de R.

# Para el caso particular de nuestro proyecto tenemos

# distancia.txt  que contiene información general de las 68 estaciones (nombre de
# de la estación, serial, estaciones adyacentes, distancias entre las estaciones). Por medio
# de este archivo  formamos los nombres de las estaciones y posteriormente cargamos las bases de
# datos en el directorio de trabajo. 

# Para convertir distancia.txt  en un objeto de R leemos con la función read.table() el archivo
# y se carga en un objeto, en este caso datos

datos <- read.table("distancia.txt",header = T, sep = ",", dec = ".", as.is = F, na.strings = "NA",colClasses = c("integer","character","integer","character","real","integer","character","real","integer","character","real","integer","character","real","integer","integer","character"))

# El procedimiento de como cargar las bases de datos en el directorio de trabajo se describe en 
# el algoritmo read_distancia.

# Asumiendo que las bases de datos ya  están cargadas como objetos de R tenemos:

# FUNCIÓN DE  ESTIMACIÓN DE VALORES FALTANTES GENERAL

# est.clave es la estación donde se desea  hacer la estimación 
# ady1, ady2, ady3 son estaciones adyacentes a est.clave

val_falt <- numeric(0) # guarda la posición que el valor faltante ocupa dentro de la base de datos 
val_falt <-which(is.na(est.clave$PRE))

ano.fal<-numeric(length(val_falt)) # guarda el año faltante
ano.fal<-est.clave[which(is.na(est.clave$PRE)),1]

mes.fal<-numeric(length(val_falt)) # guarda el mes faltante
mes.fal<-est.clave[which(is.na(est.clave$PRE)),2]

val_fal_estm <- numeric(length(val_falt)) # guarda las estimaciones

# m permite acceder a una posición  en particular de cualquiera de los vectores antes mencionados, 
# y puede ser utilizado con cualquier estructura de repetición.

#********* FUNCIÓN DE ESTIMACIÓN *************
estima_fal <- function(est.clave,ady1,ady2,ady3,m)  {

vec.min <- numeric(3)
vec.max <- numeric(3)

# Calcule el año inicial y final de la estación adyacente y anotelo en el vector correspondiente
vec.min[1]<-min(ady1$A)
vec.max[1]<-max(ady1$A)

vec.min[2]<-min(ady2$A)
vec.max[2]<-max(ady2$A)

vec.min[3]<-min(ady3$A)
vec.max[3]<-max(ady3$A)

minimo<-max(vec.min,min(est.clave$A))
maximo<-min(vec.max,max(est.clave$A))

# Valores de las estaciones de referencia a usar en la Predicción
	p1<-ady1[ady1$A==ano.fal[m]& ady1$M==mes.fal[m],3]
	p2<-ady2[ady2$A==ano.fal[m]& ady2$M==mes.fal[m],3]
	p3<-ady3[ady3$A==ano.fal[m]& ady3$M==mes.fal[m],3]

# Si es REGRESIÓN LINEAL

# base<-est.clave[est.clave$A >= minimo & est.clave$A <=maximo & est.clave$M == mes.fal[m],3]
#	ref1<-ady1[ady1$A >= minimo & ady1$A <= maximo & ady1$M == mes.fal[m],3]
#	ref2<-ady2[ady2$A >= minimo & ady2$A <= maximo & ady2$M == mes.fal[m],3]
#	ref3<-ady3[ady3$A >= minimo & ady3$A <= maximo & ady3$M == mes.fal[m],3]
#
#	reg.df<-data.frame(base,ref1,ref2,ref3)

#	reg.lm<-lm(base ~ ref1 + ref2 + ref3, data=reg.df, na.action=na.exclude)
#	print(summary(reg.lm))

	#Estimación
#	valor_estimado <-predict(reg.lm,list(ref1=p1,ref2=p2,ref3=p3))
#	if (valor_estimado < 0) valor_estimado <- 0   #para evitar valores de precipitación negativos

#	return(valor_estimado)
  
#Si es PROMEDIO SIMPLE

   #     valor_estimado <- (p1 + p2 + p3)/3
#	return(valor_estimado)

#Si es NORMAL DE PRECIPITACIÓN 

	min.a<-minimo
	max.a<-minimo + 30

	    base.pan<-sum(na.exclude(est.clave[(est.clave$A >= min.a & est.clave$A <= max.a)&(est.clave$M==mes.fal[m]),3]))/30
	    ady1.pan<-sum(na.exclude(ady1[(ady1$A >= min.a & ady1$A <= max.a)&(ady1$M==mes.fal[m]),3]))/30
	    ady2.pan<-sum(na.exclude(ady2[(ady2$A >= min.a & ady2$A <= max.a)&(ady2$M==mes.fal[m]),3]))/30
	    ady3.pan<-sum(na.exclude(ady3[(ady3$A >= min.a & ady3$A <= max.a)&(ady3$M==mes.fal[m]),3]))/30

	    #Estimación
	    valor_estimado <- 1/3*((base.pan/ady1.pan)*p1+(base.pan/ady2.pan)*p2+(base.pan/ady3.pan)*p3)
   	    return(valor_estimado)
}

# OBSERVACIÓN
# La función estima_fal hace la estimación con un sólo método. El método dependera del usuario
# para este caso, el método activo es Normal de Precipitación, pero se deja como practica hacer 
# la estimación con Regresión Lineal ó Promedio Simple.