Manually Deploy Cloud Architecture

This is a step by step tutorial of how to setup LAMP stack on two virtual machines on OpenStack cloud infrastructure. After configuration, we want to deploy the <a href="https://github.com/Anirban2404/phpMySQLapp.git">phpMySQLapp</a> on our servers. The two virtuall machines will be functioned as web server and database server respectively.

Creation and configuration of virtual machines on open stack can be found <a href="https://www.openstack.org/software/start/">here</a>. We will be using Ubuntu 16.04 as the OS of both virtuall machines. Note, we need to associate a floating IP address to our public web server.

Step 1: Install LAMP stack on Web Server

Use ssh to login into the Web server. We will first install the LAMP stack manually. LAMP stack includes Linux, Apache, MySQL and PHP. It's considered by many as the platform of choice for development and deployment of web applications.

Note: you need to have a non-root user with sudo privileges to continue.

To install Apache, use the following command:

1
2

$ sudo apt-get update
$ sudo apt-get install apache2

Then for MySQL, the command is very similar:

1

$ sudo apt-get install mysql-server

For PHP, we need some additional metapackages:

1

$ sudo apt-get install php libapache2-mod-php php-mcrypt php-mysql

Your server should restart the apache server; if not, use the following command:

1

$ sudo /etc/init.d/apache2 restart

You should now have the Ubuntu LAMP server works. You can check php works by running any php files in /var/www/ directory.

Step 2: Configure Apache on Web Server

Before we test our apache, the first thing we need to do is to modify our firework to allow outside access to the default web ports. For a default setting, you should have a UFW firewall configured to restrict access to your server.

During installation, Apache registers itself with UFW to provide a few application profiles. We can use these profiles to simplify the process of enabling or disabling access to Apache through our firewall.

To check the ufw application profiles, use

1

$ sudo ufw app list

You should get an output like the following:

1
2
3
4
5

Available applications:
  Apache
  Apache Full
  Apache Secure
  OpenSSH

For our purposes, we will allow incoming traffic for the Apache Full profile by typing:

1

$ sudo ufw allow 'Apache Full'

You can verify the change by:

1

$ sudo ufw status

The output should look like:

1
2
3
4
5
6
7
8

Status: active

To                         Action      From
--                         ------      ----
OpenSSH                    ALLOW       Anywhere                  
Apache Full                ALLOW       Anywhere                  
OpenSSH (v6)               ALLOW       Anywhere (v6)             
Apache Full (v6)           ALLOW       Anywhere (v6)

Now, hopefully you will have your Apache working.

Step 3 Get the App repository and Configure PHP

Use git clone to clone the repository into your front end server (Web Server).

1

$ git clone https://github.com/Anirban2404/phpMySQLapp.git

Move the directory to under /var/www/html/ to enable php script. Use tree can clearly see the structure of the directory.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

ubuntu@mywebserver:/var/www/html$ tree
.
├── admin_area
│   ├── insertbook.php
│   ├── insert_books.php
│   ├── insertmovie.php
│   └── insert_movies.php
├── books
│   ├── functions
│   │   ├── fetch.php
│   │   ├── functions.php
│   │   └── getbook.php
│   ├── home.php
│   ├── images
│   │   └── background_image.jpg
│   └── includes
│       └── bookDatabase.php
├── homePage.JPG
├── index.php
├── index.html
├── movies
│   ├── functions
│   │   ├── fetch.php
│   │   ├── functions.php
│   │   └── getmovie.php
│   ├── home.php
│   ├── images
│   │   └── background_image.jpg
│   └── includes
│       └── movieDatabase.php
├── mySqlDB
│   ├── bookDB.sql
│   └── movieDB.sql
├── README.md
└── siteImages
    ├── books.jpg
    └── movies.jpg

Note that index.html is the default page in the Apache server. You can either delete that file or configure PHP so that it will look into index.php to setup the front page.

To configure PHP, use vim or any texteditor like nano to open /etc/apache2/mods-enabled/dir.conf as:

1

$ sudo vim /etc/apache2/mods-enabled/dir.conf

Change the order of the directory index into:

1

DirectoryIndex index.php index.html index.cgi index.pl index.xhtml index.htm

Now php will look for index.php before index.html.

Step 4 Install MySQL and Configure it on Database Server

SSH into our database server. As in the front end, we use apt-get to install MySQL.

1

$ sudo apt-get install mysql-server

When installing MySQL, it will prompt for password. Please remember it and we will use it later.

We also need to clone the repository of the Web App here. We need to import the database in the repository to our MySQL. So we do git clone again:

1

$ git clone https://github.com/Anirban2404/phpMySQLapp.git

Following the instruction on the repository README, we import the database into our MySQL database by

