#' ---	
#' title: "Data Science and Predictive Analytics (UMich HS650)"	
#' subtitle: "<h2><u>Specialized Machine Learning Topics</u></h2>"	
#' author: "<h3>SOCR/MIDAS (Ivo Dinov)</h3>"	
#' date: "`r format(Sys.time(), '%B %Y')`"	
#' tags: [DSPA, SOCR, MIDAS, Big Data, Predictive Analytics] 	
#' output:	
#'   html_document:	
#'     theme: spacelab	
#'     highlight: tango	
#'     includes:	
#'       before_body: SOCR_header.html	
#'       after_body: SOCR_footer_tracker.html	
#'     toc: true	
#'     number_sections: true	
#'     toc_depth: 2	
#'     toc_float:	
#'       collapsed: false	
#'       smooth_scroll: true	
#' ---	
#' 	
#' In this chapter, we will discuss some technical details about data formats, streaming, optimization of computation, and distributed deployment of optimized learning algorithms. [Chapter 21](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/21_FunctionOptimization.html) provides additional optimization details.	
#' 	
#' The Internet of Things (IoT) leads to a paradigm shift of scientific inference - from static data interrogated in a batch or distributed environment to an on-demand service-based Cloud computing. Here, we will demonstrate how to work with specialized data, data-streams, and SQL databases, as well as develop and assess on-the-fly data modeling, classification, prediction and forecasting methods. Important examples to keep in mind throughout this chapter include high-frequency data delivered real time in hospital ICU's (e.g., [microsecond Electroencephalography signals, EEGs](https://physionet.org/physiobank/database/)), dynamically changing stock market data (e.g., [Dow Jones Industrial Average Index, DJI](http://www.marketwatch.com/investing/index/djia)), and [weather patterns](http://weather.rap.ucar.edu/surface/).	
#' 	
#' We will present (1) format conversion and working with XML, SQL, JSON, CSV, SAS and other data objects, (2) visualization of bioinformatics and network data, (3) protocols for managing, classifying and predicting outcomes from data streams, (4) strategies for optimization, improvement of computational performance, parallel (MPI) and graphics (GPU) computing, and (5) processing of very large datasets.	
#' 	
#' 	
#' 	
require("knitr")	
opts_knit$set(root.dir = "C:\\Users\\Dinov\\Desktop")	
#' 	
#' 	
#' # Working with specialized data and databases	
#' 	
#' Unlike the case-studies we saw in the previous chapters, some real world data may not always be nicely formatted, e.g., as CSV files. We must collect, arrange, wrangle, and harmonize scattered information to generate computable data objects that can be further processed by various techniques. Data wrangling and preprocessing may take over 80% of the time researchers spend interrogating complex multi-source data archives. The following procedures will enhance your skills collecting and handling heterogeneous real world data.  Multiple examples of handling long-and-wide data, messy and tidy data, and data cleaning strategies can be found in this [JSS `Tidy Data` article by Hadley Wickham](https://www.jstatsoft.org/article/view/v059i10).	
#' 	
#' ## Data format conversion	
#' 	
#' The R package `rio` imports and exports various types of file formats, e.g., tab-separated (`.tsv`), comma-separated (`.csv`), JSON (`.json`), Stata (`.dta`), SPSS (`.sav` and `.por`), Microsoft Excel (`.xls` and `.xlsx`), Weka (`.arff`), and SAS (`.sas7bdat` and `.xpt`).	
#' 	
#' `rio` provides three important functions `import()`, `export()` and `convert()`. They are intuitive, easy to understand, and efficient to execute. Take Stata (.dta) files as an example. First we can download [02_Nof1_Data.dta](https://umich.instructure.com/files/1760330/download?download_frd=1) from our [datasets folder](https://umich.instructure.com/courses/38100/files/folder/data). 	
#' 	
#' 	
# install.packages("rio")	
library(rio)	
# Download the SAS .DTA file first locally	
# Local data can be loaded by: 	
#nof1<-import("02_Nof1_Data.dta")	
# the data can also be loaded from the server remotely as well:	
nof1<-read.csv("https://umich.instructure.com/files/330385/download?download_frd=1")	
str(nof1)	
#' 	
#' 	
#' The data is automatically stored as a data frame. Note that `rio` sets `stingAsFactors=FALSE` as default.	
#' 	
#' `rio` can help us export files into any other format we chose. To do this we have to use the `export()` function. 	
#' 	
#' 	
#Sys.getenv("R_ZIPCMD", "zip")   # Get the C Zip application 	
Sys.setenv(R_ZIPCMD="E:/Ivo.dir/Ivo_Tools/ZIP/bin/zip.exe")	
Sys.getenv("R_ZIPCMD", "zip")	
export(nof1, "C:/Users/Dinov/Desktop/02_Nof1.xlsx")	
#' 	
#' 	
#' This line of code exports the *Nof1* data in `xlsx` format located in the `R` working directory. Mac users may have a problem exporting `*.xslx` files using `rio` because of a lack of a zip tool, but still can output other formats such as ".csv". An alternative strategy to save an `xlsx` file is to use package `xlsx` with default `row.name=TRUE`. 	
#' 	
#' `rio` also provides a one-step process to convert-and-save data into alternative formats. The following simple code allows us to convert and save the `02_Nof1_Data.dta` file we just downloaded into a CSV file.	
#' 	
#' 	
# convert("02_Nof1_Data.dta", "02_Nof1_Data.csv")	
convert("C:/Users/Dinov/Desktop/02_Nof1.xlsx", "C:/Users/Dinov/Desktop/02_Nof1_Data.csv")	
#' 	
#' 	
#' You can see a new CSV file pop-up in the working directory. Similar transformations are available for other data formats and types.	
#' 	
#' ## Querying data in SQL databases	
#' 	
#' The [CDC](https://www.cdc.gov) [Behavioral Risk Factor Surveillance System (BRFSS) Data, 2013-2015](https://www.cdc.gov/brfss/annual_data/annual_2015.html). 	
#' This file for the combined landline and cell phone data set was exported from SAS V9.3 in the XPT transport format. This file contains 330 variables. This format can be imported into SPSS or STATA. Please note: some of the variable labels get truncated in the process of converting to the XPT format	
#' 	
#' Be careful - this compressed (ZIP) file is over 315MB in size!	
#' 	
#' 	
# install.packages("Hmisc")	
library(Hmisc)	
	
memory.size(max=T)	
	
pathToZip <- tempfile()	
download.file("http://www.socr.umich.edu/data/DSPA/BRFSS_2013_2014_2015.zip", pathToZip)	
# let's just pull two of the 3 years of data (2013 and 2015)	
brfss_2013 <- sasxport.get(unzip(pathToZip)[1])	
brfss_2015 <- sasxport.get(unzip(pathToZip)[3])	
	
dim(brfss_2013); object.size(brfss_2013)	
# summary(brfss_2013[1:1000, 1:10])  # subsample the data	
	
# report the summaries for 	
summary(brfss_2013$has_plan)	
brfss_2013$x.race <- as.factor(brfss_2013$x.race)	
summary(brfss_2013$x.race)	
	
# clean up	
unlink(pathToZip)	
#' 	
#' 	
#' Let's try to use logistic regression to find out if self-reported race/ethnicity predicts the binary outcome of having a health care plan.	
#' 	
#' 	
brfss_2013$has_plan <- brfss_2013$hlthpln1 == 1 	
	
