# 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 來設定亂數種子，例如:

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 語言實作:

#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++ 的意思是甚麼?

- ++ 在前的優先權與 * 相同，因此關聯性為從右到左
- ++ 在後的優先權高於 ++在前與 *，因此關聯性為從左到右

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)

以下是個可能的程式碼實作:

#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)