1

$ mysql -u <username> -p <databasename> < <filename.sql>

Here, default username is root. Database names can be anything you like, but we are going to use them later. I chose 'bookstore' for my book database and 'movieDb' for my movie database. The files you need to import are in the 'mySqlDB' folder, as shown in the above tree. Argument -u is for username, -p means using password. This command will prompt you to enter your password for the database.

Now that we have import our database, we need to configure MySQL. First, use whatever texteditor you prefer to open the configuration file of MySQL. The file location may vary depend on different package. For this demonstration, the command is:

1

$ sudo vim /etc/mysql/mysql.conf.d/mysqld.cnf

In the configuration file, we need to bind our IP address to our database. Scroll down to find a line that looks like:

1

bind-address		= xxx.xxx.xx.xx

Replace the IP address with our database server's IP address. If there is a line starts with skip-networking, comment it out:

1

#skip-networking

Now save and exit the file. Restart the MySQL service:

1

$ service mysql restart

The final step is to grant the privilege of our user from other IP address. First Enter MySQL database by

1

$ mysql -u root -p

Replace root with your username if needed. In MySQL, type the following command:

1

mysql> GRANT ALL PRIVILEGES ON *.* TO 'USERNAME'@'%' IDENTIFIED BY 'PASSWORD' WITH GRANT OPTION;

Here 'USERNAME' should be your username, 'PASSWORD' should be your password. This command will grant privilege to user from any IP address. You can check that privileges are granted by:

1

SELECT * from information_schema.user_privileges;

You should see your user at the bottom of the chart. Now your database server is ready for use.

Step 5 Final Configuration on Web App

A final modification is on our Web server. SSH into it as usual. Use text editor to open /var/www/html/books/includes/bookDatabase.php and /var/www/html/movies/includes/movieDatabase.php. You need to replace the IP address there with the database server's IP address. If your password is not default password, replace it as well.

After this final step, Congratulations! You should have manually deploy a cloud architecture!

To run this app, open a browser and enter your web server's public floating IP address. You should see a page like:

<img src="images/Others/homePage.JPG">

Reference

1. http://howtoubuntu.org/how-to-install-lamp-on-ubuntu
2. https://stackoverflow.com/questions/8348506/grant-remote-access-of-mysql-database-from-any-ip-address
3. https://www.digitalocean.com/community/tutorials/how-to-install-the-apache-web-server-on-ubuntu-16-04

Research Note 10

Last week I was seriously illed so I did skip the note. Since my school starts, I can only use my leisure time on this research project, so without doubt, my progress is much slower compared to in summer. However, in last week, I have made some breakthrough and I thought it's worth noting.

C1 Smoothness Condition

Right now our main goal is to solve scatter fitting problem using tetrahedron partition in 3D space. Professor Schumaker has a outline to attempt this problem. We first need to partition the space using tetrahedron, then we need to figure out how to enforce C1 smoothness condition between each unit tetrahedron. After that, we can use the idea of minimum energy by adding a penalizing terms of the smoothness condition. If we can reach this far, we can then choose the right parameter and solve the problem.

Now the only missing key is to get the C1 smoothness condition. By C1 smoothness, it means that splines between different tetrahedron not only join on the sharing faces, but the directional partial derivative also agrees on the border. Since we are representing splines with coefficent at domain points, it can be proved that if the coefficient on the sharing faces agree, we will have C0 smoothness condition. And from that on, with one more layer, or I should say internal faces, of coefficient agree, we will one degree higher smoothness condition.

Enforcing coefficient to agree on the sharing face is easy. Because of the way we store the spline, it automatically satisfy the C0 smoothness condition. C1, on the other hand, is far more tricky. It involves lots of geometric manipulation of the coefficient tetrahedrons.

Long story short, Professor Larry came up with several functions to get the C1 smoothness condition with d = 3. The idea is to loop over all the interior face. For each face, get the left and right tetrahedrons' coefficient. Perform some geometric manipulation on them and obtain the internal layer of coefficient by rotations. Then enforce them to agree by setting up a linear system of equation.

After he show us the code, I spent a few days studying it. And in today's afternoon I combined his functions with my getindextetra and other helpers. With some struggles, I believe I have extended the functions to give C1 smoothness condition for any value of d. So far these functions pass all the test cases I have. This means we have obtained the final missing piece to solve scatter data fitting problems.

Probably in next week, we will put all the pieces together and finally obtain applicable solver.

New Beginning

By the end of summer, my research project officially ended. In last week, cooperating with Shiying Li, PhD student in Math, we came up with a Tetrahedron adaptive program. I will summarize our work in last week in this note.

T3adaptive

