Searching for chemical substances based on ATC codes from WHO search engine.

The concept of how to grep search engine results with R had been one of my lingering thoughts. Because sometimes it's just feels easier if I could just grep search engine results than to find a dataset or complete database of a certain thing.

So when at work I had to use ATC codes to get the substances out or the pharmaceutical(?) use. So while no-one was looking I naturally had to devise the following script. This snippet is more or less done, a part of a imaginary whole if you will. It lacks the functionality to search for partial matches. So you can't just look for all the "B01" codes. Probably doing that extension next when I have the time. Also you can't use this for substance to ATC code. But for now this one way ticket, the basic concept as a snippet. Naturally comments about making it more readable and/or efficient are always welcome.

	require(gdata)
	require("RCurl")
	require("XML")
	AtcSearch <- function(AtcCodes){
	#Function inside a function is bad practice...
	#DummySearch searches for one atc code
	DummySearch <- function(Code){
	#A link to WHO search
	htmlLink <- paste("http://www.whocc.no/atc_ddd_index/?code=", Code, sep="")
	doc <- htmlParse(getURL(htmlLink))
	Drug <- xpathSApply(doc, "//td", xmlValue)
	#No idea why these three are needed
	Drug <- gsub("Â", "", Drug)
	Drug <- gsub("\\s", " ", Drug)
	Drug <- trim(Drug) #gdata just for this...
	#if no drug is found, returning NULL
	if(sum(Drug == Code) == 0 | length(Drug) == 0) return()
	atcLevels <- xpathSApply(doc, "//b", xmlValue)
	out <- matrix(Drug[which(Drug == Code):length(Drug)], ncol= which(Drug == Code)-1, byrow=TRUE)
	#Basic handling of multiple DDDs and Units
	if(nrow(out) > 1) out <- t(apply(out, 2, function(x) paste(x[x != ""], collapse="/")))
	out <- data.frame(out, t(atcLevels))
	colnames(out) <- c(Drug[1:6], "Anatomical", "Therapeutic", "Subgroup", "ChemicalSubgroup")
	out
	}
	#lapply to get all the codes.
	ret <- lapply(AtcCodes, DummySearch)
	do.call("rbind", ret)
	}
	#Some examples
	AtcSearch("N01AX15")
	AtcSearch(c("B01AA01", "N01AX15"))
	#Not perfect for example, the who database isn't complete.
	AtcSearch("QN01AX99")

view raw AtcSearch.R hosted with ❤ by GitHub

A Merging Test Bench

As requested here's the packed data and a test bench you can test your own merging function ideas and replicate my results (hopefully). If you want the plots you can use the end part of scripts in part1 part2.

The data is a bunch of super secret Eve Online killmails. The first part of the script handles the json to data.frame transformation. After that I introduce the R core merging in a for-loop, data.table variant, reduce function variants, and finally the fast.merging. There's actually alot more data then we are using in the bench. Totalling 166106 killmails, merging only 3000 of those. That is because as noted in previous posts the merge for-looping seems to work in exponential time so going much higher then 3000 lines seems really excessive. Even now it took over an hour to run the whole bench on intel i5 3570k (not overclocked). So if you want faster results choose 1000ish frames and set rbenchmark replications to 1. With low frame count (under 100) the results are fairly similiar.

