x <- 0:6
print(x)
print(as.logical(x))
print(as.complex(x))

"""
Explanation: R Programming
<a name="top"></a>
Vahid Mirjalili, Data Scientist

R Data Types
Atomic classes
Object attributes
Create Vectors
Explicit Coercion
Matrices
Lists
Factors
Missing Values
Data Frmaes

Names


Subsetting R objects

Subsetting lists
Partial matching

Removing missing values


Reading and Writing Data

Reading data
Writing data
Reading large datasets
Reading and writing in textual format

Connection interfaces


Flow Control


Functions


Special argument "..."


Scopting Rules


Dates and Times


Format Strings to dates


Apply Functions

Shortcut apply functions

tapply() and split()


Random Number Generation


Gaussian Noise


R Profiler


Using Rprof()


Optimization in R


<a name="sec1_datatypes"></a>
1. R Data Types
top
<a name="atomic_classes"></a>
Atomic classes in R:

charcater
numeric
integer (L)
complex
logical (TRUE/FALSE)

<b> Special numbers </b> Inf and NaN: 1/0 = Inf  and 1/Inf = 0;  0/0=NaN
<a name="attribs"></a>
Object attributes:

names, dimnames
dimensions (for matrices and arrays)
class
length
Other user defined attributes

<b> attributes() </b> to get list of attributes of an object
<a name="cr_vecs"></a>
Creating vectors

vector() :  to create an empty vector
vector("numeric", length=10)
c('a', 'b')  : to concatenate multiple objects into a single vector
Mixing objects from different classes: c(1.5, 'a') : coercion

<a name="coercion"></a>
Explicit coercion

as.character()
as.numeric()
as.logical()
as.complex()
End of explanation
"""


m <- matrix(nrow=2, ncol=3)
print(m)
attributes(m)

# Method1: create a matrix by matrix()
m1 <- matrix(1:6, nrow=2, ncol=3)
print(m1)
# Method2: change the dimension of a vector
m2 <- 1:6
dim(m2) = c(2,3)
print(m2)

x <- 1:4
y <- 6:9
cbind(x,y)
rbind(x,y)

"""
Explanation: <a name="matrices"></a>
Matrices
To create a vector:

create a matrox by matrix() function
change the dimension of a vector dim()
binding vectors and matruces: cbind() and rbind()
End of explanation
"""


x <- list(23, "abc", 1+5i, TRUE)
print(x)

"""
Explanation: <a name="lists"></a>
Lists
Lists are special vectors that can store elements from different classes

Create a list by list()
End of explanation
"""


## create a vector of factors
xf <- factor(c("male", "female", "female", "female", "male", "female", "female", "male", "female"))
str(xf)
## print levels of xf
levels(xf)
unclass(xf)

table(xf)

"""
Explanation: <a name="factors"></a>
Factors
Factors are a special data type to represent categorical variables, such as gender ('male', 'female'), patterm classes ('1', '2') ..

Create factors by factor()
Access the labels: levels()
Make a table view of number of items in each factor: table()

Factor levels are ordered by default in alphabetic order. In order to directly specify the order: factor(vec, levels=c("male", "female"))
End of explanation
"""


x <- c(1,2,3, NA, 6, NaN, 8)
is.na(x)
is.nan(x)

"""
Explanation: <a name="missing_vals"></a>
Missing values
Missing values are denoted by NA or NaN
NaN is for undefined mathematical operations such as 0/0. NaN is a sub-type of NA

is.na() and is.nan() to test whether vector elements are missing
End of explanation
"""


df <- data.frame(foo=11:14, bar=c(T, T, F, T))
df
nrow(df)
ncol(df)

"""
Explanation: <a name="dataframes"></a>
Data Frames
To store tabular data. One of the most useful things in R.
Data frames can store different object classes in different columns, unlike matrices that have to have the same elements class.

special attributes: row.names, col.names
Rading tabular data from files: read.table() and read.csv()
Convert data.frame to a matrix: data.matrox()    (coercion may happen)
End of explanation
"""


