Difference between revisions of "OpasnetUtils/Drafts"

library(OpasnetUtils)

############ Shuffles columns of a data.frame so that they match a pre-defined correlation matrix

#### THIS SHOULD BE UPDATED FOR OVARIABLES AS WELL: SHUFFLING ACROSS Iter WITH
#### CORRELATION MATRIX ACROSS DEFINED INDICES AND THEIR LOCATIONS. OTHER INDICES ARE
#### KEPT UNCHANGED, SO THE SHUFFLING HAS TO HAPPEN WITHIN EACH UNIQUE LOCATION COMBINATION.

correlvar <- function(
  vars, # multivariable object to be correlated.
  Sigma # covariance matrix wanted.
) {
  
  # Method from http://www.r-bloggers.com/easily-generate-correlated-variables-from-any-distribution-without-copulas/
  require(MASS)
  mu <- rep(0,ncol(vars))
  rawvars <- as.data.frame(mvrnorm(n = nrow(vars), mu = mu, Sigma = Sigma))
  out <- as.data.frame(
    lapply(
      1:ncol(vars), 
      FUN = function(i, vars, rawvars) {
        pvars <- rank(rawvars[[i]], ties.method = "random")
        tmp <- sort(vars[[i]]) # Make sure you start with ordered data.
        tmp <- tmp[pvars] # Order based on correlated ranks
        return(tmp)
      },
      vars = vars,
      rawvars = rawvars
    )
  )
  colnames(out) <- colnames(vars)
  return(out)
}

##################### Forgets decisions so that decision indices will be recreated.

forgetDecisions <- function() {
	for(i in ls(envir = openv)) {
		if("dec_check" %in% names(openv[[i]])) openv[[i]]$dec_check <- FALSE
	}
	return(cat("Decisions were forgotten.\n"))
}

#################### Perus-ggplot ovariableille

oggplot <- function (
	ova, 		# ovariable to be plotted
	x, 		# Index for x axis
	weight = NULL, 	# Index for y axis (default: result column)
	fill = NULL, 	# Index for colour code
	base_size = NULL, # Base size for graph font (object BS or 18 is used if not given)
	turnx = FALSE, 	# Turn x axis labels vertically?
	binwidth = NULL	# Width of bins
) 
{
	if(is.null(base_size)) if(exists("BS")) base_size <- BS else base_size <- 18
	if(exists("translate")) {
		ova <- translate(ova)
		x <- translate(x)
		weight <- translate(weight)
		fill <- translate(fill)
	}
	if (is.null(weight)) 
		weight <- paste(ova@name, "Result", sep = "")
	if("Iter" %in% colnames(ova@output)) 
		ova <- oapply(ova, cols = "Iter", FUN = mean)
	plo <- ggplot(ova@output, aes_string(x = x, fill = fill, weight = weight)) +
		theme_gray(base_size = base_size)  
	if(any(ova@output[[weight]] >= 0)) {
		plo <- plo + geom_bar(
			data = subset(ova@output, ova@output[[weight]] >= 0),
			position = "stack",
			binwidth = binwidth
		)
	}
	if(any(ova@output[[weight]] < 0)) {
		plo <- plo + geom_bar(
			data = subset(ova@output, ova@output[[weight]] < 0), 
			position = "stack",
			binwidth = binwidth
		) + geom_hline(aes(yintercept = 0))
	}

	if (turnx) 
		plo <- plo + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
	return(plo)
}

################## Sähkön hinta tunneittain

#price <- opbase.data(ident="op_en7353")
#temperature <- opbase.data("op_en6315.2014_5_2015")
#temperature$Date <- substr(temperature$Date, 0, 11)
#price$Date <- substr(price$Date, 0, 11)
#mon <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
#for (i in mon) {
# price$Date <- gsub(i, which(mon == i), as.character(price$Date))
#}
#for (i in mon) {
# temperature$Date <- gsub(i, which(mon == i), as.character(temperature$Date))
#}
#price$Hours <- substr(price$Hours, 0, 2)
#price$Hours <- paste(price$Hours, ":00:00", sep="")
#temperature$Time <- paste(temperature$Time, ":00", sep="")
#as.character(temperature$Result)
#as.numeric(temperature$Result)
#cut(temperature$Result, breaks = c(-21, -18, -15, -12, -9, -6, -3, 0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30), #include.lowest=TRUE)
#DateTime <- as.POSIXct(paste(temperature$Date, temperature$Time), format="%Y-%m-%d %H:%M:%S")
#DateHours <- as.POSIXct(paste(price$Date, price$Hours), format="%Y-%m-%d %H:%M:%S")

