kashiの日記

動機

プログラム

#ifndef MATPLOTLIB_HPP
#define MATPLOTLIB_HPP

#include <cstdio>

class matplotlib {
	FILE *p;

	public:

	bool open() {
		p = popen("python -c 'import code; import os; import sys; sys.stdout = sys.stderr = open(os.devnull, \"w\"); code.InteractiveConsole().interact()'", "w");
		if (p == NULL) return false;
		send_command("import matplotlib.pyplot as plt");
		send_command("import matplotlib.patches as patches");
		send_command("fig, ax = plt.subplots()");
		send_command("plt.show(block=False)");
		return true;
	}

	bool close() {
		send_command("plt.close()");
		send_command("quit()");
		if (pclose(p) == -1) return false;
		return true;
	}

	void screen(double x1, double y1, double x2, double y2, bool EqualAspect = false) const {
		if (EqualAspect == true) {
			fprintf(p, "ax.set_aspect('equal')\n");
		}
		fprintf(p, "plt.xlim([%.17f,%.17f])\n", x1, x2);
		fprintf(p, "plt.ylim([%.17f,%.17f])\n", y1, y2);
		fprintf(p, "plt.draw()\n");
		fflush(p);
	}

	void line(double x1, double y1, double x2, double y2, const char *color = "blue", const char *opt = "") const {
		fprintf(p, "ax.draw_artist(ax.plot([%.17f,%.17f],[%.17f,%.17f], color='%s', %s)[0])\n", x1, x2, y1, y2, color, opt);
		fprintf(p, "fig.canvas.blit(ax.bbox)\n");
		fprintf(p, "fig.canvas.flush_events()\n");
		fflush(p);
	}

	void point(double x, double y, const char *color = "blue", const char *opt = "") const {
		fprintf(p, "ax.draw_artist(ax.scatter([%.17f],[%.17f], color='%s', %s))\n", x, y, color, opt);
		fprintf(p, "fig.canvas.blit(ax.bbox)\n");
		fprintf(p, "fig.canvas.flush_events()\n");
		fflush(p);
	}

	void rect(double x1, double y1, double x2, double y2, const char *edgecolor = "blue", const char *facecolor = NULL, const char *opt = "") const {
		
		if (facecolor == NULL) {
			fprintf(p, "ax.draw_artist(ax.add_patch(patches.Rectangle(xy=(%.17f,%.17f), width=%.17f, height=%.17f, fill=False, edgecolor='%s', %s)))\n", x1, y1, x2-x1, y2-y1, edgecolor, opt);
		} else {
			fprintf(p, "ax.draw_artist(ax.add_patch(patches.Rectangle(xy=(%.17f,%.17f), width=%.17f, height=%.17f, fill=True, edgecolor='%s', facecolor='%s', %s)))\n", x1, y1, x2-x1, y2-y1, edgecolor, facecolor, opt);
		}
		fprintf(p, "fig.canvas.blit(ax.bbox)\n");
		fprintf(p, "fig.canvas.flush_events()\n");
		fflush(p);
	}

	void ellipse(double cx, double cy, double rx, double ry, const char *edgecolor = "blue", const char *facecolor = NULL, const char *opt = "") const {
		if (facecolor == NULL) {
			fprintf(p, "ax.draw_artist(ax.add_patch(patches.Ellipse(xy=(%.17f,%.17f), width=%.17f, height=%.17f, fill=False, edgecolor='%s', %s)))\n", cx, cy, rx*2, ry*2, edgecolor, opt);
		} else {
			fprintf(p, "ax.draw_artist(ax.add_patch(patches.Ellipse(xy=(%.17f,%.17f), width=%.17f, height=%.17f, fill=true, edgecolor='%s', facecolor='%s', %s)))\n", cx, cy, rx*2, ry*2, edgecolor, facecolor, opt);
		}
		fprintf(p, "fig.canvas.blit(ax.bbox)\n");
		fprintf(p, "fig.canvas.flush_events()\n");
		fflush(p);
	}

	void circle(double cx, double cy, double r, const char *edgecolor = "blue", const char *facecolor = NULL, const char *opt = "") const {
		ellipse(cx, cy, r, r, edgecolor, facecolor, opt);
	}

	void polygon(double *x, double *y, int n, const char *edgecolor = "blue", const char *facecolor = NULL, const char *opt = "") const {
		int i;

		fprintf(p, "ax.draw_artist(ax.add_patch(patches.Polygon((");
		for (i=0; i<n; i++) {
			fprintf(p, "(%.17f,%.17f),", x[i], y[i]);
		}
		
		if (facecolor == NULL) {
			fprintf(p, "), fill=False, edgecolor='%s', %s)))\n", edgecolor, opt);
		} else {
			fprintf(p, "), fill=True, edgecolor='%s', facecolor='%s', %s)))\n", edgecolor, facecolor, opt);
		}
		fprintf(p, "fig.canvas.blit(ax.bbox)\n");
		fprintf(p, "fig.canvas.flush_events()\n");
		fflush(p);
	}

	void save(const char *filename) const {
		fprintf(p, "plt.savefig('%s')\n", filename);
		fflush(p);
	}

	void send_command(const char *s) const {
		fprintf(p, "%s\n", s);
		fflush(p);
	}

	void clear() const {
		send_command("plt.clf()");
	}
};

#endif // MATPLOTLIB_HPP

#include "matplotlib.hpp"

