%%% ts.rfp - find an nearly optimal solution of the
%%% traveling salesman problem with a genetic algorithm

%%% the algorithm uses a population of fixed size with
%%% diploid genomes, each representing a tour through
%%% all cities. The phenotype of an individual is
%%% represented by the shorter tour. Each new generation
%%% is formed by random selection and mixing of parents
%%% and their genomes.

%%% HOW TO USE:
%%% 1. the program needs an ENORMOUS amount of memory,
%%%    allocate as many memory cells as you can
%%% 2. type test(pln,p,q,r), where pl is a map of
%%%    coordinate pairs, n is an integer, p,q,r values
%%%    between 0 and 1;
%%%    n is the population size,
%%%    p the probability for each individual that a
%%%      mutation takes place
%%%    q the probability for a crossover during a mating
%%%    r determines the probability for the selection
%%%      of the better of an individuals chromosomes
%%%      for mating
%%% 3. the length of the shortest route, the route
%%%    itself and the average route length of the
%%%    whole population is printed; results will
%%%    be successively generated by 'more` requests


%%%%%%%% test function

test(Plan,Pop_size,Mut_rate,Cross_rate,Better_rate) :-
        Map .= generate_distmap(Plan),
        Len .= /(len(Plan),2)
        &  ts(init_pop(Pop_size,Len,Map,init_list(1,Len)),
              Map,Mut_rate,Cross_rate,Better_rate).


plan1() :& [
    03.5,00.0,12.5,00.0,00.0,02.0,06.0,02.5,04.0,03.0,
    09.5,04.0,04.5,05.0,02.5,05.5,07.0,05.5,04.0,06.0,
    01.0,06.5,07.5,08.0,12.0,07.0,02.0,07.5,03.5,07.5,
    06.5,07.5,07.0,09.0,05.5,09.5,11.0,10.0,09.9,11.9,
    02.0,11.5,07.5,12.0,10.5,13.0,02.5,13.5,07.0,15.0,
    11.5,15.0,03.5,16.5,03.0,18.5,05.5,18.5,10.5,20.5].

plan2() :& [
    0.0,0.0,1.0,0.0,2.0,0.0,3.0,0.0,
    0.0,1.0,1.0,1.0,2.0,1.0,3.0,1.0,
    0.0,2.0,1.0,2.0,2.0,2.0,3.0,2.0,
    0.0,3.0,1.0,3.0,2.0,3.0,3.0,3.0].


%%%%%%%%  main function

%% generates subsequently new generations and prints
%% the shortest tour, its length and the average
%% tour length of the population

ts(Pop,_,_,_,_) :-
    indiv[BestLength,BestRoute,_,_]
      .= select_best(Pop,nth(Pop,1))
    & tup(BestLength,BestRoute,/(sum(Pop),len(Pop))).
ts(Pop,Map,Mut_rate,Cross_rate,Better_rate) :-
    New_pop .= next_generation(Pop,Map,Mut_rate,
                               Cross_rate,Better_rate)
    & ts(New_pop,Map,Mut_rate,Cross_rate,Better_rate).


%%%%%%%% activation function for a passive structure

indiv(L1,R1,L2,R2) :& indiv[L1,R1,L2,R2].


%%%%%%%% initialization routines

%% generates from a N-element sequence of
%% coordinate pairs a N*N tableau of distances
%% between each coordinate tuple using
%% euclidean metrics

generate_distmap(Tab) :& generate_distmap1(Tab,Tab).

generate_distmap1([],_) :& [].
generate_distmap1([X,Y|Rest],T) :&
        tup(generate_distmap2(X,Y,T) |
            generate_distmap1(Rest,T)).

generate_distmap2(X,Y,[]) :& [].
generate_distmap2(X,Y,[X1,Y1|Rest]) :&
        tup(sqrt(+(*(-(X,X1),-(X,X1)),*(-(Y,Y1),-(Y,Y1))))
            | generate_distmap2(X,Y,Rest)).


%% initialize a whole population of size Indnum and
%% genome length of Indsize

%% structure of an individual: indiv[L1,Chr1,L2,Chr2]
%% L1,L2 are the lengths of the tours represented by
%% the chromosomes Chr1 and Chr2. The shorter tour
%% is always at the first place.