ograph <- function( # Määritellään yleisfunktio peruskuvaajan piirtämiseen.
	ovariable, 
	x, 
	y = character(), 
	type = character(), 
	other = character(),
	fill = NA, 
	...
) {
	cat("This function ograph is depreciated. Use oggplot instead.\n")

	if(class(ovariable) == "ovariable")  {
		if(nrow(ovariable@output) == 0) ovariable <- EvalOutput(ovariable)
		data <- ovariable@output
		title <- ovariable@name
		if(length(y) == 0) y <- paste(title, "Result", sep = "")
	} else {
		data <- ovariable
		title <- character()
		if(length(y) == 0) y <- "Result"
	}
	if(length(type) == 0) {
		if("Iter" %in% colnames(data)) type <- geom_boxplot() else type <- geom_bar(stat = "identity")
	}
	out <- ggplot(data, aes_string(x = x, y = y, fill = fill)) # Määritellään kuvan sarakkeet
	out <- out + type
	out <- out + theme_grey(base_size=24) # Fontin kokoa suurennetaan
	out <- out + labs(
			title	= title,
			y = paste(unique(data[[paste(title, "Yksikkö", sep = "")]]), sep = "", collapse = ", ")
		)
	out <- out + theme(axis.text.x = element_text(angle = 90, hjust = 1)) # X-akselin tekstit käännetään niin että mahtuvat
	if(length(other) != 0) out <- out + other
	return(out)
}

# fillna takes a data.frame and fills the cells with NA with each level in that column.
# fillna was updated in OpasnetUtils and therefore removed from here.


## collapsemarg is a placeholder for a better functionality within CollapseMarginals.
## It takes an ovariable, and summarises all indices in cols using tapply and a user-defined function.
## However, you can also use function "pick" to select locations defined in a list picks found in indices cols.
## Function "unkeep" simply drops the unkept indices without any other operation.
## The output is an ovariable with the same name as the input.
## This was first created for [[:op_fi:Radonin terveysvaikutukset]]

collapsemarg <- function(variable, cols, fun = "sum", picks = list(), ...) { # cols is a character vector, while probs is a list
	out <- dropall(variable@output)
	marginals <- colnames(out)[variable@marginal]

	if(tolower(fun) == "unkeep") { # The function must be a string, otherwise this row will fail.
		out <- out[!colnames(out) %in% cols]
	} else {
		if(tolower(fun) == "pick") {
			for(i in cols) {
				out <- out[out[[i]] %in% picks[[match(i, cols)]] , ]
			}
			cols <- "" # Those locations that were picked are still marginals.
		} else {
			margtemp <- colnames(out)[colnames(out) %in% marginals & !colnames(out) %in% cols] 
			# You must leave at least one index.

			out <- as.data.frame(as.table(tapply(result(variable), out[margtemp], fun)))
			out <- out[!is.na(out$Freq) , ]

			colnames(out)[colnames(out) == "Freq"] <- ifelse(
				length(variable@name) == 0, 
				"Result", 
				paste(variable@name, "Result", sep = "")
			)
		}
	}

	variable@output <- out
	variable@marginal <- colnames(out) %in% marginals & ! colnames(out) %in% cols
	
	return(variable)
}

