Tuesday, July 9, 2013

Why sometimes use DATA Step instead of PROC SQL

A question
There is an interesting thread about the efficiency of PROC SQL on Mitbbs.com.
Both datasets have1 and have2 have few variables, but many observations. Very simple program, but without any result after running it for a whole afternoon.
proc sql;
create table want as
select *
from have1 as a left have2 as b
on b.serie_beg <= a.serie <= b.serie_end;
quit;
If I replace left join with inner join, the result is the same. Why it takes so long?
Four types of joins in PROC SQL
To understand this question, we may have to go back into PROC SQL. A technical paper introduces that like other more generalized database systems, there are four kinds of joins to combine two data sets (tables) in the SQL procedure of SAS: Cartesian join, Hash join, Merge join and Indexed join. Using the command _method at the proc sql statement can let the log of SAS tell which join the current join is based on. For example, sqxjsl stands for Cartesian join; sqxjhsh is for Hash join; sqxjm is for Merge join; sqxjndx is for Indexed join.
  • Cartesian join
Cartesian join uses two nested loops to join two tables. In other SQL systems, Cartesian join is often referred as Nested loop join. Suppose one table has N elements and the other has M elements, the complexity of this algorithm is O(N*M ), and is a bilateral merge that costs more resources. It can be shown in an example below.
************(0)Separate as 2 datasets***************;
data male female;
set sashelp.class;
if sex = 'F' then output female;
else output male;
run;

************(1) Cartesian join ***************;
proc sql _method;
create table friends1 as
select a.name as name1, b.name as name2, a.age as age1, b.age as age2
from male as a left join female as b on a.age > b.age
;quit;
The _method option outlines the execution plan of the join. Here this Cartesian join can't be replaced by more advanced join, and SAS honestly indicates that it is the only approach available.
NOTE: The execution of this query involves performing one or more Cartesian product joins that can not be optimized.

NOTE: SQL execution methods chosen are:

sqxcrta
sqxjsl
sqxsrc( WORK.FEMALE(alias = B) )
sqxsrc( WORK.MALE(alias = A) )
  • Merge join
Merge join implements a sorting for each of the tables and merges them afterward. The time spent is the summation between the sorting and the scanning, and therefor the complexity is O(Nlog(N) + Mlog(M)).
************(2) Merge join ***************;
proc sql _method;
create table friends2 as
select a.name as name1, b.name as name2, a.age as age1, b.age as age2
from male as a left join female as b on a.age = b.age
;quit;
Clearly this join contains two sqxsort steps before the following sqxjm step. The logic here is quite similar to the traditional DATA Step's sort-sort-merge.
NOTE: SQL execution methods chosen are:

sqxcrta
sqxjm
sqxsort
sqxsrc( WORK.FEMALE(alias = B) )
sqxsort
sqxsrc( WORK.MALE(alias = A) )
  • Hash join
Hash join loads the tables (if they are small enough to fit memory) into memory as hashs table and then do the scanning. The in-memory process is fast and the big O is O(N+M ).
************(3) Hash join ***************;
proc sql _method;
create table friends3 as
select a.name as name1, b.name as name2, a.age as age1, b.age as age2
from male as a inner join female as b on a.age = b.age
;quit;
Usually such an inner join initiates the Hash join.
NOTE: SQL execution methods chosen are:

sqxcrta
sqxjhsh
sqxsrc( WORK.MALE(alias = A) )
sqxsrc( WORK.FEMALE(alias = B) )
  • Indexed join
For relative large data sets that are frequently read, building an index is a good idea. An index is a basically a B-tree structure that translates to O(log(N)). If we have indexes for both lookup table and base table, the complexity will be like Max[O(Mlog(N)), O(Nlog(M))]. In the demo below, I create two indexes for a lookup table and a simulated data set.
************(4) Indexed join ***************;
data simu;
do number = 1 to 1e6;
output;
end;
run;

