(* Skrypt demonstrujący aproksymację funkcji  i możliwość REGRESJI SYMBOLICZNEJ algorytmem genetyczno-ewolucyjnym.
A. Odrzywolek, odrzywolek@th.if.uj.edu.pl
UWAGI OGÓLNE:
1. Program napisany i działa w MATHEMATICE 7.0.1
2. Program obciąża Mathematice do granic mozliwości, i nie sprawdza prawie żadnych mozliwych błędów. Należy spodziewać się regularnych zawieszeń się kernela, Visual C++ library i błędów spowodowanych brakiem pamięci. Nie odpowiadam za skutki uboczne! Zapisywane jest tylko aktualne
najlepsze przybliżone rozwiązanie problemu, tzw. ,,apex predator''. Skutkiem zawieszenia się kernala jest więc ,,masowe wymieranie''. Ponowny restart powoduje wygenerowanie całej populacji, prócz ,,apex predatora'' od nowa losowo.
*)

(*
TODO:
-wyeliminowanie wieszania się, o ile nie jest to wina błędów w kernelu
-optymalizacja szybkości działania
-przepisanie np. w C
-równoleglizacja ewolucji
-rozwiązanie podproblemu identyfikacji stałych wymiernych, algebraicznych itp np. Sqrt[1+Exp[3 Pi/2]] na podstawie przyblizonej wartości numerycznej (algorytm genetyczny lub przeszukiwanie wszystkich możliwych drzew typu TreeForm[Sqrt[1+Exp[3 Pi/2]]] )
*)

(* 
Lista prymitywnych funkcji 2 zmiennych ("aminokwasów", "peptydów") z których składamy docelowe funkcje aproksymujace
("białka", "polipeptydy").
Podobnie jak w genetyce, wybór i liczba tych funkcji jest problemem wyglądającym na nierozwiązywalny. Z jednej strony aby wygenerować dowolną funkcję ciągłą n zmiennych wystarczy dodawanie, mnożenie przez liczbę RZECZYWISTĄ oraz JEDNA dowolnie wybrana nieliniowa funkcja jednej zmiennej (tw. A. Gorbana)
i ich składanie. Z drugiej strony wygenerowanie nawet prostych funkcji za pomocą np. tylko eksponenty i logarytmu, aczkolwiek możliwe, jest karkołomną procedurą i prowadzi do bardzo długich wyrażeń; podobnie jest z liczbami wymiernymi. Dodatkowo używanie liczb rzeczywistych (przestępnych) nie jest możliwe w komputerze, poza pojedynczymi przypadkami jak Pi, E itp. gdyż jest ich nieprzeliczalnie wiele.
Dla przypomnienia, cała złożoność życia na Ziemi pochodzi od składania 20 peptydów, więc mam nadzieję  że (być może) podobnie jest z funkcjami. Przyjęcie tej hipotezy (częściowo podpartej liczbą f. jednej zmiennej w
tw. Kołmogorowa o składaniu f. n-zmiennych z f. jednej zmiennej, równej 2n+2) nadal jednak nic nie mówi JAKIE to miałyby być funkcje. Pozostaje liczyć na to, że skuteczność procedury w dużej mierze jest niezależna od tego CO składamy, oraz że funkcje (Sin, Exp, itd.) i operacje (Plus, Times, Power) które my z powodów historycznych nazywamy ,,elementarnymi'' nie są jakoś specjalnie patologiczne do tych celów, lub wręcz są bardzo dobre. Za tym ostatnim poglądem przemawia mnogość problemów rozwiązana za pomocą tychże funkcji w 19 i 20 wieku.
Lista poniżej pochodzi z obserwacji:
a) struktury wyrazen zapisywanych przez Mathematice (zob. np. a/(1+Exp[-b x]) //TreeForm ) (P,F,H)
b) skuteczności f. sigmoidalnych, ułamków łańcychowych i f. wymiernych w aproksymacji (K, A, W)
c) konieczności tworzenia stałych wymiernych (J,M)
d) działania suwaka logarytmicznego (B,E)
e) dopuszczenia rozw. periodycznych (S)
Nie pozostaje nic innego jak eksperymentowanie.
*)