MyPointKML <- function( # The function creates a KML fille from a SpatialPointsDataFrame object.
	obj = NULL, # Spatial object with the data. A SpatialPointsDataFrame.
	kmlname = "", # Name of the KML fille (does this show on the map?)
	kmldescription = "", # Description of the KML fille (does this show on the map?)
	name = NULL, # Name for each datapoint (vector with the same length as data in obj).
	description = "", # Descrtion of each datapoint (vector with the same length as data in obj).
	icon = "http://maps.google.com/mapfiles/kml/pal4/icon24.png", # Icon shown on pin (?)
	col=NULL # I don't know what this does.
) {

cat("This function MyPointKML is depreciated. Use google.point_kml in OpasnetUtilsExt instead.\n")

    if (is.null(obj)) 
        return(list(header = c("<?xml version=\"1.0\" encoding=\"UTF-8\"?>", 
            "<kml xmlns=\"http://earth.google.com/kml/2.2\">", 
            "<Document>", paste("<name>", kmlname, "</name>", 
                sep = ""), paste("<description><![CDATA[", kmldescription, 
                "]]></description>", sep = "")), footer = c("</Document>", 
            "</kml>")))
    if (class(obj) != "SpatialPointsDataFrame") 
        stop("obj must be of class 'SpatialPointsDataFrame' [package 'sp']")
    if (is.null(name)) {
        name = c()
        for (i in 1:nrow(obj)) name <- append(name, paste("site", 
            i))
    }
    if (length(name) < nrow(obj)) {
        if (length(name) > 1) 
            warning("kmlPoints: length(name) does not match nrow(obj). The first name will be replicated.")
        name <- rep(name, nrow(obj))
    }
    if (length(description) < nrow(obj)) {
        if (length(description) > 1) 
            warning("kmlPoints: length(description) does not match nrow(obj). The first description will be replicated.")
        description <- rep(description, nrow(obj))
    }
    if (length(icon) < nrow(obj)) {
        if (length(icon) > 1) 
            warning("kmlPoints: length(icon) does not match nrow(obj). Only the first one will be used.")
        icon <- icon[1]
    }

	# This is some kind of a colour definition
	
	col2kmlcolor <- function(col) 
		paste(rev(sapply(
			col2rgb(col, TRUE), 
			function(x) sprintf("%02x", x))
		), collapse = "")
	
    kml <- kmlStyle <- ""
    
	# Create the KML fille.
	
	kmlHeader <- c("<?xml version=\"1.0\" encoding=\"UTF-8\"?>","<kml xmlns=\"http://earth.google.com/kml/2.2\">", "<Document>")
    kmlFooter <- c("</Document>", "</kml>")
    
	# Create rows to the KML fille from data in obj.
	
    for (i in 1:nrow(obj)) {
        point <- obj[i, ]
        pt_style <- paste("#style", ifelse(length(icon) == 1, 1, i), sep = "")
        kml <- append(kml, "<Placemark>")
        kml <- append(kml, paste(
		"  <description><![CDATA[",
		name[i],
		": ",
		description[i], 
		"]]></description>", 
		sep = ""
	))
	#kml <- append(kml, "<Style><IconStyle>")
	#kml <- append(kml, paste("<color>", col2kmlcolor(col[i]), "</color>", sep =""))
	#kml <- append(kml, paste("  <Icon><href>", icon, "</href></Icon>", sep = ""))
	#kml <- append(kml, "<scale>0.300000</scale>")
	#kml <- append(kml, "</IconStyle></Style>")
        kml <- append(kml, "  <Point>")
        kml <- append(kml, "    <coordinates>")
        kml <- append(kml, paste(point@coords[1], point@coords[2], sep = ","))
        kml <- append(kml, "    </coordinates>")
        kml <- append(kml, "  </Point>")
        kml <- append(kml, "</Placemark>")
    }
    
    return(paste(paste(c(kmlHeader, kmlStyle, kml, kmlFooter), sep = "", collapse = "\n"), collapse="\n", sep = ""))
}

ova2spat <- function( # This function converts an ovariable into a SpatialPointsDataFrame.
	ovariable, # An evaluated ovariable that has coordinate indices.
	coords, # The names of the coordinate indices as a character vector, first x then y.
	proj4string # Projection identifier or specification as character string. See http://spatialreference.org/
) {
	temp <- ovariable@output
	
	# Transform coordinates into numeric format.
	
	for(i in coords) { 
		if(is(temp[[i]], "factor"))    temp[[i]] <- levels(temp[[i]])[temp[[i]]]
		if(is(temp[[i]], "character")) temp[[i]] <- as.numeric(temp[[i]])
	}
	
	# Define the coordinate points first, then add other ovariable output to it.
	
	sp <- SpatialPoints(temp[coords], CRS(proj4string))
	out <- SpatialPointsDataFrame(sp, temp[!colnames(temp) %in% coords])

	#Transform the projection to longitude-latitude system.
	
	epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
	out <- spTransform(out,epsg4326String)

	return(out)
}

# MyRmap is a function for creating static Google maps as png.
# It is based on MyMap function without the "file destination" parameter
# Requires RgoogleMaps package

