R/data_preparation.R

Defines functions geologger.sampler.create.arrays create.proposal make.processed.light.object make.result.list correct.hours create.calibration make_likelihood_correction_function suggest.irrad.boundaries get.calibration.parameters plot_slopes get.calib.param logger.template.calibration logger.template.calibrarion.internal get.Irradiance process.twilights make.prerun.object find.stationary.location make.calibration plot_slopes_by_location

Documented in find.stationary.location make.calibration make.prerun.object plot_slopes_by_location

# data_preparation functions

#' plots log of observed versus expected slope by time for a known location
#'
#' The function calculates and plots calibration slopes for sunsets and sunrises for every day of the tracking period, based on the assumption that the tag remained in the same (calibration) location all the time.
#'
#' @param Proc.data processed data object generated by \code{\link{get.tags.data}}
#' @param location vector with longitude and latitude of calibration location (degrees).
#' @param log.light.borders numeric vector with length of 2 for minimum and maximum log(light) levels to use. Default value 'auto', will take these values from the Proc.data object.
#' @param log.irrad.borders numeric vector with length of 2 for minimum and maximum log(irradiance) values to use. Default value 'auto', will take these values from the Proc.data object.
#' @param ylim the y limits of the plot. The default value, NULL, indicates that the range of the finite values to be plotted should be used.
#' @param xlim the x limits of the plot. The default value, NULL, otherwise can be POSIXct or character in a form readable by \code{\link{as.POSIXct}}.
#' @details The plot of calibration slopes is used for finding start and end dates of a calibration period (the time period, during which the tag remained in the calibration location with coordinates (x,y)). During the calibration period, the calibration slopes vary little both, between the twilight events (sunrises and sunsets) and in time. When the tag changes location, the slopes for sunrises and sunsets start to deviate. There may potentially be several calibration periods for the same location (if the bird returned to the same location several times). The boundaries (start and end dates) of each of these periods are captured visually. If there were more than one calibration location, the procedure is repeated, once for each location. 
#' All the obtained calibration periods can be entered in a data frame 'Calibration.periods', for further analysis. Each line of the data frame contains start and end dates (if applicable) of the calibration period and geographic coordinates of the location.
#' @return 'NULL'
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' Proc.data<-get.tags.data(File)
#' plot_slopes_by_location(Proc.data=Proc.data, location=c(5.43, 52.93))
#' abline(v=as.POSIXct("2013-08-20", tz='GMT')) # end of first calibration period
#' abline(v=as.POSIXct("2014-05-05", tz='GMT')) # start of the second calibration period
#'
#' @author Eldar Rakhimberdiev
#' @export plot_slopes_by_location
plot_slopes_by_location<-function(Proc.data, location, log.light.borders='auto', log.irrad.borders='auto', ylim=NULL, xlim=NULL) {
   old.par <- graphics::par(no.readonly = TRUE) 
   Calibration.period<-data.frame(
         calibration.start=as.POSIXct("1900-01-01", tz='GMT'),
		 calibration.stop=as.POSIXct("2050-01-01", tz='GMT'),
		 lon=location[1], lat=location[2])
   if (log.light.borders[1]=='auto') log.light.borders<-Proc.data$log.light.borders
   if (log.irrad.borders[1]=='auto') log.irrad.borders<-Proc.data$log.irrad.borders
   calibration.parameters<-suppressWarnings(get.calibration.parameters(Calibration.period,
         Proc.data, model.ageing=FALSE, 
		 log.light.borders=log.light.borders,
		 log.irrad.borders=log.irrad.borders, 
		 plot.each = FALSE, plot.final = FALSE, suggest.irrad.borders=FALSE))
   tmp<-suppressWarnings(plot_slopes(calibration.parameters$All.slopes, ylim=ylim, xlim=xlim))
   #graphics::par(old.par)
   return(tmp)
}