x <- 1:3
names(x) <- c("foo", "bar", "norf")
x
names(x)

"""
Explanation: <a name="names"></a>
Names
R objects can have a name. Useful for self describing data, so each vector can be accesed by an index or by its object names

access and change the names by names()
names in list: x <- list(a=1, b=2, c=3)
names in matrices: dimnames()
End of explanation
"""


## subsetting vectors: numeric or logical index
x <- c("a", "b", "c", "d", "h")
x[1:3]
x[x > "c"]

"""
Explanation: <a name="subset"></a>
2. Subsetting R objects
Top
To extract subsets of data (objects).

subsetting with [ ]   (given numeric index or logical index)
subsetting with double bracket for lists: [[ ]] (only to extract a single element of a list, meaning no vector inside double bracket)
subsetting with $ to extract elements of lists and data.frames by their names
End of explanation
"""


## subsetting matrices
m1[1,2]
m2[,2]

m2[,2, drop=F]

"""
Explanation: For matrices, you can acess rows and columns by [row index, column index]
<b> By default, subsetting drops the column and row index, to avoid that use m[1,2, drop=F] </b>
This way, the output would still be a matrix.
End of explanation
"""


xl <- list(foo=1:5, bar=0.7, norf="ABcDE")
xl[[1]]
xl[["bar"]]
xl["bar"]
xl$norf

"""
Explanation: <a name="subset_lists"></a>
Subsetting lists
Using double bracket to subset lists.
Only extract one element of a list at time. Therefore, no vector inside double bracket. (x[[1:3]] ==> error)
End of explanation
"""


xl[c(1,3)]

"""
Explanation: Use single bracket [ ] To extract multiple elements of a list:
End of explanation
"""


xnest <- list(a=list(10, 12, 13, 14), b=c(4:7))
# to extract 13:
xnest[["a"]][[3]]

"""
Explanation: Nested list
End of explanation
"""


xp <- list(abc=c(1:8), fgh=c(T, F, F, T))
xp$a
xp[["f"]]
xp[["f", exact=F]]

"""
Explanation: <a name="partial_match"></a>
Partial Matching
Partial matching works for subsetting objects with [[ ]] and $
End of explanation
"""


x <- c(1,2,NA, 5, NaN, 8, NaN, 9)
x[!is.na(x)]

y <- c(3, 5, 6, NA, NA, 8, 9, 10)
complete.cases(x,y)

"""
Explanation: <a name="remove_nas"></a>
Removing missing values

In order to subset non-missing values onlt, use !is.na()
Complete cases of a data.frame, a matrix, or multiple vectors complete.cases()
End of explanation
"""


str(file)

con <- gzfile("data/hw1.csv.gz", "r")
xgz <- readLines(con, 10) ## read the first 1 lines
xgz[1]

## read from a web-page
con <- url("http://www.jhsph.edu", "r")
x <- readLines(con)
head(x)

"""
Explanation: <a name="read_write_data"></a>
3. Reading and Writing Data
Top
<a name="read_data"></a>
Reading data

To read in tabular data, use read.table() or read.csv()
readLines() to read lines of a text file
source() for reading R codes (inverse of dump())
dget() for reading R codes (inverse of dput())
load() for reading a saved workspace
unserialize() for reading single R objects in binary form

Attributes of read.table()

file: name of the file
header: logical (T/F) indicating if the file has a header or not
sep: a string that indicates how columns are separated
colClasses: a character vector indicating the class of each column in the data file
nrows: the number of rows in the dataset
comment.char: a character string indicating the comment character used in the file
skip: the number of lines to be skipped from the begining
stringAsFactors: logical, should character variables be coded as factors?

<a name="write_data"></a>
Writing data

write.table()
writeLines()
dump()
dput()
save()
serialize()

<a name="read_large"></a>
Reading large datasets
For efficiently reading in large datasets:

first identify column classes
ignore the comment character if you know there is no comment in the file: comment.char=""

You can automaucally specify column classes by reaading the first 10 rows of the file.
example:
init10 &lt;- read.table("large_dataset.txt", nrows=10)
classes &lt;- sapply(init10, class)
datAll &lt;- read.table("large_dataset.txt", colClasses=classes)

If you know the number of rows, it is more memory efficient to set nrows as well

<a name="textual_format"></a>
Reading and writing in textual format
Textual format has the advantage of readablity and storing element class information.

dput() for a single object
dunpy() for multiple objects

```
example for dupt:
y <- data.frame(a=c(1,2), b=c("HH", "WW"))
dput(y, file="data_y.R")
new.y <- dget("data_t.R")
```
```
example for dput: multiple objects
x <- "foo"
dump(c("x", "y"), file="data_xy.R"))
rm(x,y)
source("data_xy.R")
```
<a name="connection"></a>
Connection interfaces
Data can be read from other sources by creating a connection.

file(): opens a connection to a file
gzfile(): opens a connection to a gz compressed file
bzfile(): opens a connection to a bz compressed file
url(): opens a connection to a web-page
End of explanation
"""