MyRmap <- function (
	shp, # a spatial data object
	plotvar, # Name of the column that has the values to be illustrated on the map
	pch = 19, # Shape of the point (19: circle)
	cex = 0.3, # Size of the point
	legend_title = "", # Title of the legend
	legend_position = "topleft",
	numbins = 8, # Number of colour bins in graph
	center, # center of the map
	size = c(640, 480), # size of the map. This produces the right dimensions in Opasnet.
	MINIMUMSIZE = FALSE, 
	RETURNIMAGE = TRUE, 
	GRAYSCALE = FALSE, 
	NEWMAP = TRUE, 
	zoom, 
	verbose = 1, 
	...
) {
	plotvar	<- shp[[plotvar]]
	plotclr	<- brewer.pal(numbins, "Spectral")
	classes	<- classIntervals(plotvar, numbins, style = "quantile")
	colcode	<- findColours(classes, plotclr)
	latR	<- shp@coords[ , 2]
	lonR	<- shp@coords[ , 1]

	#get the bounding box:

	bb <- qbbox(lat = latR, lon = lonR)

	if (missing(zoom)) 
		zoom <- min(MaxZoom(latR, lonR, size))
	if (missing(center)) {
		lat.center <- mean(latR)
		lon.center <- mean(lonR)
	}
	else {
		lat.center <- center[1]
		lon.center <- center[2]
	}
	if (MINIMUMSIZE) {
		ll <- LatLon2XY(latR[1], lonR[1], zoom) # I think the latR and lonR are used here differently than how they
		ur <- LatLon2XY(latR[2], lonR[2], zoom) # are used elsewhere. Thus, if MINIMUMSIZE = TRUE, you may see problems.
		cr <- LatLon2XY(lat.center, lon.center, zoom)
		ll.Rcoords <- Tile2R(ll, cr)
		ur.Rcoords <- Tile2R(ur, cr)
		if (verbose > 1) {
			cat("ll:")
			print(ll)
			print(ll.Rcoords)
			cat("ur:")
			print(ur)
			print(ur.Rcoords)
			cat("cr:")
			print(cr)
		}
		size[1] <- 2 * max(c(ceiling(abs(ll.Rcoords$X)), ceiling(abs(ur.Rcoords$X)))) + 1
		size[2] <- 2 * max(c(ceiling(abs(ll.Rcoords$Y)), ceiling(abs(ur.Rcoords$Y)))) + 1

		if (verbose) cat("new size: ", size, "\n")
	}

	MyMap <- GetMap(
		center = c(lat.center, lon.center),
		zoom = zoom, 
		size = size, 
		RETURNIMAGE = RETURNIMAGE, 
		GRAYSCALE = GRAYSCALE, 
		verbose = verbose, 
		...
	)

	PlotOnStaticMap(MyMap) # Plot an empty map.

	PlotOnStaticMap( # Plot the data points on the map.
		MyMap,
		lat = latR,
		lon = lonR,
		pch = pch,
		cex = cex,
		col = colcode,
		add = T
	)

	legend( # Plot the legend on the map.
		legend_position, 
		legend	= names(attr(colcode, "table")),
		title	= legend_title, 
		fill	= attr(colcode, "palette"), 
		cex	= 1.0, 
		bty	= "y",
		bg	= "white"
	)
}

MyPlotKML <- function(
	shp, # a SpatialPointDataFrame object.
	result = "Result", # The name of  result column in shp.
	rasterization = TRUE, # Whether to rasterize the data or not.
	ncols = 32, # Number or columns in the raster.
	nrows = 32, # Number of rows in the raster.
	fun = mean # function to aggregate data points to the raster.
) {	

cat("Consider merging this function MyPolotKML with google.show_raster_on_maps in OpasnetUtilsExt.\n")

	if(rasterization) {
		#Create blank raster
		rast		<- raster()

		#Set raster extent to that of point data
		extent(rast)	<-extent(shp)

		#Choose number of columns and rows
		ncol(rast)	<- ncols
		nrow(rast)	<- nrows

		#Rasterize point data
		rast2	<- rasterize(shp, rast, shp[[result]], fun = fun)

	}

	start <- 0 # min(shp[[result]])
	end <- max(shp[[result]])
	steps <- approx(c(start,end),n=6)$y
	colors <- rev(rainbow(length(steps), start=0, end=0.50))

	# Create the colorstrip below the map.

	par(mfrow=c(6,1), mar=c(3,1,0,1), cex = 1.5)

	colorstrip <- function(colors, labels)
	{
		count <- length(colors)
		image(
			matrix(1:count, count, 1), 
			col = colors, 
			ylab = "", 
			axes = FALSE
		)
		axis(1,approx(c(0, 1), n=length(labels))$y, labels)
	}

	colorstrip(colors, steps)

	#Plot data

	google.show_raster_on_maps(rast2, col = colors, style = "height:500px;")
}


# Merge all but show_bins largest bins of indices cols to 'Other'.