#' Creates a calibration object, further used for calculation of coordinates in the \code{\link{run.particle.filter}}.
#' 
#' Function estimates all necessary parameters from the calibration data logged in a known location or locations.
#' 
#' @param Proc.data processed data object generated by \code{\link{get.tags.data}}
#' @param Calibration.periods a data frame containing start and end dates of all the calibration periods (POSIXct) and geographic coordinates of the corresponding calibration locations.
#' @param model.ageing if set to TRUE, accounts for the tag ageing (with opacification of its transparent shell of a light sensor), resulting into decreasing sensitivity of the device. This option is useful only if there were several calibration periods or if calibration period was very long (~ longer than a month).
#' @param plot.each Do you want every twilight to be plotted while processing
#' @param plot.final Do you want final calibration graph to be plotted. On the graph you can see all the observed versus expected light levels. All slopes should be similar.
#' @param likelihood.correction will estimate correction of likelihood for the current calibration parameters. Highly recommended not to be change from 'auto'. In this case FLightR will switch it to FALSE in case tag saved data on 10 minutes or longer period.
#' @param fixed.logSlope these are mean (1) and SD (2) for distribution of slopes. Should normally be estimated from the data (and thus default is c(NA, NA)). Change any of these two finite values if you want them to be predetermined and not estimated from the calibration data.
#' @param suggest.irrad.borders experimental parameter! If set to TRUE function will try to find the best values for the log.irrad.borders
#' @param return.slopes if true function will return estimated individual twilight slopes.
#' @return calibration object to be uses in the \code{\link{make.prerun.object}}
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' Proc.data<-get.tags.data(File, end.date=as.POSIXct('2013-08-20', tz='GMT'))
#' Calibration.periods<-data.frame(
#'        calibration.start=NA,
#'        calibration.stop=as.POSIXct("2013-08-20", tz='GMT'),
#'        lon=5.43, lat=52.93) 
#'        #use c() also for the geographic coordinates, if you have more than one calibration location
#'        # (e. g.,  lon=c(5.43, 6.00), lat=c(52.93,52.94))
#' print(Calibration.periods)
#' 
#' # NB Below likelihood.correction is set to FALSE for fast run! 
#' # Leave it as default 'auto' for real examples
#' Calibration<-make.calibration(Proc.data, Calibration.periods, likelihood.correction=FALSE)
#' 
#' @author Eldar Rakhimberdiev
#' @export
make.calibration<-function(Proc.data, Calibration.periods, model.ageing=FALSE, plot.each=FALSE, plot.final=FALSE, likelihood.correction='auto', fixed.logSlope=c(NA,NA), suggest.irrad.borders=FALSE, return.slopes=FALSE) {
   Calibration.periods[,1]<-as.POSIXct(Calibration.periods[,1], tz='GMT')
   Calibration.periods[,2]<-as.POSIXct(Calibration.periods[,2], tz='GMT')
   Calibration.periods$calibration.start[is.na(Calibration.periods$calibration.start)]<-as.POSIXct("1900-01-01", tz='GMT')
   Calibration.periods$calibration.stop[is.na(Calibration.periods$calibration.stop)]<-as.POSIXct("2100-01-01", tz='GMT')
   for (i in 1:nrow(Calibration.periods)) {
     if (Calibration.periods[i,1]>=Calibration.periods[i,2]) stop('Calibration start in some period is later than calibration end')
	 if (i>1) if ( Calibration.periods[i,1]<=Calibration.periods[(i-1),2]) stop('Calibration periods overlap')
   }
   calibration.parameters<-suppressWarnings(get.calibration.parameters(
         Calibration.periods, Proc.data,
          model.ageing=model.ageing,
		  log.light.borders=Proc.data$log.light.borders,
		  log.irrad.borders=Proc.data$log.irrad.borders,
		  plot.each= plot.each, plot.final= plot.final,
		  suggest.irrad.borders=suggest.irrad.borders))
   # here I check whether tag measures with 10 minutes or more inteval and if so change likelihood correction to FALSE
   if (likelihood.correction=='auto') {
      if (Proc.data$saving.period>=550) {
	     likelihood.correction<-FALSE
         message('likelihood correction switched to FALSE\n')
	  } else {
	  likelihood.correction<-TRUE
	  message('likelihood correction switched to TRUE\n')
	  }
   }
   Proc.data$log.irrad.borders=calibration.parameters$log.irrad.borders
   if (length(calibration.parameters$calib_outliers)>0) {
      Proc.data$FLightR.data$twilights$excluded[which(sapply(Proc.data$FLightR.data$twilights$datetime,
      FUN=function(x) min(abs(calibration.parameters$calib_outliers-as.numeric(x))))<Proc.data$saving.period*24)]<-1
      Proc.data_tmp<-process.twilights(Proc.data$FLightR.data$Data, 
         Proc.data$FLightR.data$twilights[Proc.data$FLightR.data$twilights$excluded==0,],
         measurement.period=Proc.data$measurement.period, saving.period=Proc.data$saving.period,
         impute.on.boundaries=Proc.data$impute.on.boundaries)
         calibration.parameters<-suppressWarnings(get.calibration.parameters(Calibration.periods, Proc.data_tmp, 
         model.ageing=model.ageing, log.light.borders=Proc.data$log.light.borders,
         log.irrad.borders=Proc.data$log.irrad.borders,
		 suggest.irrad.borders=suggest.irrad.borders))
         plot_slopes(calibration.parameters$All.slopes)
		 Proc.data$Twilight.time.mat.dusk<-Proc.data_tmp$Twilight.time.mat.dusk
		 Proc.data$Twilight.time.mat.dawn<-Proc.data_tmp$Twilight.time.mat.dawn
		 Proc.data$Twilight.log.light.mat.dusk<-Proc.data_tmp$Twilight.log.light.mat.dusk
		 Proc.data$Twilight.log.light.mat.dawn<-Proc.data_tmp$Twilight.log.light.mat.dawn
   }
   # second time search for outliers..
   if (length(calibration.parameters$calib_outliers)>0) {
      Proc.data$FLightR.data$twilights$excluded[which(sapply(Proc.data$FLightR.data$twilights$datetime,
      FUN=function(x) min(abs(calibration.parameters$calib_outliers-as.numeric(x))))<Proc.data$saving.period*24)]<-1
      Proc.data_tmp<-process.twilights(Proc.data$FLightR.data$Data, 
         Proc.data$FLightR.data$twilights[Proc.data$FLightR.data$twilights$excluded==0,],
         measurement.period=Proc.data$measurement.period, saving.period=Proc.data$saving.period,
         impute.on.boundaries=Proc.data$impute.on.boundaries)
         calibration.parameters<-suppressWarnings(get.calibration.parameters(
		 Calibration.periods, Proc.data_tmp, 
         model.ageing=model.ageing, log.light.borders=Proc.data$log.light.borders,
         log.irrad.borders=Proc.data$log.irrad.borders,
		 suggest.irrad.borders=suggest.irrad.borders))
         plot_slopes(calibration.parameters$All.slopes)
		 Proc.data$Twilight.time.mat.dusk<-Proc.data_tmp$Twilight.time.mat.dusk
		 Proc.data$Twilight.time.mat.dawn<-Proc.data_tmp$Twilight.time.mat.dawn
		 Proc.data$Twilight.log.light.mat.dusk<-Proc.data_tmp$Twilight.log.light.mat.dusk
		 Proc.data$Twilight.log.light.mat.dawn<-Proc.data_tmp$Twilight.log.light.mat.dawn
   }   

   if (!is.na(fixed.logSlope[1])) calibration.parameters$All.slopes$Parameters$LogSlope[1]<-fixed.logSlope[1]
   if (!is.na(fixed.logSlope[2])) calibration.parameters$All.slopes$Parameters$LogSlope[2]<-fixed.logSlope[2]
   
   # clean the LogSigma values not to use it under the fixed slope scenario
   if (!is.na(fixed.logSlope[1])) calibration.parameters$All.slopes$Parameters$LogSigma<-c(NA,NA)
   
   Calibration<-create.calibration(calibration.parameters$All.slopes,
                 Proc.data,
				 Proc.data$FLightR.data,
				 log.light.borders=Proc.data$log.light.borders, log.irrad.borders=Proc.data$log.irrad.borders,
				 ageing.model=calibration.parameters$ageing.model,
				 location=NA, likelihood.correction=likelihood.correction)
   Calibration$Calibration.periods<-Calibration.periods
   if (return.slopes) Calibration$Calibration.Slopes<-calibration.parameters$All.slopes$Slopes
   return(Calibration)
   }

 
#' find unknown calibration location
#'
#' Functions attempts to find a location where 
#' The function attempts to find a location for a time period assuming animal was not moving. Does not work well will shaded data!
#'
#' @param Proc.data processed data object generated by \code{\link{get.tags.data}}
#' @param calibration.start POSIXct time when stationary period started
#' @param calibration.stop POSIXct time when stationary period ended
#' @param plot plots every iteration
#' @param initial.coords location vector with initial values for location (longitude and latitude). Should be close (+-2000 km from the real location)
#' @param print.optimization do you want every optimization iteration to be printed? If TRUE - Lon, Lat, calibration mean and calibration sd are being printed. Optimization tries to minimize the latter.
#' @param reltol tolerance for optimization, see \code{\link{optim}} for more details
#' @details The idea behind the function is that it tries to minimize variance between slopes for the whole period by optimizing location. It can be seen as an extension of Hill-Ekstrom calibration idea.
#' @return vector with coordinates - longitude and latitude. 
#' @examples
#' #this example takes about 15 minutes to run
#' \dontrun{
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' Proc.data<-get.tags.data(File)
#' plot_slopes_by_location(Proc.data=Proc.data, location=c(5.43, 52.93))
#' abline(v=as.POSIXct("2013-08-20", tz='GMT')) # end of first calibration period
#' abline(v=as.POSIXct("2014-05-05", tz='GMT')) # start of the second calibration period
#' Location<-find.stationary.location(Proc.data, '2013-07-20', '2013-08-20', initial.coords=c(10, 50))
#' }
#' @author Eldar Rakhimberdiev
#' @export		
find.stationary.location<-function(Proc.data, calibration.start,  calibration.stop, plot=TRUE, initial.coords=NULL, print.optimization=TRUE, reltol=1e-4) {
   if (is.null(initial.coords)) stop('current function vesrion requires some inital coordinates to start search, they should not be very close but within few thousand km!')
   ll_function<-function(initial.coords, Proc.data, calibration.start, calibration.stop, plot=TRUE, stage=1) {
   sink("tmp")
        Calibration.period<-data.frame(
        calibration.start=as.POSIXct(calibration.start, tz='GMT'),
        calibration.stop=as.POSIXct(calibration.stop, tz='GMT'),
        lon=initial.coords[1], lat=initial.coords[2])
        #use c() also for the geographic coordinates, if you have more than one calibration location
        # (e. g.,  lon=c(5.43, 6.00), lat=c(52.93,52.94))
       log.light.borders<-Proc.data$log.light.borders
       log.irrad.borders<-Proc.data$log.irrad.borders
       calibration.parameters<-invisible(suppressWarnings(get.calibration.parameters(Calibration.period,
       Proc.data, model.ageing=FALSE, 
	   log.light.borders=log.light.borders,
	   log.irrad.borders=log.irrad.borders, 
       plot.each = FALSE, plot.final = FALSE,
	   suggest.irrad.borders=FALSE)))
	   Twilights_total<-calibration.parameters$Twilights_total
	if (length(calibration.parameters$calib_outliers)>0) {
      Proc.data$FLightR.data$twilights$excluded[which(sapply(Proc.data$FLightR.data$twilights$datetime,
      FUN=function(x) min(abs(calibration.parameters$calib_outliers-as.numeric(x))))<Proc.data$saving.period*24)]<-1
      Proc.data_tmp<-process.twilights(Proc.data$FLightR.data$Data, 
         Proc.data$FLightR.data$twilights[Proc.data$FLightR.data$twilights$excluded==0,],
         measurement.period=Proc.data$measurement.period, saving.period=Proc.data$saving.period,
         impute.on.boundaries=Proc.data$impute.on.boundaries)
      calibration.parameters<-suppressWarnings(get.calibration.parameters(Calibration.period,
      Proc.data_tmp, model.ageing=FALSE, 
	   log.light.borders=log.light.borders,
	   log.irrad.borders=log.irrad.borders, 
       plot.each = FALSE, plot.final = FALSE, 
	   suggest.irrad.borders=FALSE))
	}   
    if (plot) plot_slopes(calibration.parameters$All.slopes)
	     suppressWarnings(sink())
	     percent_excluded<-1-(sum(is.finite(calibration.parameters$All.slopes$Slopes$logSlope))/Twilights_total)

	if (stage==1) {
     	   if (print.optimization) message(paste(initial.coords[1], initial.coords[2], calibration.parameters$All.slopes$Parameters$LogSlope[1], calibration.parameters$All.slopes$Parameters$LogSlope[2], percent_excluded), '\n')
		   #print(table(calibration.parameters$All.slopes$Slopes$Type))
	   if (length(table(calibration.parameters$All.slopes$Slopes$Type))==1) {
		   warning('only_one_twilight_type_left!\n')
		   return(sum(is.finite(calibration.parameters$All.slopes$Slopes$logSlope)))
	   } else {
          return(calibration.parameters$All.slopes$Parameters$LogSlope[2]^2+percent_excluded^2)
		  }
	} else {
	    Dat<-calibration.parameters$All.slopes$Slopes[is.finite(calibration.parameters$All.slopes$Slopes$logSlope),]
	    if (length(table(Dat$Type))==1 | min(table(Dat$Type))<=2) {
	     	   warning('only_one_twilight_type_left!\n')
		       Val<-nrow(Dat)
	   } else {
		   #Dat<-subset(calibration.parameters$All.slopes$Slopes, is.finite(logSlope), select=c(logSlope, Time, Type))
		   #print(Dat)
     	   Val<-log(1/(stats::anova(stats::lm(logSlope~Time+I(Time^2)+Type, data=Dat),stats::lm(logSlope~Time, data=Dat))[2,6])) +percent_excluded^2 #-log(1/(percent_excluded+0.001)
		}
     	if (print.optimization) message(paste(initial.coords[1], initial.coords[2], calibration.parameters$All.slopes$Parameters$LogSlope[1], calibration.parameters$All.slopes$Parameters$LogSlope[2]), Val,  percent_excluded, '\n')
        return(Val)
	   }
   }

   message('stage 1...\n')
   tryCatch(Res<-stats::optim(initial.coords, fn=ll_function, Proc.data=Proc.data, calibration.start=calibration.start, calibration.stop=calibration.stop, plot=plot, control=list(reltol=1e-2)), finally=try(suppressWarnings(sink())))
   message('stage 2...\n')
   
   tryCatch(Res<-stats::optim(Res$par, fn=ll_function, Proc.data=Proc.data, calibration.start=calibration.start, calibration.stop=calibration.stop, plot=plot, stage=2,control=list(reltol=reltol)), finally=try(suppressWarnings(sink())))
   
   return(Res$par)
 }

 
   