S[x_,y_]:=Sin[x]+Cos[y]
A[x_,y_]:=x/(1+y)
J[x_,y_]:=1
P[x_,y_]:= x+y
M[x_,y_]:=x-y
W[x_,y_]:=x/(1+y^2)
B[x_,y_]:=Exp[x]+Exp[-y]
F[x_,y_]:=x*y
G[x_,y_]:=x/y
H[x_,y_]:=x^y
K[x_,y_]:=y/(1+Exp[-x])
L[x_,y_]:=Log[x,y]

(* Lista zmiennych niezależnych, np. dla problemu r. rożniczkowego zwyczajnego jest to jedna zmienna *)

truevars={X}

(* Lista parametrów których wartości będą optymalizowane numerycznymi metodami optmalizacji lokalnej, jeżeli jest ich zbyt dużo (>7) to łatwiej znajdujemy dopasowanie o zadanej dokładności, ale tracimy możliwość regresji symbolicznej, jeżeli jest ich zbyt mało (w szczególności zero) to opieramy się na czystej regresji symbolicznej (parametry liczbowe powstają wyłącznie z funkcji J poprzez jej wielokrotne składanie). Nie ma wtedy mozliwości znalezienia współczynników nie będących liczbami wymiernymi lub algebraicznymi. Zbieżność procedury jest też wtedy bardzo mała. 3 parametry to kompromis pomiędzy czystą regresją symboliczną a  zwykłym fitem wieloparametrowym, który zwykle pasuje do wszystkiego.
*)

pars={a,b,c}

(* ostatecznie wszystkie zmienne (niezależne i parametry) są podstawiane do funkcji prymitywnych *)

vars=Join[truevars,pars]

(* Lista funkcji uzywanych w regresji symbolicznej. Dla prostoty wszystkie są f. dwóch zmiennych, aczkolwiek wydaje się, że wystarczy jedna f. dwóch zmiennych. Jeżeli z góry wiemy, że pewne funkcje nie są nam potrzebne, to ułatwimy sobie pracę wyrzucjąc je z listy poniżej.
*)

funs={P,M,B,F,H,J,S}

(* Ostatecznie ,,DNA'' składa się z ,,kodonów'' w postaci zmiennych niezależnych, parametrów oraz funkcji prymitywnych. Listę kodonów interpretujemy jako wyrażenie w RPN (Odwrotna Notacja Polska). Np:
{X,X,A} oznacza A[X,X], {X,X,A,X,B} oznacza B[A[X,X],X]. Najłatwiej sprawdzić to funkcją  translate[] poniżej bez zdefiniowanych funkcji prymitywnych. NIE KAŻDY ciąg ,,kodonów'' odpowiada sensownej funkcji. Dlatego procedura transkrypt[] naprawia dowolny ciąg tak aby odpowiadał pewnej funkcji. Nie ma jednoznacznej metody takiej naprawy, jest ona dokonywana losowo, ale tak aby zachować sensowną składnię. Po ,,transkrypcji'' otrzymujemy listę ,,kodonów'' będącą poprawnym wyrażeniem RPN, które określamy jako ,,RNA''.
*)