If you are going to test are the results from different functions the same, note that the order of rows and columns will be different. I'm eagerly waiting for you results and ideas.

	#Using 64bit R 2.15.3
	#Packages we are going to use
	library(rjson)
	library(data.table)
	library(rbenchmark)
	#set your WD where your unpackaged file is, or add path.
	#Reading only partial data for this test, you can naturally do as you wish
	eve = readLines("kms-2013-01-06-s.json", n=3001)
	#json => data.frames (data lines)
	eve = matrix(eve, ncol=1)
	eve.list = apply(eve, 1, list)
	#First and last lines are useless
	eve.list = eve.list[-1]
	#eve.list = eve.list[-length(eve.list)] #didn't read the last line in this case
	#Handling json...
	eve.list = lapply(eve.list, function(x){
	x = unlist(x)
	return(fromJSON(x))
	})
	#Handling data.frame
	eve.data.frame = function(eve.list){
	framo = unlist(eve.list)
	dups = unique(names(framo)[duplicated(names(framo))])
	if(length(dups) > 0){
	for(i in seq(dups)){
	bol = names(framo) %in% dups[i]
	framo[unique(dups[i])] = paste(framo[bol], collapse= ",")
	}
	framo = framo[!duplicated(names(framo))]
	}
	return(as.matrix(t(framo)))
	}
	#This takes awhile if you read the whole data
	eve.frames = lapply(eve.list, eve.data.frame)
	#For example
	eve.frames[[1]]
	eve.frames[[100]]
	#The amount of variables in each data line
	N.variables = sapply(eve.frames, ncol)
	table(N.variables)
	#Converting data.frames to data.tables
	eve.data.table = lapply(eve.frames, as.data.table)
	#Is this correct?
	#The functions used
	core.forLoop <- function(data.list) {
	complete = data.list[[1]]
	for(i in 2:length(data.list)){
	part = data.list[[i]]
	complete = merge(complete, part,
	all=TRUE,
	sort=FALSE)
	}
	return(complete)
	}
	data.table.forLoop <- function(data.list){
	complete = data.list[[1]]
	for(i in 2:length(data.list)){
	part = data.list[[i]]
	complete = merge(complete, part,
	by=intersect( names(complete), names(part) ),
	all=TRUE,
	sort=FALSE)
	}
	return(complete)
	}
	core.Reduce.merge = function(data.list){
	Reduce(function(a, b) merge(a, b, all=TRUE, sort=FALSE), data.list)
	}
	data.table.Reduce.merge = function(data.list){
	Reduce(function(a, b) merge(a, b, by=intersect(names(a), names(b)), all=TRUE, sort=FALSE), data.list)
	}
	fast.merging = function(data.list, nparts){
	if(!is.list(data.list)) stop("data.list isn't a list")
	while(length(data.list) != 1){ #Loop until everything is merged
	if(length(data.list) > nparts){
	starts = seq(1, length(data.list), nparts)
	ends = seq(nparts, length(data.list), nparts) #starts and ends are of equal size if length(data.list) divides nparts.
	if(length(ends) < length(starts)) ends = c(ends, length(data.list)) #making sure things are even
	sections = matrix(c(starts, ends), ncol=2, byrow=FALSE)
	sections = apply(sections, 1, list)
	}else{
	sections = list(c(1, length(data.list)))
	}
	#We have the standard way inside lapply
	data.list = lapply(sections, function(x, data.list){
	if(is.list(x)) x = x[[1]]
	#the standard way starts ->
	part = data.list[[x[1]]]
	for(i in x[1]:x[2]){
	part = merge(part, data.list[[i]], all=TRUE, sort=FALSE)
	}
	#<- standard way ends
	return(part)
	}, data.list = data.list)
	}
	return(data.list[[1]]) #returning the merged data frame
	}
	#This will take awhile, go to lunch or something like that.
	benchmark(core.forLoop(eve.frames),
	data.table.forLoop(eve.data.table),
	core.Reduce.merge(eve.frames),
	data.table.Reduce.merge(eve.data.table),
	fast.merging(eve.frames, 10),
	replications=5)
	# test replications elapsed relative user.self sys.self user.child sys.child
	#1 core.forLoop(eve.frames) 5 1398.29 4.932 1393.98 0.06 NA NA
	#3 core.Reduce.merge(eve.frames) 5 1405.77 4.958 1401.02 0.31 NA NA
	#2 data.table.forLoop(eve.data.table) 5 797.76 2.814 794.61 0.09 NA NA
	#4 data.table.Reduce.merge(eve.data.table) 5 801.93 2.829 799.91 0.16 NA NA
	#5 fast.merging(eve.frames, 10) 5 283.51 1.000 282.95 0.00 NA NA
	#If you are checking are going to check are the results the same, remember to sort your rows and columns
	#basically data[order(data[, 1]), order(colnames(data))] for all data.

view raw a_merging_bench.R hosted with ❤ by GitHub

RMark: data.table merge vs core merge