proc sql _method;
create index number on simu(number)
;
create index age on male(age)
;
create table friends4 as
select a.name as name1, a.age as age1, b.number
from male as a inner join simu as b on a.age = b.number
;quit;
Hereby the Indexed join has been automatically activated.
NOTE: SQL execution methods chosen are:

sqxcrta
sqxjndx
sqxsrc( WORK.MALE(alias = A) )
sqxsrc( WORK.SIMU(alias = B) )
Alternative solutions to the question
The good thing for PROC SQL is that the execution plan manager always chooses the optimized join plan. The bad thing is that the user is not allowed to specify any particular join. As for the question at the beginning, the Cartesian joins are obviously chosen to realize the left joins. If the data sets are quite large, the waiting seems to be painful. However, I can imagine some other ways in SAS to apply faster algorithms and avoid the Cartesian joins if not necessary.
  • DATA Step Merge
The sort-sort-merge forces a merge join, which decreases the complexity from O(N*M ) to O(Nlog(N) + Mlog(M)).
  • PROC FORMAT
If the lookup table is not very large, a user defined format with PROC FORMATis a unique way in SAS to construct the hash table in memory to perform lookup or scanning, which further decreases the complexity from O(N*M ) to O(N+M ).

Thursday, June 27, 2013

The complexity for Fibonacci numbers in SAS

Fibonacci numbers have a lot of applications. There have been a lot of interesting researches regarding Fibonacci numbers since 800 years ago. Mathematically they are the numbers from a sequence, which is defined as F_n = F_{n-1} + F_{n-2}\hspace{5}with\hspace{5}F_0 =0, \hspace{2} F_1=1.

Experiment and Result

Fibonacci numbers can be derived from either the original definition that is a typical example of recursion, or a mathematically simplified form that is called closed-from expression. In SAS, a user-defined Fib1 function in PROC FCMP can realize the recursion expression, and another Fib2 function is created to express the closed-form expression. I insert both functions to a loop from 10 to 40 with step size of 5 and record the real time costed.
According to the figure above, both algorithms have no actually difference while N is relative small. When N becomes 30 or bigger, the recursion expression function spends a lot of system time. When N is even greater than 40, SAS enters a freezing status and I have to break from the command manually. Thus, the curve for the recursion expression seems to fit an exponential relationship fairly well. On the contrary, the algorithm of closed-from expression takes almost constant time for the execution,

Complexity

  • Recursion expression
    T(n) = T(n-1) + T(n-2) + O(1), which may imply O(2^{n-1}) + O(2^{n-2}) + O(1) = O(2^n).
  • Closed-form expression
    The equation is fairly straightforward and described as F_n= {\phi^n-(-\phi)^{-n}\over\sqrt{5}}\hspace{2}with\hspace{2} \phi={(1+\sqrt{5})\over2}.  Here \phi is the well-known golden ratio. Its complexity is a constant, which suggests O(n) = O(n-1) = O(1).

Conclusion

Mathematics sometimes significantly decreases the complexity of the computation, in this case, from O(2^n) to O(1). It is also useful for us to estimate how much resources are needed for specific purpose.

****(1) Bulid user-defined functions *****************;
proc fcmp outlib = work.test.fibonacci;
* Create the 1st function based on recursion;
function fib1(n);
if n<0 then return(0);
else if n = 1 or n = 2 then return(1);
else return(fib1(n-1) + fib1(n-2));
endsub;
* Create the 2nd function based on closed form;
function fib2(n);
phi = (1 + sqrt(5))/2;
return(round((phi**n - (1-phi)**n) / sqrt(5)));
endsub;
quit;

****(2) Run the tests sequentially********************;
options cmplib = work.test fullstimer;
%macro test();
* Create datasets for verification
%do j = 10 %to 40 %by 5;
data _fib1_&j;
do i = 1 to &j;
x = fib1(i);
output;
end;
%put Recursion method at &j;
run;
data _fib2_&j;
do i = 1 to &j;
x = fib2(i);
output;
end;
%put Closed form method at &j;
run;
%end;
%mend;

