Когда Диофант встречается с генетикой

Как-то давно, когда я изучал Python, я делал одну забавную вещь, которая иллюстрирует одну замечательную идею, привнесенную в программирование из биологии. Но, к сожалению, я потерял исходный код этой небольшой программки да и не планировал когда-либо публиковать эксперимент, и поэтому как всегда, в этой статье я попытаюсь начать с нуля и показать результат в своем любимом языке программирования…В этой статье я покажу один свой крайне любительский эксперимент в области генетических алгоритмов и решения некоторого типа уравнений в целых числах. Также представленная здесь реализация не является полной и даже сколько-нибудь надежной, поскольку представлена будет здесь скорее как учебный пример, нежели реальная концепция.

Для начала я расскажу, что собой представляют генетические алгоритмы.

Генетические алгоритмы — это целый класс алгоритмов эвристического типа, которые используют в своей реализации случайный подбор, мутацию и комбинирование. Такого рода алгоритмы используются в задачах оптимизации и моделирования, где точность и строгость результата не играет значительной роли, а простой перебор очень трудоемок и неэффективен. Работа любого генетического алгоритма состоит из нескольких этапов: подготовка популяции (т.е генерация некоторого количества «решений» задачи), селекция (т.е аналог естественного отбора), выбор родителей (т.е выбор тех решений из полученного множества решений, к которым в дальнейшем будут применяться остальные операции), скрещивание (т.е рекомбинирование случайным образом двух выбранных «решений» из популяции) и мутации (т.е случайное изменение некоторого выбранного «решения»). Множество случайно сгенерированных  «решений» далее для удобства мы будем именовать популяцией, само решение — особью. При этом, каждая особь будет иметь некоторый набор параметров «решения» (которые мы и будем искать с помощью алгоритма), который мы обычно называется хромосомой. Важным элементом для генетического алгоритма является способ оценки того, насколько подходит «решение» для конечного потребителя алгоритма, и этот способ заключается в оценке особей через их хромосомы с помощью некоторого признака. Этот признак мы будем называть функцией приспособленности, которая обычно представляет собой некоторое выражение, которое использует хромосому как набор параметров и превращает их в некоторый количественный признак.

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

Диофантово уравнение — это уравнение в целых числах, которое представляет собой некоторый многочлен с целочисленными коэффициентами и также целочисленными решениями. Обычно, такие уравнения представляются, если мне не изменяет память, в приведенном виде, когда число за знаком равенства переносится в левую часть, и тогда набор решений обращает уравнение в ноль.

Примечание: я не претендую на математическую строгость такого описания, поэтому всех дотошных и заинтересованных в вопросе, прошу обратиться к подходящим источникам. Начать можно с формального изложения сути дифоантовых уравнений здесь и оттуда же подчерпнуть набор ссылок для дальнейшего знакомства с проблемой.

Уравнение, которое мы будем решать, выглядит следующим образом:

a + 2 * b + 3 * c + 4 * d = 30

Соответственно, особью будет некоторый набор параметров, которые будут подставляться в наше уравнение. Хромосомой (лучше все-таки сказать, что геномом) будет набор параметров a,b,c,d; которые для некоторых операций в генетическом алгоритме удобнее представить в виде динамического алгоритма.

Сама особь может быть описана с помощью такого класса:

import std.stdio;
import std.random;
import std.string;

class Individual
{
	private
	{
	   int[] genome;
	}

	this(int a, int b, int c, int d)
	{
		genome ~= a;
		genome ~= b;
		genome ~= c;
		genome ~= d;
	}

	int[] get()
	{
		return genome;
	}

	void set(int[] genome)
	{
		this.genome = genome;
	}

	int fitness()
	{
		return (genome[0] + 2 * genome[1] + 3 * genome[2] + 4 * genome[3] - 30);
	}
}

В классе Individual (особь) кроме обычных методов set/get, созданных для соблюдения инкапсуляции, присутствует метод fitness. Этот метод служит воплощением функции приспособленности, которая в нашем случае означает насколько близок некоторый набор параметров, заключенный в особи, решению диофантового уравнения в приведенном виде: чем ближе к нулю результат функции приспособленности, тем ближе набор чисел к решению уравнения.

Теперь опишем две функции, которые помогут нам с реализацией основных концепций эволюционных алгоритмов — мутацией и кроссинговером. Мутация обозначает случайное изменение хромосомы, которое проявляется в случайном изменении одного или нескольких ее составляющих, а кроссинговер — это своеобразный «перехлест» хромосом, при котором происходит обмен участками хромосом между двумя хромосомами (именно для упрощения этой задачи мы используем внутри особи динамический массив). Но, в нашем случае, кроссинговер происходит не между двумя хромосомами, как это обычно бывает в генетике, а между хромосомами двух различных организмов, и поэтому сходство с биологией здесь самое отдаленное.