truncateIndex <- function( # Truncates indices to contain only the largest index bins.
		obj, # ovariable to use.
		cols, # names of the columns to truncate.
		bins = rep(10, length(cols)), # Number of bins to show, including Others. Smallest locations will be lumped to bin "Other".
		sum_others = TRUE # Should "Other" be summed to maintain marginal status
) {
	if(nrow(obj@output) == 0) stop("Ovariable ", obj@name, " not evaluated.\n")
	test <- oapply(abs(obj), INDEX = cols, sum, na.rm = TRUE)
	if(length(cols) > 1 & length(bins) == 1) bins <- rep(bins, length(cols))
	
	for(i in 1:length(cols))
	{
		test2 <- oapply(test, INDEX = cols[i], sum)
		test2@output <- test2@output[result(test2) > 0 , ]
		
		temp <- as.factor(obj@output[[cols[i]]])
		location_weight_order <- order(result(test2), decreasing = TRUE)
		keeps <- test2@output[[cols[i]]][location_weight_order[0:min(bins[i] - 1, nrow(test2@output))]]
		levels(temp)[!levels(temp) %in% keeps] <- "Other"
		temp <- factor(temp, levels = c(levels(temp)[levels(temp) != "Other"], "Other"))
		obj@output[[cols[i]]] <- temp
	}
	# After changing some locations to "Other", sum along indices to avoid problems
	if(sum_others) {
		obj <- oapply(obj, cols = "", FUN = sum, na.rm = TRUE)
	}
	return(obj)
}

findrest <- function (X, cols, total = 1)
{
	# findrest input is an ovariable that can be integrated over indices cols to result in total.
	# Often it is used with uncertain fractions. One (often the largest) fraction is 
	# omitted or NA, and it is replaced by whatever is missing from the total.
	if (nrow(X@output) == 0) 
		X <- EvalOutput(X)
	rescol <- paste(X@name, "Result", sep = "")
	marginals <- colnames(X@output)[X@marginal]

	# temp is the amount that is still missing from the total.
	temp <- total - oapply(X, cols = cols, FUN = function(x) sum(x, na.rm = TRUE))

	# Remove old rescol because it would cause trouble later.
	if(rescol != paste(temp@name, "Result", sep = "")) temp <- unkeep(temp, cols = rescol) 
	colnames(temp@output)[colnames(temp@output) == paste(temp@name, "Result", sep = "")] <- "tempResult" 
	temp@name <- "temp" # This is to make sure that merge works.

	out <- merge(X, temp)@output

	# Replace missing values with values from temp.
	out[[rescol]] <- ifelse(is.na(out[[rescol]]), out$tempResult, out[[rescol]]) # Result comes from temp
	out$tempResult <- NULL
	X@output <- out
	X@marginal <- colnames(X@output) %in% marginals
	return(X)
}

