のんびりしているエンジニアの日記

ソフトウェアなどのエンジニア的な何かを書きます。

boostの行列演算ライブラリのublasの拡張ライブラリ作っています。

Sponsored Links

皆さんこんにちは
お元気ですか。私は元気です。

今日はublasについてです。ublasのvectorとかmatrixとか正直、結構不便です。
うん、絶対値を求めたベクトルを返す関数すらないし……。ということで今、これらを解決したいなーなんて思っているうちに拡張用のライブラリを作っています。

何か作ったほうがいいかなという処理があれば、教えてもらいたいなと思います。

実装している内容

  1. 絶対値
  2. Log
  3. Sqrt
  4. Power

現在はこの程度しかありません。何か思いついたら増やそうと思います。
因みに今のところこんな感じでございます。

ソースコード

ライブラリ
#pragma once
#include <iostream>
#include <vector>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/triangular.hpp> // 三角行列用のヘッダ.前進消去,後退代入に必要
#include <boost/numeric/ublas/lu.hpp>       // LU分解,前進消去,後退代入用のヘッダ
#define NDEBUG

using namespace boost::numeric;

namespace ublas_func{
	/**
	各要素一つ一つに対して、適用する
	関数です。 
	*/
	class Function{
	protected:
		double number; //番号
	public:
		virtual double calc(double element) = 0;
		void setNumber(double num){
			number = num;
		}
	};

	class AbsFunction:public Function{
	public:
		virtual double calc(double element){
			return std::abs(element);
		}
	};

	class LogFunction:public Function{
	public:
		virtual double calc(double element){
			return std::log(element);
		}
	};

	class Log10Function:public Function{
	public:
		virtual double calc(double element){
			return std::log10(element);
		}
	};

	class SqrtFunction:public Function{
	public:
		virtual double calc(double element){
			return std::sqrt(element);
		}
	};

	class PowerFunction:public Function{
	public:
		virtual double calc(double element){
			return std::pow(element, number);
		}
	};

	class ExpFunction:public Function{
	public:
		virtual double calc(double element){
			return std::exp(element);
		}
	};

	class SigmoidFunction:public Function{
	public:
		virtual double calc(double element){
			return 1.0 / (1.0 + exp(-element));
		}
	};

	class SoftmaxFunction{
	public:
		boost::numeric::ublas::vector<double> calc(const boost::numeric::ublas::vector<double> &ublas_vector){
			double maximum = 0.0;
			for(int i = 0; i < ublas_vector.size(); i++){
				if(ublas_vector[i] > maximum){
					maximum = ublas_vector[i];
				}
			}

			double total = 0.0;
			for(int i = 0; i < ublas_vector.size(); i++){
				total += exp(ublas_vector[i] - maximum);
			}

			boost::numeric::ublas::vector<double> result(ublas_vector.size());
			for(int i = 0; i < ublas_vector.size(); i++){
				result[i] = exp(ublas_vector[i] - maximum) / total;
			}

			return result;
		}
	};

	/** 関数を作る */
	template<typename T>
	BOOST_UBLAS_INLINE T abs(const T &ublas_class){
		Function *fc= new AbsFunction();
		return apply_function_to_element(ublas_class,fc);
	}

	template<typename T>
	BOOST_UBLAS_INLINE T log(const T &ublas_class){
		Function *fc = new LogFunction();
		return apply_function_to_element(ublas_class, fc);
	}

	template<typename T>
	BOOST_UBLAS_INLINE T log10(const T &ublas_class){
		Function *fc = new Log10Function();
		return apply_function_to_element(ublas_class, fc);
	}

	template<typename T>
	BOOST_UBLAS_INLINE T sqrt(const T &ublas_class){
		Function *fc = new SqrtFunction();
		return apply_function_to_element(ublas_class, fc);
	}

	template<typename T>
	BOOST_UBLAS_INLINE T pow(const T &ublas_class,double number){
		Function *fc = new PowerFunction();
		fc->setNumber(number);
		return apply_function_to_element(ublas_class, fc);
	}

	template<typename T>
	BOOST_UBLAS_INLINE T exp(const T &ublas_class){
		Function *fc = new ExpFunction();
		return apply_function_to_element(ublas_class, fc);
	}

	template<typename T>
	BOOST_UBLAS_INLINE T sigmoid(const T &ublas_class){
		Function *fc = new SigmoidFunction();
		return apply_function_to_element(ublas_class, fc);
	}

	template<typename T>
	boost::numeric::ublas::matrix<T> softmax(const boost::numeric::ublas::matrix<T> &ublas_class){
		boost::numeric::ublas::matrix<double> result(ublas_class.size1(),ublas_class.size2());
		SoftmaxFunction sf;

		for(int i = 0; i < result.size1(); i++){
			row(result,i) = sf.calc(row(ublas_class,i));
		}
		return result;
	}

	/** vector全てに適用 */
	template<typename T>
	BOOST_UBLAS_INLINE boost::numeric::ublas::vector<T> apply_function_to_element(const boost::numeric::ublas::vector<T> &ublas_vector,Function *fc){
		int size = ublas_vector.size();
		boost::numeric::ublas::vector<T> result(size);
		for(int i = 0; i < size; i++){
			result[i] = fc->calc(ublas_vector[i]);
		}
		return result;
	}