init_pop(0,_,_,_) :& [].
init_pop(N,Len,Map,Ori_template) :-
        Chrom1 .= randomize_list(Ori_template,Len),
        Chrom2 .= randomize_list(Ori_template,Len)
        & tup(make_order(
                indiv(route_length(Chrom1,Map),Chrom1,
                      route_length(Chrom2,Map),Chrom2))
             | init_pop(1-(N),Len,Map,Ori_template)).


%% randomizing a list

randomize_list([],_) :& [].
randomize_list(List,N) :-
        NN .= 1+(random(N)),
        Elem .= nth(List,NN),
        NewList .= kill_nth(List,NN)
        & tup(Elem|randomize_list(NewList,1-(N))).


%%%%%%%% generator of successive populations

%% saves the fittest member of the older generation
%% ("elitist variant" - but take it not too seriously)
%% and generates N - 1 new members.
%% this process involves all mechanisms like mutation
%% and crossover

next_generation(Population,Map,Mut_rate,
                Cross_rate,Better_rate) :-
        H .= select_best(Population,nth(Population,1)),
        Temp_pop .= mutate(Population,Map,Mut_rate),
        Probability_range .= probability_range(Temp_pop),
        indiv[_,Test,_,_] .= H
        & tup(H|mate(Temp_pop,
                     Map,
                     len(Temp_pop),
                     Cross_rate,
                     Better_rate,
                     Probability_range,
                     len(Test))).


%% selects for each individual of a population randomly
%% a mutation operator and applicates it on the
%% genome, respectively on the route. the resulting,
%% possibly permutated, routes are merged to a new
%% population

mutate([],_,_) :& [].
mutate([First|Rest],Map,Rate) :-
        P .= random(1.0),
        <(P,Rate)
        & tup(mutate1(First,Map) |
              mutate(Rest,Map,Rate)).
mutate([First|Rest],Map,Rate) :&
        tup(First|mutate(Rest,Map,Rate)).

mutate1(indiv[_,Better,L,Worse],Map) :-
        =(random(2),1),
        Select .= random(10),
        NewBetter .= mechanism[Select](Better)
        & make_order(indiv(route_length(NewBetter,Map),
                           NewBetter,L,Worse)).

mutate1(indiv[L,Better,_,Worse],Map) :-
    Select .= random(10),
    NewWorse .= mechanism[Select](Worse)
    & make_order(indiv(L,Better,
                 route_length(NewWorse,Map),NewWorse)).


%% performs the reproduction cycle of the population.
%% The new population is generated by successive
%% selection of pairs of individuals and exchanging
%% the genome sequences. As an additional mutation
%% operator the crossover mechanism is invoked.

mate(_,_,1,_,_,_,_) :& [].
mate(Population,Map,N,Cross_rate,Better_rate,
     Probability_range,Length) :-
      P1 .= select_parent(Population,Probability_range),
      P2 .= select_parent(Population,Probability_range),
      indiv[_,R1,_,R2] .= mix(P1,P2,Better_rate),
      [R1new,R2new]
         .= cross_over([R1,R2],Cross_rate,Length),
      LR1 .= route_length(R1new,Map),
      LR2 .= route_length(R2new,Map)
      & tup(indiv[LR1,R1new,LR2,R2new] |
            mate(Population,Map,1-(N),Cross_rate,
                 Better_rate,Probability_range,Length)).


%% selects subsequences out of two routes at the
%% same position and exchanges them
cross_over([P1,P2],Rate,Length) :-
        P .= random(1.0),
        <(P,Rate)
        & cross_over1(P1,P2,1+(random(Length)),
                            1+(random(Length))).
cross_over([P1,P2],_,_) :& [P1,P2].

cross_over1(P1,P2,Pos1,Pos2) :-
        <(Pos2,Pos1)
        & cross_over1(P1,P2,Pos2,Pos1).
cross_over1(P1,P2,Pos1,Pos2) :-
        Sub1 .= get_sublist(P1,Pos1,Pos2),
        Sub2 .= get_sublist(P2,Pos1,Pos2)
        & tup(map_subseq(P2,Sub1,Pos1),map_subseq(P1,Sub2,Pos1)).


%% returns the sum of all probability weights;
%% to select a certain individual out of the population
%% each is assigned to an interval. the size of the
%% intervall is 100 times the reciprocal of the length
%% of the phenotypical route. all intervals can
%% be thought subsequential ordered on the real axis
%% beginning at 0. the larger the interval the more
%% likely the random number generator will
%% generate a value within the interval so that
%% the individual will be selected.
%% annotation: possibly the reciprocal is not the
%% most perfect weighting, there may be weights
%% that prefer the fitter individuals, but this would
%% be unnecessarily complicated

