真正使用 mathjs 並結合 es6 tag template

mathjs 是一款 javascript 上的數學庫，我認為最大功能是用來算矩陣。後來才發現他真正的使用方法，用的好可以取代 matlab 的簡單矩陣運算。同時也發現其實他算是 es6 tag template 語法可以應用的地方。

老酒重開

其實蠻久以前就用過這東西，但那時根本沒發現他的威力，感覺除了算反矩陣，沒有其它特別突出，也不是很好用。

現在修研究所的平差課，同學一怒之下用 jupyter 寫平差作業，我也效法用 js 寫，於是又把 mathjs 給翻出來用。

長算式的問題

本來只是記得他有矩陣功能，平差會用到。但越寫越不對，覺得寫到後面語法真是醜到爆炸，尤其要寫一長串矩陣、反矩陣、相加相乘。因為 js 沒有算符重載，所以要加矩陣得要呼叫 mathjs 自己的函數，變成一大串前綴 m 表達式混在一起。

// 最小二乘平差解
// X = inv(A' * P * A) * A' * P * L
const X = MathJs.multiply(
    MathJs.inv(
        MathJs.multiply(
            MathJs.transpose(A), P, A
        )
    ),
    MathJs.transpose(A),
    P, L
)

上面那串還好一點，但這次作業還有約制，完整應該寫成 inv(A'*P*A + 1/s*C'*C) * (A'*P*L + 1/s*C'*W) ，會再多一層 MathJs.add() ，到最後寫到受不了，去翻官網發現他們有實作簡單的 parser！

javascript 中的仿 matlab parser

他們的 parser 蠻直覺的，就是能讓你用像 matlab 的語法寫矩陣運算，同時還有提供 scope，就是很古典，用 object 當 hash 當變數用。

const design = [[1,2],
                [3,4]]
const weight = MathJs.diag([1/3, 1/4, 1/2])
const observation = MathJs.transpose([[7,6,4]])

const X = MathJs.eval("inv(A' * P * A) * A' * P * L", {
    A: design,
    P: weight, 
    L: observation
})

在 eval 中 mathjs 會解析語法，並把 * + 之類的算符換成 MathJs.multiply MathJs.add ，同時還提供其它函數，像 ' 就會被翻譯成 MathJs.transpose ， inv 就是 MathJs.inv 。 MathJs 也附帶了不少函數，像 sin cos log，甚至也能定義自己的函數，然後放在 scope 裡傳進去。

最後我連定義矩陣也直接用 eval 做；畢竟原本要定義矩陣要先寫成二維陣列，一大堆方括號逗號寫起來太麻煩了。雖然不能像 matlab 一樣直接用空白、換行分隔欄列，不能省去逗號和方括號、分號有點可惜。

const design = MathJs.eval(`
    [1, 3, 5;
     2, 4, 6]
`)

es6 template tag hack

後來我想到，其實可以自己定義語法，自己用空白、換行分隔後做成矩陣，又想到可以結合 es6 template tag，寫起來更直觀。

MathJs.spaceMatrix = function (template, ...argument) {
    if (argument.length > 1) {
        throw new Error('variable interpolation not support')
    }
    const string = template[0]
    const array = string
          .split(/\n/g)
          .map(row => row.trim())
          .filter(Boolean)
          .map(row => row
               .split(/\s+/g)
               .map(Number))
    return MathJs.matrix(array)
}
const $sm = MathJs.spaceMatrix.bind(MathJs)

const design = $sm`
    1 2 3
    4 5 6
`

後來又想到，與其用那麼蹩腳的 object 當 scope，不如直接用 tag template 的變數內插直接傳變數進去。這本來就是 es6 tag template 適用的範圍。而且 parser mathjs 都幫我寫好了，我只要包裝成能 tag template 的結構就能用了。

MathJs.template = function (template, ...variable) {
    const prefix = '_variable'
    const scope = {}
    const expression = template.reduce((previousExpression, nextExpression, index) => {
        if (index == 0) return previousExpression + nextExpression
        else {
            const variableCount = index - 1
            const variableName = prefix + variableCount
            scope[variableName] = variable[variableCount]

            return previousExpression + variableName + nextExpression
        }
    }, '')
    return this.eval(expression, scope)
}

基本上思路就是重用原始的 eval 語法，所以就是把 tag template 中每個原本的變數位置，換成變數名字，並把該 scope 中的那個名字定義為該變數的值。這是最簡單只透過 MathJs.eval 與 parser 溝通的定義方式。

const $math = MathJs.template.bind(MathJs)

const A = $sm`
    2 3 4
    5 6 7
`
const L = $sm`
    1
    2
`
const P = $math`diag([1/3,1/2])`
const X = $math`
    inv(${A}' * ${P} * ${A}) * ${A} * ${P} * ${L}
`

這樣就直覺多了，像 perl 或 bash 在雙引號內部可以放其它字串變數的寫法。唯一問題我想是 js 的內插寫法太醜了，需要用大括號括起來，造成有點難以閱讀。也許這就是 tag template 很少人用、失敗的原因吧。

但事後想想，搞不好原本的寫法還比較好讀，上面是算式比較簡單的情況，如果複雜的話就會難讀了；雖然還是比原本只能用前綴寫法好。另外 mathjs 中可以相乘可以省略 * 單純用空隔，算是較貼近數學的寫法。

const X = $math`
    inv(${normalConstrain}) * (
        ${designMatrix}' ${p} ${observation} +
        1/${s} ${designConstrain}' ${valueConstrain}
    )
`

我有向 mathjs 提了這個 es6 tag template 功能的 issue ，也在 issue 中貼了上面那段範例程式，不知道會不會收。（放一天沒人回了，感覺是不會。最近提 issue 都被放置，難道是新手光環用完了，英文太爛被發現？） mathjs 專案有點大，底下一堆目錄腳本 require 來 require 去的，找起來有夠麻煩，等最近忙完再看要不要幫他改丟 pull request。

mathjs 其它功能

mathjs 其實可以不和外界溝，直接寫成多行；在裡面也能賦值變數，會變成 scope 物件的 key value。只是我這次作業沒有這樣寫，可能下次可以試試看。

矩陣維度問題

最後來抱怨一下 mathjs 一個不太喜歡的地方，向量和矩陣不能相乘；因為在 mathjs 中向量不是矩陣，只有 寬度 1 的矩陣 和 高度 1 的矩陣 ，沒有 行向量 也沒有 列向量 。

同時，如果把內積 a · b 視為 a' * b ，因為矩陣就是矩陣，只會變成高度寬度 1 的矩陣，他還是矩陣，不會是純量。而高寬 1 的矩陣基本上是不能乘大部份矩陣的，所以也造成了不少問題。

例如如果要定義列向量，就必須寫成寬度為一的矩陣，以我的習慣會寫成高度為一再轉置，可以少寫一點方括號。或其實也能用分號幫矩陣換行，只是如果是高度一的二維矩陣還是得寫二個方括號。

const scope = {}
MathJs.eval(`
    L = [[1,2,3]]'
    L2 = [[1],
          [2],
          [3]]
    L3 = [1; 2; 3]
    L4 = [1, 3, 4;]  # this cause syntax error
`, scope)
scope.L, scope.L2 // matrix

在 matlab 中純量就是高寬為 1 的矩陣，也是長度為 1 的向量；雖然向量還是分行列向量，但向量也就是寬度或高度為 1 的矩陣。也就沒有這麼多雷。