In previous week, I have already finish this new refinement function named refinetetra3 which will choose the largest boundary face and split it in 3 to get three new sub-tetrahedrons. Later in fall, I plan to finalize it so that it can split internal faces.

The remaining puzzles are interpolation functions, grid error evaluations and an error function to choose the tetrahedron to further refine. Luckily, these functions are just extensions from the 2D cases. And I have gathered enough experiences on them in 2D box project and 3D cube project. Therefore, it only took me 1 whole day to finish them.

After we have gathered all the tools, it's finally the time to put everything up together.

In the first part of our big program, we need to initialize the big partition. Other than reading from files like we did before, I also wrote a small initialization function to create Delaunay partition on 3D box of any size.

Then we use our good old friend, list function tetralist to get all the partition information we need.

Suppose we are given a target function, now that we have all the parts ready, we performed a first interpolation on the unrefine setting, which gives us a long coefficient vector $c$. Then we can start our main refinement process.

In each iteration, we did the following

1. Based on the current coefficient, we calculate the error scaled by the volume in each tetrahedron using the error function. By error we mean the max error inside the tetrahedron.
2. Then we sort the error and find the largest 20% and store their indices in a list.
3. For every tetrahedron in the list, if it's on the boundary, we use refine3 to refine one of its boundary faces; if it's a interior tetrahedron, we perform a Alfeld refinement using refine4.

We predefined the number of iteration. If in an iteration we have the max error of all tetrahedron less than a very small number, say $10^{-5}$, we exit the loop.

After the loop, we use grid error evaluation to estimate the error and plot the final refinement.

In this refinement program, we can get exact result for any trivariate polynomial given that the degree of basis Bernstein polynomial is big enough, which is what the approximation theory shows.

Since we still haven't found a nice way to displace 3D functions, I don't have a nice plot to show here.

New Start

As I mentioned, my research has officially over, but there are still lots of works that need to be done. So far no tetrahedron spline library exists and more functions are needed.

I have decided to continue working with Professor Schumaker on tetrahedron spline library in the fall semester. And I also plan to write a honor thesis on this topic to satisfy the requirement of graduating with honor in Math.

Since I have not enrolled in independent study or formally work as a RA, I can only use my spare time on this project. I will write a research note once every month instead of every week.

Basic Neural Network

​ In last few months, I have taken 3 online courses on coursera and udemy about machine learning. Now I am going on a serious of deep learning course by <a href = "https://lazyprogrammer.me/">lazy programmer </a>. I do love this series a lot and I want to post some of my note on my blog.

Forward Propagation

It's quite a work to type the derivation and theory part in latex so I choose to just scratch the surface of them.

The first important function in the course is the prediction of a neural network using forward propagation. We create a 3 groups classification problem and predict the label use randomly initialized weight and bias.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

import numpy as np 
import matplotlib.pyplot as plt 

Nclass = 500 

X1 = np.random.randn(Nclass,2) + np.array([0,-2])
X2 = np.random.randn(Nclass,2) + np.array([2,2])
X3 = np.random.randn(Nclass,2) + np.array([-2,2])
X = np.vstack([X1, X2, X3])

Y = np.array([0]*Nclass + [1]*Nclass + [2]*Nclass)

plt.scatter(X[:,0], X[:,1], c=Y, s = 100, alpha = .5)
# plt.show()	

D = 2
M = 3
K = 3

W1 = np.random.randn(D,M)
b1 = np.random.randn(M)
W2 = np.random.randn(M,K)
b2 = np.random.randn(K)

def forward(X,W1,b1,W2,b2):
	Z = 1 / (1 + np.exp(-X.dot(W1) - b1))
	A = Z.dot(W2) + b2
	expA = np.exp(A)
	Y = expA / expA.sum(axis = 1, keepdims = True)
	return Y

def classification_rate(Y, P):
	n_correct = 0
	n_total = 0
	for i in range(len(Y)):
		n_total +=1 
		if Y[i] == P[i]:
			n_correct += 1
	return float(n_correct) / n_total

P_Y_given_X = forward(X,W1,b1,W2,b2)
P = np.argmax(P_Y_given_X, axis=1)

assert(len(P) == len(Y))

print("Classification rate for randomly chosen weights:", classification_rate(Y,P))

The code is pretty straightforward.

Back propagation

The back propagation however, need some math derivation for the formula.

I really like the video "The WRONG way to learn Backpropagation" in this course. People do not understand bp by just taking the difference between target and prediction and claim that multiplying this $\Delta$ with weight is somehow the derivative of the cost function. We find this by strict math derivation and observe the pattern then build this algorithm. The key is to expose the pattern yourself using gradient decent and substitute the repeating pattern.

Some Hardcore Math

Suppose we are working on a multiclass classification problem with $K$ target classes.