probability_range([]) :& 0.
probability_range([indiv[L,_,_,_]|Rest]) :&
        +(probability_range(Rest),/(100,L)).


%% selects an individual out of the population using
%% the interval returned by probability_range
%% (see there)

select_parent(Pop,Range) :- P .= random(Range)
        & select_parent1(Pop,P,0).

select_parent1([I],_,_) :& I.
select_parent1([indiv[L1,R1,X,Y]|_],P,Offset) :-
        <=(P,+(Offset,/(100,L1))) & indiv[L1,R1,X,Y].
select_parent1([indiv[L,_,_,_]|Rest],P,Offset) :&
        select_parent1(Rest,P,+(Offset,/(100,L))).


%% mixes randomly the genome pair of two individuals
%% means: two lists, each with two elements, are mixed

mix(indiv[L1,R1,_,_],indiv[L2,R2,_,_],Better_rate) :-
        <=(random(1.0),Better_rate)
        & make_order(indiv[L1,R1,L2,R2]).
mix(indiv[L1,P1,_,_],P2,_) :- =(1,random(2))
        & make_order(mix1([L1,P1],P2)).
mix(indiv[_,_,L1,P1],P2,_) :&
        make_order(mix1([L1,P1],P2)).

mix1([L1,P1],indiv[L2,P2,_,_]) :- =(1,random(2))
        & indiv[L1,P1,L2,P2].
mix1([L1,P1],indiv[_,_,L2,P2]) :& indiv[L1,P1,L2,P2].


%%%%%%%% common utility functions

%% mutation operators
%% operator 0: reversion of a whole sequence
%% operator 1: swapping of two ids
%% operator 2: placing an id at an another position
%% operator 3: placing a whole subsequence at
%%     an another position
%% operator N: catching of operator numbers bigger
%%     than 2 and mapping them to the first two

mechanism[0](Route) :-
        Length .= len(Route),
        Pos1 .= 1+(random(Length)),
        Pos2 .= 1+(random(Length))
        & reverse_sublist(Route,Pos1,Pos2).
mechanism[0](Route) :& Route.

mechanism[1](Route) :-
        Length .= len(Route),
        Pos1 .= 1+(random(Length)),
        Pos2 .= 1+(random(Length))
        & swap_elements(Route,Pos1,Pos2).
mechanism[1](Route) :& Route.

mechanism[2](Route) :-
        Length .= len(Route),
        Pos1 .= 1+(random(Length)),
        Pos2 .= 1+(random(Length))
        & insertion(Route,Pos1,Pos2).
mechanism[2](Route) :& Route.

mechanism[3](Route) :-
        Length .= len(Route),
        P1 .= +(random(1-(Length)),2),
        P2 .= +(random(1-(Length)),2),
        Pos1 .= min(P1,P2),
        Pos2 .= max(P1,P2),
        Head .= get_sublist(Route,1,1-(Pos1)),
        Tail .= get_sublist(Route,Pos2,Length),
        Middle .= get_sublist(Route,Pos1,1-(Pos2)),
        HT .= uni(Head,Tail),
        L .= len(HT),
        Pos .= +(random(1-(L)),2),
        H .= get_sublist(HT,1,1-(Pos)),
        T .= get_sublist(HT,Pos,L)
        & uni(H,uni(Middle,T)).
mechanism[3](Route) :& Route.

mechanism[N](Route) :- NN .= mod(N,4)
        & mechanism[NN](Route).


%% changes the order within a genome in the way
%% that the shorter route is the first element;
%% the shorter route represents also the phenotype
%% therefore only one of the genomes has to be
%% tested to determine the fitness of an indivdual

make_order(indiv[L1,R1,L2,R2]) :-
        <(L2,L1)
        & indiv[L2,R2,L1,R1].
make_order(X) :& X.


%% computes the length of a route from a given
%% table which is a 2-dimensual array of
%% real-numbers whose entries are the distances
%% from one city to another

route_length([First|Rest],Map) :&
        route_length1(First,[First|Rest],Map).

route_length1(F,[A],Map) :& nth(nth(Map,A),F).
route_length1(F,[A,B|Rest],Map) :&
        +(nth(nth(Map,A),B),route_length1(F,[B|Rest],Map)).


%% selects the individual which represents
%% the shortest route out of the whole population

