Add files via upload

code/chapter2/BML.R

@@ -0,0 +1,62 @@
+library(R6)
+
+BML = R6Class(
+  "BML",
+  public = list(
+    # alpha is the parameter of the uniform distribution to control particle distribution's density
+    # m*n is the dimension of the lattice
+    alpha = NULL,
+    m = NULL,
+    n = NULL,
+    lattice = NULL,
+    initialize = function(alpha, m, n) {
+      self$alpha = alpha
+      self$m = m
+      self$n = n
+      self$initialize_lattice()
+    },
+    initialize_lattice = function() {
+      # 0 -> empty site
+      # 1 -> blue particle
+      # 2 -> red particle
+      u = runif(self$m * self$n)
+      # the usage of L is to make sure the elements in particles are of type integer;
+      # otherwise they would be created as double
+      particles = rep(0L, self$m * self$n)
+      # https://en.wikipedia.org/wiki/Inverse_transform_sampling
+      particles[(u > self$alpha) & (u <= (self$alpha + 1.0) / 2)] = 1L
+      particles[u > (self$alpha + 1.0) / 2] = 2L
+      self$lattice = array(particles, c(self$m, self$n))
+    },
+    odd_step = function() {
+      blue.index = which(self$lattice == 1L, arr.ind = TRUE)
+      # make a copy of the index
+      blue.up.index = blue.index
+      # blue particles move 1 site up
+      blue.up.index[, 1] = blue.index[, 1] - 1L
+      # periodic boundary condition
+      blue.up.index[blue.up.index[, 1] == 0L, 1] = self$m
+      # find which moves are feasible
+      blue.movable = self$lattice[blue.up.index] == 0L
+      # move blue particles one site up
+      # drop=FALSE prevents the 2D array degenerates to 1D array
+      self$lattice[blue.up.index[blue.movable, , drop = FALSE]] = 1L
+      self$lattice[blue.index[blue.movable, , drop = FALSE]] = 0L
+    },
+    even_step = function() {
+      red.index = which(self$lattice == 2L, arr.ind = TRUE)
+      # make a copy of the index
+      red.right.index = red.index
+      # red particles move 1 site right
+      red.right.index[, 2] = red.index[, 2] + 1L
+      # periodic boundary condition
+      red.right.index[red.right.index[, 2] == (self$n + 1L), 2] = 1
+      # find which moves are feasible
+      red.movable = self$lattice[red.right.index] == 0L
+      # move red particles one site right
+      self$lattice[red.right.index[red.movable, , drop = FALSE]] = 2L
+      self$lattice[red.index[red.movable, , drop = FALSE]] = 0L
+    }
+  )
+)
+

code/chapter2/BML.py

@@ -0,0 +1,37 @@
+import numpy as np
+
+class BML:
+    def __init__(self, alpha, m, n):
+        self.alpha = alpha
+        self.shape = (m, n)
+        self.initialize_lattice()
+
+    def initialize_lattice(self):
+        u = np.random.uniform(0.0, 1.0, self.shape)
+        # instead of using default list, we use np.array to create the lattice
+        self.lattice = np.zeros_like(u, dtype=int)
+        # the parentheses below can't be ignored
+        self.lattice[(u > self.alpha) & (u <= (1.0+self.alpha)/2.0)] = 1
+        self.lattice[u > (self.alpha+1.0)/2.0] = 2
+
+    def odd_step(self):
+        # please note that np.where returns a tuple which is immutable
+        blue_index = np.where(self.lattice == 1)
+        blue_index_i = blue_index[0] - 1
+        blue_index_i[blue_index_i < 0] = self.shape[0]-1
+        blue_movable = self.lattice[(blue_index_i, blue_index[1])] == 0
+        self.lattice[(blue_index_i[blue_movable],
+                      blue_index[1][blue_movable])] = 1
+        self.lattice[(blue_index[0][blue_movable],
+                      blue_index[1][blue_movable])] = 0
+
+    def even_step(self):
+        red_index = np.where(self.lattice == 2)
+        red_index_j = red_index[1] + 1
+        red_index_j[red_index_j == self.shape[1]] = 0
+        red_movable = self.lattice[(red_index[0], red_index_j)] == 0
+        self.lattice[(red_index[0][red_movable],
+                      red_index_j[red_movable])] = 2
+        self.lattice[(red_index[0][red_movable],
+                      red_index[1][red_movable])] = 0
+

code/chapter2/benchmark.R