system.time(	
  gml1 <- glm(has_plan ~ as.factor(x.race), data=brfss_2013,	
              family=binomial)	
)   # report execution time	
summary(gml1)	
#' 	
#' 	
#' Next, we'll examine the [odds](http://wiki.socr.umich.edu/index.php/SMHS_OR_RR) (rather the log odds ration, LOR) of having a health care plan (HCP) by race (R). The LORs are calculated for two array dimensions, separately for each *race* level (presence of *health care plan* (HCP) is binary, whereas *race* (R) has 9 levels, $R1, R2, ..., R9$). For example, the odds ratio of having a HCP for $R1:R2$ is:	
#' 	
#' $$ OR(R1:R2) =  \frac{\frac{P \left( HCP \mid R1 \right)}{1 - P \left( HCP \mid R1 \right)}}{\frac{P \left( HCP \mid R2 \right)}{1 - P \left( HCP \mid R2 \right)}} .$$	
#' 	
#load the vcd package to compute the LOR	
library("vcd")	
	
lor_HCP_by_R <- loddsratio(has_plan ~ as.factor(x.race), data = brfss_2013)	
lor_HCP_by_R	
	
#' 	
#' 	
#' Now, let's see an example of querying a database containing structured relational collection of data records. A *query* is a machine instruction (typically represented as text) sent by a user to remote database requesting a specific database operation (e.g., search or summary). One database communication protocol relies on SQL (Structured query language). MySQL is an instance of a database management system that supports SQL communication that many web applications utilize, e.g., YouTube, Flickr, Wikipedia, biological databases like GO, ensembl, etc. Below is an example of an SQL query using the package `RMySQL`. An alternative way to interface an SQL database is using the package `RODBC`,	
#' 	
#' 	
# install.packages("DBI"); install.packages("RMySQL")	
# install.packages("RODBC"); library(RODBC)	
library(DBI)	
library(RMySQL)	
	
ucscGenomeConn <- dbConnect(MySQL(),	
                user='genome',	
                dbname='hg38',	
                host='genome-mysql.cse.ucsc.edu')	
	
result <- dbGetQuery(ucscGenomeConn,"show databases;"); 	
	
# List the DB tables	
allTables <- dbListTables(ucscGenomeConn); length(allTables)	
	
# Get dimensions of a table, read and report the head	
dbListFields(ucscGenomeConn, "affyU133Plus2")	
affyData <- dbReadTable(ucscGenomeConn, "affyU133Plus2"); head(affyData)	
	
# Select a subset, fetch the data, and report the quantiles	
subsetQuery <- dbSendQuery(ucscGenomeConn, "select * from affyU133Plus2 where misMatches between 1 and 3")	
affySmall <- fetch(subsetQuery); quantile(affySmall$misMatches)	
	