Код процедур мутации и кроссинговера:

void mutate(ref Individual lhs)
{
	auto A = lhs.get;

	auto rnd = Random(unpredictableSeed);
	auto mutateIndex = uniform(0, 3, rnd);
	A[mutateIndex] = uniform(0, 30, rnd);
	lhs.set(A);
}

void cross(ref Individual lhs, ref Individual rhs)
{
	auto A = lhs.get;
	auto B = rhs.get;

	auto rnd = Random(unpredictableSeed);
	auto crossIndex = uniform(0, 3, rnd);

	if (crossIndex == (A.length - 1))
	{
		auto tmp = A[crossIndex];
		A[crossIndex] = B[crossIndex];
		B[crossIndex] = tmp;
		lhs.set(A);
		rhs.set(B);
	}
	else
	{
		auto tmp = A[crossIndex..$];
		A[crossIndex..$] = B[crossIndex..$];
		B[crossIndex..$] = tmp;
		lhs.set(A);
		rhs.set(B);
	}
}

А теперь код процедуры скрещивания, которая генерирует новую особь со случайным набором генов, полученных из двух отдельно взятых особей:

Individual generate(ref Individual lhs, ref Individual rhs)
{
    int[] genomes;
    int[] newGenome;
    
    genomes ~= lhs.get;
    genomes ~= rhs.get;
    
    auto rnd = Random(unpredictableSeed);
    
    foreach (i; 0..4)
    {
        newGenome ~= genomes[uniform(0, genomes.length, rnd)];
    }
    
    return new Individual(newGenome[0], newGenome[1], newGenome[2], newGenome[3]);  
}

Опишем теперь саму популяцию и процессы отбора внутри популяции:

class Population
{
	private
	{
		Individual[] population;
	}

	this(int numberOfIndividuals)
	{
		auto rnd = Random(unpredictableSeed);

		foreach (e; 0..numberOfIndividuals)
		{
			population ~= new Individual(
				uniform(0, 30, rnd),
				uniform(0, 30, rnd),
				uniform(0, 30, rnd),
				uniform(0, 30, rnd)		
			);
		}
	}

	void selection(float mutationCoefficient)
	{
		auto rnd = Random(unpredictableSeed);
		auto numberOfSelections = uniform(0, population.length, rnd);
		
		foreach (e; 0..numberOfSelections)
		{
			auto indexA = uniform(0, population.length, rnd);
			auto indexB = uniform(0, population.length, rnd);
			auto lhs = population[indexA];
			auto rhs = population[indexB];
			cross(lhs, rhs);
		}
		
		foreach (e; 0..numberOfSelections / 2)
		{
            auto indexA = uniform(0, population.length, rnd);
			auto indexB = uniform(0, population.length, rnd);
			auto lhs = population[indexA];
			auto rhs = population[indexB];
			population ~= generate(lhs, rhs);
		}

		foreach (e; 0..cast(int) (numberOfSelections * mutationCoefficient))
		{
			auto indexA = uniform(0, population.length, rnd);
			mutate(population[indexA]);
		}
	}

	auto get()
	{
		return population;
	}
}

Для того, чтобы упростить саму процедуру поиска, мы ограничим диапазон значений для каждого из параметров хромосомы особи отрезком значений [0,30]. Это ограничение используется в конструкторе класса популяции при формировании популяции из случайных особей: генераторы псевдослучайных чисел уже настроены и дают числа в указанном диапазоне, а единственным параметром популяции является количество особе в ней, которое нужно сгенерировать. В методе под названием selection выполняются процедуры, которые лежат в основе естественного отбора — мутации и кроссинговер, которые применяются к случайно выбранным особям. Также, обе процедуры несколько разделены, чтобы обеспечить наилучшее перемешивание внутри популяции и отделение процедуры мутации от кроссинговера позволяет параметризовать саму мутацию внутри популяции с помощью коэффициента мутации. Данный коэффициент устанавливает в популяции долю особей с мутацией, и его можно подобрать по своему вкусу. Помимо этого, мы постоянно дополняем популяцию новыми особями с помощью процедуры generate, которая добавляет новых особей в количестве равном половине размера популяции.

Испытаем код с помощью unittest:

unittest
{
	import std.algorithm;
	auto a = new Individual(1,2,3,4);
	auto b = new Individual(0,3,6,1);
	writeln(a.fitness);
	writeln(b.fitness);
	cross(a,b);
	mutate(a);
	writeln(a.fitness);
	writeln(b.fitness);

	auto c = new Population(10000);
	c.selection(0.20);
	writeln(c.get.filter!(a => a.fitness == 0).map!(a => a.get));
}