int main()
{
	matplotlib g;

	// initialize
	g.open();

	// set drawing range
	g.screen(0, 0, 10, 10);
	// aspect ratio = 1
	// g.screen(0, 0, 10, 10, true);

	g.line(1,1,3,4);
	g.line(1,2,3,5, "red");

	g.rect(4,6,5,8);
	g.rect(6,6,7,8, "green");
	g.rect(8,6,9,8, "black", "red");

	g.point(4,2);
	g.ellipse(8,2,2,1);
	g.circle(2,8,2);

	double xs[5] = {6., 7., 6., 5., 5.};
	double ys[5] = {4., 5., 6., 6., 5.};
	g.polygon(xs, ys, 5, "black", "yellow");

	g.line(1,3,3,6, "green", "alpha=0.2");
	g.line(1,4,4,5, "black", "alpha=0.5, linestyle='--'");
	g.point(4, 3, "red", "s=100");

	g.save("test.pdf");

	getchar();

	// finish drawing
	g.close();
}

基本

#include "matplotlib.hpp"

int main()
{
	matplotlib g;

	g.open();
	g.close();
}

描画範囲

	g.screen(0, 0, 10, 10);

	g.screen(0, 0, 10, 10, true);

線分

	g.line(1,1,3,4);

	g.line(1,2,3,5, "red");

長方形

	g.rect(4,6,5,8);

	g.rect(6,6,7,8, "green");

	g.rect(8,6,9,8, "black", "red");

点

	g.point(4,2);

楕円

	g.ellipse(8,2,2,1);

円

	g.circle(2,8,2);

多角形

	double xs[5] = {6., 7., 6., 5., 5.};
	double ys[5] = {4., 5., 6., 6., 5.};
	g.polygon(xs, ys, 5, "black", "yellow");

ファイルにセーブ

	g.save("test.pdf");

画面消去

追加オプション

	g.line(1,3,3,6, "green", "alpha=0.2");
	g.line(1,4,4,5, "black", "alpha=0.5, linestyle='--'");
	g.point(4, 3, "red", "s=100");

• alpha=0.2として薄く表示
• alpha=0.5として、更に線種を変更
• 点を大きく

技術的な細かいこと

• 最初にpopenしてるところ、単にpythonインタプリタを呼ぶのでなく妙な引数がついてますが、自分の環境ではpythonインタプリタをそのまま呼んだ場合標準入力を一行ずつ読んでくれなかったので、小細工しました。もっとスマートな方法募集中です。

• 普通に書くと新しく何かを描画する毎に全体を再描画してしまいgnuplotと同様遅くなってしまったので、blittingとかいう技を使って再描画を抑制しているつもりです。しかし、matplotlibの構造を全く知らずに適当に実験しながら作ったので、とんでもなく的外れなことをしているかもしれません。

• C++とlibpythonをリンクしてpythonインタプリタ自身をプログラム中に取り込んでしまう方法でも同じことが出来ます。しかし、コンパイル時のオプションが必要だったり気軽さが失われてしまうのと、描画の重さに比べてパイプのオーバーヘッドはあまり大きく無さそうなので、気軽なパイプで作ってみました。

オマケ

#include <boost/numeric/ublas/vector.hpp>
#include "matplotlib.hpp"

namespace ub = boost::numeric::ublas;

template <class T, class F>
void rk(F f, ub::vector<T>& init, T start, T end) {

	ub::vector<T> x, dx1, dx2, dx3, dx4, x1, x2, x3;

	T t = start;
	T dt = end - start;

	x = init;
	dx1 = f(x, t) * dt;
	x1 = x + dx1 / 2.;
	dx2 = f(x1, t + dt/2.) * dt;
	x2 = x + dx2 / 2.;
	dx3 = f(x2, t + dt/2.) * dt;
	x3 = x + dx3;
	dx4 = f(x3, t + dt) * dt;

	init = x + dx1 / 6. + dx2 / 3. + dx3 / 3. + dx4 / 6.;
}

struct Lorenz {
	template <class T> ub::vector<T> operator() (const ub::vector<T>& x, T t){
		ub::vector<T> y(3);

		y(0) = 10. * ( x(1) - x(0) );
		y(1) = 28. * x(0) - x(1) - x(0) * x(2);
		y(2) = (-8./3.) * x(2) + x(0) * x(1);

		return y;
	}
};

int main()
{
	int i, j;
	ub::vector<double> x, x_new;
	double t, t_new, dt;
	const char *colors[3] = {"blue", "red", "green"};

	matplotlib g;

	g.open();

	g.screen(0., -30., 10., 50.);

	x.resize(3);
	x(0) = 15.; x(1) = 15.; x(2) = 36.;
	t = 0.;
	dt = pow(2., -5);

	for (i=0; i < (1 << 5) * 10; i++) {
		x_new = x;
		t_new = t + dt;
		rk(Lorenz(), x_new, t, t_new);
		for (j=0; j<3; j++) {
			g.line(t, x[j], t_new, x_new[j], colors[j]);
		}
		x = x_new;
		t = t_new;
	}

	getchar();
	g.close();
}

おわりに

2019年4月10日追記

• コンパイル時にマクロKV_USE_TPFMAを定義すると、twoproductの計算にfma命令を使うようになる。
• dka.hpp (Durand Kerner Aberth法) の重解時の安定性が向上。
• 精度保証付き数値計算とkvライブラリの概要を説明するスライドを掲載