timing <- function(dat, timecol = NA, weeks = 6, tz = "EET") 
# timing converts character or numeric inputs into times using a few default formats
{
	if(is.data.frame(dat)) # Turn dat into a data.frame in all cases.
	{
		timesall <- (dat[timecol])
	} else {
		timesall <- data.frame(Timecol = dat)
		timecol <- "Timecol"
	}

	temprow <- character() # A whole row of timesall collapsed into a string.
	
	for(i in 1:nrow(timesall))
	{
		temprow[i] <- tolower(paste(t(timesall)[ , i], collapse = ""))
	}

	weekdays <- data.frame(
		Name = c("su", "ma", "ti", "ke", "to", "pe", "la", 
		"sun", "mon", "tue", "wed", "thu", "fri", "sat",
		"sön", "mon", "tis", "ons", "tor", "fre", "lör"),
		Number = rep(0:6, times = 3)
	)
	
	# day is the number of the weekday in case of repeating events.

	day <- rep(NA, nrow(timesall))
	for(i in 1:nrow(weekdays)) 
	{
		temp1 <- grepl(weekdays$Name[i], temprow) # Is any weekday mentioned in temprow?
		day <- ifelse(is.na(day) & temp1, weekdays$Number[i], day) # Find the weekday number.
	}

	# repall are repeating events, x are non-repeating events.
	
	repall <- timesall[!is.na(day) , ]
	xall <- timesall[is.na(day) , ]
	starting <- NA
	
	########### First, change non-repeating timedates.
	# Change timedate into POSIXct assuming formats 15.3.2013 or 2013-03-15 and 15:24 or 15.24.
	# Note! 13 and 2013 mean 1.1.2013 and 3.2013 means 1.3.2013
		
	if(nrow(xall) > 0)
	{
		xout <- data.frame(Datrow = (1:nrow(timesall))[is.na(day)]) # Row numbers from dat.
			
		for(j in colnames(xall))
		{
			x <- xall[[j]]
			if(is.factor(x)) x <- levels(x)[x]
			
			x <- ifelse(grepl("^[0-9][0-9]$", x), paste("20", x, sep = ""), x)
			x <- ifelse(grepl("^[0-9][0-9][0-9][0-9]$", x), paste("1.1.", x, sep = ""), x)
			x <- ifelse(grepl("^[0-9].[0-9][0-9][0-9][0-9]$", x), paste("1.", x, sep = ""), x)
			x <- ifelse(grepl("^[0-9][0-9].[0-9][0-9][0-9][0-9]$", x), paste("1.", x, sep = ""), x)

			temp <- x
			x <- as.POSIXct(temp, format = "%d.%m.%Y %H:%M", tz = tz)
			x <- ifelse(!is.na(x), x, as.POSIXct(temp, format = "%d.%m.%Y %H.%M", tz = tz))
			x <- ifelse(!is.na(x), x, as.POSIXct(temp, format = "%Y-%m-%d %H:%M", tz = tz))
			x <- ifelse(!is.na(x), x, as.POSIXct(temp, format = "%Y-%m-%d %H.%M", tz = tz))
			x <- ifelse(!is.na(x), x, as.POSIXct(temp, format = "%d.%m.%Y", tz = tz))
			x <- ifelse(!is.na(x), x, as.POSIXct(temp, format = "%d.%m.%y", tz = tz))

			xout <- cbind(xout, X = as.POSIXct(x, origin = "1970-01-01", tz = tz)) 
			# Intermediate values turn into numeric, therefore turned back to POSIX.

			starting <- min(c(starting, xout$X), na.rm = TRUE)
			
			colnames(xout)[colnames(xout) == "X"] <- j
		}
	}
	
	############## Then, change weekly repeating timedates.
	
	# repall must have format ma 9:00 or TUESDAY 8.24. Weekday is case-insensitive and can be abbrevieated.
	# The start and end times are assumed to be on the same day. The name of day can be on either column.

	if(nrow(repall) > 0)
	{
		
		for(j in colnames(repall))
		{
			
			reptime <- gsub("[[:alpha:] ]", "", repall[[j]]) # Remove alphabets and spaces.
			reptime <- gsub("\\.", ":", reptime)
			reptime <- strsplit(reptime, split = ":")

			temp2 <- numeric()
			for(i in 1:length(reptime)) 
			{
				temp2[i] <- as.numeric(reptime[[i]][1]) * 3600 + as.numeric(reptime[[i]][2]) * 60
			}

			reptime <- temp2

			if(is.na(starting)) starting <- paste(format(Sys.Date(), format = "%Y"), "-01-01", sep = "")
			starting <-  as.POSIXlt(starting, origin = "1970-01-01") # First day of year.
			starting$mday <- starting$mday - as.numeric(format(starting, format = "%w")) + day[!is.na(day)] # Previous Sunday plus weekdaynumber.

			reps <- data.frame()
			temp3 <- starting 

			for(i in 1:weeks)
			{
				reps <- rbind(reps, data.frame(
					Datrow = (1:nrow(timesall))[!is.na(day)],
					X = as.POSIXct(temp3, origin = "1970-01-01", tz = tz) + reptime
				))
				temp3$mday <- temp3$mday + 7 # Make a weekly event.
			}
			colnames(reps)[colnames(reps) == "X"] <- j

			if(j == colnames(repall)[1]) repout <- reps else repout <- cbind(repout, reps[j])
		}
	}
	out <- rbind(xout, repout)
	
	if(is.data.frame(dat))
	{
		dat$Datrow <- 1:nrow(dat)
		dat <- dat[!colnames(dat) %in% timecol]
		out <- merge(dat, out)
		out <- out[colnames(out) != "Datrow"]
	} else {
		out <- out[[timecol]]
	}
	
	return(out)
}


# Funktio makeTimeline ottaa tapahtumalistauksen ja rakentaa siitä aikajanan. Parametrit:
# event  =  data.framena tapahtumalistaus, joka sisältää ainakin alku- ja loppuajan (Alku, Loppu) ja
# tapahtumatiedon sekä mahdollisesti toiston ja keston (Toisto = aikaväli päivinä, 
# Kesto = viimeinen tapahtuma-aika).
# timeformat = jos TRUE, oletetaan POSIX-muotoiseksi ja muutetaan operointia varten sekunneiksi.
# Jos FALSE, oletetaan reaaliluvuksi ja operoidaan suoraan luvuilla.