@@ -0,0 +1,15 @@
+library(microbenchmark)
+source('binary_search.R')
+source('find_pos.R')
+
+v=1:10000
+
+# call the find_pos 100 times; 
+# each time we randomly select an integer as the target value 
+
+# for-loop solution
+set.seed(2019)
+print(microbenchmark(find_pos(v,sample(10000,1)),times=1000))
+# binary-search solution
+set.seed(2019)
+print(microbenchmark(binary_search(v,sample(10000,1)),times=1000))

code/chapter2/benchmark.py

@@ -0,0 +1,21 @@
+from binary_search import binary_search
+from find_pos import find_pos
+import timeit
+import random
+
+v=list(range(1,10001))
+
+def test_for_loop(n):
+  random.seed(2019)
+  for _ in range(n):
+    find_pos(v,random.randint(1,10000))
+
+def test_bs(n):
+  random.seed(2019)
+  for _ in range(n):
+    binary_search(v,random.randint(1,10000))
+
+# for-loop solution
+print(timeit.timeit('test_for_loop(1000)',setup='from __main__ import test_for_loop',number=1))
+# binary_search solution
+print(timeit.timeit('test_bs(1000)',setup='from __main__ import test_bs',number=1))

code/chapter2/binary_search.R

@@ -0,0 +1,17 @@
+binary_search=function(v,x){
+  if (x>v[length(v)]){return(NULL)}
+  start = 1
+  end = length(v)
+  while (start<end){
+    mid = (start+end) %/% 2 # %/% is the floor division operator
+    if (v[mid]>=x){
+      end = mid
+    }else{
+      start = mid+1
+    }
+  }
+  return(start)
+}
+
+
+

code/chapter2/binary_search.py

@@ -0,0 +1,10 @@
+def binary_search(v,x):
+  if x>v[-1]: return
+  start,end = 0,len(v)-1
+  while start<end:
+    mid = (start+end)//2
+    if v[mid]>=x:
+      end = mid
+    else:
+      start = mid+1
+  return start

code/chapter2/binary_search_buggy.R

@@ -0,0 +1,20 @@
+binary_search_buggy=function(v,x){
+  start = 1
+  end = length(v)
+  while (start<end){
+    mid = (start+end) %/% 2 # %/% is the floor division operator
+    if (v[mid]>=x){
+      end = mid
+    }else{
+      start = mid+1
+    }
+  }
+  return(start)
+}
+v=c(1,2,5,10)
+print(binary_search_buggy(v,-1))
+print(binary_search_buggy(v,5))
+print(binary_search_buggy(v,11))
+
+
+

code/chapter2/binary_search_buggy.py

@@ -0,0 +1,15 @@
+def binary_search_buggy(v,x):
+  start,end = 0,len(v)-1
+  while start<end:
+    mid = (start+end)//2  # // is the floor division operator
+    if v[mid]>=x:
+      end = mid
+    else:
+      start = mid+1
+  return start
+
+v=[1,2,5,10]
+print(binary_search_buggy(v,-1))
+print(binary_search_buggy(v,5))
+print(binary_search_buggy(v,11))
+

code/chapter2/binary_search_buggy_debug.R

@@ -0,0 +1,19 @@
+binary_search_buggy=function(v,x){
+  browser()
+  start = 1
+  end = length(v)
+  while (start<end){
+    mid = (start+end) %/% 2 # %/% is the floor division operator
+    if (v[mid]>=x){
+      end = mid
+    }else{
+      start = mid+1
+    }
+  }
+  return(start)
+}
+v=c(1,2,5,10)
+print(binary_search_buggy(v,11))
+
+
+

code/chapter2/binary_search_buggy_debug.py

@@ -0,0 +1,15 @@
+from pdb import set_trace
+def binary_search_buggy(v,x):
+  set_trace()
+  start,end = 0,len(v)-1
+  while start<end:
+    mid = (start+end)//2
+    if v[mid]>=x:
+      end = mid
+    else:
+      start = mid+1
+  return start
+
+v=[1,2,5,10]
+print(binary_search_buggy(v,11))
+

code/chapter2/find_pos.R

@@ -0,0 +1,11 @@
+find_pos=function(v,x){
+  for (i in 1:length(v)){
+    if (v[i]>=x){
+      return(i)
+    }
+  }
+}
+v=c(1,2,5,10)
+print(find_pos(v,-1))
+print(find_pos(v,4))
+print(find_pos(v,11))

code/chapter2/find_pos.py

