В этой статье мы расскажем про один из интересных на наш взгляд экспериментов, которые мы давно планировали, но не могли реализовать, поскольку оригинальный код написан был на Perl, а разбираться с алгоритмом не было времени.
Автор оригинального кода довольно известен и для его блога eax.me, мы как-то делали гостевую статью про FPGA, и именно его идею мы и хотели воспроизвести в D.
Эксперимент, а именно так можно это назвать, который будет описываться здесь связан с генетическими алгоритмами…
Про генетические алгоритмы, мы уже рассказывали в статье и там все было завязано на уравнениях.
Эта статья не будет исключением, однако, мы займемся несколько иной задачей: будем аппроксимировать функции полиномами некоторой степени. Если на этом моменте, вы решили статью закрыть, то не спешите: все описанное в этой статье не выходит за пределы обычной школьной математики и не требует от вас каких-либо познаний в области высшей математики.
Давайте сразу поясним, что мы будем сегодня делать: в этой статье, мы будем искать некую формулу, которая будет описывать некоторую зависимость между одним набором данных и другим набором. Эта формула будет представлять собой полином, который принимает произвольную величину из первого набора и дает значение, которое почти совпадает с некоторым другим значением из второго набора. Если вам непонятно, то мы будем пытаться построить зависимость между набором точек X и точек Y c помощью полинома произвольный степени, которую можно выбрать.
Полином – это многочлен, то есть выражение, следующего вида:
x + a0 * x + a1 * (x ^^ 2) + a2 * (x ^^ 3) + ... + an * (x ^^ n - 1)
где a0,…,an – произвольные коэффициенты, а n – это максимальная степень полинома (многочлена).
Т.е. по сути, мы будем искать приближение некоторой функции Y = F(X) с помощью многочлена степени n.
Поскольку, степени меняются последовательно и есть даже доказательства того, что любую функциональную зависимость можно приблизить полиномом некоторой степени, то наша задача состоит лишь в том, чтобы для выбранной степени n полинома подобрать коэффициенты.
Существует множество разных методов, как найти эти коэффициенты, но в основном, как мы поняли все сводится к решению некоторой системы линейных уравнений, с учетом того, что данных может быть очень много, то такие подходы – это ресурсозатратные вычисления. Кроме того, порой не требуется высокая точность приближения, но важно получить хотя бы примерный результат.
И вот тут, можно найти своего рода компромиссное решение – использовать методы поиска, которые не ориентированы на высокую точность, и по сути своей являются эвристическими. Одним из таких методов является генетический алгоритм. Оригинальную идею, как уже говорилось ранее, реализовал автор блога eax.me Александр, который предложил подбирать эти коэффициенты с помощью генетического алгоритма.
Алгоритм, который предложил он, нам показался довольно простым и логичным, однако, как говорилось ранее, у нас были трудности понимания некоторых деталей (поскольку скрипт был написан на Perl, который никто из авторов LightHouse Software никогда не использовал) да и просто хотелось бы иметь скрипт, который работает также. Насчет последнего стоит объяснить: несмотря на казалось бы узкую специализацию такого алгоритма, он удивительно полезен во многих областях от просто оптимизации до поиска разных зависимостей между полученными с приборов данными.
Алгоритм, в общих чертах, работает следующим образом: считывается некоторый CSV-файл, в котором записано множество точек в виде x;y и загружается в некоторую табличку, которая будет применяться далее в обработке. Далее, на основании предварительно заданных констант, таких как степень полинома, максимальное отклонение полинома от аппроксимируемых значений и некоторых других, создается начальное поколение особей, к которому далее применяются методы генетического алгоритма. Особь представляет собой небольшую структуру, которая содержит в себе массив, длины равной максимальной степени многочлена и значение некоторой функции приспособленности. Функция приспособленности, в данном случае, это максимальное отклонение значения полинома (коэффициенты которого и содержит особь) от табличных значений, и это обычная разница между ходом “графиков” полинома и исходной функциональной зависимости. Далее, после формирования первоначального поколения, применяется основной алгоритм, который включает в себя скрещивание и мутацию. Скрещивание происходит просто, как взятие двух некоторых особей и генерации новой с помощью простого попарного усреднения всех значений полинома, хранимых обеими особями. Мутация осуществляется на основании взятия коэффициентов, которые хранит особь, и умножения их на некое число с последующим сложением с текущими значениями. Внутри алгоритма происходит постоянное применение скрещивания и мутаций, после чего происходит подсчет функций приспособленности для получившейся популяции (т.е набора особей) и сортировка всей популяции на основании функции приспособленности (именно для ее вычисления и нужна таблица значений). Чем функция приспособленности меньше, тем меньше отклонение от исходной функции, и следовательно, лучше результат. Для улучшения результатов некоторая часть популяции отмирает и эта та часть, у которой значение функции приспособленности больше всего. После проведения всех этих манипуляций цикл повторяется, и как правило один такой цикл называется поколением.
В оригинальном алгоритме, цикл повторяется до тех пор, пока не будет достигнуто минимальное возможное отклонение, которое устанавливается отдельной константой, однако, мы внесли в алгоритм ряд своих улучшений. Мы добавили ряд констант, которые управляют тем, в каких пределах будут генерироваться коэффициенты, а также сделали возможность установки необходимого количества поколений, которые должен пройти алгоритм в дополнение к изначальному условию минимизации отклонения от графика. Еще мы добавили возможность в алгоритме указать долю мутаций, а также ряд случайных перемешиваний самой популяции, чтобы улучшить разнообразие особей и увеличить сходимость алгоритма.
Несмотря на ряд изменений, сам алгоритм по сути остался без изменений и в нем по прежнему можно многое улучшить.
А сама реализация на D выглядит так:
import std.algorithm;
import std.conv;
import std.file;
import std.math;
import std.random;
import std.stdio;
import std.string;
enum MAXIMUM_POLYNOMIAL_POWER = 5;
enum NUMBER_OF_BEST_INDIVIDUALS = 32;
enum POINT_SELECTION_FREQUENCY = 1;
enum double MINIMAL_POLYNOMIAL_DEVIATION = 0.0000005;
enum MINIMUM_X = -4.0;
enum MAXIMUM_X = 4.0;
auto loadCSV(string filename)
{
double[double] points;
File file = File(filename);
string line;
// счетчик точек
int pointsCounter;
while ((line = file.readln.strip) != "")
{
if (++pointsCounter == POINT_SELECTION_FREQUENCY)
{
pointsCounter = 0;
auto xy = line.split(";");
if (xy.length == 2)
{
auto x = to!double(xy[0]);
auto y = to!double(xy[1]);
points[x] = y;
}
}
}
file.close;
return points;
}
struct Individual
{
double[] data;
double fitness;
}
/// Сгенерировать нулевое поколение
auto firstGeneration()
{
Individual[] population;
auto rng = Random(unpredictableSeed);
foreach (e; 0..NUMBER_OF_BEST_INDIVIDUALS)
{
double[] data;
foreach (_; 0..MAXIMUM_POLYNOMIAL_POWER+1)
{
auto p = uniform(MINIMUM_X, MAXIMUM_X, rng);
if (p.isNaN)
{
p = 0.0;
}
data ~= p;
}
population ~= Individual(data, double.max);
}
return population;
}
/// Вычислить функцию приспособленности
double calculateFitness(ref const double[] data, ref double[double] table)
{
double maxDelta = 0.0f;
foreach (x, y; table)
{
double delta = abs(y - calculatePolynomial(data, x));
maxDelta = delta > maxDelta ? delta : maxDelta;
}
return maxDelta;
}
/// Вычислить значение полинома в некоторой точке
double calculatePolynomial(ref const double[] data, double x)
{
double result = 0.0;
foreach (i, e; data)
{
result += e * (x ^^ i);
}
return result;
}
/// Запустить генетический алгоритм
auto runGeneticAlgorithm(ref double[double] pointsTable, double mutationRate = 0.5, size_t maximumGenerations = 150)
{
assert(mutationRate < 1.0 && mutationRate > 0.0);
// начальная популяция
Individual[] population = firstGeneration();
// номер поколения
size_t generationNumber = 0;
auto rng = Random(unpredictableSeed);
while (population[0].fitness >= MINIMAL_POLYNOMIAL_DEVIATION)
{
for (size_t i = 0; i < NUMBER_OF_BEST_INDIVIDUALS - 1; ++i)
{
for (size_t j = i + 1; j < NUMBER_OF_BEST_INDIVIDUALS; ++j)
{
double[] tmp;
// скрещивание
for (size_t k = 0; k < MAXIMUM_POLYNOMIAL_POWER + 1; k++)
{
tmp ~= (population[i].data[k] + population[j].data[k]) * 0.5;
}
// мутация
if (rng.dice(1.0 - mutationRate, mutationRate))
{
foreach (ref e; tmp)
{
e += (e + uniform(MINIMUM_X, MAXIMUM_X, rng));
}
}
population ~= Individual(tmp);
}
}
// вымирание предков
population = population[NUMBER_OF_BEST_INDIVIDUALS..$];
// перемешивание
population.randomShuffle(rng);
// вычисляем функции приспособленности
foreach (ref e; population)
{
e.fitness = calculateFitness(e.data, pointsTable);
}
// сортировка по значению функции приспособленности
sort!((a, b) => a.fitness < b.fitness)(population);
// отбираем лучших
population = population[0..NUMBER_OF_BEST_INDIVIDUALS];
// перемешивание
population.randomShuffle(rng);
// выводим номер поколения и значение целевой функции у "лучшего" потомка
writef("%d;%f;%f;", generationNumber, population[0].fitness, population[$-1].fitness);
// выводим формулу для D
string formula = "";
foreach (i, e; population[0].data)
{
if (i == 0)
{
formula ~= format(`%.05f`, e);
}
else
{
auto w = format(`%.05f * (x ^^ %d)`, e, i);
formula ~= (e >= 0) ? (" + " ~ w) : w.replace("-", " - ");
}
}
// пишем формулу
write(formula);
writeln;
//
if (generationNumber == maximumGenerations )
{
return population[0];
}
// увеличиваем номер поколения
generationNumber++;
}
return population[0];
}
void main()
{
auto w = loadCSV(`test.csv`);
auto p = runGeneticAlgorithm(w, 0.76, 2048);
}Как видите, сам алгоритм несложен и состоит из нескольких довольно простых функций, почти каждая из которых прокомментирована. Сам алгоритм требует для своей работы заранее подготовленный CSV-файл, который мы делали следующим скриптом:
import std.math;
import std.stdio;
void main()
{
File file;
file.open(`test.csv`, `w`);
for (float i = 0.0f; i < 2 * PI; i += 0.000001)
{
file.writeln(i, `;`, sin(i));
}
file.close;
}Здесь дан пример для синуса, но вы можете заменить его вычисление на свою функцию, а также можете изменить шаги вычисления и диапазон начальных значений (также, можно просто отказаться от вычислений и подать свой, заранее подготовленный файл).
В функцию runGeneticAlgorithm подаются заранее подготовленная таблица значений (делается функцией loadCSV), количество поколений и доля мутации в популяции, а результатом выполнения ее является набор формул, которые последовательно выводятся в стандартный вывод вместе с номером поколения и значением функции приспособленности.
Сам формат вывода данного алгоритма выглядит следующим образом:
номер поколения;функция приспособленности лучшего потомка;функция приспособленности худшего потомка;формула полинома для лучшего потомка
А вот пример формулы полинома, который был нами найден для синуса на промежутке от 0 до 2 * PI:
auto d = 0.00994 + 0.86994 * (x ^^ 1) + 0.28182 * (x ^^ 2) - 0.39955 * (x ^^ 3) + 0.08809 * (x ^^ 4) - 0.00559 * (x ^^ 5);
И для более наглядного представления того, что получилось можно сравнить два графика, которые мы сделали через онлайн-сервис Mathway:
Красный график – это исходная функция синуса, а синий график – это график того полинома, чтио мы нашли. Чем ближе график полинома проходит от графика самой функциональной зависимости, тем точнее значение полученной формулы для значений из файла!
Это не единственный пример и вы можете сами испытать приведенный в статье код и найти свои примеры аппроксимации известных функций – им все это с помощью всего лишь сложений и умножений. Конечно, можно многое можно улучшить, к примеру, найти способы выхода алгоритма из стационарного состояния или улучшить вычисление полинома через схему Горнера.
На этом все, и если у вас найдутся идеи, где это применить или как это улучшить, то пишите их в комментариях.