#' combines data, calibration and sets up priors
#' 
#' This function is one step before \code{\link{run.particle.filter}}. It combines data, calibration, spatial extent and movement priors and estimates spatial likelihoods that used later in the particle filter.
#' @param Proc.data Processed data object created by \code{\link{get.tags.data}}.
#' @param Grid Spatial grid created by \code{\link{make.grid}}.
#' @param start release location (lat, lon).
#' @param end end of the track location. Will use \code{start} by default. Use NA in case of unknown end point.
#' @param Calibration Calibration object created by \code{\link{make.calibration}}.
#' @param threads number of parallel threads to use. default is -1, which means FLightR will use all available threads except 1. Value 1 will force sequential evaluation
#' @param Decision prior for migration probability values from 0 to 1 are allowed
#' @param Direction Direction prior for direction of migration (in degrees) with 0 pointing to the North
#' @param Kappa concentration parameter for vonMises distribution, 0 means uniform or even distribution. Will set some prior for direction for all the track, so is not recommended to be changed
#' @param M.mean Prior for mean distance travelled between consecutive twilights, km
#' @param M.sd Prior for sd of distance travelled between consecutive twilights, the higher the value is the wider is the the distribution
#' @param likelihood.correction Should likelihood correction estimated during \code{\link{make.calibration}} run be used?
#' @return Object to be uses in the \code{\link{run.particle.filter}}
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' # to run example fast we will cut the real data file by 2013 Aug 20
#' Proc.data<-get.tags.data(File, end.date=as.POSIXct('2013-07-02', tz='GMT'))
#' Calibration.periods<-data.frame(
#'        calibration.start=NA,
#'        calibration.stop=as.POSIXct("2013-08-20", tz='GMT'),
#'        lon=5.43, lat=52.93) 
#'        #use c() also for the geographic coordinates, if you have more than one calibration location
#'        # (e. g.,  lon=c(5.43, 6.00), lat=c(52.93,52.94))
#' print(Calibration.periods)
#' 
#' # NB Below likelihood.correction is set to FALSE for fast run! 
#' # Leave it as default TRUE for real examples
#' Calibration<-make.calibration(Proc.data, Calibration.periods, likelihood.correction=FALSE)
#' 
#' Grid<-make.grid(left=0, bottom=50, right=10, top=56,
#'   distance.from.land.allowed.to.use=c(-Inf, Inf),
#'   distance.from.land.allowed.to.stay=c(-Inf, Inf))
#'
#' all.in<-make.prerun.object(Proc.data, Grid, start=c(5.43, 52.93),
#'                              Calibration=Calibration, threads=2)
#'
#' @author Eldar Rakhimberdiev
#' @export
make.prerun.object<-function(Proc.data, Grid, start, end=start, Calibration, threads=-1, Decision=0.05, Direction=0, Kappa=0, M.mean=300, M.sd=500, likelihood.correction=TRUE) {
   if (length(Decision)>1) stop("Decision has to have length of 1, to sepcify it per twilight, change it in the result object of this function")
   if (length(Direction)>1) stop("Direction has to have length of 1, to sepcify it per twilight, change it in the result object of this function")
   if (length(Kappa)>1) stop("Kappa has to have length of 1, to sepcify it per twilight, change it in the result object of this function")
   if (length(M.mean)>1) stop("M.mean has to have length of 1, to sepcify it per twilight, change it in the result object of this function")
   if (length(M.sd)>1) stop("M.sd has to have length of 1, to sepcify it per twilight, change it in the result object of this function")
   if (any(is.na(Proc.data$FLightR.data$twilights$datetime))) {
      warning('NAs found and deleted in Proc.data$FLightR.data$twilights$datetime\n')
	  Proc.data$FLightR.data$twilights<- stats::na.omit(Proc.data$FLightR.data$twilights)
   }
   Processed.light<-make.processed.light.object(Proc.data$FLightR.data)

   Index.tab<-create.proposal(Processed.light, start=start, Grid=Grid)
   Index.tab$Decision<-Decision # prob of migration
   Index.tab$Direction<- Direction # direction 0 - North
   Index.tab$Kappa<-Kappa # distr concentration 0 means even
   Index.tab$M.mean<- M.mean # distance mu
   Index.tab$M.sd<- M.sd # distance sd

   all.in<-geologger.sampler.create.arrays(Index.tab, Grid, start=start, stop=end)

   all.in$Calibration<-Calibration
   all.in$Data<-Proc.data$FLightR.data

   if (threads!=1 ) {
      Possible.threads<-parallel::detectCores()
      if (threads<=0) Threads<-max(Possible.threads+threads, 1)
      if (threads>0) Threads<-min(Possible.threads,threads)
   } else {
      Threads<-1
   }
   Phys.Mat<-get.Phys.Mat.parallel(all.in, Proc.data$Twilight.time.mat.dusk,
      Proc.data$Twilight.log.light.mat.dusk,
      Proc.data$Twilight.time.mat.dawn,
      Proc.data$Twilight.log.light.mat.dawn,
      threads=Threads, calibration=all.in$Calibration, likelihood.correction=likelihood.correction)
   all.in$Spatial$Phys.Mat<-Phys.Mat
   return(all.in)
}



   
   