transkrypt[DNA_]:=Module[{position,n,dna},
dna=DNA;
If[Length[dna]==0,AppendTo[dna,RandomChoice[vars]]];
If[Length[dna]==1,AppendTo[dna,RandomChoice[vars]];AppendTo[dna,RandomChoice[vars]]];
If[Length[dna]==2,AppendTo[dna,RandomChoice[funs]]];
While[Length[dna]>2&&Not[MemberQ[vars,dna[[1]]]&&MemberQ[vars,dna[[2]]]],dna=Drop[dna,1]] ;
If[Length[dna]==2&&MemberQ[vars,First[dna]]&&MemberQ[funs,Last[dna]],PrependTo[dna,RandomChoice[vars]]];
If[Length[dna]==2&&MemberQ[funs,First[dna]]&&MemberQ[vars,Last[dna]],PrependTo[dna,RandomChoice[vars]];PrependTo[dna,RandomChoice[vars]];AppendTo[dna,RandomChoice[funs]]];
If[Length[dna]==2&&MemberQ[funs,First[dna]]&&MemberQ[funs,Last[dna]],PrependTo[dna,RandomChoice[vars]];PrependTo[dna,RandomChoice[vars]];PrependTo[dna,RandomChoice[vars]]];
position=2;
While[position<Length[dna]&&Sum[If[MemberQ[vars,dna[[n]]],1,-1],{n,1,position}]>1,position++] ;
dna=Take[dna,position];
While[Sum[If[MemberQ[vars,dna[[n]]],1,-1],{n,1,Length[dna]}]>1, AppendTo[dna,RandomChoice[funs]]] ;
Append[dna,RandomChoice[vars]]
]

(* ,,RNA'' jest gotowe do ,,translacji'' czyli przetłumaczenia na funkcję. Odbywa się to prostym interpreterem RPN poniżej. Wynikiem działania funkcji translate[] jest pewna funkcja zmiennych niezależnych i parametrów *)

translate[RNA_]:=
Module[{stos={},p,rna},
rna=RNA;
While[Length[rna]>0,
If[
MemberQ[vars,First[rna]],
AppendTo[stos, First[rna]];rna=Delete[rna,1],
p=First[rna][stos[[-2]],stos[[-1]]];
stos=Delete[stos,-1];
stos=Delete[stos,-1];
AppendTo[stos,p];
rna=Delete[rna,1]]
];
First[stos]
]

(* Później będziemy potrzebować jeszcze procedury replikacji funkcji. Początkowo ,,RNA'' miało być wklejane w ,,DNA'', ale ostatecznie zwracamy po prostu ,,RNA'' , które staje się od tego momentu ,,DNA'' funkcji potomnej. Zabezpiecza to przed nadmiernym (faktycznie usuwa w całości) rozrostem tzw. niekodujacego lub inaczej ,,śmieciowego'' DNA. To jest jeszcze istotny punkt do przemyślenia, gdyż w przyrodzie występuje ,,śmieciowy'' DNA, aczkolwiek jego rola w skuteczności działania mechanizmów genetycznych jest mętna. Konieczność wycięcia wynikła w praktyce, gdyż funkcje broniły się przed wpływem mutacji poprzez rozrost śmieciowego DNA do list o rozmiarach przekraczających milion,
podczas gdy sekwencje kodujące miały długość 3-10. Wydaje się, że gdyby częstość mutacji rosła z rozmiarem DNA ten problem można by rozwiązać, ale nie próbowałem.
*)

replicate[DNA_,RNA_]:=
If[LongestCommonSubsequence[DNA,Most[RNA]]==Most[RNA],RNA,RNA,RNA]

(* ,,DNA'' podlega mutacjom. Funkcja poniżej losowo podmienia jeden kodon. *)

mutant[DNA_]:=ReplacePart[DNA,RandomInteger[{1,Length[DNA]}]->RandomChoice[Union[vars,funs]]]

(* Funkcja poniżej ,,krzyżuje'' dwa DNA, rozrywając i sklejając je w losowych miejscach *)