This is the third post concerning fast merging in R, first here and second here. This time we are going to look at how the merge function from data.table package works in our case, requested by Uwe Block.

As a reminder the first post concerns doing a merging scheme using the lapply function. After that in the second post we looked how it translates to a parallel solution for even more speed.

As a reminder the problem is the following: given N amount of files with a random amount of variables in each file. How do we merge these files into a complete data within reasonable amount of time? For example you can consider that the medical instruments records it's each run (patient) to a different file. With only the values measured by the instrument are recorded to the file.

In practice I'm using Eve Online killmails. So far I have around 160k mails. Hence the need for a fast solution.

Personally I haven't really used data.table, but it looks good enough. So I recommend checking it out. But in a sense nothing's really changed on my code part. I did a fast merging function based on data.table merge called the lady, but won't give the code because it isn't fast enough (spoiler alert). Also the data.table merge is a generic function. Meaning if you are using data.tables in a merge it automaticaly uses the data.table merge function. But enough it this prattle. Here's the results.



data.table merge is faster then the R core merge function if you are merging a for-loop solution. The problem still is that the data.table merging in a for-loop works in exponential time. Which doesn't translate well for big data. The single core fast merging scheme still beats the data.table merge in a for-loop. The biggest wonder is why when we combine the data.table merge with parallel fast merging scheme we get a slower outcome. Naturally I though we could squeeze few seconds out with low frame count, but seems that it's actually few minutes worse. So far I'm running dry in practical ideas and reasons why the lady seems so fat. Maybe it doesn't cleane up it's diet ? (memory management?) If you have some insight please comment.