##############
# functions below were developed before FLightR 0.3.8
# process.data
process.twilights<-function(All.p, Filtered_tw, measurement.period=60, saving.period=NULL, impute.on.boundaries=FALSE) {
# this function just prepares data for the next steps..
##########
## Dusk
# processing Dusk
Dusk.all<-Filtered_tw$datetime[Filtered_tw$type==2]
Twilight.index.mat.dusk<-sapply(which(All.p$gmt %in% Dusk.all & All.p$type==2), FUN=function(x) (x-24):(x+24))
Twilight.index.mat.dusk<-apply(Twilight.index.mat.dusk, c(1,2), FUN=function(x) ifelse (x>0, x, NA))
Max.Index<-nrow(All.p)
Twilight.index.mat.dusk<-apply(Twilight.index.mat.dusk, c(1,2), FUN=function(x) ifelse (x>Max.Index, NA, x))

Twilight.time.mat.dusk<-apply(Twilight.index.mat.dusk, c(1,2), FUN=function(x) as.numeric(All.p$gmt[x]))
Twilight.time.mat.dusk<-apply(Twilight.time.mat.dusk, c(1,2), FUN=function(x) ifelse(is.finite(x), x, 0))
Twilight.log.light.mat.dusk<-apply(Twilight.index.mat.dusk, c(1,2), FUN=function(x) log(All.p$light[x]))
#Twilight.log.light.mat.dusk<-apply(Twilight.index.mat.dusk, c(1,2), FUN=function(x) All.p$light[x])
Twilight.log.light.mat.dusk<-apply(Twilight.log.light.mat.dusk, c(1,2), FUN=function(x) ifelse(is.finite(x), x, -1))


Twilight.time.mat.dusk<-Twilight.time.mat.dusk-(saving.period-measurement.period)

# processing Dawn
Dawn.all<-Filtered_tw$datetime[Filtered_tw$type==1]

Twilight.index.mat.dawn<-sapply(which(All.p$gmt %in% Dawn.all & All.p$type==1), FUN=function(x) (x-24):(x+24))
Twilight.index.mat.dawn<-apply(Twilight.index.mat.dawn, c(1,2), FUN=function(x) ifelse (x>0, x, NA))
Max.Index<-nrow(All.p)
Twilight.index.mat.dawn<-apply(Twilight.index.mat.dawn, c(1,2), FUN=function(x) ifelse (x>Max.Index, NA, x))

Twilight.time.mat.dawn<-apply(Twilight.index.mat.dawn, c(1,2), FUN=function(x) as.numeric(All.p$gmt[x]))
Twilight.time.mat.dawn<-apply(Twilight.time.mat.dawn, c(1,2), FUN=function(x) ifelse(is.finite(x), x, 0))
Twilight.log.light.mat.dawn<-apply(Twilight.index.mat.dawn, c(1,2), FUN=function(x) log(All.p$light[x]))
#Twilight.log.light.mat.dawn<-apply(Twilight.index.mat.dawn, c(1,2), FUN=function(x) All.p$light[x])
Twilight.log.light.mat.dawn<-apply(Twilight.log.light.mat.dawn, c(1,2), FUN=function(x) ifelse(is.finite(x), x, -1))

Res<-list(Twilight.time.mat.dusk=Twilight.time.mat.dusk, Twilight.log.light.mat.dusk=Twilight.log.light.mat.dusk, Twilight.time.mat.dawn=Twilight.time.mat.dawn,  Twilight.log.light.mat.dawn=Twilight.log.light.mat.dawn, measurement.period=measurement.period, saving.period=saving.period, impute.on.boundaries=impute.on.boundaries)
return(Res)
}




get.Irradiance<-function(alpha, r=6378, s=6.9, intigeo.template.correction=FALSE) {
	# function from Ekstrom 2007
	erf <- function(x) 2 * stats::pnorm(x * sqrt(2)) - 1
	## (see Abramowitz and Stegun 29.2.29)
	## and the so-called 'complementary error function'
	erfc <- function(x) 2 * stats::pnorm(x * sqrt(2), lower = FALSE)
	## and the inverses
	#erfinv <- function (x) qnorm((1 + x)/2)/sqrt(2)
	#erfcinv <- function (x) qnorm(x/2, lower = FALSE)/sqrt(2)
	
	Res<-alpha
	u<-sqrt(r/(2*s))*sin(alpha)
	Res[(u<=0)]<-(exp(-u^2)/(1+erf(-u)))[(u<=0)]
	Res[(u>0)]<-(exp(-u^2)/(erfc(u)))[(u>0)]
	
	if (intigeo.template.correction) Res[Res>0.001]=exp(log(Res[Res>0.001])*1.279+log(Res[Res>0.001])^2*0.091)
	
	return(Res)
}


logger.template.calibrarion.internal<-function( Twilight.time.mat.Calib.dawn, Twilight.log.light.mat.Calib.dawn, Twilight.time.mat.Calib.dusk, Twilight.log.light.mat.Calib.dusk, positions=NA, plot.each=TRUE, plot.final=TRUE,log.light.borders=NA,  log.irrad.borders=c(-8, 1.3), impute.on.boundaries=FALSE) {

   oldpar <- graphics::par(no.readonly = TRUE)    
   on.exit(graphics::par(oldpar))   
	# =================
	# in this function I'll add a new lnorm calibration...
	#
	# now I'll try to do them both..
	if (!is.list(positions)) {
		if (length (positions) !=2) stop("the positions should be a list with $dawn  and $dusk dataframes with pairs of coordiantes OR just one pair\n")
		dawn=matrix(ncol=2, nrow=dim(Twilight.time.mat.Calib.dawn)[2])
		dawn[,1]<-positions[1]
		dawn[,2]<-positions[2]
		#
		dusk=matrix(ncol=2, nrow=dim(Twilight.time.mat.Calib.dusk)[2])
		dusk[,1]<-positions[1]
		dusk[,2]<-positions[2]
		#
		positions=list(dawn=dawn, dusk=dusk)
	} 
	Twilight.time.mat.Calib.dawn<-Twilight.time.mat.Calib.dawn[-25,]
	Twilight.log.light.mat.Calib.dawn<-Twilight.log.light.mat.Calib.dawn[-25,]
	Twilight.time.mat.Calib.dusk<-Twilight.time.mat.Calib.dusk[-25,]
	Twilight.log.light.mat.Calib.dusk<-Twilight.log.light.mat.Calib.dusk[-25,]
	# let's try to create Calib.data.all first!!
	Calib.data.dawn<-data.frame()
	if (plot.each) graphics::par(ask=FALSE)
	for (dawn in 1:dim(Twilight.time.mat.Calib.dawn)[2]) {
    message("\r checking dawn", dawn)
		#Twilight.solar.vector<-solar(as.POSIXct(Twilight.time.mat.Calib.dawn[, dawn], tz="gmt", origin="1970-01-01"))
		Data<-check.boundaries(positions$dawn[dawn,], Twilight.solar.vector=NULL,  Twilight.log.light.vector = Twilight.log.light.mat.Calib.dawn[,dawn], plot=plot.each, verbose=FALSE,  log.light.borders=log.light.borders, log.irrad.borders=log.irrad.borders, dusk=FALSE, Twilight.time.vector=Twilight.time.mat.Calib.dawn[, dawn], impute.on.boundaries=impute.on.boundaries)
		if (length(Data)==0) {
		message("\r dawn", dawn, "was excluded from the calibration")
		} else {
		Calib.data.dawn<-rbind(Calib.data.dawn, cbind(LogLight=Data[,1], LogIrrad=Data[,2], Day=dawn, Time=Data[,3], Elevs=Data[,4]))	
		}		
	}
   if (nrow(Calib.data.dawn)>0) {Calib.data.dawn$type<-"Dawn"}

	Calib.data.dusk<-data.frame()
	for (dusk in 1:dim(Twilight.time.mat.Calib.dusk)[2]) {
    message("\r checking dusk", dusk )

		#Twilight.solar.vector<-solar(as.POSIXct(Twilight.time.mat.Calib.dusk[, dusk], tz="gmt", origin="1970-01-01"))
		Data<-check.boundaries(positions$dusk[dusk,], Twilight.solar.vector=NULL,  Twilight.log.light.vector=Twilight.log.light.mat.Calib.dusk[,dusk], plot=plot.each, verbose=FALSE,  log.light.borders=log.light.borders, log.irrad.borders=log.irrad.borders, dusk=TRUE,  Twilight.time.vector=Twilight.time.mat.Calib.dusk[, dusk], impute.on.boundaries=impute.on.boundaries)
		if (length(Data)==0) {
		message("\r dusk", dusk, "was excluded from the calibration")
		} else {
		Calib.data.dusk<-rbind(Calib.data.dusk, cbind(LogLight=Data[,1], LogIrrad=Data[,2], Day=dim(Twilight.time.mat.Calib.dawn)[2]+dusk, Time=Data[,3], Elevs=Data[,4]))
		}
	}
	if (nrow(Calib.data.dusk)>0) {Calib.data.dusk$type<-"Dusk"}
	message('\r')
	if (nrow(Calib.data.dusk)==0 & nrow(Calib.data.dawn)==0) stop('Error: None of twilights could happen in the location... Something went wrong at the previous stages')
	Calib.data.all<-rbind(Calib.data.dawn, Calib.data.dusk)
	if (plot.final) {
		graphics::plot(LogLight~LogIrrad,data=Calib.data.all, type="n")
		for (Day in unique(Calib.data.all$Day)) {
		graphics::lines(LogLight~LogIrrad, data=Calib.data.all[Calib.data.all$Day==Day,], type="b", pch="+", col=grDevices::rainbow(length(unique(Calib.data.all$Day)))[Day])
		}
	}
	
	Calib.data.all$fDay<-as.factor(Calib.data.all$Day)
	message('\r\n')
	message('\r\n')
	return(Calib.data.all)	
}

