カタベログ

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

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

Python 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つぐらい課題とにらめっこして作りたいけど・・・まずは今出来ているところまでで提出できる形にしないとね．