As requested, here's the bench.

	#Using 64bit R 2.15.3
	#Packages we are going to use
	#download the packed data here: https://www.dropbox.com/s/c1rtcryp8p6500w/kms-2013-01-06-s.json.bz2
	library(rjson)
	library(data.table)
	library(snow)
	library(rbenchmark)
	#set your WD where your unpackaged file is, or add path.
	#Reading only partial data for this test, you can naturally do as you wish
	eve = readLines("kms-2013-01-06-s.json", n=80000)
	#json => data.frames (data lines)
	eve = matrix(eve, ncol=1)
	eve.list = apply(eve, 1, list)
	#First and last lines are useless
	eve.list = eve.list[-1]
	#eve.list = eve.list[-length(eve.list)] #didn't read the last line in this case
	#Handling json...
	eve.list = lapply(eve.list, function(x){
	x = unlist(x)
	return(fromJSON(x))
	})
	#Handling data.frame
	eve.data.frame = function(eve.list){
	framo = unlist(eve.list)
	dups = unique(names(framo)[duplicated(names(framo))])
	if(length(dups) > 0){
	for(i in seq(dups)){
	bol = names(framo) %in% dups[i]
	framo[unique(dups[i])] = paste(framo[bol], collapse= ",")
	}
	framo = framo[!duplicated(names(framo))]
	}
	return(as.matrix(t(framo)))
	}
	#This takes awhile if you read the whole data
	eve.frames = lapply(eve.list, eve.data.frame)[1:60000]
	#data.frames to data.tables
	eve.frames2 = lapply(eve.frames, as.data.table)
	#unique(sapply(eve.frames2, class)) #data.frame data.tables
	##
	#Functions
	##
	#Standard FM
	fast.merging = function(data.list, nparts){
	if(!is.list(data.list)) stop("data.list isn't a list")
	while(length(data.list) != 1){ #Loop until everything is merged
	if(length(data.list) > nparts){
	starts = seq(1, length(data.list), nparts)
	ends = seq(nparts, length(data.list), nparts) #starts and ends are of equal size if length(data.list) divides nparts.
	if(length(ends) < length(starts)) ends = c(ends, length(data.list)) #making sure things are even
	sections = matrix(c(starts, ends), ncol=2, byrow=FALSE)
	sections = apply(sections, 1, list)
	}else{
	sections = list(c(1, length(data.list)))
	}
	#We have the standard way inside lapply
	data.list = lapply(sections, function(x, data.list){
	if(is.list(x)) x = x[[1]]
	#the standard way starts ->
	part = data.list[[x[1]]]
	for(i in x[1]:x[2]){
	part = merge(part, data.list[[i]], all=TRUE, sort=FALSE)
	}
	#<- standard way ends
	return(part)
	}, data.list = data.list)
	}
	return(data.list[[1]]) #returning the merged data frame
	}
	#MC FM
	mc.fast.merging = function(data.list, nparts, cluster){
	if(!is.list(data.list)) stop("data.list isn't a list")
	while(length(data.list) != 1){ #Loop until everything is merged
	if(length(data.list) > nparts){
	starts = seq(1, length(data.list), nparts)
	ends = seq(nparts, length(data.list), nparts) #starts and ends are of equal size if length(data.list) divides nparts.
	if(length(ends) < length(starts)) ends = c(ends, length(data.list)) #making sure things are even
	sections = matrix(c(starts, ends), ncol=2, byrow=FALSE)
	sections = apply(sections, 1, list)
	}else{
	sections = list(c(1, length(data.list)))
	}
	if(length(sections) !=1){
	data.list = parLapply(cluster, sections, function(x, data.list){
	if(is.list(x)) x = x[[1]]
	#the standard way starts ->
	part = data.list[[x[1]]]
	for(i in x[1]:x[2]){
	part = merge(part, data.list[[i]], all=TRUE, sort=FALSE)
	}
	#<- standard way ends
	return(part)
	}, data.list = data.list)
	}else{
	data.list = lapply(sections, function(x, data.list){
	if(is.list(x)) x = x[[1]]
	part = data.list[[x[1]]]
	for(i in x[1]:x[2]){
	part = merge(part, data.list[[i]], all=TRUE, sort=FALSE)
	}
	return(part)
	}, data.list = data.list)
	}
	}
	return(data.list[[1]]) #returning the merged data frame
	}
	#Here's the lady
	LadyMergAlot = function(data.list, nparts, cluster){
	if(!is.list(data.list)) stop("data.list isn't a list")
	while(length(data.list) != 1){ #Loop until everything is merged
	if(length(data.list) > nparts){
	starts = seq(1, length(data.list), nparts)
	ends = seq(nparts, length(data.list), nparts) #starts and ends are of equal size if length(data.list) divides nparts.
	if(length(ends) < length(starts)) ends = c(ends, length(data.list)) #making sure things are even
	sections = matrix(c(starts, ends), ncol=2, byrow=FALSE)
	sections = apply(sections, 1, list)
	}else{
	sections = list(c(1, length(data.list)))
	}
	if(length(sections) !=1){
	data.list = parLapply(cluster, sections, function(x, data.list){
	if(is.list(x)) x = x[[1]]
	#the standard way starts ->
	part = data.list[[x[1]]]
	for(i in x[1]:x[2]){
	part = merge(part, data.list[[i]],
	by=intersect(names(part), names(data.list[[i]])),
	all=TRUE, sort=FALSE)
	}
	#<- standard way ends
	return(part)
	}, data.list = data.list)
	}else{
	data.list = lapply(sections, function(x, data.list){
	if(is.list(x)) x = x[[1]]
	part = data.list[[x[1]]]
	for(i in x[1]:x[2]){
	part = merge(part, data.list[[i]],
	by=intersect(names(part), names(data.list[[i]])),
	all=TRUE, sort=FALSE)
	}
	return(part)
	}, data.list = data.list)
	}
	}
	return(data.list[[1]]) #returning the merged data frame
	}
	##
	#Here starts the bench.
	##
	sekvenssi = seq(2, 5002, 100)
	sekvenssi2 = seq(2, 50002, 1000)
	#Without data.table
	time1 = rep(NA, length(sekvenssi))
	for(i in 1:length(sekvenssi)){
	complete = eve.frames[[1]]
	time1[i] = system.time(
	for(j in 2:sekvenssi[i]){
	complete = merge(complete, eve.frames[[j]], all=TRUE, sort=FALSE)
	}
	)[3]
	print(paste(round(time1[i],2), nrow(complete),i))
	}
	write(time1, "time1.dat")
	#With data.table
	time2 = rep(NA, length(sekvenssi))
	for(i in 1:length(sekvenssi)){
	complete = eve.frames2[[1]]
	time2[i] = system.time(
	for(j in 2:sekvenssi[i]){
	complete = merge(complete, eve.frames2[[j]], by=intersect(names(complete), names(eve.frames2[[j]])),all=TRUE, sort=FALSE)
	}
	)[3]
	print(paste(round(time2[i],2), nrow(complete),i))
	}
	write(time2, "time2.dat")
	#FM with data.frames
	time3 = rep(NA, length(sekvenssi2))
	for(i in 1:length(sekvenssi2)){
	complete = eve.frames[1:sekvenssi2[i]]
	time3[i] = system.time(fast.merging(complete, 10))[3]
	print(c(i, time3[i]))
	if(time3[i] > 400) stop()
	}
	write(time3, "time3.dat")
	#MCFM with data.frames
	cl = makeCluster(4, type="SOCK")
	time4 = rep(NA, length(sekvenssi2))
	for(i in seq(sekvenssi2)){
	complete = eve.frames[1:sekvenssi2[i]]
	time4[i] = system.time(mc.fast.merging(complete, 10, cl))[3]
	print(c(i, time4[i]))
	if(time4[i] > 400) stop()
	}
	stopCluster(cl)
	write(time4, "time4.dat")
	#MCFM with data.tables
	cl = makeCluster(4, type="SOCK")
	time5 = rep(NA, length(sekvenssi2))
	for(i in seq(sekvenssi2)){
	complete = eve.frames2[1:sekvenssi2[i]]
	time5[i] = system.time(LadyMergAlot(complete, 10, cl))[3]
	print(c(i, time5[i]))
	if(time5[i] > 400) stop()
	}
	stopCluster(cl)
	write(time5, "time5.dat")
	time1 = scan("time1.dat")
	time2 = scan("time2.dat")
	time3 = scan("time3.dat")
	time4 = scan("time4.dat")
	time5 = scan("time5.dat")
	png("mark3.png", width=960)
	plot(sekvenssi, time1/60, type="l", xlim=c(0,50000), ylim=c(0,5),
	ylab="Elapsed time (minutes)", xlab="Frames merged")
	points(sekvenssi, time2/60, type="l", col="blue")
	points(sekvenssi2, time3/60, type="l", col="red")
	points(sekvenssi2, time4/60, type="l", col="pink")
	points(sekvenssi2, time5/60, type="l", col="green")
	legend("bottomright", c("R standard", "data.table", "FM 1 cores", "FM 4 cores","The Lady"),
	lty=c(1,1,1,1), col=c("black", "blue", "red", "pink", "green"))
	dev.off()