# ok now we need template.calibration.function..
logger.template.calibration<-function(Twilight.time.mat.Calib.dawn, Twilight.log.light.mat.Calib.dawn, Twilight.time.mat.Calib.dusk, Twilight.log.light.mat.Calib.dusk, time.shift=0, positions, log.light.borders=NA,  log.irrad.borders=c(-15, 50), plot.each=TRUE, plot.final=TRUE, impute.on.boundaries=FALSE) {
	
	Calibration.original<-logger.template.calibrarion.internal(Twilight.time.mat.Calib.dawn+time.shift, Twilight.log.light.mat.Calib.dawn, Twilight.time.mat.Calib.dusk+time.shift, Twilight.log.light.mat.Calib.dusk, positions=positions, log.light.borders=log.light.borders,  log.irrad.borders= log.irrad.borders, plot.each=plot.each, plot.final=plot.final,impute.on.boundaries=impute.on.boundaries)

	return(Calibration.original)
	}

	

get.calib.param<-function(Calib.data.all, plot=FALSE, calibration.type=NULL) {

if (is.null(calibration.type)) calibration.type="parametric.slope"
message("calibration method used:", calibration.type, "\n")

Calib.data.all$fTwilight<-Calib.data.all$fDay

#============================
# we should do some outlier tests...
cur.data<-Calib.data.all
p.lm1<-stats::lm(LogLight~fTwilight/LogIrrad,data=cur.data)
cur.slope<-matrix(stats::coef(p.lm1),ncol=2) # foo[,2] - slope
cur.slope<-as.data.frame(cur.slope)
Diag<-diag(stats::vcov(p.lm1))
Diag.slopes<-Diag[which(sapply(strsplit(names(Diag), ":"), length)==2)]

cur.slope$sd[!is.na(cur.slope[,2])]<-sqrt(Diag.slopes)
cur.slope$type<-stats::aggregate(cur.data[,"type"],by=list(Day=cur.data$fTwilight),FUN=function(x) unique(x))[,2]
cur.slope$time<-stats::aggregate(cur.data[,"Time"],by=list(Day=cur.data$fTwilight),FUN=function(x) x[1])[,2]

cur.slope$Duration<-stats::aggregate(cur.data[,"Time"],by=list(Day=cur.data$fTwilight),FUN=function(x) max(x)-min(x))[,2]
names(cur.slope)[2]<-"slope"

#===========================
# ok I've decided to make separate lm
Slopes<-c()
Slopes.sd<-c()
Intercept=c()
Sigma=c()
Type=c()
Time=c()
Elevs<-c()
Day<-c()
for (i in (unique(cur.data$fTwilight))) {
# lm
#plot(LogLight~LogIrrad, data=cur.data[cur.data$fTwilight==i,])
Data<-cur.data[cur.data$fTwilight==i,]
Data<-Data[Data$LogLight>0,]

Lm<-stats::lm(LogLight~LogIrrad,data=Data)
Cur_slope_mean<-stats::coef(Lm)[2]
Cur_slope_sd<-sqrt(stats::vcov(Lm)[4])


if (calibration.type=="nonparametric.slope") {
	# Here I want to physically estimate slopes

	All_slopes_cur<-diff(Data$LogLight)/diff(Data$LogIrrad)
	All_slopes_mean<-mean(All_slopes_cur)
	All_slopes_sd<-stats::sd(All_slopes_cur)
	# exclude outliers
	All_slopes_diff<-abs(All_slopes_cur-All_slopes_mean)
	Slopes_outliers<-which(All_slopes_diff>3*All_slopes_sd)
	if (length(Slopes_outliers)>0) {
		All_slopes_cur<-All_slopes_cur[-Slopes_outliers]
		All_slopes_mean<-mean(All_slopes_cur)
		All_slopes_sd<-stats::sd(All_slopes_cur)
	}
	Cur_slope_mean<-All_slopes_mean
	Cur_slope_sd<-All_slopes_sd
}

Slopes<-c(Slopes, Cur_slope_mean)
Slopes.sd<-c(Slopes.sd, Cur_slope_sd)
Intercept=c(Intercept, stats::coef(Lm)[1])
Sigma<-c(Sigma, summary(Lm)$sigma)
Type=c(Type, cur.data$type[cur.data$fTwilight==i][1])
Time<-c(Time,ifelse(cur.data$type[cur.data$fTwilight==i][1]=="Dusk", max(cur.data$Time[cur.data$fTwilight==i]), min(cur.data$Time[cur.data$fTwilight==i])))
Elevs=c(Elevs, mean(cur.data$Elevs[cur.data$fTwilight==i], na.rm=TRUE))
Day=c(Day, cur.data$Day[cur.data$fTwilight==i][1])
# nonparametric slope estimation

}
#plot(cur.slope$slope~Slopes)
#plot(cur.slope$sd~Slopes.sd)

# ok, we now want to replace old slopes with new ones.
# but first we want to add 0.5 instead of NA - at least this is how it works now inside the functions and how it was done in a simulations
cur.slope$slope<-Slopes
cur.slope$Intercept<-Intercept
cur.slope$Sigma<-Sigma
#Slopes.sd[is.na(Slopes.sd)]<-sd.fast
cur.slope$sd<-Slopes.sd

names(cur.slope)[2]<-"slope"
#=====================
if (plot) graphics::hist(log(cur.slope$slope))


cur.slope.int<-cur.slope[is.finite(log(cur.slope$slope)),]

Parameters<-list(Intercept=c(mean(cur.slope.int$Intercept, na.rm=TRUE), stats::sd(cur.slope.int$Intercept, na.rm=TRUE)), LogSlope=c(mean(log(cur.slope.int$slope), na.rm=TRUE), stats::sd(log(cur.slope.int$slope), na.rm=TRUE)), LogSigma=c(mean(log(cur.slope.int$Sigma[!is.na(cur.slope.int$sd)])), stats::sd(log(cur.slope.int$Sigma[!is.na(cur.slope.int$sd)]))), mean.of.individual.slope.sigma=mean(cur.slope.int$sd, na.rm=TRUE), calibration.type=calibration.type)

#cur.slope$time<-aggregate(cur.data[,"Time"],by=list(Day=cur.data$fTwilight),FUN=function(x) x[1])[,2]
#cur.slope$time<-aggregate(cur.data[,"Time"],by=list(Day=cur.data$fTwilight),FUN=mean)[,2]

Res<-list(Parameters=Parameters, Slopes=data.frame(Slope=cur.slope$slope, Time=Time, Intercept=cur.slope$Intercept, Sigma=cur.slope$Sigma, Slopes.sd=Slopes.sd, Type=Type, Elevs=Elevs, Day=Day))
return(Res) 
}