makeTimeline <- function(event, timeformat = TRUE) { 
	
	# eventiin luodaan Toisto ja Kesto, jos jompikumpi puuttuu.
	if(!all(c("Toisto", "Kesto") %in% colnames(event)))
	{
		event$Toisto <- NA
		event$Kesto <- NA
	}
	
	if(timeformat)
	{
		for(m in c("Alku", "Loppu", "Kesto"))
		{
			event[[m]] <- as.double(as.POSIXct(event[[m]])) # Muutetaan aika sekunneiksi.
		}
		event$Toisto <- event$Toisto * 3600 * 24 # Muutetaan Toisto päivistä sekunneiksi.
	}
	
	# Jos on puuttuvaa tietoa kestosta tai toistosta, korvataan inertillä datalla.
	test <- event$Toisto == 0 | event$Kesto == 0 | is.na(event$Toisto) | is.na(event$Kesto)
	event$Toisto <- ifelse(test, 1, event$Toisto)
	event$Kesto <- ifelse(test, event$Alku, event$Kesto)

	# Luodaan aikajanan alkupiste. Eventrow otetaan aikanaan tyhjältä riviltä.

	timeline <- data.frame(Time = min(event$Alku), Eventrow = nrow(event)+1) 
	
	for(i in 1:nrow(event)) { # Toista jokaiselle havaintoriville
		times <- 0:floor((event$Kesto[i] - event$Alku[i]) / event$Toisto[i]) # Toistojen määrä

		# Toista tapahtumaa Kestoon asti.
		temp <- data.frame(
			Time = event$Alku[i] + times * event$Toisto[i], 
			End = event$Loppu[i] + times * event$Toisto[i], 
			EventrowStart = i
		)
		
		timeline <- merge(timeline, temp[ , c("Time", "EventrowStart")], all = TRUE) # Lisätään tapahtumat aikajanaan.
		colnames(temp) <- c("Remove", "Time", "EventrowEnd") # Muutetaan otsikot, koska nyt halutaan mergata loppuhetket aikajanaan.
		timeline <- merge(timeline, temp[, c("Time", "EventrowEnd")], all = TRUE)
		for(j in 2:nrow(timeline)) {
		# Tämä luuppi käy aikajanan läpi ja täydentää tapahtumat.
		
			# Ensin kaikkiin uusiin aikapisteisiin jatketaan sitä aiempaa toimintaa, joka oli menossa edellisessä pisteessä.
			if(is.na(timeline$Eventrow[j])) 
			{
				timeline$Eventrow[j] <- timeline$Eventrow[j-1]
			}

			# Sitten jatketaan uutta toimintaa niihin uusiin pisteisiin, jotka eivät ole loppupisteitä.
			if(is.na(timeline$EventrowStart[j]) & is.na(timeline$EventrowEnd[j])) 
			{
				timeline$EventrowStart[j] <- timeline$EventrowStart[j-1]
			}
		}
			
		# Jos uutta toimintaa on olemassa, statuksena käytetään sitä, muutoin statusta eli aiempaa toimintaa.
		timeline$Eventrow <- ifelse(!is.na(timeline$EventrowStart), timeline$EventrowStart, timeline$Eventrow)
		
		# Leikataan turhat sarakkeet pois ja siirrytään seuraavalle event-riville
		timeline <- timeline[ , c("Time", "Eventrow")] 
	}

	event <- rbind(event, rep(NA, ncol(event))) # Lisää rivi loppuhetkeä varten
	event$Alku[nrow(event)] <- max(timeline$Time)
	event$Eventrow <- row(event)[ , 1]

	timeline <- merge(timeline, event) # Yhdistetään Statuksen eli eventin rivinumeron avulla.
	timeline <- timeline[!colnames(timeline) %in% c("Eventrow", "Alku", "Loppu", "Toisto", "Kesto")]
	timeline <- timeline[order(timeline$Time) , ]

	if(timeformat) timeline$Time <- as.POSIXct(timeline$Time, origin = "1970-01-01")
	
	return(timeline)
}

# Calculate the cumulative impact of the events on building stock to given years