	/** matrix全てに適用 */
	template<typename T>
	BOOST_UBLAS_INLINE boost::numeric::ublas::matrix<T> apply_function_to_element(const boost::numeric::ublas::matrix<T> &ublas_matrix,Function *fc){
		boost::numeric::ublas::matrix<double> result(ublas_matrix.size1(),ublas_matrix.size2());
		for(int i = 0; i < ublas_matrix.size1(); i++){
			for(int j = 0; j < ublas_matrix.size2(); j++){
				result(i,j) = fc->calc(ublas_matrix(i,j));
			}
		}
		return result;
	}

	/**
	ここから以下は、列ごと行ごとに実行するプログラムです。
	Ex:列ごとに平均を取ったベクトルを返す、
	*/
	class FunctionPerVector{
	public:
		virtual double calc(const boost::numeric::ublas::vector<double> &ublas_vector) = 0;
	};

	class MeanFunctionPerVector : public FunctionPerVector{
	public:
		virtual double calc(const boost::numeric::ublas::vector<double> &ublas_vector){
			return sum(ublas_vector) / (double)ublas_vector.size();
		};
	};

	/** 統計 axisは列か行か 0:行 1:列 */
	template<typename T>
	boost::numeric::ublas::vector<T> apply_to_function_per_vector(const boost::numeric::ublas::matrix<T> &ublas_matrix,FunctionPerVector *fcpv,int axis){
		boost::numeric::ublas::vector<double> result;
		if (axis == 0){
			int n_rows = ublas_matrix.size1();
			result.resize(n_rows);
			for(int i = 0; i < n_rows; i++){
				double r = fcpv->calc(row(ublas_matrix,i));
				result[i] = r;
			}
		}else{
			int n_cols = ublas_matrix.size2();
			result.resize(n_cols);
			for(int i = 0; i < n_cols; i++){
				double r = fcpv->calc(column(ublas_matrix,i));
				result[i] = r;
			}
		}
		return result;
	}

	template<typename T>
	boost::numeric::ublas::vector<T> mean_by_rows(boost::numeric::ublas::matrix<T> &ublas_matrix){
		FunctionPerVector *fcpv = new MeanFunctionPerVector();
		return apply_to_function_per_vector(ublas_matrix, fcpv, 0);
	}

	template<typename T>
	boost::numeric::ublas::vector<T> mean_by_column(boost::numeric::ublas::matrix<T> &ublas_matrix){
		FunctionPerVector *fcpv = new MeanFunctionPerVector();
		return apply_to_function_per_vector(ublas_matrix, fcpv, 1);
	}

	/** -1 ~ 1までで作成するmatrix */
	static boost::numeric::ublas::matrix<double> generateMatrixRandom(int n_row, int n_col){
		boost::numeric::ublas::matrix<double> result(n_row,n_col);
		for(int i = 0; i < n_row; i++){
			for(int j = 0; j < n_col; j++){
				result(i,j) = (double)rand() / (RAND_MAX / 2.0) - 1.0;
			}
		}
		return result;
	}

	/** -1 ~ 1までで作成するvector */
	static boost::numeric::ublas::vector<double> generateVectorRandom(int size){
		boost::numeric::ublas::vector<double> result(size);
		for(int i = 0; i < result.size(); i++){
			result[i] = (double)rand() / (RAND_MAX / 2.0) - 1.0;
		}
		return result;
	}
};
実験動作
#include "UblasLibrary.hpp" //上記ライブラリ
int main(int argc, char const *argv[])
{
	boost::numeric::ublas::vector<double> v(3);
    boost::numeric::ublas::matrix<double> m(3,3);

    v(0) = 1; v(1) = -2; v(2) = 3;

    m(0,0) = 1; m(0,1) = -2; m(0,2) = 3;
    m(1,0) = 4; m(1,1) = -5; m(1,2) = 6;
    m(2,0) = 7; m(2,1) = -8; m(2,2) = 9;

    std::cout << "absのテスト" << std::endl;
    std::cout << ublas_func::abs(v) << std::endl;
    std::cout << ublas_func::abs(m) << std::endl;

    std::cout << "logのテスト" << std::endl;
    std::cout << ublas_func::log(v) << std::endl;
    std::cout << ublas_func::log(m) << std::endl;

    std::cout << "powのテスト" << std::endl;
    std::cout << ublas_func::pow(v,2.0) << std::endl;
    std::cout << ublas_func::pow(m,2.0) << std::endl;

    std::cout << "sqrtのテスト" << std::endl;
    std::cout << ublas_func::sqrt(v) << std::endl;
    std::cout << ublas_func::sqrt(m) << std::endl;
	return 0;
}

出力

absのテスト
[3](1,2,3)
[3,3]((1,2,3),(4,5,6),(7,8,9))
logのテスト
[3](0,nan,1.09861)
[3,3]((0,nan,1.09861),(1.38629,nan,1.79176),(1.94591,nan,2.19722))
powのテスト
[3](1,4,9)
[3,3]((1,4,9),(16,25,36),(49,64,81))
sqrtのテスト
[3](1,nan,1.73205)
[3,3]((1,nan,1.73205),(2,nan,2.44949),(2.64575,nan,3))

nanは負の値なので、ある意味正常ですが、exceptionとか投げるようにしておくべきかな…。少々設計に困っています。
逆行列とか擬似逆行列とかは追加しても損はないかな。ちょっと考えておきます。