«Техническая» часть АРСС моделирования

В предыдущей статье я уделил внимание исключительно алгоритму и «физике» метода. В этой статье будет описано немного про получение результатов и их представлении.

0. Содержание:

1. Программирование как оно было
2. Представление результатов

1. Программирование как оно было:

Целью были результаты, потому ничего странного не было в том, что бы выбрать готовый продукт который уже умеет АРСС. Выбор был сделан в пользу oldCronos. Причинами, сподвигнувшими меня на такой выбор были:

• Распространение под GNU GPL что развязывает руки, в известных пределах, конечно же.
• C++ как основной язык.
• Использование GSL для сложным мат. методов — методы Левенберга-Марквардта, Гаусса-Жордана etc.
• Неплохая архитектура — извлечь в независимый модуль всю модель не составило труда.

Для всех пар из множества была посчитана функция максимального правдоподобия.

Все расчёты проводились на маленьком «кластере» из 16 ПК на базе AMD Athlon 64 X2 4400+ c 2 Гб оперативной памяти. Впрочем, это и всё что понадобилось.

Для 120 дней (а это в общей сложности 480 точек для чистых данных и 120 для прореженных) и для 9 существенных точках на гринвичском меридиане имеем 800 моделей. Расчёты заняли существенное время — 1 неделя и 3 дня.

Как же всё это было сделано, для начала подключили все нужные header-файлы

#include <boost/mpi.hpp>
namespace mpi = boost::mpi;
#include <boost/serialization/access.hpp>

#include <setjmp.h>		// jmp_buf, longjmp, setjmp

#include <gsl/gsl_errno.h>	// GSL: gsl_set_error_handler

// oldCronos headers
#include "ARMAModel.hpp"
#include "TimeSeries.hpp"

boost.serialization нужен для передачи через boost.mpi собственноручно написанного простенького класса, который будет хранить порядки модели и значение.

class task {
private:
	friend class boost::serialization::access;
	
	template<class Archive>
	void serialize(Archive &ar, const unsigned int version) {
		ar &p;
		ar &q;
		ar &llh;
	}
	
public:
	int p;
	int q;
	double llh;
	
	task() {}
	explicit task(int p, int q, double llh)
		: p(p), q(q), llh(llh)
	{}
};

Следующим важным моментом стало задание своего обработчика ошибок для GSL. Это необходимо чтобы GSL не делал abort() при расходимости методов оптимизации — в частности, используемого ЛМ метода.

jmp_buf pos;
int p;
int q;
bool error;

void my_handler(const char *reason, const char *file, int line, int gsl_errno)
{
    std::cout << "t" << p << " " << q << " # " << file << ":" << line << ' ' << reason << std::endl;
    error = true;
    longjmp(pos, 1);
}

auto old_handler = gsl_set_error_handler(&my_handler);	// setting own gsl_error handler
…
setjmp(pos);
…
gsl_set_error_handler(old_handler);	// returning default gsl error handler

Создать контейнер с набором не составляет проблем, потому сразу отправляем задания на удалённые машины и ждём ответ. И снова отправляем.

Полный код можно увидеть здесь.

Концом заданий является специально сформированный объект класса task, где p или q равны -1. В результате из файлы примерно такого содержимого, получаем такое.

2. Представление результатов

С помощью элементарного кода на LaTeX и используя пакет pgfplots можно получить красивые картинки (1 и 2):

documentclass{standalone}

usepackage{pgfplots}
pgfplotsset{compat=newest}

pgfplotsset{every axis/.append style={
	font=tiny,
	tiny,
	tick style={ultra thin}}}

begin{document}
	begin{tikzpicture}
		begin{axis}[height=65mm
			,width=65mm,
			scatter/@pre marker code/.code={
				pgfmathparse{100-33.3*ln(pgfplotspointmetatransformed+1)}
				letopacity=pgfmathresult
				defmarkopts{mark=*,color=red!opacity,fill=red!opacity}
				expandafterscopeexpandafter[markopts]}
			,scatter/@post marker code/.code={endscope}
			,scatter src=explicit
			,xlabel={$p$}
			,ylabel={$q$}]

			addplot[only marks,scatter] table[x index=0, y index=1,meta index=2] {<data>};
		end{axis}
	end{tikzpicture}
end{document}

С помощью ImageMagick, или, к примеру, этого сервиса были получены картинки в формате png.

Рис. 1: Максимальное правдоподобие для 45°S, 0°.

Значение для пары сначала линейно отображается на шкалу от 0 до 1000, а затем элементарно преобразовывается формулой .

Рис. 2: Максимальное правдоподобие для 0°, 0°. Прореженные данные.

А с помощью простенького приложения на Qt можно посмотреть IONEX формат и сохранить изображение на память.

А с помощью LaTeX это всё можно приукрасить:

documentclass{standalone}

usepackage{pgfplots}

begin{document}
	begin{tikzpicture}
		begin{axis}[height=45mm
				,width=45mm
				,scale only axis
				,enlargelimits=false
				,axis on top
				,point meta min=0
				,point meta max=6.5e100
				,xtick={-120,0,120}
				,ytick={-60,0,60}
				,xlabel={$lambda,,{}^circ$}
				,ylabel={$varphi,,{}^circ$}
				,xlabel shift={-1mm}
				,ylabel shift={-5mm}]

			addplot graphics[xmin=-180,xmax=180,ymin=-87.5,ymax=87.5] <image>};
		end{axis}
	end{tikzpicture}
end{document}

И в результате получим pdf с картинкой, осями и подписями.

Автор: m0nhawk