timepoints <- function(
	# Function timepoints takes an event list and turns that into existing crosscutting situations at 
	# timepoints defined by years. The output will have index Time.
	# In other words, this will integrate over obstime at specified timepoints.
	X, # X must be an ovariable with a column of the same name as obstime.
	obstime, # obstime must be a single-column data.frame of observation times. 
	sumtimecol = TRUE # Should the timecol be summed up?
	# obstime and timecol may be numeric (by coercion) or POSIXt.
) {
	timecol <- colnames(obstime)
	marginals <- colnames(X@output)[X@marginal]

	# tapply (and therefore possibly oapply) changes continuous indices to factors! Must change back by hand.
	if("factor" %in% c(class(obstime[[timecol]]), class(X@output[[timecol]]))) {
		X@output[[timecol]] <- as.numeric(as.character(X@output[[timecol]]))
		obstime[[timecol]] <- as.numeric(as.character(obstime[[timecol]]))
	}

	out <- data.frame()
	if(sumtimecol) by <- setdiff(marginals, timecol) else by <- marginals

	for(i in obstime[[timecol]]) {
		temp <- X@output[X@output[[timecol]] <= i , ]
		if(nrow(temp) > 0) {
			temp <- aggregate(
				temp[paste(X@name, "Result", sep = "")],
				by = temp[colnames(temp) %in% by],
				FUN = sum
			)
		}
		if(nrow(temp) > 0) out <- rbind(out, data.frame(Time = i, temp))
	}
	X@output <- out
	X@marginal <- colnames(out) %in% c("Time", marginals) # Add Time to marginal

	return(X)
}

# ana2ova takes in variable tables from Analytica (produced with Copy Table). The output is a data.frame in long format.
ana2ova <- function(dat) 
{
	i <- 1 # Number of indices
	cols <- character() # Names of indices
	locs <- list() # Vectors of locations for each index.
	out <- data.frame() # output
	repeat{ # Find out cols and locs from dat.
		temp <- strsplit(dat[i + 1], split = "\t")[[1]]
		cols[i] <- temp[1]
		if(length(temp) == 1) break
		locs[[i]] <- temp[2:length(temp)]
		i <- i + 1
	}
	locs[[i]] <- NA # The innermost index should have its place in locs even if the locations are not yet known.

	inner <- (1:length(dat))[dat %in% cols[i]] # Places where the innermost index is mentioned (the data starts from the next row).
	if(!all(dat[inner - 1] == dat[i]) & i > 1) stop("Structure incorrect\n")
	if(length(inner) == 1) # len is the length of the innermost index.
	{
		len <- length(dat) - i - 1
	} else {
		len <- inner[2] - inner[1] - i
	}

	for(k in inner) # Go through each table with the inner and second most inner indices.
	{
		ins <- strsplit(dat[(k + 1):(k + len)], split = "\t") # Make a data.frame.
		ins <- data.frame(matrix(unlist(ins), nrow = len, byrow = TRUE))

		if(i == 1) cn <- "value" else cn <- locs[[i - 1]]
		colnames(ins) <- c(cols[i], cn) # Give location names to columns
		if(i > 1) ins <- melt(ins, id.vars = cols[i], variable.name = cols[i - 1]) # If several columns.
		if(i > 2) # If there are several tables.
		{
			uplocs <- data.frame(Removethis = 0)
			for(l in 1:(i - 2))
			{
				uplocs[[cols[l]]] <- strsplit(dat[k - i + l], split = "\t")[[1]][2]
			}
			out <- rbind(out, cbind(uplocs, ins)) # Collects all tables into a data.frame.
		} else {
			out <- ins
		}
	}
	out$Removethis <- NULL
	colnames(out)[colnames(out) == "value"] <- "Result"
	return(out)
}

orbind2 <- function( # Like orbind but the value is an ovariable.
	o1, # ovariable whose slots are used in the value.
	o2, # ovariable
	use_fillna = FALSE, # Do we use fillna to fill in the NA values in indices?
	warn = "" # What warning is given if fillna is used?
) {
	x <- unkeep(o1 * 1, prevresults = TRUE, sources = TRUE)
	y <- unkeep(o2 * 1, prevresults = TRUE, sources = TRUE)
	xmarg <- colnames(x@output)[x@marginal]
	ymarg <- colnames(y@output)[y@marginal]
	for(i in xmarg) x@output[[i]] <- as.factor(x@output[[i]])
	for(i in ymarg) y@output[[i]] <- as.factor(y@output[[i]])
	out <- o1
	out@output <- orbind(x, y)

	if(use_fillna) {
		b <- character()
		for(i in colnames(out@output)[out@marginal]) {if(any(is.na(out@output[[i]]))) b <- c(b, i)}
		if(length(b) > 0) {
		out@output <- fillna(out@output, b)
		warning(warn, "\nMissing values had to be filled by function fillna in indices: ", b, "\n")
		}
	}

	colnames(out@output)[colnames(out@output) == "Result"] <- paste(o1@name, "Result", sep = "")
	out@marginal <- colnames(out@output) %in% c(xmarg, ymarg)

	return(out)
}

rm(wiki_username)
objects.store(list = ls())
cat("All objects in the global namespace were stored:", ls(), "\n")