# Get repeat mask	
bedFile <- 'repUCSC.bed'	
df <- dbSendQuery(ucscGenomeConn,'select genoName,genoStart,genoEnd,repName,swScore, strand,	
                  repClass, repFamily from rmsk') %>%	
        dbFetch(n=-1) %>%	
        mutate(genoName = str_replace(genoName,'chr','')) %>%	
        tbl_df %>%	
        write_tsv(bedFile,col_names=F)	
message('written ', bedFile)	
	
# Once done, close the connection	
dbDisconnect(ucscGenomeConn)	
#' 	
#' 	
#' To complete the above database SQL commands, it requires access to the remote UCSC SQL Genome server and user-specific credentials. The example below can be done by all users, as it relies only on local services.	
#' 	
#' 	
# install.packages("RSQLite")	
library("RSQLite")	
	
# generate an empty DB and store it in RAM	
myConnection <- dbConnect(RSQLite::SQLite(), ":memory:")	
myConnection	
dbListTables(myConnection)	
	
# Add tables to the local SQL DB	
data(USArrests); dbWriteTable(myConnection, "USArrests", USArrests)	
dbWriteTable(myConnection, "brfss_2013", brfss_2013)	
dbWriteTable(myConnection, "brfss_2015", brfss_2015)	
	
# Check again the DB content	
dbListFields(myConnection, "brfss_2013")	
dbListTables(myConnection); 	
	
# Retrieve the entire DB table (for the smaller USArrests table)	
dbGetQuery(myConnection, "SELECT * FROM USArrests")	
	
# Retrieve just the average of one feature	
myQuery <- dbGetQuery(myConnection, "SELECT avg(Assault) FROM USArrests"); myQuery	
	
myQuery <- dbGetQuery(myConnection, "SELECT avg(Assault) FROM USArrests GROUP BY UrbanPop"); myQuery	
	
# Or do it in batches (for the much larger brfss_2013 and brfss_2015 tables)	
myQuery <- dbGetQuery(myConnection, "SELECT * FROM brfss_2013")	
	
# extract data in chunks of 2 rows, note: dbGetQuery vs. dbSendQuery	
# myQuery <- dbSendQuery(myConnection, "SELECT * FROM brfss_2013")	
# fetch2 <- dbFetch(myQuery, n = 2); fetch2	
# do we have other cases in the DB remaining?	
# extract all remaining data	
# fetchRemaining <- dbFetch(myQuery, n = -1);fetchRemaining	
# We should have all data in DB now	
# dbHasCompleted(myQuery)	
	
# compute the average (poorhlth) grouping by Insurance (hlthpln1)	
# Try some alternatives: numadult nummen numwomen genhlth physhlth menthlth poorhlth hlthpln1	
myQuery1_13 <- dbGetQuery(myConnection, "SELECT avg(poorhlth) FROM brfss_2013 GROUP BY hlthpln1"); myQuery1_13	
	
# Compare 2013 vs. 2015: Health grouping by Insurance	
myQuery1_15 <- dbGetQuery(myConnection, "SELECT avg(poorhlth) FROM brfss_2015 GROUP BY hlthpln1"); myQuery1_15	
	
myQuery1_13 - myQuery1_15	
	
# reset the DB query	
# dbClearResult(myQuery)	
	
# clean up	
dbDisconnect(myConnection)	
#' 	
#' 	
#' ## Real Random Number Generation	
#' 	
#' We are already familiar with (pseudo) random number generation (e.g., `rnorm(100, 10, 4)` or `runif(100, 10,20)`), which generate *algorithmically* computer values subject to specified distributions. There are also web-services, e.g., [random.org](http://random.org), that can provide *true random* numbers based on atmospheric noise, rather than using a pseudo random number generation protocol. Below is one example of generating a total of 300 numbers arranged in 3 columns, each of 100 rows of random integers (in decimal format) between 100 and 200.	
#' 	
#' 	
siteURL <- "http://random.org/integers/"  #   base URL	
shortQuery<-"num=300&min=100&max=200&col=3&base=10&format=plain&rnd=new"	
completeQuery <- paste(siteURL, shortQuery, sep="?")  # concat url and submit query string	
rngNumbers <- read.table(file=completeQuery)        # and read the data	
rngNumbers      	
#' 	
#' 	
#' ## Downloading the complete text of web pages	
#' 	
#' `RCurl` package provides an amazing tool for extracting and scraping information from websites. Let's install it and extract information from a SOCR website.	
#' 	
#' 	
# install.packages("RCurl")	
library(RCurl)	
web<-getURL("http://wiki.socr.umich.edu/index.php/SOCR_Data", followlocation = TRUE)	
str(web, nchar.max = 200)	
#' 	
#' 	
#' The `web` object looks incomprehensible. This is because most websites are wrapped in XML/HTML hypertext or include JSON formatted meta-data. `RCurl` deals with special HTML tags and website meta-data.	
#' 	
#' To deal with the web pages only, `httr` package would be a better choice than `RCurl`. It returns a list that makes much more sense.	
#' 	
#' 	
#install.packages("httr")	
library(httr)	
web<-GET("http://wiki.socr.umich.edu/index.php/SOCR_Data")	
str(web[1:3])	
#' 	
#' 	
#' ## Reading and writing XML with the `XML` package	
#' 	
#' A combination of the `RCurl` and the `XML` packages could help us extract only the plain text in our desired webpages. This would be very helpful to get information from heavy text based websites. 	
#' 	
#' 	
web<-getURL("http://wiki.socr.umich.edu/index.php/SOCR_Data", followlocation = TRUE)	
#install.packages("XML")	
library(XML)	
web.parsed<-htmlParse(web, asText = T)	
plain.text<-xpathSApply(web.parsed, "//p", xmlValue)	
cat(paste(plain.text, collapse = "\n"))	
#' 	
#' 	
#' Here we extracted all plain text between the starting and ending *paragraph* HTML tags, `<p>` and `</p>`. 	
#' 	
#' More information about [extracting text from XML/HTML to text via XPath is available here](http://www.r-bloggers.com/htmltotext-extracting-text-from-html-via-xpath).	
#' 	
#' ## Web-page Data Scraping	
#' 	
#' The process that extracting data from complete web pages and storing it in structured data format is called `scraping`. However, before starting a data scrap from a website, we need to understand the underlying HTML structure for that specific website. Also, we have to check the terms of that website to make sure that scraping from this site is allowed.	
#' 	
#' The R package `rvest` is a very good place to start "harvesting" data from websites.	
#' 	
#' To start with, we use `read_html()` to store SOCR data website into a `xmlnode` object. 	
#' 	
#' 	
library(rvest)	
SOCR<-read_html("http://wiki.socr.umich.edu/index.php/SOCR_Data")	
SOCR	
#' 	
#' 	
#' From the summary structure of `SOCR`, we can discover that there are two important hypertext section markups `<head>` and `<body>`. Also, notice that the SOCR data website uses `<title>` and `</title>` tags to separate title in the `<head>` section. Let's use `html_node()` to extract title information based on this knowledge. 	
#' 	
#' 	
SOCR %>% html_node("head title") %>% html_text()	
#' 	
#' 	
#' Here we used `%>%` operator, or pipe, to connect two functions. The above line of code creates a chain of functions to operate on the `SOCR` object. The first function in the chain `html_node()` extracts the `title` from `head` section. Then, `html_text()` translates HTML formatted hypertext into English. [More on `R` piping can be found in the `magrittr` package](https://cran.r-project.org/web/packages/magrittr/vignettes/magrittr.html).	
#' 	
#' Another function, `rvest::html_nodes()` can be very helpful in scraping. Similar to `html_node()`, `html_nodes()` can help us extract multiple nodes in an `xmlnode` object. Assume that we want to obtain the meta elements (usually page description, keywords, author of the document, last modified, and other metadata) from the SOCR data website. We apply `html_nodes()` to `SOCR` object for lines start with `<meta` in `<head>` section. It is optional to use `html_attrs()`(extracts attributes, text and tag name from html) to make texts prettier.	
#' 	
#' 	
meta<-SOCR %>% html_nodes("head meta") %>% html_attrs()	
meta	
#' 	
#' 	
#' ## Parsing JSON from web APIs	
#' 	
#' Application Programming Interfaces (APIs) allow web-accessible functions to communicate with each other. Today most API is stored in JSON (JavaScript Object Notation) format.	
#' 	
#' JSON represents a plain text format used for web applications, data structures or objects. Online JSON objects could be retrieved by packages like `RCurl` and `httr`. Let's see a JSON formatted dataset first. We can use [02_Nof1_Data.json](https://umich.instructure.com/files/1760327/download?download_frd=1) in the class file as an example. 	
#' 	
#' 	
library(httr)	
nof1<-GET("https://umich.instructure.com/files/1760327/download?download_frd=1")	
nof1	
#' 	
#' 	
#' We can see that JSON objects are very simple. The data structure is organized using hierarchies marked by square brackets. Each piece of information is formatted as a `{key:value}` pair.	
#' 	
#' The package `jsonlite` is a very useful tool to import online JSON formatted datasets into data frame directly. Its syntax is very straight forward. 	
#' 	
#' 	
#install.packages("jsonlite")	
library(jsonlite)	
nof1_lite<-fromJSON("https://umich.instructure.com/files/1760327/download?download_frd=1")	
class(nof1_lite)	
#' 	
#' 	
#' ## Reading and writing Microsoft Excel spreadsheets using XLSX	
#' 	
#' We can transfer a *xlsx* dataset into CSV and use `read.csv()` to load this kind of dataset. However, R provides an alternative `read.xlsx()` function in package `xlsx` to simplify this process. Take our `02_Nof1_Data.xls` data in the class file as an example. We need to download the file first.	
#' 	
#' 	
# install.packages("xlsx")	
library(xlsx)	
nof1<-read.xlsx("C:/Users/Dinov/Desktop/02_Nof1.xlsx", 1)	
str(nof1)	
#' 	
#' 	
#' The last argument, `1`, stands for the first excel sheet, as any excel file may include a large number of tables in it. Also, we can download the `xls` or `xlsx` file into our R working directory so that it is easier to find file path.	
#' 	
#' # Working with domain-specific data	
#' 	
#' Powerful Machine-Learning methods have already been applied in many fields. Some of them are very specialized and require unique approaches to address their characteristics.	
#' 	
#' ## Working with bioinformatics data	
#' 	
#' Genetic data are stored in widely varying formats and usually have more feature variables than observations. They could have 1,000 columns and only 200 rows. One of the commonly used pre-processng steps for such datasets is *variable selection*. We will talk about this in [Chapter 16]( http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/16_FeatureSelection.html).	
#' 	
#' The Bioconductor project created powerful R functionality (packages and tools) for analyzing genomic data, see [Bioconductor for more detailed information](http://www.bioconductor.org).	
#' 	
#' ## Visualizing network data	
#' 	
#' Social network data and graph datasets describe the relations between nodes (vertices) using connections (links or edges) joining the node objects. Assume we have *N* objects, we can have $N*(N-1)$ directed links establishing paired associations between the nodes. Let's use an example with *N=4* to demonstrate a simple graph potentially modeling the following linkage table.	
#' 	
#' objects| 1  | 2  | 3  | 4	
#' -------|----|----|----|---	
#' 1|.....| $1\rightarrow 2$|$1\rightarrow 3$|$1\rightarrow 4$	
#' 2|$2\rightarrow 1$|.....|$2\rightarrow 3$|$2\rightarrow 4$	
#' 3|$3\rightarrow 1$|$3\rightarrow 2$|.....|$3\rightarrow 4$      	
#' 4|$4\rightarrow 1$|$4\rightarrow 2$|$4\rightarrow 3$|.....	
#' 	
#' If we change the $a\rightarrow b$ to an indicator variable (0 or 1) capturing whether we have an edge connecting a pair of nodes, then we get the graph *adjacency matrix*.	
#' 	
#' Edge lists provide an alternative way to represent network connections. Every line in the list contains a connection between two nodes (objects).	
#' 	
#' Vertex|Vertex	
#' ------|------	
#' 1 |2	
#' 1 |3	
#' 2 |3	
#' 	
#' The above edge list is listing three network connections: object 1 is linked to object 2; object 1 is linked to object 3; and object 2 is linked to object 3. Note that edge lists can represent both *directed* as well as *undirected* networks or graphs. 	
#' 	
#' We can imagine that if *N* is very large, e.g., social networks, the data representation and analysis may be resource intense (memory or computation). In R, we have multiple packages that can deal with social network data. One user-friendly example is provided using the `igraph` package. First, let's build a toy example and visualize it using this package.	
#' 	
#' 	
#install.packages("igraph")	
library(igraph)	
g<-graph(c(1, 2, 1, 3, 2, 3, 3, 4), n=10)	
plot(g)	
#' 	
#' 	
#' Here `c(1, 2, 1, 3, 2, 3, 3, 4)` is an edge list with 4 rows and `n=10` means we have 10 nodes (objects) in total. The small arrows in the graph shows us directed network connections. We might notice that 5-10 nodes are scattered out in the graph. This is because they are not included in the edge list, so there are no network connections between them and the rest of the network.  	
#' 	
#' Now let's examine the [co-appearance network of Facebook circles](https://snap.stanford.edu/data/egonets-Facebook.html). The data contains anonymized `circles` (friends lists) from Facebook collected from survey participants using [a Facebook app](https://www.facebook.com/apps/application.php?id=201704403232744). The dataset only includes edges (circles, 88, 234) connecting pairs of nodes (users, 4, 039) in the ego networks.	
#' 	
#' The values on the connections represent the number of links/edges within a circle. We have a huge edge-list made of scrambled Facebook user IDs. Let's load this dataset into R first. The data is stored in a text file. Unlike CSV files, text files in table format need to be imported using `read.table()`. We are using `header=F` option to let R know that we don't have a header in the text file that contains only tab-separated node pairs (indicating the social connections, edges, between Facebook users).	
#' 	
#' 	
soc.net.data<-read.table("https://umich.instructure.com/files/2854431/download?download_frd=1", sep=" ", header=F)	
head(soc.net.data)	
#' 	
#' 	
#' Now the data is stored in a data frame. To make this dataset ready for `igraph` processing and visualization, we need to convert `soc.net.data` into a matrix object.	
#' 	
#' 	
soc.net.data.mat <- as.matrix(soc.net.data, ncol=2)	
#' 	
#' 	
#' By using `ncol=2`, we made a matrix with two columns. The data is now ready and we can apply `graph.edgelist()`.	
#' 	
#' 	
# remove the first 347 edges (to wipe out the degenerate "0" node)	
graph_m<-graph.edgelist(soc.net.data.mat[-c(0:347), ], directed = F)	
#' 	
#' 	
#' Before we display the social network graph we may want to examine our model first.	
#' 	
#' 	
summary(graph_m)	
#' 	
#' 	
#' This is an extremely brief yet informative summary. The first line `U--- 4038 87887` includes potentially four letters and two numbers. The first letter could be `U` or `D` indicating undirected or directed edges. A second letter `N` would mean that the objects set has a "name" attribute. A third letter is for weighted (`W`) graph. Since we didn't add weight in our analysis the third letter is empty ("`-`"). A forth character is an indicator for bipartite graphs (whose vertices can be divided into `two disjoint sets` and (that is, represent independent sets where each vertex from one set connects to one vertex in the other set). The two numbers following the 4 letters represent the `number of nodes` and the `number of edges`, respectively. Now let's render the graph. 	
#' 	
#' 	
plot(graph_m)	
#' 	
#' 	
#' This graph is very complicated. We can still see that some words are surrounded by more nodes than others. To obtain such information we can use `degree()` function which list the number of edges for each node.	
#' 	
#' 	
degree(graph_m)	
#' 	
#' 	
#' Skimming the table we can find that `the 107-th` user has as many as 1,044 connections, which makes the user a *highly-connected hub*. Likely, this node may have higher social relevance.	
#' 	
#' Some edges might be more important than other edges because they serve as a bridge to link a cloud of nodes. To compare their importance, we can use the betweenness centrality measurement. *Betweenness centrality* measures centrality in a network. High centrality for a specific node indicates influence. `betweenness()` can help us to calculate this measurement.	
#' 	
#' 	
betweenness(graph_m)	
#' 	
#' 	
#' Again, `the 107-th` node has the highest betweenness centrality ($3.556221e+06$).	
#' 	
#' We can try another example using [SOCR hierarchical data, which is also available for dynamic exploration as a tree graph](http://socr.ucla.edu/SOCR_HT_ResourceJavaViewer.html). Let's read its JSON data source using the `jsonlite` package.	
#' 	
#' 	
tree.json<-fromJSON("http://socr.ucla.edu/SOCR_HyperTree.json", simplifyDataFrame = FALSE)	
#' 	
#' 	
#' This generates a `list` object representing the hierarchical structure of the network. Note that this is quite different from edge list. There is one root node, its sub nodes are called *children nodes*, and the terminal notes are call *leaf nodes*. Instead of presenting the relationship between nodes in pairs, this hierarchical structure captures the level for each node. To draw the social network graph, we need to convert it as a `Node` object. We can utilize `as.Node()` function in `data.tree` package to do so.	
#' 	
#' 	
# install.packages("data.tree")	
library(data.tree)	
tree.graph<-as.Node(tree.json, mode = "explicit")	
#' 	
#' 	
#' Here we use `mode="explicit"` option to allow "children" nodes to have their own "children" nodes. Now, the `tree.json` object has been separated into four different node structures - `"About SOCR", "SOCR Resources", "Get Started", ` and `"SOCR Wiki"`. Let's plot the first one using `igraph` package.	
#' 	
#' 	
plot(as.igraph(tree.graph$`About SOCR`), edge.arrow.size=5, edge.label.font=0.05)	
#' 	
#' 	
#' In this graph, `"About SOCR"` that located at the center is the root.	
#' 	
#' # Data Streaming	
#' 	
#' The proliferation of Cloud services and the emergence of modern technology in all aspects of human experiences leads to a tsunami of data much of which is steamed real-time. The interrogation of such voluminous data is an increasingly important area of research. *Data streams* are ordered, often unbounded sequences of data points created continuously by a data generator. All of the data mining, interrogation and forecasting methods we discuss here are also applicable to data streams.	
#' 	
#' ## Definition	
#' Mathematically, a *data stream* in an ordered sequence of data points	
#' $$Y = \{y_1, y_2, y_3, \cdots, y_t, \cdots \},$$	
#' where the (time) index, $t$, reflects the order of the observation/record, which may be single numbers, simple vectors in multidimensional space, or objects, e.g., [structured Ann Arbor Weather (JSON)](http://weather.rap.ucar.edu/surface/index.php?metarIds=KARB) and [its corresponding structured form](http://weather.rap.ucar.edu/surface/index.php?metarIds=KARB&std_trans=translated). Some streaming data is *streamed* because it's too large to be downloaded shotgun style and some is *streamed* because it's continually generated and serviced. This presents the potential problem of dealing with data streams that may be unlimited.	
#' 	
#' **Notes**:	
#' 	
#' * *Data sources*: Real or synthetic stream data can be used. Random simulation streams may be created by `rstream`. Real stream data may be piped from financial data providers, the WHO, World Bank, NCAR and other sources.	
#' * *Inference Techniques*: Many of the data interrogation techniques we  have seem can be employed for dynamic stream data, e.g., `factas`, for PCA, `rEMM` and `birch` for clustering, etc. Clustering and classification methods capable of processing data streams have been developed, e.g., *Very Fast Decision Trees* (VFDT), *time window-based Online Information Network* (OLIN), *On-demand Classification*, and the *APRIORI* streaming algorithm.	
#' * *Cloud distributed computing*: Hadoop2/HadoopStreaming, SPARK, Storm3/RStorm provide an environments to expand batch/script-based R tools to the Cloud.	
#' 	
#' ## The `stream` package	
#' 	
#' The `R` `stream` package provides data stream mining algorithms using `fpc`, `clue`, `cluster`, `clusterGeneration`, `MASS`, and `proxy` packages. In addition, the package `streamMOA` provides an `rJava` interface to the Java-based data stream clustering algorithms available in the *Massive Online Analysis* (MOA) framework for stream classification, regression and clustering.	
#' 	
#' If you need a deeper exposure to data streaming in R, we recommend you go over the [stream vignettes](https://cran.r-project.org/web/packages/stream/vignettes/stream.pdf).	
#' 	
#' ## Synthetic example - random Gaussian stream	
#' 	
#' This example shows the creation and loading of a *mixture of 5 random 2D Gaussians*, centers at (*x_coords*, *y_coords*) with paired correlations *rho_corr*, representing a simulated data stream.	
#' 	
#' * Generate the stream:	
#' 	
#' 	
# install.packages("stream")	
library("stream") 	
x_coords <- c(0.2,0.3, 0.5, 0.8, 0.9)	
y_coords <- c(0.8,0.3, 0.7, 0.1, 0.5)	
p_weight <- c(0.1, 0.9, 0.5, 0.4, 0.3) # A vector of probabilities that determines the likelihood of generated a data point from a particular cluster	
set.seed(12345)	
stream_5G <- DSD_Gaussians(k = 5, d = 2, mu=cbind(x_coords, y_coords), p=p_weight)	
#' 	
#' 	
#' * K-Means clustering	
#' 	
#' We will now try [k-means](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/12_kMeans_Clustering.html) and [density-based data stream clustering algorithm, D-Stream](https://pdfs.semanticscholar.org/1c6e/cca7bd2f03b55233159cb7c0095a14c4b4c3.pdf), where micro-clusters are formed by grid cells of size *gridsize* with density of a grid cell (Cm) is least 1.2 times the average cell density. The model is updated with the next 500 data points from the stream.	
#' 	
#' 	
dstream <- DSC_DStream(gridsize = .1, Cm = 1.2) 	
update(dstream,  stream_5G,  n  =  500)	
#' 	
#' 	
#' First, let's run the [k-means clustering](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/12_kMeans_Clustering.html) with $k=5$ clusters and plot the resulting micro and macro clusters	
#' 	
#' 	
kmc <- DSC_Kmeans(k = 5)	
recluster(kmc,   dstream)	
plot(kmc, stream_5G, type = "both", xlab="X-axis", ylab="Y-axis")	
#' 	
#' 	
#' In this clustering plot, *micro-clusters are shown as circles* and *macro-clusters are shown as crosses* and their sizes represent the corresponding cluster weight estimates.	
#' 	
#' Next try the density-based data stream clustering algorithm D-Stream. Prior to updating the model with the next 1,000 data points from the stream, we specify the grid cells as micro-clusters, grid cell size (gridsize=0.1), and a micro-cluster (Cm=1.2) that specifies the density of a grid cell as a multiple of the average cell density.	
#' 	
#' 	
dstream <- DSC_DStream(gridsize = 0.1, Cm = 1.2)	
update(dstream,  stream_5G,  n=1000)	
#' 	
#' 	
#' We can re-cluster the data using k-means with 5 clusters and plot the resulting *micro* and *macro* clusters.	
#' 	
#' 	
km_G5 <- DSC_Kmeans(k = 5) 	
recluster(km_G5,   dstream)	
plot(km_G5, stream_5G, type = "both")	
#' 	
#' 	
#' Note the subtle changes in the clustering results between `kmc` and `km_G5`.	
#' 	
#' ## Sources of Data Streams	
#' 	
#' ### Static structure streams	
#' 	
#' -	*DSD_BarsAndGaussians* generates two uniformly filled rectangular and two Gaussians clusters with different density.	
#' -	*DSD_Gaussians* generates randomly placed static clusters with random multivariate Gaussian distributions.	
#' -	*DSD_mlbenchData* provides streaming access to machine learning benchmark data sets found in the `mlbench` package.	
#' -	*DSD_mlbenchGenerator* interfaces the generators for artificial data sets defined in the mlbench package.	
#' -	*DSD_Target* generates a ball in circle data set.	
#' -	*DSD_UniformNoise* generates uniform noise in a d-dimensional (hyper) cube.	
#' 	
#' ### Concept drift streams	
#' 	
#' -	*DSD_Benchmark* provides a collection of simple benchmark problems including splitting and joining clusters, and changes in density or size, which can be used as a comprehensive  benchmark  set for algorithm comparison.	
#' -	*DSD_MG* is a generator to specify complex data streams with concept drift. The shape as well as the behavior of each cluster over time can be specified using keyframes.	
#' -	*DSD_RandomRBFGeneratorEvents* generates streams using radial base functions with noise. Clusters move, merge and split.	
#' 	
#' ### Real data streams	
#' 	
#' -	*DSD_Memory* provides a streaming interface to static, matrix-like data (e.g., a data frame, a matrix) in memory which represent a fixed portion of a data stream. Matrix-like objects also include large objects potentially stored on disk like `ff::ffdf`.	
#' -	*DSD_ReadCSV* reads data line by line in text format from a file or an open connection and makes it available in a streaming fashion. This way data that is larger than the available main memory can be processed.	
#' -	*DSD_ReadDB* provides an interface to an open result set from a SQL query to a relational database.	
#' 	
#' ## Printing, plotting and saving streams	
#' 	
#' For `DSD` objects, some basic stream functions include `print()`, `plot()` and `write_stream()`. These can save part of a data stream to disk. `DSD_Memory` and `DSD_ReadCSV` objects also include member functions like `reset_stream()` to reset the position in the stream to its beginning.	
#' 	
#' to request a new batch of data points from the stream we use `get_points()`.  This chooses a *random cluster* (based on the probability weights in `p\_weight`) and a point is drawn from the multivariate Gaussian distribution ($mean=mu, covariance\ matrix=\Sigma$) of that cluster. Below, we pull $n = 10$ new data points from the stream.	
#' 	
#' 	
new_p  <-  get_points(stream_5G,  n  =  10)	
new_p	
	
new_p  <-  get_points(stream_5G,  n  =  100,  class  =  TRUE)	
head(new_p, n = 20)	
	
plot(stream_5G,  n  =  700,  method  =  "pc")	
#' 	
#' 	
#' Note that if you add *noise* to your stream, e.g., `stream_Noise <- DSD_Gaussians(k = 5, d = 4, noise = .1, p = c(0.1, 0.5, 0.3, 0.9, 0.1))`, then the noise points won't be part of any clusters will have an `NA` class label.	
#' 	
#' ## Stream animation	
#' 	
#' Clusters can be animated over time by `animate_data()`. Use `reset_stream()` to start the animation at the beginning of the stream and note that this method is **not implemented** for streams of class `DSD_Gaussians`, `DSD_R`, `DSD_data.frame`, and `DSD`. We'll create a new `DSD_Benchmark` data stream.	
#' 	
#' 	
set.seed(12345)	
stream_Bench <- DSD_Benchmark(1) 	
stream_Bench	
#' 	
#' 	
#' 	
library("animation")	
reset_stream(stream_Bench)	
animate_data(stream_Bench, n=10000, horizon=100, xlim = c(0, 1), ylim = c(0, 1))	
#' 	
#' 	
#' This benchmark generator creates two 2D clusters moving in 2D. One moves from *top-left* to *bottom-right*, the other from *bottom-left* to *top-right*. Then they meet at the center of the domain, the 2 clusters overlap and then split again.	
#' 	
#' Concept drift in the stream can be depicted by requesting ($10$) times $300$ data points from the stream and animating the plot. Fast-forwarding the stream can be accomplished by requesting, but ignoring, ($2000$) points in between the ($10$) plots.	
#' 	
#' 	
for(i in 1:10) {	
   plot(stream_Bench, 300, xlim = c(0, 1), ylim = c(0, 1))	
   tmp <- get_points(stream_Bench, n = 2000)	
}	
	
reset_stream(stream_Bench)	
animate_data(stream_Bench, n=8000, horizon = 120,	xlim=c(0, 1), ylim=c(0, 1))	
	
# Animations can be saved as HTML or GIF	
#saveHTML(ani.replay(), htmlfile = "stream_Bench_Animation.html")	
#saveGIF(ani.replay())	
#' 	
#' 	
#' Streams can also be saved locally by `write_stream(stream_Bench,  "dataStreamSaved.csv",  n =  100,  sep=",")` and loaded back in `R` by `DSD_ReadCSV()`.	
#' 	
#' ## Case-Study: SOCR Knee Pain Data	
#' 	
#' These data represent the $X$ and $Y$ spatial knee-pain locations for over $8,000$ patients, along with *labels* about the knee $F$ront, $B$ack, $L$eft and $R$ight. Let's try to read the [SOCR Knee Pain Datasest](http://wiki.socr.umich.edu/index.php/SOCR_Data_KneePainData_041409) as a stream. 	
#' 	
#' 	
library("XML"); library("xml2"); library("rvest")	
	
wiki_url <- read_html("http://wiki.socr.umich.edu/index.php/SOCR_Data_KneePainData_041409")	
html_nodes(wiki_url, "#content")	
	
kneeRawData <- html_table(html_nodes(wiki_url, "table")[[2]])	
normalize<-function(x){	
  return((x-min(x))/(max(x)-min(x)))	
}	
kneeRawData_df <- as.data.frame(cbind(normalize(kneeRawData$x), normalize(kneeRawData$Y), as.factor(kneeRawData$View)))	
colnames(kneeRawData_df) <- c("X", "Y", "Label")	
# randomize the rows of the DF as RF, RB, LF and LB labels of classes are sequential	
set.seed(1234)	
kneeRawData_df <- kneeRawData_df[sample(nrow(kneeRawData_df)), ]	
summary(kneeRawData_df)	
# View(kneeRawData_df)	
#' 	
#' 	
#' We can use the `DSD::DSD_Memory` class to get a stream interface for matrix or data frame objects, like the Knee pain location dataset. The number of true clusters $k=4$ in this dataset.	
#' 	
#' 	
# use data.frame to create a stream (3rd column contains the label assignment)	
kneeDF <- data.frame(x=kneeRawData_df[,1], y=kneeRawData_df[,2], 	
  class=as.factor(kneeRawData_df[,3]))	
head(kneeDF)  	
	
streamKnee <- DSD_Memory(kneeDF[,c("x", "y")], class=kneeDF[,"class"], loop=T)	
streamKnee	
	
# Each time we get a point from *streamKnee*, the stream pointer moves to the next position (row) in the data.	
get_points(streamKnee, n=10)	
streamKnee	
	
# Stream pointer is in position 11 now	
	
# We can redirect the current position of the stream pointer by:	
reset_stream(streamKnee,  pos  =  200)	
get_points(streamKnee, n=10)	
streamKnee	
#' 	
#' 	
#' ## Data Stream clustering and classification (DSC)	
#' 	
#' Let's demonstrate clustering using `DSC_DStream`, which assigns points to cells in a grid. First, initialize the clustering, as an empty cluster and then use the `update()` function to implicitly alter the mutable `DSC` object.	
#' 	
#' 	
dsc_streamKnee <- DSC_DStream(gridsize = 0.1, Cm = 0.4, attraction=T) 	
dsc_streamKnee	
# stream::update	
reset_stream(streamKnee,  pos  =  1)	
update(dsc_streamKnee, streamKnee,  n  =  500) 	
dsc_streamKnee	
head(get_centers(dsc_streamKnee))	
	
plot(dsc_streamKnee,  streamKnee, xlim=c(0,1), ylim=c(0,1))	
# plot(dsc_streamKnee,  streamKnee, grid = TRUE)	
# Micro-clusters are plotted in red on top of gray stream data points	
# The size of the micro-clusters indicates their weight - it's proportional to the number of data points represented by each micro-cluster. 	
# Micro-clusters are shown as dense grid cells (density is coded with gray values).	
#' 	
#' 	
#' Next, we can use K-means clustering.	
#' 	
#' 	
kMeans_Knee <- DSC_Kmeans(k = 5) # choose 4-5 clusters as we have 4 knee labels	
recluster(kMeans_Knee, dsc_streamKnee)	
plot(kMeans_Knee, streamKnee, type = "both")	
	
animate_data(streamKnee, n=1000, horizon=100, xlim = c(0, 1), ylim = c(0, 1))	
	
# purity <- animate_cluster(kMeans_Knee, streamKnee, n=2500, type="both", xlim=c(0,1), ylim=c(-,1), evaluationMeasure="purity", horizon=10)	
	
animate_cluster(kMeans_Knee, streamKnee,  horizon = 100,  n = 5000, measure = "purity", plot.args = list(xlim = c(0, 1), ylim = c(0, 1)))	
#' 	
#' 	
#' ##	Evaluation of data stream clustering	
#' 	
#' 	
# Synthetic Gaussian example	
# stream <- DSD_Gaussians(k = 3, d = 2, noise = .05) 	
# dstream <- DSC_DStream(gridsize = .1)	
# update(dstream,  stream,  n  =  2000)	
# evaluate(dstream,  stream,  n  =  100)	
	
evaluate(dsc_streamKnee, streamKnee, measure = c("crand", "SSQ", "silhouette"), n = 100, type = c("auto", "micro", "macro"), assign = "micro", assignmentMethod = c("auto", "model", "nn"), noise = c("class", "exclude"))	
	
clusterEval  <- evaluate_cluster(dsc_streamKnee, streamKnee, measure = c("numMicroClusters", "purity"), n = 5000, horizon = 100) 	
head(clusterEval)	
	
plot(clusterEval[ , "points"], clusterEval[ , "purity"], type = "l", ylab  =  "Avg Purity",  xlab  =  "Points")	
	
animate_cluster(dsc_streamKnee, streamKnee,  horizon  =  100,  n  =  5000, measure = "purity", plot.args = list(xlim = c(0, 1), ylim = c(0, 1)))	
#' 	
#' 	
#' The `dsc_streamKnee` is the evaluated clustering, where $n$ is the data points are taken from `streamKnee`. The evaluation `measure` can be specified as a vector of character strings. Points are assigned to clusters in `dsc_streamKnee` using `get_assignment()` and can be used to assess the quality of the classification. By default, points are assigned to *micro-clusters*, or can be assigned to *macro-cluster* centers by `assign = "macro"`. ALso, new points can be assigned to clusters by the rule used in the clustering algorithm by `assignmentMethod = "model"` or using nearest-neighbor assignment (`nn`). 	
#' 	
#' # Optimization and improving the computational performance	
#' 	
#' Here and in previous chapters, e.g., [Chapter 14](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/14_ImprovingModelPerformance.html), we notice that R may sometimes be slow and memory-inefficient. These problems may be severe, especially for datasets with millions of records or when using complex functions. There are [packages for processing large datasets and memory optimization](https://rstudio-pubs-static.s3.amazonaws.com/72295_692737b667614d369bd87cb0f51c9a4b.html) -- `bigmemory`, `biganalytics`, `bigtabulate`, etc.	
#' 	
#' ## Generalizing tabular data structures with `dplyr`	
#' 	
#' We have also seen long execution times when running processes that ingest, store or manipulate huge `data.frame` objects. The `dplyr` package, created by Hadley Wickham and Romain Francoi, provides a faster route to manage such large datasets in R. It creates an object called `tbl`, similar to `data.frame`, which has an in-memory column-like structure. R reads these objects a lot faster than data frames.	
#' 	
#' To make a `tbl` object we can either convert an existing data frame to `tbl` or connect to an external database. Converting from data frame to `tbl` is quite easy. All we need to do is call the function `as.tbl()`.	
#' 	
#' 	
#install.packages("dplyr")	
library(dplyr)	
nof1_tbl<-as.tbl(nof1)	
nof1_tbl	
#' 	
#' 	
#' This looks like a normal data frame. If you are using R Studio by viewing the `nof1_tbl` you can see the same output as `nof1`.	
#' 	
#' ## Making data frames faster with data.table	
#' 	
#' Similar to `tbl`, the `data.table` package provides another alternative to data frame object representation. `data.table` objects are processed in R much faster compared to standard data frames. Also, all of the functions that can accept data frame could be applied to `data.table` objects as well. The function `fread()` is able to read a local CSV file directly into a `data.table`.	
#' 	
#' 	
#install.packages("data.table")	
library(data.table)	
nof1<-fread("C:/Users/Dinov/Desktop/02_Nof1_Data.csv")	
#' 	
#' 	
#' Another amazing property of `data.table` is that we can use subscripts to access a specific location in the dataset just like `dataset[row, column]`. It also allows the selection of rows with Boolean expression and direct application of functions to those selected rows. Note that column names can be used to call the specific column in `data.table`, whereas with data frames, we have to use the `dataset$columnName` syntax).	
#' 	
#' 	
nof1[ID==1, mean(PhyAct)]	
#' 	
#' 	
#' This useful functionality can also help us run complex operations with only a few lines of code. One of the drawbacks of using `data.table` objects is that they are still limited by the available system memory.	
#' 	
#' ## Creating disk-based data frames with `ff`	
#' 	
#' The `ff` (fast-files) package allows us to overcome the RAM limitations of finite system memory. For example, it helps with operating datasets with billion rows. `ff` creates objects in `ffdf` formats, which is like a map that points to a location of the data on a disk. However, this makes `ffdf` objects inapplicable for most R functions. The only way to address this problem is to break the huge dataset into small chunks. After processing a batch of these small chunks, we have to combine the results to reconstruct the complete output. This strategy is relevant in parallel computing, which will be discussed in detail in the next section. First, let's download one of the [large datasets in our datasets archive, UQ_VitalSignsData_Case04.csv](https://umich.instructure.com/files/366335/download?download_frd=1).	
#' 	
#' 	
# install.packages("ff")	
library(ff)	
# vitalsigns<-read.csv.ffdf(file="UQ_VitalSignsData_Case04.csv", header=T)	
vitalsigns<-read.csv.ffdf(file="https://umich.instructure.com/files/366335/download?download_frd=1", header=T)	
#' 	
#' 	
#' As mentioned earlier, we cannot apply functions directly on this object. 	
#' 	
#' 	
mean(vitalsigns$Pulse)	
#' 	
#' 	
#' For basic calculations in datasets we can download another package `ffbase`. This allows operations on `ffdf` objects using simple tasks like: mathematical operations, query functions, summary statistics and bigger regression models using packages like `biglm`, which will be mentioned later in this chapter.	
#' 	
#' 	
# install.packages("ffbase")	
library(ffbase)	
mean(vitalsigns$Pulse)	
#' 	
#' 	
#' ## Using massive matrices with `bigmemory`	
#' 	
#' The previously introduced packages include alternatives to `data.frames`. FOr instance, the `bigmemory` package creates alternative objects to 2D matrices (second-order tensors). It can store huge datasets and can be divided into small chunks that can be converted to data frames. However, we cannot directly apply machine learning methods on this types of objects. More [detailed information about the `bigmemory` package is available online](http://www.bigmemory.org).	
#' 	
#' # Parallel computing	
#' 	
#' In previous chapters, we saw various machine learning techniques applied as serial computing tasks. The traditional protocol involves: First, applying *function 1* to our raw data. Then, using the output from *function 1* as an input to *function 2*. This process in iterated for a series of functions. Finally, we have the terminal output generated by the last function. This serial or linear computing method is straight forward but time consuming and perhaps sub-optimal.	
#' 	
#' Now we introduce a more efficient way of computing - *parallel computing*, which provides a mechanism to deal with different tasks at the same time and combine the outputs for all of processes to get the final answer faster. However, parallel algorithms may require special conditions and cannot be applied to all problems. If two tasks have to be run in a specific order, this problem cannot be parallelized.	
#' 	
#' ## Measuring execution time	
#' 	
#' To measure how much time can be saved for different methods, we can use function `system.time()`.	
#' 	
#' 	
system.time(mean(vitalsigns$Pulse))	
#' 	
#' 	
#' This means calculating the mean of `Pulse` column in the `vitalsigns` dataset takes 0.001 seconds. These values will vary between computers, operating systems, and states of operations.	
#' 	
#' ## Parallel processing with multiple cores	
#' 	
#' We will introduce two packages for parallel computing `multicore` and `snow` (their core components are included in the package `parallel`). They both have a different way of multitasking. However, to run these packages, you need to have a relatively modern multicore computer. Let's check how many cores your computer has. This function `parallel::detectCores()` provides this functionality. `parallel` is a base package, so there is no need to install it prior to using it.	
#' 	
#' 	
library(parallel)	
detectCores()	
#' 	
#' 	
#' So, there are eight (8) cores in my computer. I will be able to run up to 6-8 parallel jobs on this computer.	
#' 	
#' The `multicore` package simply uses the multitasking capabilities of the *kernel*, the computer's operating system, to "fork" additional R sessions that share the same memory. Imagine that we open several R sessions in parallel and let each of them does part of the work. Now, let's examine how this can save time when running complex protocols or dealing with large datasets. To start with, we can use the `mclapply()` function, which is similar to `lapply()`, which applies functions to a vector and returns a vector of lists. Instead of applying functions to vectors `mcapply()` divides the complete computational task and delegates portions of it to each available core. We will apply a simple, yet time consuming, task-generating random numbers for demonstrating this procedure. Also, we can use the `system.time()` to track the time differences.	
#' 	
#' 	
set.seed(123)	
system.time(c1<-rnorm(10000000))	
	
# Note the multi core calls may not work on Windows, but will work on Linux/Mac.	
#This shows a 2-core and 4-vore invocations	
# system.time(c2<-unlist(mclapply(1:2, function(x){rnorm(5000000)}, mc.cores = 2)))	
# system.time(c4<-unlist(mclapply(1:4, function(x){rnorm(2500000)}, mc.cores = 4)))     	
	
# And here is a Windows (single core invocation)	
system.time(c2<-unlist(mclapply(1:2, function(x){rnorm(5000000)}, mc.cores = 1)))	
#' 	
#' 	
#' The `unlist()` is used at the end to combine results from different cores into a single vector. Each line of code creates 10,000,000 random numbers. `c1` is a regular R command, which used longest time. `c2` used two cores to finish the task (each core handle 5,000,000 numbers) and used less time than the first one. `c4` used all four cores to finish the task and successfully reduce the time again. We can see that when we use more cores the time is significantly reduced.	
#' 	
#' The `snow` package allows parallel computing on multicore multiprocessor machines or a network of multiple machines. It might be more difficult to use but it's also certainly more flexible. First we can set how many cores we want to use via `makeCluster()` function.	
#' 	
#' 	
# install.packages("snow")	
library(snow)	
cl<-makeCluster(2)	
#' 	
#' 	
#' This call might cause your computer to pop up a message warning about access though the firewall. To do the same task we can use `parLapply()` function in the `snow` package. Note that we have to call the object we created with the previous `makeCluster()` function.	
#' 	
#' 	
system.time(c2<-unlist(parLapply(cl, c(5000000, 5000000), function(x) {rnorm(x)})))	
#' 	
#' 	
#' While using `parLapply()`, we have to specify the matrix and the function that will be applied to this matrix. Remember to stop the cluster we made after completing the task, to release back the system resources.	
#' 	
#' 	
stopCluster(cl)	
#' 	
#' 	
#' ## Parallelization using `foreach` and `doParallel`	
#' 	
#' The `foreach` package provides another option of parallel computing. It relies on a loop-like process basically applying a specified function for each item in the set, which again is somewhat similar to `apply()`, `lapply()` and other regular functions. The interesting thing is that these loops can be computed in parallel saving substantial amounts of time. The `foreach` package alone cannot provide parallel computing. We have to combine it with other packages like `doParallel`. Let's reexamine the task of creating a vector of 10,000,000 random numbers. First, register the 4 compute cores using `registerDoParallel()`.	
#' 	
#' 	
# install.packages("doParallel")	
library(doParallel)	
cl<-makeCluster(4)	
registerDoParallel(cl)	
#' 	
#' 	
#' Then we can examine the time saving `foreach` command. 	
#' 	
#' 	
#install.packages("foreach")	
library(foreach)	
system.time(c4<-foreach(i=1:4, .combine = 'c')	
            %dopar% rnorm(2500000))	
#' 	
#' 	
#' Here we used four items (each item runs on a separate core), `.combine=c` allows `foreach` to combine the results with the parameter `c()` generating the aggregate result vector.	
#' 	
#' Also, don't forget to close the `doParallel` by registering the sequential backend.	
#' 	
#' 	
unregister<-registerDoSEQ()	
#' 	
#' 	
#' ## GPU computing	
#' 	
#' Modern computers have graphics cards, GPU (Graphics Processing Unit), that consists of thousands of cores, however they are very specialized, unlike the standard CPU chip. If we can use this feature for parallel computing, we may reach amazing performance improvements, at the cost of complicating the processing algorithms and increasing the constraints on the data format. Specific disadvantages of GPU computing include relying on a proprietary manufacturer (e.g., NVidia) frameworks and Complete Unified Device Architecture (CUDA) programming language. CUDA allows programming of GPU instructions into a common computing language. This [paper provides one example of using GPU computation to improve significantly the performance of advanced neuroimaging and brain mapping processing of multidimensional data](http://dx.doi.org/10.1016/j.cmpb.2010.10.013).	
#' 	
#' The R package `gputools` is created for parallel computing using NVidia CUDA. Detailed [GPU computing in R information is available online](https://cran.r-project.org/web/packages/gputools/gputools.pdf).	
#' 	
#' # Deploying optimized learning algorithms	
#' 	
#' As we mentioned earlier, some tasks can be parallelized easier than others. In real word situations, we can pick the algorithms that lend themselves well to parallelization. Some of the R packages that allow parallel computing using ML algorithms are listed below.	
#' 	
#' ## Building bigger regression models with `biglm`	
#' 	
#' `biglm` allows training regression models with data from SQL databases or large data chunks obtained from the `ff` package. The output is similar to teh standard `lm()` function that builds linear models. However, `biglm` operates efficiently on massive datasets.	
#' 	
#' ## Growing bigger and faster random forests with `bigrf`	
#' 	
#' The `bigrf` package can be used to train random forests combining the `foreach` and `doParallel` packages. In [Chapter 14](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/14_ImprovingModelPerformance.html), we presented random forests as machine learners ensembling multiple tree learners. With parallel computing, we can split the task of creating thousands of trees into smaller tasks that can be outsourced to each available compute core. We only need to combine the results at the end. Then, we will obtain the exact same output in a relatively shorter amount of time.	
#' 	
#' ## Training and evaluation models in parallel with `caret`	
#' 	
#' Combining the `caret` package with `foreach` and we can obtain a powerful method to deal with time-consuming tasks like building a random forest learner. Utilizing the same example we presented in [Chapter 14(http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/14_ImprovingModelPerformance.html), we can see the time difference of utilizing the `foreach` package.	
#' 	
#' 	
qol<-read.csv("https://umich.instructure.com/files/481332/download?download_frd=1")	
qol<-qol[!qol$CHARLSONSCORE==-9 , -c(1, 2)]	
qol$CHARLSONSCORE<-as.factor(qol$CHARLSONSCORE)	
library(caret)	
ctrl<-trainControl(method="cv", number=10)	
grid_rf<-expand.grid(.mtry=c(2, 4, 8, 16))	
#' 	
#' 	
#library(caret)	
system.time(m_rf <- train(CHARLSONSCORE ~ ., data = qol, method = "rf", 	
metric = "Kappa", trControl = ctrl, tuneGrid = grid_rf))	
#' 	
#' 	
#' It took more than a minute to finish this task in standard execution model purely relying on the regular `caret` function. Below, this same model training completes much faster using parallelization (less than half the time) compared to the standard call above.	
#' 	
#' 	
set.seed(123)	
cl<-makeCluster(4)	
registerDoParallel(cl)	
getDoParWorkers()	
system.time(m_rf <- train(CHARLSONSCORE ~ ., data = qol, method = "rf", 	
metric = "Kappa", trControl = ctrl, tuneGrid = grid_rf))	
unregister<-registerDoSEQ()	
#' 	
#' 	
#' # Practice problem	
#' 	
#' Try to analyze [the co-appearance network in the novel "Les Miserablese"](https://umich.instructure.com/files/330389/download?download_frd=1). The data contains the weighted network of co-appearances of characters in Victor Hugo's novel "Les Miserables". Nodes represent characters as indicated by the labels and edges connect any pair of characters that appear in the same chapter of the book.  The values on the edges are the number of such co-appearances. 	
#' 	
#' > miserablese<-read.table("https://umich.instructure.com/files/330389/download?download_frd=1", sep="", header=F)	
#' > head(miserablese)	
#' 	
#' Also, try to interrogate some of the larger datasets we have using alternative parallel computing and big data analytics.