x <- 5
y <- if(x>3) {
   10
}else{
   0
}

print(y)

## three ways to iterate through a vector:
x <- c("a", "b", "c")

for (i in 1:length(x)) {
   print(x[i])
}

for (i in seq_along(x))  print(x[i])

for (letter in x)  print(letter)

"""
Explanation: <a name="flow_control"></a> Top
4. Flow Control

if, else, else if
for
while
repeat an infinite loop, (to stop use break)
break break the execution of a loop
next skip one iteration
return exit a function
End of explanation
"""


args(lm)

"""
Explanation: <a name="functions"></a>
5. Functions
Top
Functions are first class objects.
We can have nested functions, functions inside other functions.
Function arguments can have default values. Function arguments can be mapped positionally, or by name.
The arguments that don't have default values, are required when function is called.
Note: Split your functions so that each function does a single task.
End of explanation
"""


myplot <- function(x, y, mytype="l", ...) {
   plot(x, y, type=mytype)
}

"""
Explanation: <a name="special_arg"></a>
special argument "..."
Three dots indicates variable number of arguments.
It is usually used when extending another function:
End of explanation
"""


args(paste)

args(cat)

paste("a", "b", sep=":")
# partial matching doesn't work after ...
paste("a", "b", se=":")

"""
Explanation: ... is also useful when the number of input arguments is not known in advance:
End of explanation
"""


search()

make.Gaussian <- function(mu, sigma) {
    gauss <- function(x) {
       exp(-((x-mu)/sigma)^2)
    }
    gauss
}

mygauss_0 <- make.Gaussian(mu=0, sigma=1)
mygauss_1 <- make.Gaussian(mu=1, sigma=1)
mygauss_0

library(IRdisplay)

x <- seq(-3, 5, by=0.1)

png("figs/mygauss.png")
par(mar=c(4.7, 4.7, 0.6, 0.6))
plot(x, mygauss_0(x), type="l", col=rgb(0.2, 0.4, 0.7, 0.5), lwd=6, xlab="x", ylab="Gaussian Distn.", cex.axis=1.5, cex.lab=1.5)
lines(x, mygauss_1(x), col=rgb(0.7, 0.3, 0.2, 0.5), lwd=6)
dev.off()
PNG("figs/mygauss.png")

"""
Explanation: <a name="scoping"></a>
6. Scoping Rules
Top
R binds a value to every symbol.
If you define a new objects with a name that is used previously in R:
R uses lexical scoping

R first searches for a free variable in the environment in which the function is defined
If no matched found, it searches in the parent environment

The global enfironemnt is the workspace.
When a free variable is found, R searches through the search space, from the workspace (global environment) until the symbol is found or it hits the emoty environment.
End of explanation
"""


t <- as.Date("2014-05-11")
print(t)
unclass(t)  ## number of days since 1970-01-01


x <- Sys.time()
print(x)
p <- as.POSIXlt(x)
names(unclass(p))
p$hour