plot_slopes<-function(all.slopes, ylim=NULL, xlim=NULL) {
   oldpar <- graphics::par(no.readonly = TRUE)   
   on.exit(graphics::par(oldpar))    

all.slopes$Slopes<-all.slopes$Slopes[all.slopes$Slopes$Slope>0,]
if (is.null(xlim)) {
   graphics::plot(log(all.slopes$Slopes$Slope)~as.POSIXct(all.slopes$Slopes$Time, tz="GMT", origin="1970-01-01"), type="n", main="red - dawn, black - dusk", xlab="time", ylab="log(Slope)", ylim=ylim, las=1)
} else {
   graphics::plot(log(all.slopes$Slopes$Slope)~ as.POSIXct(all.slopes$Slopes$Time, tz="GMT", origin="1970-01-01"), type="n", main="red - dawn, black - dusk", xlab="time", ylab="log(Slope)", ylim=ylim, las=1, xlim=  as.POSIXct(xlim, tz='GMT'))

}

graphics::lines(log(Slope)~ as.POSIXct(Time, tz="GMT", origin="1970-01-01"), data=all.slopes$Slopes[all.slopes$Slopes$Type=="Dusk",])
graphics::points(log(Slope)~ as.POSIXct(Time, tz="GMT", origin="1970-01-01"), data=all.slopes$Slopes[all.slopes$Slopes$Type=="Dusk",], pch="+")
graphics::points(log(Slope)~ as.POSIXct(Time, tz="GMT", origin="1970-01-01"), data=all.slopes$Slopes[all.slopes$Slopes$Type=="Dawn",], pch="+", col="red")
graphics::lines(log(Slope)~ as.POSIXct(Time, tz="GMT", origin="1970-01-01"), data=all.slopes$Slopes[all.slopes$Slopes$Type=="Dawn",], col="red")
#invisible()
return(NULL)
}


get.calibration.parameters<-function(Calibration.periods, Proc.data, model.ageing=TRUE, log.light.borders=log(c(2, 63)),  log.irrad.borders=c(-7,1.5), plot.each=FALSE, plot.final=TRUE, calibration.type=NULL, suggest.irrad.borders=FALSE) {

# ageing.model this option should be used only in case there are two periods of calibration or there is one but very long.
if (nrow(Calibration.periods)==1 & model.ageing) warning("you have only one calibration period ageing estimation is unreliable and should be turned to F!!!")

Calibration.periods.int<-Calibration.periods
Calibration.periods.int$calibration.start<-as.numeric(Calibration.periods.int$calibration.start)
Calibration.periods.int$calibration.stop<-as.numeric(Calibration.periods.int$calibration.stop)


Dusk.calib.days<-c()
Dawn.calib.days<-c()
Positions=list(dawn=c(), dusk=c())
#sum(apply(Calibration.periods.int, 1, FUN=function(x) x[2]-x[1]))/3600/24
for (i in 1:nrow(Calibration.periods)) {
   Dusk.calib.days.cur<-NULL
   Dusk.calib.days.cur<-which(Proc.data$Twilight.time.mat.dusk[25,]>Calibration.periods.int$calibration.start[i] & Proc.data$Twilight.time.mat.dusk[25,] < Calibration.periods.int$calibration.stop[i])
   if (length(Dusk.calib.days.cur) > 0) {
      Dusk.calib.days<-unique(c(Dusk.calib.days, Dusk.calib.days.cur))
      dusk.cur=matrix(ncol=2, nrow=length(Dusk.calib.days.cur))
      dusk.cur[,1]<-Calibration.periods.int$lon[i]
      dusk.cur[,2]<-Calibration.periods.int$lat[i]
      Positions$dusk<-rbind(Positions$dusk, dusk.cur)
   }
   Dawn.calib.days.cur<-NULL
   Dawn.calib.days.cur<-which(Proc.data$Twilight.time.mat.dawn[25,]>Calibration.periods.int$calibration.start[i] & Proc.data$Twilight.time.mat.dawn[25,] < Calibration.periods.int$calibration.stop[i] )
   if (length(Dawn.calib.days.cur) > 0) {
      Dawn.calib.days<-unique(c(Dawn.calib.days,Dawn.calib.days.cur ))
      dawn.cur=matrix(ncol=2, nrow=length(Dawn.calib.days.cur))
      dawn.cur[,1]<-Calibration.periods.int$lon[i]
      dawn.cur[,2]<-Calibration.periods.int$lat[i]
      Positions$dawn<-rbind(Positions$dawn, dawn.cur)
   }
}
Twilight.time.mat.Calib.dusk<-Proc.data$Twilight.time.mat.dusk[,Dusk.calib.days]
Twilight.log.light.mat.Calib.dusk<-Proc.data$Twilight.log.light.mat.dusk[,Dusk.calib.days]

## Dawn

Twilight.time.mat.Calib.dawn<-Proc.data$Twilight.time.mat.dawn[,Dawn.calib.days]
Twilight.log.light.mat.Calib.dawn<-Proc.data$Twilight.log.light.mat.dawn[,Dawn.calib.days ]

Twilights_total<-length(c(Dusk.calib.days,Dawn.calib.days))
#-------------------#
  if (suggest.irrad.borders) {
     Calib.data.all<-logger.template.calibration(Twilight.time.mat.Calib.dawn, Twilight.log.light.mat.Calib.dawn, Twilight.time.mat.Calib.dusk, Twilight.log.light.mat.Calib.dusk, positions=Positions, log.light.borders=log.light.borders,  log.irrad.borders=c(-10,10), plot.each=plot.each, plot.final=plot.final, impute.on.boundaries=Proc.data$impute.on.boundaries)

     log.irrad.borders<-suggest.irrad.boundaries(Calib.data.all) 
     warning('!!! log.irrad.borders were updated to new values:', log.irrad.borders[1], log.irrad.borders[2], '\n')
  }

Calib.data.all<-logger.template.calibration(Twilight.time.mat.Calib.dawn, Twilight.log.light.mat.Calib.dawn, Twilight.time.mat.Calib.dusk, Twilight.log.light.mat.Calib.dusk, positions=Positions, log.light.borders=log.light.borders,  log.irrad.borders=log.irrad.borders, plot.each=plot.each, plot.final=plot.final, impute.on.boundaries=Proc.data$impute.on.boundaries)

All.slopes<-suppressWarnings(get.calib.param(Calib.data.all, plot=FALSE, calibration.type=calibration.type))

All.slopes$Slopes$logSlope<-suppressWarnings(log(All.slopes$Slopes$Slope))

if (model.ageing) {
All.slopes.int<-All.slopes
All.slopes.int$Slopes<-All.slopes.int$Slopes[is.finite(All.slopes.int$Slopes$logSlope),]
Time.start=min(c(Twilight.time.mat.Calib.dusk, Twilight.time.mat.Calib.dawn))
Model=stats::lm(logSlope~I(Time-Time.start), data=All.slopes.int$Slopes)
Model$Time.start<-Time.start
calib_outliers<-All.slopes.int$Slopes$Time[which(abs(stats::residuals(Model))>3*stats::sd(stats::residuals(Model))*sqrt(length(stats::residuals(Model))))]
Res<-list(calib_outliers=calib_outliers, ageing.model=Model, All.slopes=All.slopes,log.irrad.borders=log.irrad.borders)
} else {
All.slopes.int<-All.slopes

All.slopes.int$Slopes<-All.slopes.int$Slopes[is.finite(All.slopes.int$Slopes$logSlope),]

calib_outliers<-All.slopes.int$Slopes$Time[which(abs(All.slopes.int$Slopes$logSlope-mean(All.slopes.int$Slopes$logSlope, na.rm=TRUE))>3*stats::sd(All.slopes.int$Slopes$logSlope, na.rm=TRUE))]

Res<-list(calib_outliers=calib_outliers, All.slopes=All.slopes, log.irrad.borders=log.irrad.borders, Twilights_total=Twilights_total)
}
return(Res)
}