select_best([],Best) :- ! & Best.
select_best([indiv[L1,R1,L2,R2]|Rest],
            indiv[L3,R3,L4,R4]) :-
        <(L1,L3)
        ! & select_best(Rest,indiv[L1,R1,L2,R2]).
select_best([_|Rest],Best) :- ! & select_best(Rest,Best).


%% returns the sum of all route lengths in a population

sum([]) :& 0.
sum([indiv[L,_,_,_]|Rest]) :& +(L,sum(Rest)).


%% embeds a sublist in a list and eliminates duplicates

map_subseq(L,[],_) :& L.
map_subseq(L,[First|Rest],Pos) :-
        NewPos .= get_pos(L,First)
        & map_subseq(swap_elements(L,Pos,NewPos),Rest,1+(Pos)).


%%%%%%%% general utility functions

%% reversion of a sublist from element no. Pos1 to Pos2

reverse_sublist(List,Pos1,Pos2) :- >(Pos1,Pos2)
        & reverse_sublist1(List,Pos2,Pos1).
reverse_sublist(List,Pos1,Pos2) :&
        reverse_sublist1(List,Pos1,Pos2).

reverse_sublist1(List,1,1) :& List.
reverse_sublist1(List,1,Pos2) :&
        reverse_sublist2(List,[],Pos2).
reverse_sublist1([First|Rest],Pos1,Pos2) :&
   tup(First|reverse_sublist1(Rest,1-(Pos1),1-(Pos2))).

reverse_sublist2(R,List,0) :& uni(List,R).
reverse_sublist2([First|Rest],L,N) :&
        reverse_sublist2(Rest,[First|L],1-(N)).

%% returns the union of two lists

uni([],L2) :& L2.
uni([First|Rest],L2) :& tup(First|uni(Rest,L2)).


%% returns nth element of a list

nth([First|_],1) :& First.
nth([_|Rest],N) :& nth(Rest,1-(N)).


%% returns length of a list

len([]) :& 0.
len([_|Rest]) :& +(len(Rest),1).


%% initializes list of length Begin - End
%% with integer values Begin ... End

init_list(End,End) :& [End].
init_list(Begin,End) :&
        tup(Begin | init_list(1+(Begin),End)).


%% replaces element of a list at position Pos
%% with Elem

replace_elem([_|Rest],1,Elem) :& [Elem|Rest].
replace_elem([First|Rest],Pos,Elem) :&
        tup(First|replace_elem(Rest,1-(Pos),Elem)).


%% swaps two elements of a list
%% at positions Pos1 and Pos2

swap_elements(List,Pos1,Pos2) :-
        A .= nth(List,Pos1),
        B .= nth(List,Pos2),
        L .= replace_elem(List,Pos1,B)
        & replace_elem(L,Pos2,A).

%% removes element at position Pos1 and
%% inserts it at position Pos2

insertion(List,Pos1,Pos2) :-
        >(Pos2,Pos1)
        & insertion1(List,Pos1,Pos2).

insertion(List,Pos1,Pos2) :& insertion2(List,Pos1,Pos2).

insertion1(List,Pos1,Pos2) :-
        Elem .= nth(List,Pos1),
        NewList .= kill_nth(List,Pos1)
        & insert(NewList,Elem,1-(Pos2)).

insertion2(List,Pos1,Pos2) :-
        Elem .= nth(List,Pos1),
        NewList .= kill_nth(List,Pos1)
        & insert(NewList,Elem,Pos2).

%% returns list without nth element

kill_nth([_|Rest],1) :& Rest.
kill_nth([First|Rest],N) :&
        tup(First|kill_nth(Rest,1-(N))).


%% returns subsequence of a list from
%% position Pos1 to pos. Pos2

get_sublist(L,1,Pos2) :& get_sublist1(L,Pos2).
get_sublist([First|Rest],Pos1,Pos2) :&
        get_sublist(Rest,1-(Pos1),1-(Pos2)).

get_sublist1([First|_],1) :& [First].
get_sublist1([First|Rest],Pos2) :&
        tup(First | get_sublist1(Rest,1-(Pos2))).


%% returns position of an element in a list

get_pos([First|_],First) :& 1.
get_pos([_|Rest],E) :& 1+(get_pos(Rest,E)).


%% inserts element into list at nth position

insert([First|Rest],Elem,1) :& [Elem,First|Rest].
insert([First|Rest],Elem,N) :&
        tup(First|insert(Rest,Elem,1-(N))).