"""
Explanation: <a name="date_time"></a>
7. Dates and Times
Top

dates are represented by Date class
times are represented by POSIXct or POSIXlt classes
POSIXct is a very large integer; It uses a data frame to store time
POSIXlt store other information, such as day of month, day of year, ...

internal reference time is 1970-01-01


weekdays() what day of the week that time/date is

months() gives the month name

quarters() gives the quarter number


objects of the same class can be used in mathematical operations
End of explanation
"""


dstring <- c("January 01 2012, 10:40", "March 05 2013, 11:20")
strptime(dstring, "%B %d %Y, %H:%M")

## the formats can be found in help page ?strptime

"""
Explanation: <a name="format_strptime"></a>
Format Strings to dates
To convert a string in an arbitrary format to a Date format:
End of explanation
"""


## Example for lapply()
x <- list(a=1:20, b=rnorm(10, 2))
lapply(x, mean)
sapply(x, mean)

## Example for apply()
x <- matrix(rnorm(50, mean=1), nrow=10, ncol=5)
apply(x, 2, mean)
apply(x, 1, mean)

"""
Explanation: <a name="apply_functions"></a>
8. Apply Functions
Top
Apply functions can be used to apply a function over elements of a vector, list or data frame.
They make the codes shorter, and faster through vectorized operations.

lapply()  loop over elements of a list, and evaluate a function on each element.
if input is not a list, it will be corced to a list (if possible)


sapply()  same as lapply, but also simplifies the result
apply()   apply a function over margins of an array
tapply()  apply a function over subsets of a vector
mapply()  multivariate version of lapply
End of explanation
"""


# Benchmark: Comparing apply() and rowMeans()
library(rbenchmark)
a <- matrix(rnorm(100000), nrow=1000)

benchmark(apply(a, 1, mean))
benchmark(rowMeans(a))

"""
Explanation: <a name="shortcut_apply"></a>
Shortcut functions
For means and sums of a matrix, these shortcut functions perform much faster on big matrices.

rowSums() equivalent to apply(x, 1, sum)
colSums() equivalent to apply(x, 2, sum)
rowMeans() equivalent to apply(x, 1, mean)
colMeans() equivalent to apply(x, 2, mean)
End of explanation
"""


## Example for tapply()
x <- c(rnorm(10, mean=3), runif(10, min=1, max=3), rnorm(10, mean=5, sd=2))
f <- gl(3, 10)
tapply(x, f, mean)

"""
Explanation: <a name="tapply_split"></a>
tapply() and split()
tapplly() can be used to apply a function on each sub-group (determied by a factor vector).
End of explanation
"""


split(x, f)

library(datasets)
s <- split(airquality, airquality$Month)

sapply(s, function(x) colMeans(x[,c("Ozone", "Solar.R", "Wind", "Temp")], na.rm=T))

"""
Explanation: split()  is not a loop function, but it can be used to split a vector based on a given factor vector.
The output is a list.
End of explanation
"""


x <- rnorm(n=10, mean=20, sd=5)
x
str(rpois)

"""
Explanation: <a name="random"></a>
9. Random Number Generation
Top

rnorm() generates normal random number variates with a given mean and SD
rnorm(n, mean=0, sd=1)
dnorm() returns the probablity density function at a point x, with given mean and CD
pnorm() returns the cumulative distribution function at point x
qnorm() quantile
rpois() generates random numbers with Poisson distribution with a given rate (lambda)
rbinom() generates binomial random numbers with a given probability
End of explanation
"""


## adding Gaussian noise to sin(x)
x <- seq(0, 2*pi, by=0.02)
y <- sapply(sin(x), function(x) rnorm(n=1, mean=sin(x), sd=0.2))
## the same as y <- sin(x) + rnorm(n=length(x), mean=0, sd=0.2)

par(mar=c(4.7, 4.7, 0.7, 0.7))
plot(x, y, pch=19, col=rgb(0.2, 0.3, 0.7, 0.6), cex=1.3, cex.lab=1.6, cex.axis=1.6)

"""
Explanation: <a name="gauss_noise"></a>
Gaussian Noise
End of explanation
"""


x <- rnorm(200)
log.mu <- 0.5 + 0.3 * x
y <- rpois(200, exp(log.mu))