* Output SAS log to text file;
proc printto log = "c:\tmp\test.log" new;
run;
%test()
proc printto;
run;

****(3) Display the results***************************;
proc datasets nolist;
delete _:;
quit;

data one;
infile "c:\tmp\test.log" truncover;
input @1 line $80.;
if index(line, 'real') > 0;
run;

data two;
set one;
length method $20.;
retain number 5;
* Ignore the time used by PROC PRINTTO;
if _n_ > 1;
if mod(_n_ , 2) = 0 then do;
method = "Recursion";
number + 5;
end;
else method = "Closed form";
time = input(substr(line, 21, 5), 5.2);
run;

proc sgplot data = two;
series x = number y = time / group = method;
yaxis label = "Time elapsed by seconds" grid;
run;

Tuesday, June 25, 2013

Use doParallel for rolling regression in R

Rolling regression for a large data set costs lots of resources. There are a few strategies to speed up this process. For example, in R, there is a rollapply function in the dynlm package. In SAS, PROC FCMP is one of the options for optimization.
Revolution Analytics' doParallel package utilizes the power of the foreach package and realizes parallel processing. It seems that the recent update version 1.0.3 has fixed a bug on Windows system. I am interested in seeing if the doParallel package would allow any efficiency boost for rolling regression, on my old Windows computer which has two cores (Intel Core 2 Duo 3.06GHz).
I first created the random vectors of 20000 for x and y, and set the rolling window size to be 20. I selected 3 scenarios: sequential processing, parallel processing with 2 cores and parallel processing with 4 cores. The parameters solved by the rolling regressions are shown in the picture above. The matrices by the three methods have no difference.
It turned out that the time cost has been significantly improved under the parallel mode. Since I actually have no more than 2 cores on this computer, the registerDoParallel(cores=4) automatically killed the redundant connections and performed the same as the cores=2 mode. Since currently most computers have multiple cores, the doParallel package has a lot of potentials in statistics.
####(1) Import libraries and set up directory#################
library(doParallel)
library(ggplot2)
setwd('C:/Google Drive/r')
getwd()

####(2) Simulate the data set ################################
# Set a few parameters for regressions
n <- 20000
beta0 <- 0.01; beta1 <- 0.4
mse = 0.9
set.seed(1234)

# Generate x and y variables
x <- rnorm(n)/100
y <- beta0 + beta1*x + rnorm(n, 0, mse)/100

####(3) Run the tests #######################################
# Set the rollling window size to be 20
windowSize <- 20

# Run sequentially
time1 <- system.time(
z1 <- foreach(i=seq(1, (n-windowSize+1), 1), .combine=rbind) %do% {
idx <- i:(i+windowSize-1)
coefficients(lm(y[idx]~x[idx]))
}
)

# Run with 2 cores registered
registerDoParallel(cores=2)
time2 <- system.time(
z2 <- foreach(i=seq(1, (n-windowSize+1), 1), .combine=rbind) %dopar% {
idx <- i:(i+windowSize-1)
coefficients(lm(y[idx]~x[idx]))
}
)

# Run with 4 cores registered
registerDoParallel(cores=4)
time3 <- system.time(
z3 <- foreach(i=seq(1, (n-windowSize+1), 1), .combine=rbind) %dopar% {
idx <- i:(i+windowSize-1)
coefficients(lm(y[idx]~x[idx]))
}
)

####(4) Display the results #################################
# Show the parameters from rolling regression
z <- data.frame(cbind(z1, 1:nrow(z1)))
names(z) <- c('beta0', 'beta1', 'exp')
ggplot(data=z, aes(x=exp, y=beta1, alpha=0.5)) + geom_line() + xlab("Number of regressions") +
geom_hline(aes(yintercept=0.4, color = 'red3')) + opts(legend.position="none")
ggplot(data=z, aes(x=exp, y=beta0, alpha=0.5)) + geom_line() + xlab("Number of regressions") +
geom_hline(aes(yintercept=0.01, color = 'red3')) + opts(legend.position="none")

