Matrix (2014/08/06 21:46)

-- [[
[[

# matrix API みんなでつくればこわくない #

* TODO
    * どのような形式でデータを保持するかの策定
    *
    * 加減算は、それぞれの行と列の数が揃ってないとできないので、そのチェック
    * 逆行列は正方行列(縦と横の数が同じ)じゃないと出てこないので、チェック
    * 逆行列は、正則じゃないと出てこないので、計算として違法なことをしたのか数学的に出せないからなのかで戻り値を変える？
    * 2次はサラスの公式使って簡略化

* 細かいこと
    * 行：たてに増えていく。 1行取り出したら横に長い。 row
    * 列：よこに増えていく。						column

-- (kssr) 3DCGなら4x4まであれば良い、んですよね？
-- (Quartz) ええ。しかし情報処理とかに使いたいときは、例えば通信の符号語長次元要りますね。
-- (Quartz) まあそれは後回しで。
-- (Noki 基本的には、n×m行列で、それぞれの成分をもった行列を生成するメソッドを基本に、
-- (Noki 各演算子に、メタテーブルで演算を定義していくのが最初の到達地点ですかね。
-- (Quartz) その辺は下に貼ったvectorAPI 参考に。
-- (Noki) ですね。もしかしたら、最初は2×2の行列のみを扱うのを作って、慣れるのもありかもしれません。
-- (Noki) 一般の行列は、ちょっち大変そう。
-- (h2)まずはnew関数および最初の行列のデータ構造かな
-- (Noki) ですかね。このAPIのnewメソッドを呼び出すたびに、新しい行列を生成(インスタンスを作る)でよさそうです。
-- (Noki) まんま、vectorAPIがそうですｶﾞ。
-- (h2)Luaのテーブルは実値じゃなくて参照だから・・・。それを忘れてプログラム書いて変なことになった経験がｗ
-- (Noki) 必ず参照渡しと言うのも、潔い言語ですよねw
-- (Quartz) ああ、私もあります。
-- (Noki) よし、僕は場が整って満足したので、今日は撤退しますw 明日も早いので…。
-- (h2)おつかれさまー
-- (Noki おつかれさまです。
-- (kssr)おつかれさまですー
-- (Quartz) お疲れ様でしたー

-- (Quartz) プライベートな変数要りますかね
-- (h2)まずは最初からmxnの行列を作るnew関数？配列いじるlocal関数がいくつかいるかな
-- (Quartz) mulScalar とかの簡単なのから実装していきますか。
-- (h2)forでまわすか、再帰使うか。for でいいよね。あとは中でtable.insert ? 行数と列数も保存しておけば計算速い？　あぁ、入力チェックがいる。
-- (k)引数でもらった行列のデータを、何も加工せずに持っておくだけってダメ？
-- (Q) テーブルなので、うっかり元のを変更すると行列の中身まで変わりますね。
-- (h2)じゃあ独自書式で？ だめだ応用がきかないｗ
-- (Q) いったん自身の内部に複製する必要があるような。内部的には多次元配列のかたちのテーブルで
-- (Q) それぞれの行列の要素をシャローコピーでいいんじゃないでしょうか
-- (h2)念のため、中身表示用のpp関数欲しいな。僕はそれから作るね。
-- (Q) 中で持ってるデータは変数なにに入れます？ (h2)なにとは？　とりあえずおまかせで。したがいます
-- (Q) あれ？ いまの構造のままだとプライベートメンバ持てない？
-- (k) 今の話の流れと関係ないですが、行列のデータ部分だけ取得出来たら良いと思う。（textutils.serialize()とか使って遊びたくなるかも）
-- (Q) 行列のオブジェクトそれぞれが関数持っててもデータがメモリ食うので、__index ありで行きましょう。privateなメンバはなしで。
-- (h)そういえば、入力データのチェック考えてないな。チェック用関数、private関数はどこに定義すればよい？　らじゃ
-- (h)ほかにどんなMatrixの関数を実装すればよいのかな。

2014/08/06
-- (Noki) ここって、今誰がいるのか、わかると便利ですね…。
--

]]
--]]




-- コードここから
-- ▼▼▼▼▼▼▼▼▼▼
local _new
-- map関数ほしくない？ map({1,2,3},{4,5,6}, plus) => {5,7,9} みたいの
local map = function(self, o, func)
  if #self ~= #o then error("Invalid arguments: different length") end
  local tmp = {}
  for i, v in ipairs(self) do table.insert(tmp, func(v, o[i])) end
  return tmp
end
local mapAdd = function(self, o) return map(self, o, function(x,y) return x+y end) end
local mapSub = function(self, o) return map(self, o, function(x,y) return x-y end) end

-- mmap({1,2,3},
-- 単位行列の作成  makeIdentity(3,2) => 3x2の単位行列
local makeIdentityMatrix(m,n)
  local tmp={}
  for i=1,n do
    local tmp2={}
    for j=1,m do
      table.insert(tmp2, ((i==j) and 1) or 0)
    end
    table.insert(tmp, tmp2)
  end
  return tmp
end
--

-- 2つの行列の縦横がそろっているかチェック(使わなかった
local eqaulColRow = function(self, o)
  return (self:col()==o:col() and self:row()==o:row())
end

-- おもしろいものつくったよ。名づけてマトリックスイテレータ  行列の中身を 横、縦の順に回します。(h お昼休みの終焉　では夜に
-- for col,row,v in makeMatrixIterator({{1,2,3},{4,5,6},{7,8,8}) do
--  print(col, ' ', row,' ',v)
-- end
local makeMatrixIterator =function(m)
  local col, row = 0,1
  local num_col, num_row =#(m[1]), #m
  return function()
    while row <= num_row do
      while col < num_col do
        col=col+1
        return col,row, m[row][col]
      end
      row = row+1
      col=0
    end
  end
end
-- matrixMap2( {{1,2,3},{4,5,6}}, (function(a,b) return (a+b) end), 1)
--  => {{2,3,4},{5,6,7}}
local matrixMap2 = function(m, func, ...)
  local tmp={}
  for col,row,v in makeMatrixIterator(m) do
    tmp[row] = tmp[row] or {}
    tmp[row][col]=func(v,...)
  end
  return tmp
end

local matrixMap = function(m, func)
  local tmp = {}
  for r, row in ipairs(m) do
    local tmp2 = {}
    for c, val in ipairs(row) do
      table.insert(tmp2,func(r,c,val))
    end
    table.insert(tmp,tmp2)
  end
  return tmp
end

local array = function(...)
  local args = {...}
  local a = {}
  if #args == 1 then
    for i=1,args[1] do
      a[i]=0
    end
  else
    local n = table.remove(args,1)
    for i=1,n do
      a[i]=array(unpack(args))
    end
  end
  return a
end

local isSquare = function(m) return m:row() == m:col() end
local inv; -- 長さnの二次元配列をn次正方行列としてみなして逆行列を返す
local cofactor;
local det = function(m) -- 長さnの二次元配列をn次正方行列としてみなして行列式を返す
  if #m == 2 then return m[1][1]*m[2][2]-m[1][2]*m[2][1] end
  local ret = 0
  for i=1,#m do
    ret = ret + confactor(m,i,1)*m[i][1]
  end
  return ret
end

local inv = function(m)
  local d = det(m)
  if d == 0 then return nil end
  local tmp = {}
  for r = 1,#m do
    local tmp2 = {}
    for c = 1,#m do
      table.insert(tmp2, cofactor(m,c,r)/d)
    end
    table.insert(tmp, tmp2)
  end
  return tmp
end

cofactor = function(m, i, j) -- 長さnの二次元配列（n次正方行列）m の i,j 余因子
  local tmp = {}
  for r = 1,#m do
    if r ~= i then
      local tmp2 = {}
      for c = 1,#m do
        if c ~= j then
          table.insert(tmp2, m[r][c])
        end
      end
      table.insert(tmp, tmp2)
    end
  end
  return -1^(i+j)*det(tmp)
end

-- validation 関数が必要だと思う　全ての行(row)の各要素が同じ数であることを確認。配列の形が直方体であること（？）
local isMatrix = function(tbl)
  if type(tbl)~="table" or #tbl==0 then return false end
  local numArg=#tbl[1]
  for i,v in ipairs(tbl) do
    if #v~=numArg then return false end
  end
  return true
end


local matrix = {
  add = function(self, o) return _new(map(self, o, mapAdd)) end,
  sub = function(self, o) return _new(map(self, o, mapSub)) end,

  mul = function( self, o ) --掛け算 mul( m1, m2 ) => m3
    local r = self:row()
    local n = self:col()
    local c = o:col()
    local tmp = {}
    for i=1,r do
      local tmp2 = {}
      for j=1,c do
        local s = 0
        for k=1,n do
          s = s+self._a[i][k]*o._a[k][j]
        end
        table.insert(tmp2,s)
      end
      table.insert(tmp,tmp2)
    end
    return _new(tmp)
  end,
  muls = function( self, n )　--(行列 x 数)、または(数 x 行列)  muls( m1, n ) => m2
  local tmp = {}
    for i,v in ipairs(self._a) do
      local tmp2 = {}
      for j,w in ipairs(v) do
        table.insert(tmp2, n*w)
      end
      table.insert(tmp, tmp2)
    end
    return _new(tmp)
  end,
  det = function( self ) -- 行列式? det( m ) => num, 1次元はスカラー、2次元たすきがけ、3次元はサラス、4次って？
    -- 3 次以上は余因子行列を用いて算出します
    if not isSquare(self) then return nil end
    return det(self._a)
  end,
  inv = function( self ) -- 逆行列 inv(m1) => m2
    if not isSquare(self) then return nil end
    return _new(inv(self._a))
  end,
  toArray = function( self )
    local tmp = {}
    for i,v in ipairs(self._a) do
      table.insert(tmp, {unpack(v)})
    end
    return tmp
  end,
  clone = function( self )
    return _new( self._a )
  end,
  tostring = function( self )  --一瞬 textutils.serialize使おうとおもったけど、行列としてみづらいｗ
    local tmp=""; local sep=" " -- じゃあ、末端のスペースそのままでもいいね？（邪悪顔 ありりｗ もち
	for i,row in ipairs( self ) do -- 入れ子が深い関数には罪悪感を感じる
      for j,col in ipairs( row ) do -- 表示用セパレータ（横）はカンマ？スペース？ 一般にはスペースのようです ok
        tmp=tmp..col..sep --(k)桁をそろえるように機能追加していい？ (h2)おねがい。一瞬考えたけど、全部の要素を一旦みるのめんｄ
        -- (k)あとでtostring2をつくっておきます。とかいいつつ今日はもう寝ます。お先です。(h)新しいのできたらこの関数を消して上書きしてください。
      end
      tmp=tmp..(i~=#self and "\n" or "")
    end
    return tmp
  end,
  round = function(self, nTolerance) -- 各要素の小数値を丸める。nToleranceは公差。round(m)=>全ての要素を整数に丸めた行列
    local nTolerance = nTolerance or 1.0
    local _round = function(num)
      return math.floor( num + ((nTolerance*0.5)/nTolerance) * nTolerance)
    end
    return _new( matrixMap2(self, _round) )
  end,
  row = function(self) return #self end,
  col = function(self) return #(self[1]) end -- 各行データの長さは同じだって信じてる
}


-- debug 用
-- debug用の行列表示用ヘルパー pretty print matrix
local ppm = function(m) print(matrix.tostring(m)) end

--[[
-- 和と差の検証用コード(h2 ではおやすみ　ノシ
a={{1,2,1},{1,3,5},{1,8,3} }
b={{3,2,1},{3,2,1},{3,2,1} }
c={{1.533,9,1.0001},{3.499,900.3,0.01}}

ppm(a)
ppm(b)
ppm(a:add(b))
ppm(a:sub(b))
ppm(c)
ppm(round(c))

--]]





local mmetatable = {
	__index = matrix,
	__add = matrix.add,
	__sub = matrix.sub,
	__mul = function( s, o ) return type(o) == "table" and s:mul(o) or s:muls(o) end,
	__unm = function( m ) return v:mul(-1) end,
	__tostring = matrix.tostring,
}
-- ex. new({{1,1,1},{1,1,1},{1,1,1}}) h)こちらのほうがいいね
-- unpackで一回分解しているから違うテーブルになるはずだけど
-- (Q) ええ、これでいいと思います
_new = function(...)
  local args = {...}
  local m = {
    _a = {},
  }
  local tbl
  if type(args[1]) == "table" then
    if args[1].toArray then
      tbl = args[1]:toArray()
    else tbl = args[1] end
    for i,v in ipairs(tbl) do
      table.insert(m._a, {unpack(v)})
    end
  else

  end
  setmetatable( m, mmetatable )
  return m
end

new = _new

-- ▲▲▲▲▲▲▲▲▲▲
-- コードここまで

-- 参考用 vector API

local vector = {
	add = function( self, o )
		return vector.new(
			self.x + o.x,
			self.y + o.y,
			self.z + o.z
		)
	end,
	sub = function( self, o )
		return vector.new(
			self.x - o.x,
			self.y - o.y,
			self.z - o.z
		)
	end,
	mul = function( self, m )
		return vector.new(
			self.x * m,
			self.y * m,
			self.z * m
		)
	end,
	dot = function( self, o )
		return self.x*o.x + self.y*o.y + self.z*o.z
	end,
	cross = function( self, o )
		return vector.new(
			self.y*o.z - self.z*o.y,
			self.z*o.x - self.x*o.z,
			self.x*o.y - self.y*o.x
		)
	end,
	length = function( self )
		return math.sqrt( self.x*self.x + self.y*self.y + self.z*self.z )
	end,
	normalize = function( self )
		return self:mul( 1 / self:length() )
	end,
	round = function( self, nTolerance )
	    nTolerance = nTolerance or 1.0
		return vector.new(
			math.floor( (self.x + (nTolerance * 0.5)) / nTolerance ) * nTolerance,
			math.floor( (self.y + (nTolerance * 0.5)) / nTolerance ) * nTolerance,
			math.floor( (self.z + (nTolerance * 0.5)) / nTolerance ) * nTolerance
		)
	end,
	tostring = function( self )
		return self.x..","..self.y..","..self.z
	end,
}

local vmetatable = {
	__index = vector,
	__add = vector.add,
	__sub = vector.sub,
	__mul = vector.mul,
	__unm = function( v ) return v:mul(-1) end,
	__tostring = vector.tostring,
}

function new( x, y, z )
	local v = {
		x = x or 0,
		y = y or 0,
		z = z or 0
	}
	setmetatable( v, vmetatable )
	return v
end