hybryda[DNA1_,DNA2_]:=Join[Partition[DNA1,RandomInteger[{2,Length[DNA1]-2}]]//First,Partition[DNA2,RandomInteger[{2,Length[DNA2]-2}]]//Last]

(* 
Funkcja fitness[] poniżej jest kluczowa i decydująca. Szereguje ona funkcje od najlepiej radzącej sobie z zadanym problemem,
do najgorszej. Nie ma jednoznacznego sposobu na zdefiniowanie tej funkcji, a  poniżej jest tylko bardzo prosty przykład. Funkcja próbuje mianowicie znaleźć numerycznie fit do pewnych dancych liczbowych, a następnie sprawdza jakość tego fitu sumując różnice kwadratów. Dodatkowo nakładana jest ,,kara'' dla funkcji stałych i liniowych w truevars. Jest tu nadal spore pole do inwencji. Mozna próbować, np. wstawiać funkcję do równania, użyć Simplify a nastepnie kategoryzować ,,jakość'' licząc np. LeafCount czy znaki ASCII w wyniku.
Można próbować liczyć granicę w Infinity, aby preferować funkcje zachowujące  się ,,lepiej'' asymptotycznie itp. itd. Można też mieszać różne kryteria. Istotne jest zachowanie możliwości ,,wspinania się'' na ,,szczyt nieprawdopodobieństwa'' małymi kroczkami: kryteria typu TAK/NIE są złe, musimy ilościowo (najlepiej liczbami zmiennoprzecinkowymi, nie całkowitymi)
określić ,,jakość'' funkcji.  
*)

fitness[data_,fun_]:=MemoryConstrained[TimeConstrained[
Module[{f,vals,delta,ii,dopasowanie},
dopasowanie=FindFit[data, Re[fun], Transpose[{pars,RandomReal[{-2^RandomInteger[],1},Length[pars]]}],truevars];
f=If[NumericQ[fun /.dopasowanie /.  X->Sqrt[Pi]],fun /. dopasowanie,0,0];
vals=Transpose[data]//First;
delta=(Transpose[data]//Last) - Table[f ,{X,vals}];
(Sum[Re[delta[[ii]]]^2,{ii,1,Length[vals]}]//Sqrt)/(Length[vals](Max[vals]-Min[vals]))+
If[D[f,X]==0,100000,0,0]+If[D[f,X,X]==0,100000,0,0]
],60,300000],500*1024^2,400000]

(* Przykładowe równanie które chcemy rozwiązać *)

eq=D[y[X],X,X]-y'[X]+Sin[X]^2==0

(* Rozwiązanie analityczne dla niecierpliwych, algorytm nie potrzebuje oczywiście tej wiedzy i nigdzie z niej nie korzysta*)

analSOL=y/.DSolve[{eq,y[0]==0,y'[0]==1},y,X][[1]]

(* Rozwiązanie numeryczne *)

sol=y /.First[NDSolve[{eq,y[0]==0,y'[0]==1},y,{X,-100,100}]]

(* Losowo wygenerowane wartości rozwiązania powyżej, które będziemy usiłowali odtworzyć *)

goalDAT1=({#,sol[#]}&)@RandomReal[{-1,1},20]//Transpose;
goalDAT2=({#,sol[#]}&)@RandomReal[{-10,10},100]//Transpose;

goalDAT =Join[goalDAT1,goalDAT2];

(* Początkowa populacja losowych ciągów ,,DNA''. Długość ,,DNA'' oraz ilość ,,gatunków'' są losowe w sensownyc granicach.  *)

popinit=pop=newpop=Table[{,Table[RandomChoice[{RandomChoice[truevars],RandomChoice[funs],RandomChoice[pars],RandomChoice[funs]}],{ii,1,1+2^RandomInteger[{0,3}]}],},{ii,1,RandomInteger[{100,1000}]}];

(* Pierwszy z listy jest (chwilowo) uznany na najlepszy gatunek *)

apexRNA=popinit[[1,2]]//transkrypt

translate[apexRNA]

initFIT=fitness[goalDAT,%]

(* Odczytujemy ,,DNA'' z dysku jeżeli kontynuujemy rachunek *)

apex=If[FileExistsQ["apex.m"],Import["apex.m"],{initFIT,apexRNA,apexRNA}]

best=%[[1]]

oldbest=10^100.

gen=0

(* Pętla poniżej ,,ewoluuje'' w najprostszy sposób populację funkcji. Ewolucja jest dosyć prymitywna. W każdym pokoleniu eliminujemy najgorzej przystosowanych, kontrolujemy populację ,,top predatora'' i dokonujemy mutacji i krzyżówek. "Apex predator" jest nietykalny, co daje nam złudzenie monotonicznej zbieżności algorytmu. Dzięki temu jeżeli nie uda się regresja symboliczna problemu, to na pocieszenie zostaje nam fit, zwykle o imponującej dokładności i zachwujący asymptotykę.

Pętlę można przerwać w dowolnym momencie z menu Evaluation/Abort i sprawdzić wizualnie jakość fitu.
*)

While[best>10^(-16),
pop=newpop;
num=Dimensions[pop][[1]];
If[best<oldbest,
Print[
gen, "\t",num, "\t",best,"\n","APEX PREDATOR:\t",apex,"\n\t",translate[apex[[3]]]
]
];
If[best<oldbest,Export["apex.m",apex]];
gen++;
Print["GENERATION=",gen,"\tPOPULATION=",num] ;
(* transkrypcja i dodanie RNA *)
popT=Table[{pop[[ii,1]],pop[[ii,2]],pop[[ii,2]]//transkrypt},{ii,1,num}];
(* translacja RNA na funkcje *)
popL=Table[{Evaluate[translate[popT[[ii,3]]]],popT[[ii,2]],popT[[ii,3]]},{ii,1,num}];
(* obliczenie dopasowania każdego osobnika i posortowanie *)
popFIT=Table[{fitness[goalDAT,popL[[ii,1]]],popT[[ii,2]],popT[[ii,3]]},{ii,1,num}] //Sort;
(* sprawdzamy czy pojawił się nowy apex predator, jeżeli nie stary nadal rzadzi *)
apex=If[popFIT[[1,1]]<apex[[1]],popFIT[[1]],apex,apex];
oldbest=best;
best=If[apex[[1]]<oldbest,apex[[1]],oldbest,oldbest];
newgen=Table[{popFIT[[ii,1]],replicate[popFIT[[ii,2]],popFIT[[ii,3]]],},{ii,1,RandomInteger[{1,num}]}];
newgen=TakeWhile[Reverse[newgen],#[[1]]≠best&];
mutants=Table[{,mutant[newgen[[RandomInteger[{1,Dimensions[newgen][[1]]}],2]]],},{ii,1,RandomInteger[{1,20}]}];
hybrids=Table[{,hybryda[newgen[[RandomInteger[{1,Dimensions[newgen][[1]]}],2]],newgen[[RandomInteger[{1,Dimensions[newgen][[1]]}],2]]],},{ii,1,RandomInteger[{1,20}]}];
apexmutants=Table[{,mutant[apex[[2]]],},{ii,1,RandomInteger[{1,20}]}];
apexhybrids=Table[{,hybryda[newgen[[RandomInteger[{1,Dimensions[newgen][[1]]}],2]],apex[[2]]],},{ii,1,RandomInteger[{1,20}]}];
randoms=Table[{,Table[RandomChoice[{RandomChoice[truevars],RandomChoice[funs],RandomChoice[pars],RandomChoice[funs]}],{ii,1,1+2^RandomInteger[{0,4}]}],},{ii,1,RandomInteger[{1,10}]}];
newpop=Join[newgen,mutants,hybrids,apexmutants,apexhybrids,randoms];
PrependTo[newpop,apex];
]

apex[[3]]//translate

symb=% /.FindFit[goalDAT,%,pars,truevars]

fitness[goalDAT,%]

symb//TreeForm

(* W przykładach poniżej zakres zmiennej X zwiększono 2x w stosunku do obszaru fitowania *)

Plot[{symb,sol[X]},{X,-20,20},PlotStyle->{Red,Black}] (* Porównanie fitu i rozw. analitycznego *)

Plot[symb-sol[X],{X,-20,20},PlotLabel->"Błąd absolutny"]

Plot[(symb-sol[X])/sol[X],{X,-20,20}, PlotLabel->"Błąd względny"]