# Show the time elapsed
time <- data.frame(rbind(time1, time2, time3))
row.names(time) <- c('No parallel', '2 cores', '4 cores')
ggplot(data=time, aes(x=row.names(time), y=elapsed)) + geom_bar(aes(fill=row.names(time))) +
ylab("Total time elapsed by seconds") + xlab("")

Thursday, June 13, 2013

The US airports with most flight routes

The summer just begins and a lot of people start to experience flight delay. I am also interested in seeing which cites in the US are well connected by flights. The OpenFlights project provides free flight information. Then I used R to download the data and loaded them into a SQLite database, since I try to keep them as the persistent data. RSQLite facilities R to join and query tables in the database. Finally the flight routes and the airports were visualized by ggplot2 and ggmap.
RankIATA codeCityArriving flight routes
1ORDChicago494
2LAXLos Angeles438
3DENDenver401
4JFKNew York363
5ATLAtlanta339
6DFWDallas-Fort Worth281
7SFOSan Francisco259
8IAHHouston244
9MIAMiami244
10EWRNewark242
Wiki says that Atlanta is the busiest airport in the US according to total passenger boardings. However, from the number of the incoming flight routs, it only ranks 5th, following Chicago, Los Angeles, Denver and New York. Possibly Atlanta is the hub mostly for passengers to do connection. If somebody really loves air traveling, Chicago(with ORD) and New York(with both JFK and EWR) are the two most convenient cities to stay with, because they have the most options.
This post is inspired by one post on the blog Data Science and R
# Import libraries and set up directory
library(ggmap)
library(RSQLite)
setwd("C:/Google Drive/Codes")

# Read data directly from URLs
airport <- read.csv("http://openflights.svn.sourceforge.net/viewvc/openflights/openflights/data/airports.dat", header = F)
route <- read.csv("http://openflights.svn.sourceforge.net/viewvc/openflights/openflights/data/routes.dat", header = F)

# Remove the airports without IATA codes and rename the variables
airport <- airport[airport$V5!='', c('V3', 'V4', 'V5','V7','V8','V9')]
colnames(airport) <- c("City", "Country", "IATA", "lantitude", "longitude", "altitude")
route <- route[c('V3', 'V5')]
colnames(route) <- c("Departure", "Arrival")

# Store data to SQLite database
conn <- dbConnect("SQLite", dbname = "air.db")
dbSendQuery(conn, "drop table if exists airport;")
dbWriteTable(conn, "airport", airport)
dbSendQuery(conn, "drop table if exists route;")
dbWriteTable(conn, "route", route)
dbDisconnect(conn)

# Manipulate data in SQLite database
conn <- dbConnect("SQLite", dbname = "air.db")
sqlcmd01 <- dbSendQuery(conn, "
select a.type, a.city as iata, a.frequency, b.city, b.country, b.lantitude, b.longitude
from (select 'depart' as type, departure as city, count(departure) as frequency
from route
group by departure
order by frequency desc, type) as a
left join airport as b on a.city = b.iata
order by frequency desc
;")
top <- fetch(sqlcmd01, n = -
1)
sqlcmd02 <- dbSendQuery(conn, "
select route.rowid as id, route.departure as point, airport.lantitude as lantitude, airport.longitude as longitude
from route left join airport on route.departure = airport.iata
union
select route.rowid as id, route.arrival as point, airport.lantitude as lantitude, airport.longitude as longitude
from route left join airport on route.arrival = airport.iata
order by id
;")
combine <- fetch(sqlcmd02, n = -
1)
dbDisconnect(conn)

# Draw the flight routes and the airports on Google map
ggmap(get_googlemap(center = 'us', zoom =
4, maptype = 'roadmap'), extent = 'device') +
geom_line(data = combine, aes(x = longitude, y = lantitude, group = id), size =
0.1,
alpha =
0.05,color = 'red4') +
geom_point(data = top, aes(x = longitude, y = lantitude, size = frequency), colour = "blue", alpha =
0.3) +
scale_size(range=c(
0,15))