suggest.irrad.boundaries<-function(Calib.data.all, leave.out=0.01) {

   Model<-stats::lm(LogLight~LogIrrad + fDay, data=Calib.data.all)

   MinLight<-round(min(Calib.data.all$LogLight),1)
   Left.border<- -(stats::coef(Model)[1] + stats::quantile(stats::coef(Model)[3:length(stats::coef(Model))], leave.out) - MinLight)*stats::coef(Model)[2]

   MaxLight<-round(max(Calib.data.all$LogLight),1)
   Right.border<-(MaxLight- (stats::coef(Model)[1] + stats::quantile(stats::coef(Model)[3:length(stats::coef(Model))], 1-leave.out))) *stats::coef(Model)[2]
   return(c(Left.border, Right.border))
}


make_likelihood_correction_function<-function(calib_log_mean, calib_log_sd, cur_sd_range, npoints=300, plot=FALSE) {
   cur_mean_range<-c(exp(calib_log_mean)-5*exp(calib_log_sd), exp(calib_log_mean)+5*exp(calib_log_sd))
   Res<-c()
   for (i in 1:npoints) {
   	  message('\r simulation:  ', round(i/npoints*100), '%', sep='')
	  message('\r')
      cur_mean<-stats::runif(300, cur_mean_range[1], cur_mean_range[2])
	  cur_sd<-stats::runif(1, cur_sd_range[1], cur_sd_range[2])
      Cur_max_real<- stats::optimize(f=function(x) 
	        mean(stats::dlnorm(stats::rnorm(200000,x , cur_sd),calib_log_mean, calib_log_sd)),
			interval=cur_mean_range, maximum = TRUE)$maximum
      Res<-rbind(Res, cbind(calib_log_mean, calib_log_sd, cur_sd, cur_mean_max=Cur_max_real))
   }
   Res<-as.data.frame(Res)
   Res$Corr <- exp(Res$calib_log_mean)-Res$cur_mean_max
   message('\r  estimating correction function...')
   message('\r')

   MvExp<-mgcv::gamm(Corr~s(cur_sd), data=Res, weights=nlme::varExp(form =~ cur_sd))
   # now we check for the outliers..
   Resid<-stats::resid(MvExp$lme, type='normalized')
   Outliers<-NULL
   Outliers<-which(abs(Resid)>3*stats::sd(Resid))
   if (length(Outliers)>0) {
      Res<-Res[-Outliers,]
      MvExp<-mgcv::gamm(Corr~s(cur_sd), data=Res, weights=nlme::varExp(form =~ cur_sd))
   }
   XX<-seq(cur_sd_range[1], cur_sd_range[2], by=0.01)
   c_fun<-stats::approxfun(mgcv::predict.gam(MvExp$gam, newdata=data.frame(cur_sd=XX)) ~XX)

   if (plot) {
      graphics::plot(Res[,4]~Res[,3], pch='+', ylab='cur_mean', xlab='cur_sd')
      graphics::lines(exp(Res$calib_log_mean[1])-mgcv::predict.gam(MvExp$gam, newdata=data.frame(cur_sd=XX))~XX, col='red')
	  }
	Out<-list(c_fun=c_fun, Res=Res)
	message('\r')
	message('\n')
	return(Out)
}

create.calibration<-function( All.slopes, Proc.data, FLightR.data, location, log.light.borders, log.irrad.borders, ageing.model=NULL, likelihood.correction=TRUE) {
# Now we create 'parameters' object that will have all the details about the calibration
Parameters<-All.slopes$Parameters # LogSlope # 
Parameters$measurement.period<-Proc.data$measurement.period 
Parameters$saving.period<-Proc.data$saving.period 
Parameters$log.light.borders<-log.light.borders # these are the boundaries in which one should use the BAS tag.. for the other types of tags they will be different.
Parameters$min.max.values<-c(min(FLightR.data$Data$light), max(FLightR.data$Data$light))
Parameters$log.irrad.borders=log.irrad.borders
Parameters$impute.on.boundaries=Proc.data$impute.on.boundaries
Parameters$location=location

Parameters$LogSlope_1_minute<-Parameters$LogSlope
 if (is.null(ageing.model)) {
time_correction_fun= eval(parse(text=paste("function (x,y) return(", Parameters$LogSlope[1], ")")))
lat_correction_fun<-function(x, y, z) return(0)
} else {
time_correction_fun= function(x, y) return(0)
lat_correction_fun<-eval(parse(text=paste("function (x,y) return(",  stats::coef(ageing.model)[1], "+", stats::coef(ageing.model)[2], "* (y-",ageing.model$Time.start," ))")))
}
c_fun=NULL
if (likelihood.correction) {
     if (is.na(All.slopes$Parameters$LogSigma[1])) {
	     cur_sd_range=c(0, 1)
	 } else {
	     cur_sd_range=c(0, max(exp(All.slopes$Parameters$LogSigma[1]+5*All.slopes$Parameters$LogSigma[2]),1))
	 }
     c_fun=make_likelihood_correction_function(Parameters$LogSlope[1], Parameters$LogSlope[2], cur_sd_range=cur_sd_range)$c_fun
}
Calibration<-list(Parameters=Parameters, time_correction_fun=time_correction_fun, lat_correction_fun=lat_correction_fun, c_fun=c_fun)

return(Calibration)
}


correct.hours<-function(datetime) {
	# this function is supposed to correct hours by adding some amount of them.
   hours <- as.numeric(format(datetime,"%H"))+as.numeric(format(datetime,"%M"))/60
   hours[as.POSIXlt(datetime, tz='GMT')$isdst==1] <-hours[as.POSIXlt(datetime, tz='GMT')$isdst==1]-1
	cor <- rep(NA, 24)
	for(i in 0:23){
		cor[i+1] <- max(abs((c(hours[1],hours)+i)%%24 - 
		            (c(hours,hours[length(hours)])+i)%%24),na.rm=TRUE)
	}
	hours <- (hours + (which.min(round(cor,2)))-1)%%24
	return(hours)
}

make.result.list<-function(Data, raw.X, raw.Y) {
	Int<-c(min(raw.Y), max(raw.Y))
	Index<-sapply(raw.X,  FUN=function(x) which.min(abs(as.numeric(Data$d$gmt-x))))
	Res<-Data$d[Index,]
	# now I want to adjust time and light
	Res$light<-stats::approx(x=Data$d$gmt, y=Data$d$light, xout=raw.X)$y
	Res$gmt<-raw.X
	Res$Hour<-raw.Y
	Result<-list(Data=Res)
	Result$Data$gmt.adj<-Result$Data$gmt
	Result$Data$gmt<- as.POSIXct(Result$Data$gmt, tz="GMT", origin="1970-01-01")
	Result$Data$gmt.adj<- as.POSIXct(Result$Data$gmt.adj, tz="GMT", origin="1970-01-01")
	return(Result)
}

make.processed.light.object<-function(FLightR.data) {
# dusk
raw.Y.dusk<-correct.hours(FLightR.data$twilights$datetime[FLightR.data$twilights$type==2 & FLightR.data$twilights$excluded==0])
raw.X.dusk<-as.numeric(FLightR.data$twilights$datetime[FLightR.data$twilights$type==2 & FLightR.data$twilights$excluded==0])
Data_tmp<-list(d=FLightR.data$Data)
Result.Dusk<-make.result.list(Data_tmp, raw.X.dusk, raw.Y.dusk)

# dusk
raw.Y.dawn<-correct.hours(FLightR.data$twilights$datetime[FLightR.data$twilights$type==1 & FLightR.data$twilights$excluded==0])
raw.X.dawn<-as.numeric(FLightR.data$twilights$datetime[FLightR.data$twilights$type==1 & FLightR.data$twilights$excluded==0])
Result.Dawn<-make.result.list(Data_tmp, raw.X.dawn, raw.Y.dawn)

# combine
processed.light<-list(Final.dusk=Result.Dusk, Final.dawn=Result.Dawn)

return(processed.light)
}