void main()
{
	
}

Результат:

В общем, у нас получилась довольно простая экспериментальная реализация генетического алгоритма, которая позволяет понять общий принцип и идею алгоритма. Но, несмотря на то, что полученная программа работоспособна, она не лишена и недостатков, среди которых: дублирование кода, отсутствие более строгих проверок (иногда приводит к падению) и эмпирический подбор необходимых параметров алгоритма (коэффициент мутации и количество новых особей) — и все это усугубляется еще и тем, что программа учебная, и написана была очень давно…

Однако, генетические алгоритмы — это неплохой опыт, особенно, если учишь новый язык и переносишь собственный алгоритм с другого языка.

В заключение, привожу свою старую реализацию на Python2:

# -*- coding: utf-8 -*-
"""
Created on Sat May 18 00:47:03 2013
 
@author: snowcat
@description : Решение Диофантова Уравнения a + 2*b + 3*c +4*d = 30
              с помощью генетического алгоритма.
 
 
"""
 
from random import randint
 
class Genetic:
    def __init__(self,gen=None,id=None):
        '''
       Инициализация отдельной особи.
       self.gen - генотип особи
       self.id - уникальный идентификатор для кроссовера
       '''
        if gen != None:
            self.gen = gen
        else:
            self.gen = [randint(1,30) for i in range(1,5)]
           
        if id != None:
            self.id = gen
        else:
            self.id = randint(0,2)
           
    def fitness(self):
        '''
       Функция приспособленности.
       '''
        a,b,c,d = self.gen
        return (a+2*b+3*c+4*d)-30
   
    def crossover(self,f):
        '''
       Кроссовер.
       Перемешивание генотипов.        
       '''
       
        # гены исходной особи
        a,b,c,d = self.gen
        mid = self.id
       
        # гены некоторой другой особи
        aa,bb,cc,dd = f.genotype()
        fid = f.show()
       
        if mid == 0:
            if mid == fid:
                return Genetic([a,bb,cc,dd])
            else:
                return Genetic([aa,b,c,d])
        elif mid == 1:
            if mid == fid:
                return Genetic([a,b,cc,dd])
            else:
                return Genetic([aa,bb,c,d])
        else:
            if mid == fid:
                return Genetic([a,b,c,dd])
            else:
                return Genetic([aa,bb,cc,d])
                   
    def genotype(self):
        '''
       Генотип особи.
       '''
        return self.gen
 
    def show(self):
        '''
       Идентификатор кроссовера.
       '''
        return self.id
       
       
class Population:
    def __init__(self,num=10):
        '''
       Создаем популяцию из num экземпляров.
       По умолчанию - 10 особей.
       '''
        self.p = [Genetic() for i in xrange(1,num+1)]
       
    def cross(self):
        '''
       Кроссовер внутри популяции.
       '''
        from random import sample
        l = len(self.p)
        # создаем две перемешанных копии популяции
        # для выполнения кроссовера
        f = sample(self.p,l)
        m = sample(self.p,l)
        z = []
        # сам кроссовер
        for i in xrange(0,len(f)):
            z.append(f[i].crossover(m[i]))
           
        # промежуточная функция для выполнения
        # "естественного отбора"
        def check(a,b):
            if a.fitness() < b.fitness():
                return a
            else:
                return b
        # популяция после отбора
        self.p = map(check,self.p,z)
       
    def mutate(self):
        '''
       Генератор мутаций внутри популяции.
       '''
        l = len(self.p)
        num = int(l*0.4213)+1 # количество мутаций
       
        for i in xrange(1,num):
            v = randint(0,l-1) # элемент популяции нужный для мутации
            p = self.p.pop(v) # запоминаем этот элемент и удаляем его из популяции
            ex = p.genotype() # выясняем его генотип
            ind = randint(0,3) # выбираем элемент генотипа для мутации
            ex.pop(ind)   # удаляем его
            ex.insert(ind,randint(1,30)) # вставляем мутацию
            self.p.append(Genetic(ex)) # вставляем мутированную особь
       
    @property
    def view(self):
        '''
       Все особи популяции.
       '''
        return self.p
 
# популяция из 50 особей                
x = Population(150)
 
# запускаем развитие в течение 15 поколений.
for i in range(1,15):
    x.cross()
    x.mutate()
   
# отсеиваем особей с генотипами, соответствующими решениям
# диофантова уравнения    
k = [i for i in x.view if i.fitness() == 0]
 
if len(k) == 0:
    print u'Решений не обнаружено ! Попробуйте еще раз :)'
else:
    for i in k:
        print i.genotype()

Добавить комментарий