view raw FM_bench3.R hosted with ❤ by GitHub

RMark: Parallel fast merging

Last time we looked at a different method of merging multiple files via fast mergin g. I highly recommend checking the fast merging function before we will translate it to a parallel solution. Parallel programming in this case means we are using more then one processor core. So even if you have a fancy 8 core machine R only uses one core if you don't use parallel functions.
In R parallel programming can be done for example with snow package. When we change our function from 1 core solution to a multicore solution we use snow packages parLapply function instead of lapply function. This change will magically make our independent loops run faster without doing any major changes to the script. As parLapply behaves almost exactly the same as lapply. So let us present the script first, the explanation will be given after that. There's alot more case handling in multicore case as we can't call our multicore solution if we only need one core.

	#Multicore (mc) requires snow package as multicore isn't on windows...
	require(snow)
	mc.fast.merging = function(data.list, nparts, cluster){
	if(!is.list(data.list)) stop("data.list isn't a list")
	while(length(data.list) != 1){ #Loop until everything is merged
	if(length(data.list) > nparts){
	starts = seq(1, length(data.list), nparts)
	ends = seq(nparts, length(data.list), nparts) #starts and ends are of equal size if length(data.list) divides nparts.
	if(length(ends) < length(starts)) ends = c(ends, length(data.list)) #making sure things are even
	sections = matrix(c(starts, ends), ncol=2, byrow=FALSE)
	sections = apply(sections, 1, list)
	}else{
	sections = list(c(1, length(data.list)))
	}
	if(length(sections) !=1){
	data.list = parLapply(cluster, sections, function(x, data.list){
	if(is.list(x)) x = x[[1]]
	#the standard way starts ->
	part = data.list[[x[1]]]
	for(i in x[1]:x[2]){
	part = merge(part, data.list[[i]], all=TRUE, sort=FALSE)
	}
	#<- standard way ends
	return(part)
	}, data.list = data.list)
	}else{
	data.list = lapply(sections, function(x, data.list){
	if(is.list(x)) x = x[[1]]
	part = data.list[[x[1]]]
	for(i in x[1]:x[2]){
	part = merge(part, data.list[[i]], all=TRUE, sort=FALSE)
	}
	return(part)
	}, data.list = data.list)
	}
	}
	return(data.list[[1]]) #returning the merged data frame
	}