Let's say we have a cost function $L = \sum_n\sum_k t_k^n \log(y_k^n)$ where $n$ denotes the sampel number and $k$ is for all the class. This is the famous cross entropy function and I have also discussed it in my previous note Beauty in Math.

Suppose we have $D$ input features and 1 hidden layer with $M$ nodes. We use $W_{DM}$ and $V_{MK}$ to denote the weight between input layer and hidden layer and weights between hidden layer and output layer respectively. For simplicity, suppose there is no bias terms. And further assume that we are using sigmoid function as activation function on hidden layer and softmax function on output layer.

Now our goal is to derive $\frac{\partial L}{\partial V_{MK}}$ and $\frac{\partial L}{\partial W_{DM}}$. First apply chain rule $$ \frac{\partial L}{\partial V_{mk}} = \sum_n \sum_{k'} t_{k'}^n \frac{1}{y_{k'}^n}\frac{\partial y_{k'}^n}{\partial V_{mk}} $$ Note here $k'$ is not the same as $k$ and $V_{mk}$ denotes one entry in weight matrix $V$.

We know that $$ y_{k'}^n = softmax(V_{K}^TZ) = \frac{\exp^{a_k}}{\sum_j \exp^{a_j}} \ a_k = V_K^T Z = \sum_m V_{mk} Z_m $$ where $Z$ denote the output at the hidden layer.

If $k = k'$ , we have $\frac{\partial y_{k'}}{\partial a_k} = y_{k'}(1-y_k)$, and else $\frac{\partial y_{k'}}{\partial a_k} = -y_{k'}y_k$. So we can rewrite them using Kronecker Delta as $$ \frac{\partial y_{k'}}{\partial a_k} = y_{k'}(\delta_{k'k} - y_k) $$ Then $$ \frac{\partial y_{k'}^n}{\partial V_{mk}} = \frac{\partial y_{k'}^n}{\partial a_{k}} \frac{\partial a_k}{\partial V_{mk}} =y_{k'}(\delta_{k'k} - y_k) Z_m $$ So $$ \begin{align} \frac{\partial L}{\partial V_{mk}} &= \sum_n \sum_{k'} t_{k'}^n \frac{1}{y_{k'}^n}y_{k'}^n(\delta_{k'k} - y_k^n) Z_m^n \ &= \sum_n \sum_{k'} (t_{k'}^n \delta_{k'k} - t_{k'}^ny_k^n)Z_m^n \ &= \sum_n (t_k^n - \sum_{k'} t_{k'}^ny_k^n )Z_m^n \ &= \sum_n (t_k^n - y_k^n )Z_m^n \end{align} $$ since $\sum_k t_k \equiv 1$ as each sample belongs to only one class.

The derivation for $\frac{\partial L}{\partial W_{DM}}$ is just some additional work. Since we have $$ Z_M = \frac{1}{1 + \exp^{-W_{MK}^TX}} $$ we can easily derive that $$ \begin{align} \frac{\partial L}{\partial W_{dm}} & = \sum_n (t_k^n - y_k^n ) \frac{\partial Z_m^n}{\partial W_{dm}} \ & = \sum_n (t_k^n - y_k^n )Z_m^n (1-Z_m^n)X_d^n \end{align} $$ The pattern here is actually already visiable. If we have one more hidden layer with same activation function, all we need to do is just to replace $X_d^n$ with the hidden value and multiply by addition factors. Because of the recursive nature of derivatives, back propagation is able to handle any deep neural network.

When putting these into code, we definitely want to vectorize them. It's pretty easy to make mistake as there are so many matrices to keep track of. One easy way I use to double check is to check the size of the matrix. Make sure their inner sizes match and the output size is the same as the target weight.

Code

In this code example, we keep working on the previous three groups problem.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87

import numpy as np 
import matplotlib.pyplot as plt 

def forward(X,W1,b1,W2,b2):
	Z = 1 / (1 + np.exp(-X.dot(W1)-b1))
	A = Z.dot(W2) + b2
	expa = np.exp(A)
	Y = expa / expa.sum(axis = 1, keepdims = True)
	return Y, Z

def classification_rate(Y,P):
	n_correct = 0
	n_total = 0
	for i in range(len(Y)):
		n_total += 1
		if Y[i] == P[i]:
			n_correct += 1

	return float(n_correct) / n_total


def cost(T, Y):
	tot = T * np.log(Y)
	return tot.sum()

def derivative_w2(Z, T, Y):
	return Z.T.dot(T - Y)

def derivative_b2(T,Y):
	return (T-Y).sum(axis=0)


def derivative_w1(X, Z, T, Y, W2):
	dZ = (T - Y).dot(W2.T) * Z * (1 - Z)
	return X.T.dot(dZ)

def derivative_b1(T, Y, W2, Z):
	return ((T-Y).dot(W2.T) * Z * (1-Z)).sum(axis=0)

def main():
	Nclass = 500
	D = 2 # input
	M = 3 # hidden layer
	K = 3 # output

	X1 = np.random.randn(Nclass,D) + np.array([0,-2])
	X2 = np.random.randn(Nclass,D) + np.array([2,2])
	X3 = np.random.randn(Nclass,D) + np.array([-2,2])
	X = np.vstack([X1, X2, X3])

	Y = np.array([0]*Nclass + [1]*Nclass + [2] * Nclass)
	N = len(Y)

	T = np.zeros((N,K))
	for i in range(N):
		T[i,Y[i]] = 1

	# plt.scatter(X[:,0], X[:,1], c=Y, s=100, alpha=.5)
	# plt.show()

	W1 = np.random.randn(D, M)
	b1 = np.random.randn(M)
	W2 = np.random.randn(M, K)
	b2 = np.random.randn(K)

	learning_rate = 10e-7
	costs = []
	for epoch in range(100000):
		output, hidden = forward(X,W1,b1,W2,b1)
		if epoch % 100 == 0:
			c = cost(T, output)
			P = np.argmax(output, axis = 1)
			r = classification_rate(Y,P)
			print("epoch", epoch, "cost", c, "classification rate:",r)
			costs.append(c)

		W2 += learning_rate *  derivative_w2(hidden, T, output)
		b2 += learning_rate *  derivative_b2(T, output)
		W1 += learning_rate *  derivative_w1(X, hidden, T, output, W2)
		b1 += learning_rate *  derivative_b1(T, output, W2, hidden)

	plt.plot(costs)
	plt.show()


if __name__ == '__main__':
	main()

Research Note 9

This week we mainly work on the tetrahedron project. Besides those little helper functions, I implemented two new main components of our project.

First is another new refine function. We have decided to renamed the previous refine function of tetrahedrons. The new one split a tetrahedron into 3 new tetrahedrons by finding the largest facet and its center point, and connect the opposite vertex to that center point to insert 3 new internal faces. In that way an old tetrahedron is splited into 3.

The implementation is quite similar to the previous function and only took me one day to finish this function.

<img src="/images/Research/Tetra/refine3.jpg" width = "500">

Here is a demonstration of the refinement.

Stochastic Walk

Another function is called stochastic walk. In the final stage of our project, we will have many points in our large tetrahedron partition, and we need to locate each point, aka, find the tetrahedron contain a given point.

Matlab has a built in function called pointLocation, which should accomplish that. However, I tested it with some simple cases and find that it fail in some. I have posted a question on the MathWork forem but receive no response so far.

So we decide to implement our own pointLocation function. The very first and foremost idea is to use the barycentric coordinate's property. If a point is located inside a tetrahedron then the barycentric coordinate of that point relative to that tetrahedron will have all positive components. However, this means we need to calculate the barycentric coordinate of a point for all the tetrahedron. This is too inefficient.

Another way is to build a K-D tree and search on the branch. However, K-D tree fails in some cases, e.g. triangulation with holes in it.

I did some research on Google and found a paper introducing some nice way to find point location. The paper is called Walk in a Triangulation by Olivier, etc. It can be found <a href="https://hal.inria.fr/inria-00072509/document">here</a>.

In the paper, there are 4 different ways to find point location. The idea is quite similar. Start from a tetrahedron, according to some condition, move to next tetrahedron and so on, and finally stop in the target tetrahedron.

I choose to implement the Stochastic Walk as it's easier to implement and we have the available helper functions. It starts from a random tetrahedron. Choose a random facet of it, and check whether the point is on the other side of the facet. If so, move to the next tetrahedron sharing that facet; if not, move to another facet of the current tetrahedron. In the mean time, it also 'remember' the last tetrahedron visited so that we won't ends up in an infinite loop in a pair of tetrahedron.

There are also some other little details, but I don't want to discuss it here.

Research Note 8

This week, I continued to work on the tetrahedron project and implement several functions.

Refine1

The refine1 function that I started last week is pretty tricky. There are many edge cases but fortunately I have finished it.

As mentioned in last note, the way this function works is to find the longest edge of the target tetrahedron and use an internal face to split that edge and also the whole tetrahedron into two smaller ones. These two smaller tetrahedron will always have same volume and share a face.

The hard part about this function is to keep track of all other list information. For example, after we insert an interior face, we create some two subedges of the longer edge which we split, and we also add two edges on the other two faces. If a neighbor tetrahedron has been refined before, then some of these new edges may already exist, so we have to first check their existence and then update or insert new information. Same process has to be applied to all the list information we keep.

Another tricky part is the orientation of all the list. For instance, we keep track of a #tetra by 6 matrix for all the edges of each tetrahedron. In order to obtain the right coefficient of spline when we later implement other functions, we need to keep an order among edges of a tetrahedron. So when new edges come, we have to maintain the order.

I don't want to explain all the details here. The final result is promising.

Refine3

On Wednesday, I met with the Professor. After our discussion, we thought a more common way to refine a tetrahedron is to split it into four smaller tetrahedron by connecting each vertex to the center of the original tetrahedron. So I start to work on a new function to accomplish that.

To be honest, this one is easier than refine1. Since all the new vertex, new edges, and new faces are inside the refined tetrahedron, I no longer need to check their existence and separate the cases. Though orientation is still a problem, with experience in refine1, I successfully implemented this function in this week.

The way we initialize a partition is always reading from files. It's a simple fact that a cube can be partitioned into 5 tetrahedron, 4 of which are identical and one in the middle, as the following diagram shows.

<img src="/images/Research/Tetra/tetra0.jpg" width="500">

After we refine the first tetrahedron using refine3 the output will be like:

<img src="/images/Research/Tetra/tetra1.jpg" width="500">

Research Note 7

This week I first put the 3d function pieces together to make a adaptive refinement code penalizing with C0 smoothness condition. The whole program is based on the 2d version that Professor Schumaker wrote.

Tadpat03d

The program is called Tadapt03d, where T for tensor product, adapt0 means adaptive fitting with C0 smoothness and 3d indicates that it's on 3d cubes.

Same as usual, the program begins with calling the list function reclist2 to get list of information of the partition. And a function handler is created, which is our target function to approximate.

It then called the interpdp03d function to interpolate the given function f with C0 smoothness condition enforced. The tsp0smooth3d function, which return the smoothness condition matrix M, is called inside the interpdp03d function. Now that we obtain the coefficient vector $c$, we start to refine.

The criterial is to calculate the error of approximation in each unit cube, and find the biggest 20% of all the cube to further refine. I wrote an error calculating function called rerr13d, which gives the errors in each cube. We then refine those cubes with large error. I would prefer to call them a single batch. After we refine all the cube in the batch, we call the interpdp03d function again to obtain the new coefficient vector $c$. Then we start another iteration.

The number of batches or level of refinement is given by user input.

We still haven't found a good way to visualize 3d functions so now I can only use error to check my program.

One thing to notice, at this point the tsp0smooth3d function is actually half-complete. As mention in last note, there are at least four cases in 3d C0 smoothness condition. I still haven't figure out the right way to put in hanging edges and hanging vertices on multiple edges. So the function actually miss some of the condition.

Tetra-mesh

Plot and readfile

On the other hand, I have started the tetrahedron project. This one will be the final project of this summer, and I believe it will be more challenging than the previous two.

As I did in the previous two projects, I started by looking for a way to visualize tetrahedrons. The function patch is still a reasonable choice. As I learn from visualizing cubes, I did basically the same — wrote a function to plot all the tetrahedrons given a matrix of points and a matrix of vertices of each tetra.

For tetrahedron, it's more like the 2d triangle case. It's hard to just initialize a bunch of tetras with few arguments. So I prefer to read in the files like Professor Schumaker did in his splinepak. I wrote a readtetra function to read in formatted tetrahedron data and get the basic information about the tetrahedron group.

List function

Now that I've got a working plot function and a read file function, I move on to the first cornerstone — list function.

Similar to both the cubelist and trilist case, I need to get all the needed information based on the very basic list I read from files. Specifically, I need a matrix E for tetra edge relation, a matrix F for tetra face relation, a matrix evertex for edge vertex relation, a matrix fvertex for face vertex relation and so on.

One tricky thing that's unique for tetrahedrons is their orientations. For rectangles and cubes, it's pretty easy to figure out how to put the edges and vertices in order. But for tetrahedron, so far I haven't get a nice and easy way to do the same.

I've tried to orient the edges in such a way that three edges of one face go first in an order that obey the right hand rule, and then the other three edges go next with same order. But later in my refine function, things become more complicated and I decide to leave such idea behind now.

It took me about two days or so to get a nice working list function.

Refine1

A tetrahedron can be refined in several ways, one into two, three or four. I started with the easiest, find the longest edge and use a interior face to split the edge and the whole tetra into two small ones.

I am still working on this function. There are some tricky place that needs to be handled. I believe by next week, I will have a nice working refine1 funcion ready.

Research Note 6

This week we debug the 2d code and keep working on 3d functions.

In 2d case, Professor Schumaker came up with a function to add C1 smoothness condition between rectangles. This function makes use of recl and rec, vectors denoting the rectangles on the left and right of each edge. Unfortunately, due to the way we develop refinerec, there is quite a few bugs in the function. So I spent rought two days on debugging the refinerec function and finally get a version that pass enormous test cases.

Interdp3d and Valtspgrid3d

On the 3d side, I develop some functions as the extensions of their 2d correspondences. Two of them are interdp3d and valtspgrid3d.

Interdp3d interpolate a given function on the whole rectangulation and return the coefficient as a 3d Matrix. Since Matlab does not support nested list structure, we look to tensor toolbox for help. A tensor is just a n dimension matrix. The toolbox include basic functionality such as addition, subtraction, scalar multiplication and scalar division. It also supports multiplication between tensor and matrix on a given direction. That's what we mainly use.

Back to my functions. Interdp3d returns a 3d tensor storing the coefficient of the 3d spline on the rectangulation. Then I use valtspgrid0 to obtain approximation of values of the given function at 3d grid points using the coefficient we get. To varify the result, I compute the norm of error terms.

We haven't found a good way to visualize the function. So far we have try to plot the isosurface. There is a nice package on MathWork forem named sliceomatic. It gives a live isosurface plot, but it still does not meet our expectation.

Other functions

Since my refinecube function only returns lists of relation among cubes, faces, edges and vertices, I need to find a way to get the other indispensible information such as the hanging vertices, subedges, etc.

So I wrote a bunch of small helper functions, like gethv for get hanging vertex, getse for get subedges, getsf for get subfaces and gethe for get hanging edges. Some of the cases are pretty tricky, for example in 3d a vertex can be hanging on multiplt edges.

Now that I prefer not to update these other important list in the refinement function, that function will definitely be less error-prone. However, in the view of efficiency, separating them out will without doubt be slower. So far it has not been a problem.

tsp0smooth3d

Now I am currently working on the tsp0smooth3d function. It enforces the C0 smoothness condition on the 3d partition. Compared to the 2d tsp0smooth I implemented half month ago, this one is much more complicated. In 2d, only two cases are in consideration, hanging vertices and subedges. On the other side, in 3d, there are at least 4 cases, hanging vertices, hanging edges, subedges and subfaces. And each cases is much harder to implement, since we have higher degree of Bernstein basis polynomial.

It's not that easy to describe the problem I ran into while implementing this code. The vital issue is that I cannot just easily visualize the work. Most of the cases need pure imagination. When 3d is involved, things always get more complicated.

I hope I can finish it in few days next week.

P.S.

One interesting and weird thing I noticed about Matlab is its meshgrid function. In 3d cases, one needs to do

​ [X,Y,Z] = meshgrid(ty,tx,tz)

to get the correct mesh in 3d, where ty, tx, tz are generated using linspace on their corresponding direction. Notice the weird order of argument. I am still wondering why these people prefers to pass the argument in this way. Just because of this ordering, I spent 2 more hours debugging my code.

Research Note 5

I was seriously ill last week so I haven't worked a lot and updated my Research Note.

This two week I keep working on the T-mesh 3D functions. By today, I have finished the cubelist and refinecube function. These two function corresponds to those in 2D case, a.k.a. reclist and refinerec. In particular, refinecube is much more complicated than refinerec so I took a completely different approach on writing it. I wrote three helper functions to update information about vertices, edges and faces respectively.

Other than the 3d functions, now that we have finish most part of 2d cases, we can finally create some sophisticated program on top of these Bernstein basis polynomial on rectangulation.

For example, in this week, Professor Schumaker had come up with an adaptive fitting program using the functions we have written last month.

The program fit a function on a rectangulation, calculate the error in each iteration and refine those rectangles that will decrease the error in greatest extent.

An output case:

​ For fitting function $$ f(x,y) = \tanh(40y-80x^2)-\tanh(40x-80y^2) $$ with a initial 5 by 5 rectangulation using Bernstein basis polynomial of degree 5 on both direction and with a 50 level adaptive fitting, we get the final rectangulation as

<img src="/images/Research/Tmesh/adapt1.jpg" width="500">

The function plot and the error plot are

<img src="/images/Research/Tmesh/adapt2.jpg" width="500">

<img src="/images/Research/Tmesh/adapt3.jpg" width="500">

Research Note 4

This week my Professor Schumaker left for a conference, and won't be back until next Friday. Before leaving, he has assigned me several tasks.

Tmesh Functions

The very imminent one is to continue work on the Tmesh functions. We need a function to perform penalize least square approximation on refined rectangulations. It consists of two major steps:

1. Getting the smoothness condition between unit rectangles.
2. Perform penalize least square approximation.

Smoothness Condition Matrix

Since we have refining our rectangles, we may have very weird partition of the original domain region. In order to obtain a smooth surface after we approximate a function, we need to calculate the smoothness condition on sharing edges and sharing vertices.

The idea is not very complicate. First, we obtain the number of smoothness condition $ns$ we need by counting sharing edges and sharing vertices. The formula also involves the dimensions on x, y directions.

Then we construct a $ns$ by $nc$ zero matrix, where $nc$ here is the number of coefficient for a spline on the rectangulation. We can use the sparse function in Matlab to save storage space, since the smoothness condition only involves the coefficients on the same edge.

To get the value at specific entry, we need to traverse all the shared edges and hanging vertices. Note the order we put the smoothness condition in each row does not matter, because when we enforce the condition, we will simply do it by minimize the product of $cM$ where $c$ is the spline coefficient and $M$ is our matrix.

As we have different degree of polynomial on different direction, I chose to first traverse all the horizontal edge and then vertical. In each iteration, I get all the coefficients of domain points on the current edge, and then solve for a linear system. This system represent the fact that the Bernstein polynomial on different side of the edge reach the same value on the domain points. In other words, this is also a interpolation problem, where I represent the coefficient of the Bernstein polynomial using the coefficient on the other side. We basically do the same on all the hanging vertices.

Finally, the function return the well-constructed matrix $M$ as output.

Penalize Least Square Approximation

This function is a little similar to the one that interpolate any given function, except for the fact that we are solving least square fitting problem accompany with smoothness condition.

I have finished the outline of the function and there is a problem that I have to wait till Professor get back. Basically I have to use the Kronecker product of Bernstein polynomial to get the linear system. I will add details about this function in next Research Note.

Refinerec

In addition to the two functions above, I also spent two whole days on debugging and fixing refinerec function.

Professor notice a significant bug on the previous version. It's pretty complicated to articulate the bug. I added 250 lines and finally fixed it. Now the whole refine function is 750 line, but I believe I could have combine some of the cases or extract out some helper functions to make the whole code more readable.

Other Directions

Other than the Tmesh functions that I need to work on, Professor pointed two direction for me in the next two weeks. First is the 3D version of the Tmesh problems. In other words, I need to extend our current work on 2D to 3D. We both know that I will be a huge amount of work, so he only expects me to first get all the geometric functions.

The other direction is 3D Tetrahedron partition. This one is even harder since we can't easily visualize what happen in 3D space. So I chose to first work on the 3D Tmesh direction.

3D Tmesh

Till this Sunday, I have finish several function we need. The path along which I work on these function is pretty straightforward. First I think we need to come up with a way to visualize the boxes in 3D. I found the patch function in Matlab and implement a simple function to display a cube. One example with label is as follow

<img src="/images/Research/Tmesh3D/reference.jpg">

Then I need a function that can plot any number of cube partition together, which calls on the simple plot function. Since we will definitely need a 3D version of reclist, which I named it cubelist, I choose to use the output from the list function to plot the cubes.

In 3D there are a lot of possible output containing information about the whole partition that we may need. Thus I need to think carefully before I move on. We have 4 different dimension of objects: vertices, edges, faces and cubes (or boxes). First and foremost, given $\textbf{[nx, ny, nz]}$ as a parameter for the initial partition, we inevitably need to output the coordinate of each vertices, and that is $\textbf{[x,y,z]}$, three outputs. Second, we need to have the relation between cubes and vertices, cubes and edges, and also cubes and faces; they are three matrices, V of size $ncube$ by 8, E of size $ncube$ by 12 and F of size $ncube$ by 6. In these matrices, we use the vertices number, edge number and face number to represent each cube.

After that, we need the relation between edges and vertices, between faces and edges, and between faces and vertices. These matrices can help us quickly get the information we need. I named them as evertex, fedge and fvertex.

Finally, I think in the future, we also need some opposite relation, which are matrices that connect vertices and cubes, connect edges and cubes, and connect faces and cubes. They help me to find the cube that are arround a given vertex, edge or face. In addition, they can also be used to find the boundary vertices, edges and faces as we put zero for nothing (or boundary). I call these matrices vstar, estar and fstar.

All these output I chose is entirely based on my own experience. Therefore, after Professor Schumaker comes back, we decide whether to keep all of them and add some more, or dump some of them.

With these output, it's far more than enough to generate a plot of the whole partition. Here is a 3 by 3 by 4 partition, where the numbers are points on each side.

<img src="/images/Research/Tmesh3D/cube1.jpg">

It's possible to labe all of the object but that will be too cumbersome and messy especially when we get to more partition.

I have already finish the cubelist function now, and I decide to move on to the hard part, the refinecube function in next week. Wish me good luck!