plot(x, y, pch=19, col=rgb(0.2, 0.3, 0.7, 0.6), cex=1.5, cex.axis=1.6, cex.lab=1.6)

"""
Explanation: Use set.seed() to make the random numbers reproducible

<a name="general_linear"></a>
Simulate random numbers from a Generalized Linear Model
Assume $y \approx Poisson(\mu)$
and $log(\mu) = \beta_0 + \beta_1 x$
End of explanation
"""


str(sample)
sample(letters, 5)
sample(1:10, size=6, replace=T)

"""
Explanation: <a name="random_sample"></a>
Random Sampling

sample() to draw a random sample from a set
replace=TRUE 
replace=FALSE
End of explanation
"""


system.time(rnorm(n=100000, mean=0, sd=1))

system.time(readLines("http://www.jhsph.edu"))

"""
Explanation: <a id="rprofiler"></a>
10. R Profiler


First design your code, unit-test it, then optimize it


Using system.time() to measure how much time it takes to perform a step

returns an object of class proc_time
user time  --> CPU time
elapsed time  --> wallclock time
End of explanation
"""


system.time({
    n <- 1000
    r <- numeric(n)
    for (i in 1:n) {
        x = rnorm(i)
        r[i] = mean(x)
    }
})

# comparing time of vectorized operations vs. for loops
n=1000000
x <- rnorm(n, mean=0, sd=1)

system.time(sum(x))
system.time({
    s = 0
    for (i in 1:n) {
        s = s + x[i]
    }
})

"""
Explanation: Measureing multiple expressions:
End of explanation
"""


Rprof(filename="data/rprof_test.dat")

for (i in 1:1000) {
    x <- seq(0, 5, by=0.01)
    y <- 1.5+x + rnorm(n=length(x), mean=0, sd=1)
    lm(y ~ x)
}

Rprof(NULL)

summaryRprof(filename="data/rprof_test.dat")

"""
Explanation: <a id="using-rprof"></a>
Using Rprof()

summaryRprof() will tabulates the output
End of explanation
"""


## define a constructor function:
make.NegLogLik <- function(data, fixed=c(FALSE, FALSE)) {
    params <- fixed
    function(p) {
        params[!fixed] <- p
        mu <- params[1]
        sigma <- params[2]
        a <- -0.5 * length(data) * log(2*pi*sigma^2)
        b <- -0.5 * sum((data-mu)^2) / (sigma^2)
        -(a + b)
    }
}

# Generate a random dataset with normal distribution
d <- rnorm(n=1000, mean=2, sd=3)
nLL <- make.NegLogLik(d)
nLL
ls(environment(nLL))

## Minimze Neative of log-likelihood
optim(c(mu=0, sigma=1), nLL)$par

"""
Explanation: <a id="r-optim"></a>
11. Optimization in R


optim(), nlm(), optimize()


Pass them a function and a vector of parameters. 

We can hold some parameters fixed, by fixed=c(F, T)

Example: maximize log-likelihood
End of explanation
"""


## generate a range of possible values for mean
x <- seq(1, 3, len=100)
## when fixing sigma, use the true value
nLL.mean <- make.NegLogLik(d, fixed=c(FALSE, 3)) 
y <- sapply(x, nLL.mean)

#windows.options(width=40, height=10)
#layout(matrix(c(1,2), 1, 2, byrow = TRUE),
#   widths=c(5,5), heights=c(1,1))



par(mar=c(4.7, 4.7, 0.7, 0.7), mfrow=c(1,2))
plot(x, -(y-min(y)), type="l", lwd=4, col="steelblue", 
     cex.axis=1.7, cex.lab=1.7, xlab=expression(mu), ylab="log-likelihood")
plot(x, exp(-(y-min(y))), type="l", lwd=4, col="steelblue", 
     cex.axis=1.7, cex.lab=1.7, xlab=expression(mu), ylab="likelihood")


"""
Explanation: <a id="plot-loglik"><a>
Plot the log-likelihood function
End of explanation
"""


