# 2020q3 Homework9 (quiz 9)
contributed by < `JimmyLiu0530` >
###### tags: `進階電腦系統理論與實作`
## 測驗1: 用 ANN 推理出 2-input multiplexer 的邏輯
[類神經網路 (ANN) 介紹](https://www.youtube.com/watch?v=aircAruvnKk)
在若干作業系統 (如 Linux, OpenBSD, MS-Windows, macOS) 中，[Address space layout randomization](https://en.wikipedia.org/wiki/Address_space_layout_randomization) (簡稱 ASLR) 是種防範記憶體損壞漏洞被利用的電腦安全技術，透過隨機放置行程關鍵資料區域的定址空間，防止攻擊者能可靠地跳躍到記憶體的特定位置。我們可利用 ASLR 來設定亂數種子，例如:
```c=
static uint64_t splitmix_x;
static void splitmix_init()
{
/* Hopefully, ASLR makes our function address random */
uintptr_t x = (uintptr_t)((void *) &splitmix_init);
struct timespec time;
clock_gettime(CLOCK_MONOTONIC, &time);
x ^= (uintptr_t) time.tv_sec;
x ^= (uintptr_t) time.tv_nsec;
splitmix_x = x;
}
```
由於函式 `splitmix_init` 的地址在 ASLR 生效後，可能會隨機得到不同的數值 (注意: 僅是「可能」，實際要看作業系統的表現)，我們可利用這特性來設定亂數種子 `splitmix_x`。
以下是一個類神經網路的 C 語言實作:
```c=
#include <assert.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
struct ann;
typedef double (*ann_act_func)(const struct ann *ctx, double a);
typedef struct ann {
int inputs, hidden_layers, hidden, outputs;
ann_act_func activation_hidden; /* used for hidden neurons */
ann_act_func activation_output; /* used for output */
int total_weights; /* weights and size of weights buffer */
int total_neurons; /* neurons + inputs and size of output buffer */
double *weight; /* All weights */
double *output; /* store input array and output of each neuron */
double *delta; /* total_neurons - inputs */
} ann_t;
#define LOOKUP_SIZE 4096
static inline double ann_act_hidden_indirect(const struct ann *ctx, double a) {
return ctx->activation_hidden(ctx, a);
}
static inline double ann_act_output_indirect(const struct ann *ctx, double a) {
return ctx->activation_output(ctx, a);
}
static const double sigmoid_dom_min = -15.0, sigmoid_dom_max = 15.0;
static double interval, lookup[LOOKUP_SIZE];
#define unlikely(x) __builtin_expect(!!(x), 0)
#define UNUSED __attribute__((unused))
static double ann_act_sigmoid(const ann_t *ctx UNUSED, double a) {
if (a < -45.0)
return 0;
if (a > 45.0)
return 1;
return 1.0 / (1 + exp(-a));
}
static void ann_init_sigmoid_lookup(const ann_t *ctx) {
const double f = (sigmoid_dom_max - sigmoid_dom_min) / LOOKUP_SIZE;
interval = LOOKUP_SIZE / (sigmoid_dom_max - sigmoid_dom_min);
for (int i = 0; i < LOOKUP_SIZE; ++i)
lookup[i] = ann_act_sigmoid(ctx, sigmoid_dom_min + f * i);
}
static double ann_act_sigmoid_cached(const ann_t *ctx UNUSED, double a) {
assert(!isnan(a));
if (a < sigmoid_dom_min)
return lookup[0];
if (a >= sigmoid_dom_max)
return lookup[LOOKUP_SIZE - 1];
size_t j = (a - sigmoid_dom_min) * interval + 0.5;
if (unlikely(j >= LOOKUP_SIZE))
return lookup[LOOKUP_SIZE - 1];
return lookup[j];
}
/* conventional generator (splitmix) developed by Steele et al. */
static uint64_t splitmix_x;
static inline uint64_t splitmix() {
splitmix_x += 0x9E3779B97F4A7C15;
uint64_t z = splitmix_x;
z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9;
z = (z ^ (z >> 27)) * 0x94D049BB133111EB;
return z ^ (z >> 31);
}
static void splitmix_init() {
/* Hopefully, ASLR makes our function address random */
uintptr_t x = (uintptr_t)((void *) &splitmix_init);
struct timespec time;
clock_gettime(CLOCK_MONOTONIC, &time);
x ^= (uintptr_t) time.tv_sec;
x ^= (uintptr_t) time.tv_nsec;
splitmix_x = x;
/* do a few randomization steps */
uintptr_t max = ((x ^ (x >> 17)) & 0x0F) + 1;
for (uintptr_t i = 0; i < max; i++)
splitmix();
}
/* uniform random numbers between 0 and 1 */
#define ANN_RAND() (((double) splitmix()) / UINT64_MAX)
/* Set weights randomly */
void ann_randomize(ann_t *ctx) {
for (int i = 0; i < ctx->total_weights; ++i)
ctx->weight[i] = ANN_RAND() - 0.5; /* weights from -0.5 to 0.5 */
}
/* Create and return a new network context */
ann_t *ann_init(int inputs, int hidden_layers, int hidden, int outputs) {
if (hidden_layers < 0)
return 0;
if (inputs < 1)
return 0;
if (outputs < 1)
return 0;
if (hidden_layers > 0 && hidden < 1)
return 0;
const int hidden_weights =
hidden_layers ? (inputs + 1) * hidden +
(hidden_layers - 1) * (hidden + 1) * hidden
: 0;
const int output_weights =
(hidden_layers ? (WWW) : (inputs + 1)) * outputs;
const int total_weights = (hidden_weights + output_weights);
const int total_neurons = (inputs + hidden * hidden_layers + outputs);
/* Allocate extra size for weights, outputs, and deltas. */
const int size =
sizeof(ann_t) + sizeof(double) * (total_weights + total_neurons +
(total_neurons - inputs));
ann_t *ret = malloc(size);
assert(ret);
ret->inputs = inputs;
ret->hidden_layers = hidden_layers;
ret->hidden = hidden;
ret->outputs = outputs;
ret->total_weights = total_weights;
ret->total_neurons = total_neurons;
ret->weight = (double *) ((char *) ret + sizeof(ann_t));
ret->output = ret->weight + ret->total_weights;
ret->delta = ret->output + ret->total_neurons;
ann_randomize(ret);
ret->activation_hidden = ann_act_sigmoid_cached;
ret->activation_output = ann_act_sigmoid_cached;
ann_init_sigmoid_lookup(ret);
return ret;
}
/* Return a new copy of network context */
ann_t *ann_copy(ann_t const *ctx) {
const int size = sizeof(ann_t) +
sizeof(double) * (ctx->total_weights + ctx->total_neurons +
(ctx->total_neurons - ctx->inputs));
ann_t *ret = malloc(size);
assert(ret);
memcpy(ret, ctx, size);
ret->weight = (double *) ((char *) ret + sizeof(ann_t));
ret->output = ret->weight + ret->total_weights;
ret->delta = ret->output + ret->total_neurons;
return ret;
}
/* Free the memory used by a network context */
void ann_free(ann_t *ctx) { free(ctx); }
/* Run the feedforward algorithm to calculate the output of the network. */
double const *ann_run(ann_t const *ctx, double const *inputs) {
double const *w = ctx->weight;
double *o = ctx->output + ctx->inputs;
double const *i = ctx->output;
/* Copy the inputs to the scratch area, where we also store each neuron's
* output, for consistency. This way the first layer isn't a special case.
*/
memcpy(ctx->output, inputs, sizeof(double) * ctx->inputs);
if (!ctx->hidden_layers) {
double *ret = o;
for (int j = 0; j < ctx->outputs; ++j) {
double sum = *w++ * (-1.0);
for (int k = 0; k < ctx->inputs; ++k)
sum += *w++ * i[k];
*o++ = ann_act_output_indirect(ctx, sum);
}
return ret;
}
/* Figure input layer */
for (int j = 0; j < ctx->hidden; ++j) {
double sum = *w++ * (-1.0);
for (int k = 0; k < ctx->inputs; ++k)
sum += *w++ * i[k];
*o++ = ann_act_hidden_indirect(ctx, sum);
}
i += ctx->inputs;
/* Figure hidden layers, if any. */
for (int h = 1; h < ctx->hidden_layers; ++h) {
for (int j = 0; j < ctx->hidden; ++j) {
double sum = *w++ * (-1.0);
for (int k = 0; k < ctx->hidden; ++k)
sum += *w++ * i[k];
*o++ = ann_act_hidden_indirect(ctx, sum);
}
i += ctx->hidden;
}
double const *ret = o;
/* Figure output layer. */
for (int j = 0; j < ctx->outputs; ++j) {
double sum = *w++ * (-1.0);
for (int k = 0; k < ctx->hidden; ++k)
sum += *w++ * i[k];
*o++ = ann_act_output_indirect(ctx, sum);
}
assert(w - ctx->weight == ctx->total_weights);
assert(o - ctx->output == ctx->total_neurons);
return ret;
}
#define LOOP_MAX 1024
int main(int argc, char *argv[]) {
splitmix_init();
/* Input and expected out data for the 2-input Multiplexer */
const double input[8][3] = {
{0, 0, 0}, {0, 0, 1}, {0, 1, 0}, {0, 1, 1},
{1, 0, 0}, {1, 0, 1}, {1, 1, 0}, {1, 1, 1},
};
const double output[8] = {
0, 0, 1, 1, 0, 1, 0, 1,
};
/* New network with 3 inputs, 1 hidden layer of 2 neurons,
* and 1 output.
*/
ann_t *ctx = ann_init(3, 1, 2, 1);
double err, last_err = LOOP_MAX;
int count = 0;
do {
++count;
if (count % LOOP_MAX == 0) {
ann_randomize(ctx);
last_err = LOOP_MAX;
}
ann_t *save = ann_copy(ctx);
/* Take a random guess at the ANN weights. */
for (int i = 0; i < ctx->total_weights; ++i)
ctx->weight[i] += RRR;
err = 0;
for (int k = 0; k < (1 << 3); k++)
err += pow(*ann_run(ctx, input[k]) - output[k], 2.0);
/* Keep these weights if they're an improvement. */
if (err < last_err) {
ann_free(save);
last_err = err;
} else {
ann_free(ctx);
ctx = save;
}
} while (err > 0.01);
printf("Finished in %d loops.\n", count);
/* Run the network and see what it predicts. */
for (int k = 0; k < (1 << 3); k++)
printf("[%1.f, %1.f, %1.f] = %1.f.\n", input[k][0], input[k][1],
input[k][2], *ann_run(ctx, input[k]));
ann_free(ctx);
return 0;
}
```
上方的程式碼嘗試由 Multiplexer 的輸入和輸出，推理出其中邏輯。Multiplexing switch 的作用:
一個具備 $2^n$ 輸入端的 Multiplexer (簡稱 MUX) 有 $n$ 個可選擇的輸入－輸出線路，可經由控制端來選擇 (select) 其中一個訊號作為輸出。
> 2-input Multiplexer Design
> 出處: [The Multiplexer](https://www.electronics-tutorials.ws/combination/comb_2.html)
該程式預期的執行輸出如下:
```
Finished in 359 loops.
[0, 0, 0] = 0.
[0, 0, 1] = 0.
[0, 1, 0] = 1.
[0, 1, 1] = 1.
[1, 0, 0] = 0.
[1, 0, 1] = 1.
[1, 1, 0] = 0.
[1, 1, 1] = 1.
```
其中 `359` 可能會替換為其他數值。請研讀程式碼，並補完程式碼，使得類神經網路得以運作。
### 程式說明
splitmix_init() 在這程式中所扮演的角色?
這個程式主要透過 feedforward algorithm 以及隨機猜測 ANN 的 weights 來模擬一個 2-input 多工器的結果。
- 函式 `ann_init`:
- 根據傳入的參數，依序為 # of input neurons, # of hidden layers, # of neurons per hidden layer, 以及 # of output neurons 來配置空間給 ANN 中的參數並初始化。
- 像是 `hidden_weights` 代表 input 跟 hidden、hidden 跟 hidden 間總共有多少個 weights (包含邊的 weights 以及點的 biases)。
- 而 `output_weights` 則是代表 hidden (或 input) 跟 output 間的 weights 總數，因此，`WWW = hidden + 1`，因為 `+1` for biases of output-layer neurons。
- 此外，在初始化的過程中，還會去呼叫函式 `ann_randomize` 來給予 ANN 中所有 weights 一個隨機的值，其中值介於 -0.5 ~ 0.5 。
- 最後，呼叫函式 `ann_init_sigmoid_lookup` 來預先計算好輸入介於 `sigmoid_dom_min` 與 `sigmoid_dom_max` 間的激勵函數值，並存於 `lookup` 表格中。
:::info
:notes: NOTE:
這邊的預先計算並不是真的把所有在這範圍內的實數都算出來，因為我們 `lookup` 表格大小有限，因此只能將這個範圍切出 `LOOKUP_SIZE` 個區間，在相同區間內的 x 都會得到相同的函數值。
:::
- 函數 `ann_run`
- line 184-193 這是 for 沒有 hidden layer 的 ANN
- 如有 hidden layer，主要分成三大類去計算他們的激勵函數值，input layer -> hidden layer、hidden layer -> hidden layer、以及 hidden layer -> output layer，但其實不管哪一類，操作都很相似。
- 舉 input layer -> hidden layer 為例:
- 對第一層 hidden layer 的每個神經元去計算他們的激勵函數值。
- 每個神經元的激勵函數的輸入為所有 input layer 的 neurons 的 feature values 乘上 weight，再加上 bias。透過呼叫函式 `ann_act_hidden_indirect`，查詢 `lookup` 表格，即可取得其激勵函數值，並存入 `*o++`
- 最後，所有 output layer 神經元的激勵函數值即為此次預測的結果，並回傳存儲存這些值的初始位置。
:::info
:notes: NOTE:
1. *o++ 的意思是甚麼?
- ++ 在前的優先權與 * 相同，因此關聯性為從右到左
- ++ 在後的優先權高於 ++在前與 *，因此關聯性為從左到右
```c=
int a[]={10, 20};
int *o = a;
printf("%d\n", *o++); /* Get 10, and o points to the next element(i.e. 20) */
printf("%d\n", (*o)++); /* Get 10, after printing, 10 will be increased by 1 */
printf("%d\n", *++o); /* Get 20 */
printf("%d\n", ++*o); /* Get 11 */
```
2. 除了 input layer 的神經元的 input 是先給定的之外，其餘神經元的激勵函數值會在計算下一層神經元的激勵函數時用到。
:::
- line 267-268 會計算此次預測的 loss function，並與上一次的預測做比較。若這次的預測比上次差，則捨棄這次的預測；否則捨棄上一次的預測。不斷猜測 ANN weights，直到預測的誤差小於等於 0.01。
### 延伸問題
#### 1. 使用 back-propagation (反向傳遞) 改寫上述程式碼，並探討激勵函數的特性
#### 2. 引入其他激勵函數，例如 ReLU，並確保針對其特性進行適度初始化
#### 3. 比照 [Tinn](https://github.com/glouw/tinn)，搭配給定的 data set，實作出能夠辨識手寫印度—阿拉伯數字的類神經網路
> - 參考 [nanonn](https://github.com/zserge/nanonn)
#### 4. 研究亂數相關文獻，並探討偽亂數產生器 (PRNG) 的實作
> - [你的程式夠「亂」嗎？](https://www.ithome.com.tw/voice/110007)
> - [亂數](https://blog.danielchen.cc/2019/01/31/random-number/)
> - [The fastest conventional random number generator that can pass Big Crush?](https://lemire.me/blog/2019/03/19/the-fastest-conventional-random-number-generator-that-can-pass-big-crush/)
> - [ARM and Intel have different performance characteristics: a case study in random number generation](https://lemire.me/blog/2019/03/20/arm-and-intel-have-different-performance-characteristics-a-case-study-in-random-number-generation/)
#### 5. 研讀 [Vectorizing random number generators for greater speed: PCG and xorshift128+](https://lemire.me/blog/2018/06/07/vectorizing-random-number-generators-for-greater-speed-pcg-and-xorshift128-avx-512-edition/)，重現實驗並思考 SIMD 加速 PRNG 的機會
#### 6. 研讀 [Bitwise Neural Networks](https://arxiv.org/abs/1602.02830) 和 [Binarized Neural Networks: Training Deep Neural Networks with Weights and Activations Constrained to +1 or -1](https://arxiv.org/abs/1601.06071)，摘錄論文重點，並試著在上述程式碼的基礎上予以實作
> - [tinier-nn](https://github.com/codekansas/tinier-nn)
-------------------------------------------------
## 測驗2: LeetCode [37. Sudoku Solver](https://leetcode.com/problems/sudoku-solver/)
數獨規則如下:
1. 每列 (row) 必須出現 1 到 9 的數字，且只能出現一次
2. 每行 (column) 必須出現 1 到 9 的數字，且只能出現一次
3. 每個 $3×3$ 構成的「宮」(即子格子，以粗框標記) 內，數字 1 到 9 必須出現，且只能出現一次
示範:
解法:
若每次挑出一空格，並展開完成該盤面的可能數字，嘗試展開所有解法路徑，可畫成一個樹 (tree) 來思考：
其中必定存在一條從 root 到 leaf 的路徑，能將其盤面填滿，並且盤面上的所有數字，能滿足數獨的所有條件規則 —— 題目假設存在唯一解。在決定空格能填什麼數字時，已將不符合規則的數字排除，以減少後續的計算量
> [參考解題思路](https://hackmd.io/@kenjin/BkTcYzRLH)
以下是個可能的程式碼實作:
```c=
#include <stdbool.h>
#include <stdint.h>
#include <stddef.h>
typedef struct {
uint16_t row[9], col[9], box[9];
} masks_t;
const uint16_t ALL_NUMS = 0x1FF;
static inline int bi(int x, int y)
{
return (y / 3) * 3 + (x / 3);
}
static bool solve(char **board, masks_t *m)
{
for (int y = 0; y < 9; ++y) {
for (int x = 0; x < 9; ++x) {
if (board[y][x] != '.')
continue;
int b = bi(x, y);
uint16_t posmoves = m->row[x] & m->col[y] & m->box[b];
while (posmoves) {
size_t n = __builtin_ctz(posmoves);
uint16_t k = ~(1U << n);
board[y][x] = '1' + n;
m->row[x] &= k, m->col[y] &= k, m->box[b] &= k;
if (solve(board, m))
return true;
m->row[x] |= ~k, m->col[y] |= ~k, m->box[b] |= ~k;
posmoves &= MOVE;
}
board[y][x] = '.';
return false;
}
}
return true;
}
void solveSudoku(char **board, int boardSize, int *boardColSize)
{
masks_t m;
for (int i = 0; i < 9; ++i)
m.row[i] = m.col[i] = m.box[i] = ALL_NUMS;
for (int y = 0; y < 9; ++y) {
for (int x = 0; x < 9; ++x) {
if (board[y][x] != '.') {
int n = board[y][x] - '1';
uint16_t k = ~(1U << n);
m.row[x] &= k, m.col[y] &= k, m.box[bi(x, y)] &= k;
}
}
}
for (int y = 0; y < 9; ++y) {
for (int x = 0; x < 9; ++x) {
if (board[y][x] != '.')
continue;
uint16_t mask = m.row[x] & m.col[y] & m.box[bi(x, y)];
if ((mask & (MMM)) == 0) {
size_t n = __builtin_ctz(mask);
uint16_t k = ~(1U << n);
board[y][x] = '1' + n;
m.row[x] &= k, m.col[y] &= k, m.box[bi(x, y)] &= k;
}
}
}
solve(board, &m);
}
```
請補完程式碼，使得程式符合 LeetCode 題目要求。
### 程式說明
先考慮函式 `solveSudoku` :
`m.row[i]`, `m.col[i]`, `m.box[i]` 這三個陣列分別用來記錄第 `i` 行, 列, 3x3 格子中用過哪些數字，並以 bit 表示。舉例來說，`m.row[0] = 0000 0000 0110 1010` 代表在第 0 行，數字 2, 4, 6, 7 已被使用過。因此，在 line 48-56 根據現有的 9x9 棋盤格去更新每行、列以及 3x3 格子中，哪些數字被用過了。此外，`bi(int x, int y)` 函式的功能在於回傳 (x, y) 格子位於哪個 3x3 方格。
至於 line 58-70 主要針對棋盤格上的空格作處理，其中 `mask` 紀錄某個 1x1 格子可以填入的數字，若只有一種數字可以填，那麼直接填入該格，並更新 `m.row[i]`, `m.col[i]`, `m.box[i]`。因此 `MMM = mask - 1`。
上述作業完成後，呼叫函式 `solve`，主要利用遞迴地去檢查每一種可能到底有沒有符合數獨的規則。一旦有任何一個 1x1 格子沒有數字可以填，則直接回傳 `false`。
因為前面已經處理過只有一種數字可填的情況了，因此剩下的空格都是至少有兩種可能。line 25-34 就是在試看看每個格子的所有可能的數字，每試一種數字都會遞迴呼叫函式 `solve` 去檢查在這種嘗試下的棋盤格是否能完成數獨。若可以，則回傳 `true`； 否則須將此種數字從 `posmoves` 中移除。因此，`MOVE = posmoves - 1`。
### 延伸問題
#### 1. 重寫不同於上述的程式碼，效能表現應優於上述
#### 2. 嘗試運用 SIMD 指令加速數獨求解
> - 參考 [Nerd Sniped: A Sudoku Story](https://t-dillon.github.io/tdoku/) 及 [simdoku](https://github.com/yotarok/simdoku)
Sign in with Wallet
Connect another wallet