カタベログ

IT技術に関するブログを書きたい.食べ物関連はInstagramをご参照の事.

MacOS(Mojave)でDatabricks koalasを利用しようとしてつまづいた

前提

  • 買ったばかりので開発系ツールがあまり入っていないMac環境
  • GithubのReadMeを読めばおおよそ回避できる

What's "Koalas"?

  • SparkのDataFrameの為にPandasライクなAPIを用意したのがKoalas
  • 米国Databricks社が開発してOSSとして公開した、Python用のフレームワーク
  • 米国Databricks社はApache Sparkを使った製品を開発・提供している(Sparkへのコミットもしている)
  • SparkではPySparkと呼ばれるPythonからSparkを扱うためのツールが用意されている(Spark自体はJVM言語のScalaで開発されている)
  • PySparkでDataFrameを用いる時はSparkが用意するものを使うのが高速
  • 一方、PythonではDataFrameを扱うときは一般的にPandasを用いる
  • PandasのDataFrameの書き方とSparkのDataFrameの書き方は異なり学習コストがかかる

つまづきポイント

Python 3.7.xでのpyarrowのバグを踏んだ

Github(https://github.com/databricks/koalas)のGet Startedに書いてるやんけ。 と言うことで、あらかじめ導入していたpyenvとvirtualenvを使って3.6.8の仮想環境を用意しよう。

If this fails to install the pyarrow dependency, you may want to try installing with Python 3.6.x, as pip install arrow does not work out of the box for 3.7 https://github.com/apache/arrow/issues/1125.

pyenvで3.6.xをインストールしようとするとCのヘッダーファイルがなかった

3.7.2入れた時にもハマりましたが、Xcodeをアップデートした時にまた元に戻った・・・のかな? http://kane-please.hatenablog.com/entry/2018/11/07/001823

cmakeの他にCythonやOpen JDKもいるやんけ

Dependencesに書いてないやん。当たり前なのか?

でも結局

In [1]: import databricks.koalas as ks

In [2]: df = ks.DataFrame({'x':range(3), 'y':['a','b','b'], 'z':['a','b','b']})
WARNING: An illegal reflective access operation has occurred
WARNING: Illegal reflective access by org.apache.spark.unsafe.Platform (file:/Users/tkm/.pyenv/versions/3.6.8/envs/mySpark/lib/python3.6/site-packages/pyspark/jars/spark-unsafe_2.12-2.4.2.jar) to method java.nio.Bits.unaligned()
WARNING: Please consider reporting this to the maintainers of org.apache.spark.unsafe.Platform
WARNING: Use --illegal-access=warn to enable warnings of further illegal reflective access operations
WARNING: All illegal access operations will be denied in a future release
19/04/26 06:52:11 WARN NativeCodeLoader: Unable to load native-hadoop library for your platform... using builtin-java classes where applicable
Using Spark's default log4j profile: org/apache/spark/log4j-defaults.properties
Setting default log level to "WARN".
To adjust logging level use sc.setLogLevel(newLevel). For SparkR, use setLogLevel(newLevel).

Docker でCentOS用意しよ…。

Raspberry Pi Type-Bに,とりあえずArch Linux入れました.

2,3週間ぐらい前にRaspberry Pi買ったんですが,奇跡的に自宅にSDカードを使うことのできる端末がなく,放置してました.
USBメモリでどうにかできないかと試行錯誤してたけど,もう諦めてリーダー買ってきました.

手順をざっくりと書くと以下の通りです.
Raspberry Piの公式ページからArchLinuxのimgファイルをダウンロードする
→SDカードをコンピュータに接続する
→diskutilやdf,mountコマンドでディスク識別子を調べる
→diskutil umount [ディスク識別子]でアンマウントする
→sudo dd if=[imgファイルのパス] of=/dev/rdisk[ディスク識別子にあった番号]
(例:disk2s1ならrdisk2)
→しばし待つと終了メッセージが現れる.
→再マウントされているので,/Volume内のSDカードのでディレクトリに入って正しく書き込まれていることを確認してdiskutil ejectする.
Raspberry PiにSDカードを挿し,電源ON.表面のLEDの内,PWR以下の三つが光ってれば無事に読み込まれている.

細かいことはぐぐって調べてください.

Windows8でBUFFALO BSHSBD04を利用できるようにする方法

FILCO Majestouch MINILA Air JP 68key 茶軸(FFBT68M/NB)を買いました.
Bluetoothのコンパクトなキーボードが欲しかったのです.
冬のボーナスも出たことなので遠慮なく買ってきた訳です.


自分の持ってるデスクトップでも利用できるようにとBurtoothアダプタを買ってきた訳ですが,なんと公式のドライバがWindows8に非対応だというのです.
ワイヤレス生活を満喫するつもりMAXだったのに・・・.


心底がっかりしながらも,藁にもすがる思いでGoogle先生に縋り付きました.
すると.価格.comの口コミにベストアンサーがありました.
その方法というのが,「東芝製ユーティリティソフトを入れることで解決できる」とのこと.
http://www.princeton.co.jp/download/ptm-ubt3s/top.html
ページ中央にあるダウンロードボタンをクリックすると,「v70005.zip」というファイルがダウンロードされる.解凍して「setup.exe」を叩けば終わり.
再起動して「新しい接続の追加」を実行すれば,あとは画面に現れる手順に従えば接続できる.


安物買いの銭失いになるところだった.

リカレントネットワークによる記憶画像の想起 [追記]

ネックであった重み計算でのループ回数の削減を試みました.今までは記憶させるパターンのN次元ベクトルをp個生成するという形でした.今回のは記憶させるパターンのi番目の画素を要素としたp次元ベクトルをN個生成します.こうすることで内積一発で重みが計算できるので1000個記憶させても2個記憶させたときと同じくらいの速度で重みが計算されます.たぶん,これ以上は内積を利用した高速化は無理だと思います.

	#直交する記憶パターンの生成
	for i in xrange(self.__size):
            tmpNDArray = np.zeros(self.__patternNum)
            tmpNDArray[0] = self.__lena[i]
            for k in xrange(1, self.__patternNum):
		tmpNDArray[k] = (1.0 if(rnd.random() > 0.5) else -1.0)
            self.__pattern.append(tmpNDArray)

	self.__weight = np.zeros((self.__size, self.__size))
	for i in xrange(self.__size):
            for j in xrange(i, self.__size):
		if i != j :
                    self.__weight[i][j] = np.dot(self.__pattern[i], (self.__pattern[j]).T) / self.__size
                    self.__weight[j][i] = self.__weight[i][j]

こんな感じで該当部分を書き換えると動くと思います.なんでこんな簡単なことを見落としていたのだろう・・・愚かですね.

リカレントネットワークによる記憶画像の想起

こんばんは.

今日は課題で作ったリカレントネットワークによる記憶画像の想起スクリプトを載せます.画像ファイルはPBM形式のものを扱ってます.課題用なのでかなり課題に特化したものなので,特にファイル操作部分に関しては各自変更が必要でしょう.

import numpy as np
import random as rnd
import copy

def sgn(num):
    return (1 if(num > 0.0) else (0 if(num == 0.0) else -1))

class AssociativeMemory:
    def __init__(self, patternNum=1):
        self.__patternNum = patternNum
        self.__pattern = []

        #ファイル操作
        tmpArray = []
        lineNum = 0
        for line in open('ファイル名','r'):
            if lineNum >= 2:
                line.rstrip('\n')
                for charNum in xrange(len(line)):
                    if ((line[charNum] == "1") or (line[charNum] == "0")):
                        tmpArray.append(1.0 if (int(line[charNum]) == 1) else -1.0)
            lineNum = lineNum + 1
        self.__size = len(tmpArray)
        self.__lena = np.array(tmpArray)
        self.__pattern.append(self.__lena)
        
        #直交する記憶パターンの生成
        for k in xrange(self.__patternNum):
            tmpNDArray = np.zeros(self.__size)
            for i in xrange(self.__size):
                tmpNDArray[i] = (1.0 if(rnd.random() > 0.5) else -1.0)
            self.__pattern.append(tmpNDArray)

        #重みの計算
        self.__weight = np.zeros((self.__size, self.__size))
        for i in xrange(self.__size):
            for j in xrange(i, self.__size):
                if i != j :
                    for p in xrange(len(self.__pattern)):
                        self.__weight[i][j] = self.__weight[i][j] + self.__pattern[p][i] * self.__pattern[p][j]
                    self.__weight[i][j] = self.__weight[i][j] / (self.__size)
                    self.__weight[j][i] = self.__weight[i][j]
                else:
                    self.__weight[i][j] = 0.0
                    self.__weight[j][i] = 0.0
        
    def recall(self, noiseRate = 0.0, step = 15):
        state = copy.deepcopy(self.__lena)
        self.__m = []

        #原画像にノイズを加えてそれを初期状態とする
        for i in xrange(len(state)):
            if noiseRate > 0.0:
                state[i] = (state[i] if(rnd.random() > noiseRate)\
                            else (1.0 if(state[i] == -1)\
                                  else -1.0))
        #初期の誤差を計算
        prodResult = np.dot(state, self.__lena.T)
        self.__m.append(prodResult / (self.__size))
        
        for t in xrange(step):
            nextState = np.zeros(self.__size)

            #想起
            for i in xrange(self.__size):
                prodResult = np.dot(self.__weight[i], state.T)
                nextState[i] = sgn(prodResult)

            #時刻tにおける誤差
            prodResult = np.dot(nextState, self.__lena.T)
            self.__m.append(prodResult / (self.__size))
            state = copy.deepcopy(nextState)

        return self.__m


ドライバはこんな感じ

#!/opt/local/bin/ python
# -*- coding: utf-8 -*-

import AssociativeMemory

test = AssociativeMemory.AssociativeMemory(2)#パターン数を自然数で
print test.recall(0.4)#入力画像として原画像にどれだけノイズを混ぜるかを指定(最高1.0)


こんな感じです.できるだけ計算量を削減させようと工夫したけど,結局,かなり遅いです.重みは対称行列なことを利用したりしたんですが,一番時間かかるのは重みの計算です.見通しを立てて,ループ削減が難しそうならC++でやる方がいいよね.


追記があります.次の記事を参照すると幸せになれます.

課題で作ったものいろいろ

こんばんは.明日提出の課題を解いたけどドキュメント化していないアカウントです.
せっかくいろいろ作ったのでなんかの役に立つかもな,もとい,俺こんなの作ったよ褒めて褒めてっという目的の基にスクリプトソースコードを晒していこうと思います.

まず,ジニ係数算出スクリプトから.
必要なものはNumPy,ScipPyです.

#!/usr/bin/env python

import sys,string
import numpy as np

inputData = []
#データをプログラムで扱うための前処理
#データにあわせて適宜書き換える
for line in open('ファイル名','r'):
    inputData.append(float(line.rstrip('\r\n')))

inputData.sort()
N = len(inputData)
x = np.array(inputData)
mean = np.mean(x)
tmp = []
for i in xrange(N):
    tmp.append(2 * i - N - 1)

val = np.array(tmp)
gini = np.dot(val, x) / (N**2 * mean)

print "Gini Index : " + str(gini)

簡単ですね.

次に平均場近似を用いた画像修復スクリプト

import sys, copy, string, math, random

q = 0.2
h = 0.5 * math.log((1-q)/q)
J = 0.4
th = 0.0001

currentError = float(sys.maxint)
reconstructRate = 0
errorList = []
baseBitString = []
degBitString = []
newDeg = []
oldDeg = []
baseDeg = [] 
base = []

#データをプログラムで扱うための前処理
#データにあわせて適宜書き換える
lineNum = 1
for line in open('ファイル名','r'):
    if lineNum >= 3:
        line.rstrip('\n')
        for charNum in xrange(len(line)):
            if ((line[charNum] == "1") or (line[charNum] == "0")):
                baseBitString.append(1 if (int(line[charNum]) == 1) else -1)
                if (random.random() < q):
                    degBitString.append((-1 if (int(line[charNum]) == 1) else 1))
                else :
                    degBitString.append(1 if (int(line[charNum]) == 1) else -1)
                
    lineNum = lineNum + 1

for i in xrange(50, len(baseBitString) - 50, 50):
    base.append(baseBitString[i:i+50])
    baseDeg.append(degBitString[i:i+50]) 

for i in xrange(len(baseDeg)):
    tmp = []
    for j in xrange(len(baseDeg[i])):
        newData = math.tanh(J * ((baseDeg[i][j+1] if ((j + 1) < len(baseDeg[i])) else 0)\
                                 + (baseDeg[i][j-1] if ((j - 1) >= 0) else 0)\
                                 + (baseDeg[i+1][j] if ((i + 1) < len(baseDeg)) else 0)\
                                 + (baseDeg[i-1][j] if ((i - 1) >= 0) else 0)) + h * baseDeg[i][j])
        tmp.append(newData)
    oldDeg.append(tmp)

print 'Base Image:'
for i in xrange(len(base)):
    print str(base[i])
    
print 'Degeneration Image :'
for i in xrange(len(baseDeg)):
    print str(baseDeg[i])

print "Error :"
while (currentError > th):
    currentError = 0.0;
    for i in xrange(len(baseDeg)):
        tmp = []
        for j in xrange(len(baseDeg[i])):
            newData = math.tanh(J * ((oldDeg[i][j+1] if ((j + 1) < len(oldDeg[i])) else 0)\
                                   + (oldDeg[i][j-1] if ((j - 1) >= 0) else 0)\
                                   + (oldDeg[i+1][j] if ((i + 1) < len(oldDeg)) else 0)\
                                   + (oldDeg[i-1][j] if ((i - 1) >= 0) else 0)) + h*baseDeg[i][j])
            currentError = currentError + abs(newData - oldDeg[i][j])
            tmp.append(newData)
        newDeg.append(tmp)
    oldDeg = copy.deepcopy(newDeg)
    errorList.append(currentError)
    print currentError
    del newDeg[:]

for i in xrange(len(oldDeg)):
    print str(oldDeg[i])
 
for i in xrange(len(oldDeg)):
    for j in xrange(len(oldDeg[i])):
        oldDeg[i][j] = (1 if (oldDeg[i][j] > 0) else (0 if (oldDeg[i][j] == 0) else -1))
        reconstructRate = reconstructRate + (1 if (abs(base[i][j] - oldDeg[i][j]) == 0) else 0)

print 'Reconstruct Image:'
for i in xrange(len(oldDeg)):
    print str(oldDeg[i])

if (float(len(base) * len(base[0])) > 0.0):
    print 'Reconstruct Rate : \n \t' + str(reconstructRate) + '/' + str(float(len(base) * len(base[0]))) + \
          ' = ' + str(reconstructRate/float(len(base) * len(base[0])))

前処理や出力が多いだけで,操作だけ見るとかなり簡単です.

とにかく,Python

  • 簡単に
  • 簡潔に
  • わかりやすく

プログラムが書ける素晴らしい言語です.ちょっとした課題ならぜひ使うべき.簡潔とわかりやすいって同じじゃね?って思ったあなた,大事なことなので二回言ったんです.
僕なんてちょっとした四則演算のためにインタプリタ起動しちゃうぐらいです(ぇ


次に,Kalmanフィルタをクラスとして実装したものです.
これはC++で書きました.なんでC++かというと,お察しください.
Boostライブラリが必要です.あと,yanoの日記で公開されている逆行列行列式を計算するコードも利用しました.勝手に拝借してごめんなさい.そしてありがとうございます.
デフォルトでboost::numeric::ublasにはないってことを今回はじめて知ったよ.
Wikipediaのカルマンフィルタの項の設定例が解けるってことは確認しましたが,Kalmanフィルタの理解が充分ではないのでそれ以外には直接使えないクソグラムかもしれません.

#pragma once

#include <cstdlib>
#include <ctime>
#include <cmath>
#include <iostream>
#include <string>
#include <vector>
#include <boost/random.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include "ublasdeterminant.hpp"
#include "ublasinvert.hpp"

using namespace std;
using namespace boost;
using namespace boost::numeric::ublas;

typedef matrix<double> dmat;

class KalmanFilter{
private:
  int accSigma, mesureSigma;
  double R;
  dmat *A,*B,*Q,*H;
  dmat KalmanGain;
public:
  KalmanFilter(dmat*,dmat*,dmat*, int const);
  KalmanFilter(dmat*,dmat*,dmat*,dmat*, double const, int const, int const);
  ~KalmanFilter();
  std::vector<dmat> Filtering(dmat*, dmat*, dmat*, dmat*, double const);
  std::vector<dmat> Filtering(dmat*, dmat*, dmat*, dmat*);
  std::vector<dmat> Correct(dmat*, dmat*, dmat*, dmat*);
  std::vector<dmat> Predict(dmat*, dmat*, double const);
  double calcBelief(dmat*, dmat*, dmat*);
};

KalmanFilter::KalmanFilter(dmat* A, dmat* B, dmat* Q, int const accSigma){
  this->A	   = A;
  this->B	   = B;
  this->Q	   = Q;
  this->accSigma    = accSigma;
}

KalmanFilter::KalmanFilter(dmat* A, dmat* B, dmat* Q, dmat* H, double const R, int const accSigma, int const mesureSigma){
  this->A	   = A;
  this->B	   = B;
  this->Q	   = Q;
  this->H	   = H;
  this->R	   = R;
  this->accSigma    = accSigma;
  this->mesureSigma = mesureSigma;
}

KalmanFilter::~KalmanFilter(){
}

std::vector<dmat> KalmanFilter::Filtering(dmat* State, dmat* Cov, dmat* Ctrl, dmat* G, double const Mesure){
  std::vector<dmat> CorrectResult = KalmanFilter::Correct(State, Cov, Ctrl, G);
  std::vector<dmat> PredictResult = KalmanFilter::Predict(&(CorrectResult.at(0)), &(CorrectResult.at(1)), Mesure);

  return PredictResult;
}

std::vector<dmat> KalmanFilter::Filtering(dmat* State, dmat* Cov, dmat* Ctrl, dmat* G){
  mt19937			mesureGen(static_cast<unsigned long>(time(0)));
  normal_distribution<> mesureDst(0, mesureSigma);
  variate_generator<mt19937&, normal_distribution<> > GaussNoise(mesureGen, mesureDst);
  double                const	Mesure = (prod(*H, *State))(0,0) + GaussNoise();

  return KalmanFilter::Filtering(State, Cov, Ctrl, G, Mesure);
}

std::vector<dmat> KalmanFilter::Correct(dmat* State, dmat* Cov, dmat* Ctrl, dmat* G){
  dmat			prioriState(State->size1(), State->size2());
  dmat			prioriCov(Cov->size1(), Cov->size2());
  std::vector<dmat>	Result;
  mt19937		accGen(static_cast<unsigned long>(time(0))); 
  normal_distribution<> accDst(0, accSigma);
  variate_generator<mt19937&, normal_distribution<> > Acceleration(accGen, accDst);

  prioriState = prod(*A, *State) + prod(*B, *Ctrl) + (*G) * Acceleration();
  prioriCov   = prod(*A, *Cov);
  prioriCov   = prod(prioriCov, trans(*A)) + *Q;
  
  Result.push_back(prioriState);
  Result.push_back(prioriCov);

  return Result;  
}

std::vector<dmat> KalmanFilter::Predict(dmat* prioriState, dmat* prioriCov, double const Mesure){
  dmat			posterioriState(prioriState->size1(), prioriState->size2());
  dmat			posterioriCov(prioriCov->size1(), prioriCov->size2());
  std::vector<dmat>	Result;
  
  KalmanGain = prod(*H, *prioriCov);
  KalmanGain = prod(*prioriCov, trans(*H)) / (prod(KalmanGain, trans(*H))(0,0) + R);
  posterioriState = *prioriState + KalmanGain * (Mesure - (prod(*H, *prioriState))(0,0));
  posterioriCov = identity_matrix<double>(posterioriCov.size1()) - prod(KalmanGain, *H);
  posterioriCov = prod(posterioriCov, *prioriCov);

  Result.push_back(posterioriState);
  Result.push_back(posterioriCov);

  return Result;
}

double KalmanFilter::calcBelief(dmat* Info, dmat* State, dmat* Cov){
  dmat invertCov;
  ublasinvert::invert(*Cov, invertCov);
  dmat Gap = prod(trans(*Info - *State), invertCov);
  return 1 / sqrt(ublasdeterminant::determinant(2.0 * M_PI * (*Cov))) * exp(-0.5 * prod(Gap, (*Info - *State))(0,0));
}

ついでに,ドライバも.

#include <cstdlib>
#include <cmath>
#include <iostream>
#include <string>
#include <vector>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include "KalmanFilter.hpp"

int main(void);

int main (){
  int const		accSigma    = 1;
  int const		mesureSigma = 10;
  double const		R	    = 10.0;
  double const		deltaT	    = 1.0;
  double const		Mesure	    = 5.0;
  dmat			A(2,2), B(2,2), G(2,1), Q(2,2), H(1,2);
  dmat			State(2,1), Cov(2,2), Ctrl(2,1);
  std::vector<dmat>	Result;

  A(0,0) = 1;A(0,1) = deltaT;
  A(1,0) = 0;A(1,1) = 1;

  B(0,0) = 0;B(0,1) = 0;
  B(1,0) = 0;B(1,1) = 0;

  G(0,0) = deltaT * deltaT * 0.5;
  G(1,0) = deltaT;
  Q = accSigma * prod(G, trans(G));
  
  H(0,0) = 1;H(0,1) = 0;

  State(0,0) = 0;
  State(1,0) = 0;

  Cov(0,0) = 1;Cov(0,1) = 0;
  Cov(1,0) = 0;Cov(1,1) = 1;

  Ctrl(0,0) = 0;
  Ctrl(1,0) = 0;

  KalmanFilter kf = KalmanFilter(&A, &B, &Q, &H, R, accSigma, mesureSigma);

  cout << "Initial State : " << State << "\nInitial Cov : " << Cov << endl;
  
  for(int i = 0; i < 5; ++i){
    if(i < 4){
      Result = kf.Filtering(&State, &Cov, &Ctrl, &G);
    }else{
      Result = kf.Filtering(&State, &Cov, &Ctrl, &G, Mesure);
    }

    State = Result.at(0);
    Cov = Result.at(1);
    cout << "Posteriori State : " << State << "\nPosteriori Cov : " << Cov << endl;
  }

  return 0;
}

予測だけ使うとかもできます.

#include <cstdlib>
#include <cmath>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <boost/lexical_cast.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include "KalmanFilter.hpp"

int main(int, char**);

int main (int argc, char **argv){
  int const		accSigma    = 1;
  int const             t = (argc == 2)?(lexical_cast<int>(argv[1])):(4);
  double const		deltaT	    = 1.0;
  dmat			A(2,2), B(2,2), G(2,1), Q(2,2);
  dmat			State(2,1), Cov(2,2), Ctrl(2,1);
  std::vector<dmat>	Result;

  A(0,0) = 1;A(0,1) = deltaT;
  A(1,0) = 0;A(1,1) = 1;

  B(0,0) = 0;B(0,1) = 0;
  B(1,0) = 0;B(1,1) = 0;

  G(0,0) = deltaT * deltaT * 0.5;
  G(1,0) = deltaT;
  Q = accSigma * prod(G, trans(G));
  
  State(0,0) = 0;
  State(1,0) = 0;

  Cov(0,0) = 1;Cov(0,1) = 0;
  Cov(1,0) = 0;Cov(1,1) = 1;

  Ctrl(0,0) = 0;
  Ctrl(1,0) = 0;

  KalmanFilter kf = KalmanFilter(&A, &B, &Q, accSigma);

  for(int i = 0; i < t; ++i){
    string fileName = "KFCorrect(k=" + lexical_cast<string>(i + 1) + ").txt";
    ofstream ofs(fileName.c_str());
    Result = kf.Correct(&State, &Cov, &Ctrl, &G);
    State = Result.at(0);
    Cov = Result.at(1);

    for(double location = -20.0; location <= 20.0; location += 0.5){
      for(double vel = 0.0; vel <= 0.0; vel += 0.5){
	dmat Info(2,1);
	Info(0,0) = location;
	Info(1,0) = vel;
	double prob = kf.calcBelief(&Info, &State, &Cov);
	
	ofs << location << "\t" << vel << "\t" << prob << endl;
      }
    }
  }

  return 0;
}

といった感じです.

とにかくC++

  • 冗長
  • 冗長
  • 冗長

ってことを痛感しましたね(
いいところもありますよ.低級に近い高級言語(shared_pointerとかC++11とかの機能使えばかなり高級かも・・・)なんでとにかく計算が早い.本当に早い.Pythonには早さが足りない!(PyPyあるけどNumPyPyがまだ使い辛い).

あと,今更だけどboost::numeric::ublasでは

ってことを知りました.余因子とかあの辺りの知識があれば自分で実装できるのかもだけど,素晴らしいコードが既に公開されていたのでそっち使っちゃった///


時間が許せば,せっかくなんでもう2つぐらい課題とにらめっこして作りたいけど・・・まずは今出来ているところまでで提出できる形にしないとね.