view raw mc_fast_merging.R hosted with ❤ by GitHub

For example if we have 100 frames to merge and a 4 core machine. We choose a section size of 25 parts (nparts=25). Which means that we have 4 sections. Section 1 will be frames from 1-25, section 2 will be frames 26-50, and so on. Now we tell the first core to merge the frames from section 1, the second core will merge frames outlined in section 2, and so on. After that is done we only have 4 frames, which can we will merge using the standard way. With this fast merging should atleast two times faster when using 2 cores. Let us see if that is the case, unfortunately I only have 4 cores.

	#The frames.
	eve.frames = lapply(eve.list, eve.data.frame)[1:40000]
	#eve.frames = eve.frames[1:3001]
	sekvenssi = seq(2, 40002, 500)[1:30]
	#sekvenssi = seq(1, 4001, 100)
	time1 = rep(NA, length(sekvenssi))
	for(i in seq(sekvenssi)){
	complete = eve.frames[1:sekvenssi[i]]
	time1[i] = system.time(fast.merging(complete, 25))[3]
	print(paste(i, round(time1[i], 2)))
	}
	write(time1, "time1.dat")
	library(snow)
	cl = makeCluster(2, type="SOCK")
	time2 = rep(NA, length(sekvenssi))
	for(i in seq(sekvenssi)){
	complete = eve.frames[1:sekvenssi[i]]
	time2[i] = system.time(mc.fast.merging(complete, 25, cl))[3]
	print(i)
	}
	stopCluster(cl)
	write(time2, "time2.dat")
	cl = makeCluster(3, type="SOCK")
	time3 = rep(NA, length(sekvenssi))
	for(i in seq(sekvenssi)){
	complete = eve.frames[1:sekvenssi[i]]
	time3[i] = system.time(mc.fast.merging(complete, 25, cl))[3]
	print(i)
	}
	stopCluster(cl)
	write(time3, "time3.dat")
	cl = makeCluster(4, type="SOCK")
	time4 = rep(NA, length(sekvenssi))
	for(i in seq(sekvenssi)){
	complete = eve.frames[1:sekvenssi[i]]
	time4[i] = system.time(mc.fast.merging(complete, 25, cl))[3]
	print(i)
	}
	stopCluster(cl)
	write(time4, "time4.dat")
	time1 = scan("time1.dat")
	time2 = scan("time2.dat")
	time3 = scan("time3.dat")
	time4 = scan("time4.dat")
	png("mark2.png")
	plot(sekvenssi, time1/60, type="l", ylab="Elapsed time (minutes)", xlab="Frames merged")
	points(sekvenssi, time2/60, type="l", col="blue")
	points(sekvenssi, time3/60, type="l", col="red")
	points(sekvenssi, time4/60, type="l", col="pink")
	legend("bottomright", c("FM 1core", "FM 2 cores", "FM 3 cores","FM 4 cores"),
	lty=c(1,1,1,1), col=c("black", "blue", "red", "pink"))
	dev.off()

view raw fm_bench_2.R hosted with ❤ by GitHub

Damn syntax highlighting

Well... It seems my syntax highlighting broke... So Github A-HOY! But don't have the time today... Curses...

R Benchmark: Merge vs fast.merging

Last summer as a trainee I ran into a bit of trouble with the merge function in R. The main problem is that it's slower then your avarage tortoise in a tard field. Also I could not use rbind because of the different lengths of the frames. In code I mean that the next scheme is really slow when you have alot of files:

complete = read.csv(yourfiles[1])

for(i in 2:length(yourfiles)){
     part = read.csv(yourfiles[i])
     complete = merge(complete, part, all=TRUE, sort=FALSE)
}

It wasn't really a problem back in the summer because I didn't have to think it that much. But now it's realy a pain as I need to merge over 160 000 lines. So because I needed to save some time I came up with the following merging algorithm:

#data.list = list of your frames you want to merge
#nparts = how many parts are merged inside apply, the optimal value is probably 
#different for every machine, but should not be a big deal.

fast.merging = function(data.list, nparts=10){

if(!is.list(data.list)) stop("data.list isn't a list") #standard error

while(length(data.list) != 1){ #Loop until everything is merged
       #define sections according to nparts
       if(length(data.list) > nparts){
            starts = seq(1, length(data.list), nparts)
            ends = seq(nparts, length(data.list), nparts) 
            #starts and ends are of equal size 
            #if length(data.list) divides nparts.
       if(length(ends) < length(starts)) ends = c(ends, length(data.list)) 
                #making sure things are even
            sections = matrix(c(starts, ends), ncol=2, byrow=FALSE)
            sections = apply(sections, 1, list)
       }else{
            sections = list(c(1, length(data.list)))
       }
 
 #We have the standard way inside lapply
 #lapply merges a length(sections) amount of nparts sections 
 
 data.list = lapply(sections, function(x, data.list){
  if(is.list(x)) x = x[[1]]
  #the standard way starts ->
  part = data.list[[x[1]]]
  for(i in x[1]:x[2]){ 
   part = merge(part, data.list[[i]], all=TRUE, sort=FALSE) 
  } 
  #<- standard way ends
  return(part)
 }, data.list = data.list)
 #rinse and repeat until all is done
 }
 return(data.list[[1]]) #returning the merged data frame
}

The scheme seems completely unintuitive at first. Because you are basically using lapply to vectorize for-loop merging. But it probably makes sense when you have a grasp how vectorizing works in R, or why you should always define the size of your variables before you use them. But enough of that, time for the actual benchmark results:



As we can see the standard for-loop merging grows exponentially in time as we merge more frames. The algorithm I used grows in linear time and didn't really depend on the nparts argument. Meaning that with nparts being 50, 25 or 10 the time difference is few seconds when we are at 4001 merges. I didn't want to go further because that for-looping get's really tedious to wait. I think we can do better then fast merging, for one we can use parallel calculations. But that's for next week (read: hopefully this weekend).

Benchmark was created with this scheme:

#Getting our lines to a list of lines
#my frames, in this case only line per data-frame.
eve.frames = lapply(eve.list, eve.data.frame) #Some Eve Online kill mails
sekvenssi = seq(1, 4001, 100) #only trying 1, 101, 201, ... 4001.


time1 = rep(NA, length(sekvenssi))
for(i in 1:length(sekvenssi)){
 complete = eve.frames[[1]]
 time1[i] = system.time(
            for(j in 2:sekvenssi[i]){
             complete = merge(complete, eve.frames[[j]], 
                                     all=TRUE, sort=FALSE)
  }
  )[3]
 print(paste(round(time1[i],2), nrow(complete),i))
}
write(time1, "time1.dat")

#backup = complete

time2 = rep(NA, length(sekvenssi))
for(i in seq(sekvenssi)){
 complete = eve.frames[1:sekvenssi[i]]
 time2[i] = system.time(fast.merging(complete, 50))[3]
 print(paste(round(time2[i],2), nrow(complete),i))
}
write(time2, "time2.dat")

time3 = rep(NA, length(sekvenssi))
for(i in seq(sekvenssi)){
 complete = eve.frames[1:sekvenssi[i]]
 time3[i] = system.time(fast.merging(complete, 25))[3]
 print(paste(round(time3[i],2), nrow(complete),i))
}
write(time3, "time3.dat")

