% Basic MATLAB script for evolving L diallelic loci gen generations 
% Demonstrates the approach to linkage equilibrium in a model of
% recombination without interference
close all; clear all;
L=5; gen=100;
r=0.2;
alpha=define_recombination_no_interference(L,zeros(2^L,1),r);
PM=zeros(2^L-1,1); PM(end)=1/2; %P_000... = P_111... = 1/2
[index1 index2 beta]=precompute(L, alpha);
Pend(1)=PM(end);
for ngen=1:gen
    PM=P_to_M(PM);
    PM=M_to_Mp(PM,index1,index2,beta);
    PM=M_to_P(PM);
    Pend(ngen+1)=PM(end);
end
plot(0:gen,Pend); hold on; plot(0:gen,ones(gen+1)/2^L,'r--');