# this create proposal is corrected for the missing loess
create.proposal<-function(processed.light, start=c(-98.7, 34.7), end=NA, Grid) {
	if (is.na(end)) end=start
All.Days<-seq(min(min(processed.light$Final.dusk$Data$gmt), min(processed.light$Final.dawn$Data$gmt)),max(max(processed.light$Final.dusk$Data$gmt), max(processed.light$Final.dawn$Data$gmt)), by="days")

## Probabilty.of.migration
All.Days.extended<-seq(min(All.Days)-86400, max(All.Days)+86400, by="days")

# ver 0.4.2
# correct twilights for polar periods..
   start_no_polar<-start
   start_no_polar[2]<-min(abs(start_no_polar[2]), 65.73) * sign(start_no_polar[2])

# these are all potential twilights that we would have for the start point...
Potential.twilights<-sort(c(maptools::sunriset(matrix(start_no_polar, nrow=1), All.Days.extended, direction="sunrise", POSIXct.out=TRUE)[,2], maptools::sunriset(matrix(start_no_polar, nrow=1), All.Days.extended, direction="sunset", POSIXct.out=TRUE)[,2]))
# now we need to add dask or dawn to it..
Sunrise<-maptools::sunriset(matrix(start_no_polar, nrow=1), Potential.twilights[1], direction="sunrise", POSIXct.out=TRUE)[,2]
Sunrises<-c(Sunrise-(3600*24), Sunrise, Sunrise+(3600*24))
Sunset<-maptools::sunriset(matrix(start_no_polar, nrow=1), Potential.twilights[1], direction="sunset", POSIXct.out=TRUE)[,2]
Sunsets<-c(Sunset-(3600*24), Sunset, Sunset+(3600*24))
First.twilight<-ifelse(which.min(c(min(abs(difftime(Sunrises,Potential.twilights[1], units="mins"))), min(abs(difftime(Sunsets,Potential.twilights[1], units="mins")))))==1, "dawn", "dusk")

Index.tab<-data.frame(Date=Potential.twilights)

if (First.twilight=="dusk") {
	Index.tab$Dusk<-rep(c(TRUE,FALSE), times=length(All.Days.extended))
} else {
	Index.tab$Dusk<-rep(c(FALSE,TRUE), times=length(All.Days.extended))
}

Index.tab$Curr.mat<-NA # this will be NA if no data and row number if there are..
Index.tab$Real.time<- as.POSIXct(NA, tz="GMT")
Index.tab$time<- as.POSIXct(NA, tz="GMT")

######################################
# !!!! this will not work if bird will move for over 12 time zones!!!
# Dusk
for (i in 1:length(processed.light$Final.dusk$Data$gmt)) {
	Row2write<-which.min(abs(difftime(processed.light$Final.dusk$Data$gmt[i], Index.tab$Date[Index.tab$Dusk==TRUE], units="mins")))
	Index.tab$time[Index.tab$Dusk==TRUE][Row2write]<-processed.light$Final.dusk$Data$gmt[i]
	Index.tab$Real.time[Index.tab$Dusk==TRUE][Row2write]<-processed.light$Final.dusk$Data$gmt.adj[i]
	Index.tab$Curr.mat[Index.tab$Dusk==TRUE][Row2write]<-i
}

# Dawn
for (i in 1:length(processed.light$Final.dawn$Data$gmt)) {
	Row2write<-which.min(abs(difftime(processed.light$Final.dawn$Data$gmt[i], Index.tab$Date[Index.tab$Dusk==FALSE], units="mins")))
	Index.tab$time[Index.tab$Dusk==FALSE][Row2write]<-processed.light$Final.dawn$Data$gmt[i]
	Index.tab$Real.time[Index.tab$Dusk==FALSE][Row2write]<-processed.light$Final.dawn$Data$gmt.adj[i]
	Index.tab$Curr.mat[Index.tab$Dusk==FALSE][Row2write]<-i
}
# cutting empty ends..
while (is.na(Index.tab$Curr.mat[1])) Index.tab<-Index.tab[-1,]
while (is.na(Index.tab$Curr.mat[nrow(Index.tab)])) Index.tab<-Index.tab[-nrow(Index.tab),]
Index.tab$Point<-NA
First.Point<-which.min(sp::spDistsN1(Grid[,1:2], start,  longlat=TRUE))
Index.tab$Point[1]<-First.Point
#
# I decided that Curr.mat is not needed anymore
#Index.tab$Main.Index[,-which(names(Index.tab$Main.Index)=="Curr.mat")]
# this need to double checked
Index.tab$yday<-as.POSIXlt(Index.tab$Date, tz="GMT")$yday

 return(Index.tab)
}



geologger.sampler.create.arrays<-function(Index.tab, Grid, start, stop=start) {
	# this function wil take in Index.table with all the proposals and will create an array for case of missed data..
	
	# the main feature is that it can account for missing data..
	Index.tab.old<-Index.tab
	Index.tab$proposal.index<-NA
	Index.tab$proposal.index[Index.tab$Dusk==TRUE]<-"Dusk"
	Index.tab$proposal.index[Index.tab$Dusk==FALSE]<-"Dawn"
	
	Missed.twilights<-which(is.na(Index.tab$Curr.mat))
	
	if (length(Missed.twilights)>0) {
	Deleted.Rows.Count<-0
		for (i in Missed.twilights) {
			i=i-Deleted.Rows.Count
			Index.tab$proposal.index[i-1]<-"Comb"
			Index.tab$Decision[i-1]<-(Index.tab$Decision[i-1]+Index.tab$Decision[i])/2 # 
			if (Index.tab$Direction[i-1] != Index.tab$Direction[i])	{
				Index.tab$Direction[i-1]<-0
				Index.tab$Kappa[i-1]<-0
			}
			Index.tab$M.mean[i-1]<-Index.tab$M.mean[i-1]+Index.tab$M.mean[i]
			Index.tab$M.sd[i-1]<-sqrt((Index.tab$M.sd[i-1])^2+(Index.tab$M.sd[i])^2)
			Index.tab<-Index.tab[-i,]
			Deleted.Rows.Count=Deleted.Rows.Count+1
		}
	}
	output<-list()
	Matrix.Index.Table<-	Index.tab[,which(names(Index.tab) %in% c("Decision", "Direction", "Kappa", "M.mean", "M.sd", "yday", "Real.time", "time", "Dusk", "Loess.se.fit", "Loess.n")) ]
	# I  also remove last line as bird was not flying after last twilight
	Matrix.Index.Table<-Matrix.Index.Table[-nrow(Matrix.Index.Table),]

	
	# main index will have the same amount of rows we have in twilight matrices without first.. 
	Main.Index<-c()
	Main.Index$Biol.Prev<-1:nrow(Matrix.Index.Table)
	Main.Index$proposal.index<-Index.tab$proposal.index[-nrow(Index.tab)]
	Main.Index<-as.data.frame(Main.Index)
	
	Indices<-list(Matrix.Index.Table=Matrix.Index.Table, Main.Index=Main.Index)
	
	output$Indices=Indices
	
	# ok now we need to add to output all other parts..
	
	output$Spatial<-list(Grid=Grid) #[,1:2]

	output$Spatial$Behav.mask<-as.integer(Grid[,3])

	output$Spatial$start.point<-which.min(sp::spDistsN1(Grid[,1:2], start,  longlat=TRUE))
    
	output$Spatial$start.location<-start
	if (!is.na(stop[1]) ) {
	   output$Spatial$stop.point<-which.min(sp::spDistsN1(Grid[,1:2], stop,  longlat=TRUE))
	   output$Spatial$stop.location<-stop
	} else {
	   output$Spatial$stop.point<-NA
	   output$Spatial$stop.location<-NA
	}
	return(output)
}

Try the FLightR package in your browser

Any scripts or data that you put into this service are public.

FLightR documentation built on July 6, 2021, 5:08 p.m.