time4 = rep(NA, length(sekvenssi))
for(i in seq(sekvenssi)){
 complete = eve.frames[1:sekvenssi[i]]
 time4[i] = system.time(fast.merging(complete, 10))[3]
 print(paste(round(time4[i],2), nrow(complete),i))
}
write(time4, "time4.dat")

#To check if they are correct try
#complete = fast.merging(complete, 10)
#backup = backup[order(backup[, 1]), order(colnames(backup))]
#complete = complete[order(complete[, 1]), order(colnames(complete))]
#sum(backup == complete, na.rm=TRUE)
#sum(is.na(backup)) 
#these should add up to the matrix size.

plot(sekvenssi, time1/60, type="l", ylab="Elapsed time (minutes)", xlab="Lines merged", main="Standard merging VS fast merging")
points(sekvenssi, time2/60, type="l", col="blue")
points(sekvenssi, time3/60, type="l", col="red")
points(sekvenssi, time4/60, type="l", col="pink")
legend("right", c("standard", "FM 50", "FM 25", "FM 10"),
        lty=c(1,1,1,1), col=c("black", "blue", "red", "pink")) 

Comment: Search and Replace

A post in R bloggers caught my attention this morning. The main idea was that how can you change objects in a string. For example given a basket of fruits we would like to change apples to bananas by using R and the ifelse funtion. There are two main solutions how to change one object into another:

#Given a basket of fruits
basket = c("apple", "banana", "lemon", "orange", "orange", "pear", "cherry")

basket = ifelse(basket == "apple", "banana" , basket)
basket[basket == "lemon"] = "pear"
#The latter is usually preferred, but sometimes the ifelse cannot be avoided

If you want to change multiple objects you would need an equal amount of ifelses. In the original post, here, R de jeu presented his function solution which removes the cumbersome writing. I recommend reading the orginal post if the idea isn't immediately clear.

The function given in the post by R de jeu is the following

decode <- function(x, search, replace, default = NULL) {

# build a nested ifelse function by recursion

 decode.fun <- function(search, replace, default = NULL)
 if (length(search) == 0L) {
  function(x) if (is.null(default)) x else rep(default, length(x))
 } else {
         function(x) ifelse(x == search[1L], replace[1L],
                            decode.fun(tail(search, -1L),
                                tail(replace, -1L),
                                default)(x))
}
      return(decode.fun(search, replace, default)(x))
}

I'm not a big fan of recursive functions as they tend to be the slowest of all possible solutions. So I decided to do my own version of the decode function without recursion.

#Same input as in the decode function
decode2 = function(x, search, replace, default=NULL){
 
 #this is if we want the "default type" output
 if(!is.null(default)){ 
    default.bol = !(x %in% search)
    x[default.bol] = default 
 }

 for(i in seq(search)){
  x[x == search[i]] = replace[i] 
 }

 return(x)
}

A standard for-looping with extra logical operations. Now it's time for some speed tests! Spoilers: my solution was around 20ish times faster, but with small vectors the time difference isn't that important.

#If it's worth doing, it's worth overdoing

picnic.basket <- c("apple", "banana",  "Lightsaber", "pineapple", "strawberry",
"lemon","joy","orange", "genelec speakers", "cat", "pear", "cherry",
"beer", "sonic screwdriver", "cat with a caption", "knife", "cheddar",
"The Book of Proofs", "evil ladies", "cheese", "The R inferno", "smoked reindeer meat",
"salmon", "Tardis", "Book of Death", "Pillows", "Blanket", "Woman")

multi.dimensional.basket = sample(picnic.basket, 1000000, replace=TRUE)
search = sample(picnic.basket, 10)
replace = sample(picnic.basket[!(picnic.basket %in% search)], 10)

#Making sure results match.
sum(decode(multi.dimensional.basket, search, replace) == 
    decode2(multi.dimensional.basket, search, replace))
sum(decode(multi.dimensional.basket, search, replace, "fig") == 
    decode2(multi.dimensional.basket, search, replace, "fig"))

system.time( decode(multi.dimensional.basket, search, replace) )
system.time( decode2(multi.dimensional.basket, search, replace) )

system.time( decode(multi.dimensional.basket, search, replace, "fig") )
system.time( decode2(multi.dimensional.basket, search, replace, "fig") )