@@ -0,0 +1,9 @@
+def find_pos(v,x):
+  for i in range(len(v)):
+    if v[i]>=x:
+      return i
+v=[1,2,5,10]
+print(find_pos(v,-1))
+print(find_pos(v,4))
+print(find_pos(v,11))
+

code/chapter2/hello_world.R

@@ -0,0 +1 @@
+print("Hello World!")

code/chapter2/hello_world.py

@@ -0,0 +1 @@
+print("Hello World!")

code/chapter2/multi_assignment.R

@@ -0,0 +1,26 @@
+`%=%` = function(left, right) {
+  # we require the RHS to be a list strictly
+  stopifnot(is.list(right))
+  # dest_env is the desitination environment enclosing the variables on LHS
+  dest_env = parent.env(environment())
+  left = substitute(left)
+
+  recursive_assign = function(left, right, dest_env) {
+    if (length(left) == 1) {
+      assign(x = deparse(left),
+             value = right,
+             envir = dest_env)
+      return()
+    }
+    if (length(left) != length(right) + 1) {
+      stop("LHS and RHS must have the same shapes")
+    }
+
+    for (i in 2:length(left)) {
+      recursive_assign(left[[i]], right[[i - 1]],  dest_env)
+    }
+  }
+
+  recursive_assign(left, right, dest_env)
+}
+

code/chapter2/vectorization_1.R

@@ -0,0 +1,17 @@
+library(microbenchmark)
+
+# generate n standard normal r.v
+rnorm_loop = function(n){
+x=rep(0,n)
+for (i in 1:n) {x[i]=rnorm(1)}
+}
+
+rnorm_vec = function(n){
+x=rnorm(n)
+}
+
+n=100
+# for loop
+print(microbenchmark(rnorm_loop(n),times=1000))
+# vectorize
+print(microbenchmark(rnorm_vec(n),times=1000))

code/chapter2/vectorization_1.py

@@ -0,0 +1,19 @@
+import timeit
+import numpy as np
+
+def rnorm_for_loop(n):
+  x=[0]*n # create a list with n 0s
+  np.random.seed(2019)
+  for _ in range(n):
+    np.random.normal(0,1,1)
+
+def rnorm_vec(n):
+  np.random.seed(2019)
+  x = np.random.normal(0,1,n)
+
+print("for loop")
+print(f'{timeit.timeit("rnorm_for_loop(100)",setup="from __main__ import rnorm_for_loop",number=1000):.6f} seconds')
+print("vectorized")
+print(f'{timeit.timeit("rnorm_vec(100)",setup="from __main__ import rnorm_vec",number=1000):.6f} seconds')
+
+

code/chapter3/gd_plot.py

@@ -0,0 +1,14 @@
+import numpy as np
+from mpl_toolkits.mplot3d import axes3d
+import matplotlib.pyplot as plt
+from matplotlib.patches import FancyArrowPatch
+
+fig = plt.figure()
+ax = fig.add_subplot(111, projection="3d")
+X, Y = np.mgrid[-1:1:30j, -1:1:30j]
+Z = X**2+Y**2 + 1
+ax.plot_surface(X, Y, Z, cmap="Greens_r", lw=0, rstride=2, cstride=2)
+ax.contour(X, Y, Z, 10, lw=3, cmap="Greens_r", linestyles="solid", offset=1)
+ax.axes.xaxis.set_ticklabels([])
+ax.axes.yaxis.set_ticklabels([])
+plt.show()

code/chapter3/linear_regression.R

@@ -0,0 +1,20 @@
+
+library(R6)
+LR = R6Class(
+  "LR",
+  public = list(
+    coef = NULL,
+    initialize = function() {
+
+    },
+    fit = function(x, y) {
+      self$qr_solver(cbind(1, x), y)
+    },
+    qr_solver = function(x, y) {
+      self$coef = qr.coef(qr(x), y)
+    },
+    predict = function(new_x) {
+      cbind(1, new_x) %*% self$coef
+    }
+  )
+)

code/chapter3/linear_regression.py

@@ -0,0 +1,17 @@
+import numpy as np
+
+
+class LR:
+    def __init__(self):
+        self.coef = None
+
+    def qr_solver(self, x, y):
+        q, r = np.linalg.qr(x)
+        p = np.dot(q.T, y)
+        return np.dot(np.linalg.inv(r), p)
+
+    def fit(self, x, y):
+        self.coef = self.qr_solver(np.hstack((np.ones((x.shape[0], 1)), x)), y)
+
+    def predict(self, x):
+        return np.dot(np.hstack((np.ones((x.shape[0], 1)), x)), self.coef)