dir='/data/www/19genomes/variants.tables/';

filename='chr1.alleles.txt';
whole=[dir filename];
fileID = fopen(whole);
%columns chromomsome, psebp, bp, no. alleles, major AF, 19 strains
raw_founder_data = textscan(fileID, '%f %f %f %f %f %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s ');
fclose(fileID);
no_founders=19;
raw_chromosome=raw_founder_data{1};
founder_pse_bp=raw_founder_data{2};
founder_bp=raw_founder_data{3};
founder_allele_total=raw_founder_data{4};
major_A_F=raw_founder_data{5};
for i=1:19
    founder_data=horzcat(founder_data, raw_founder_data{5+i});
end
for i=1:length(raw_chromosome)
    ch=num2str(raw_chromosome(i));
    chr='chr';
    chr1=[chr ch];
    chromosome{i}=chr1;
end
dir='/Net/mus/data/www/POOLING/ARABIDOPSIS/FOUNDER/MAGIC/MAGIC.';
filename='18';
whole=[dir filename];
fileID = fopen(whole);
raw_indiv_data=textscan(fileID, '%s %f %s %s');
fclose(fileID);
indiv_chr=raw_indiv_data{1};
indiv_bp=raw_indiv_data{2};
indiv_SNP=raw_indiv_data{3};

%founder matrix creation
%delete no SNP rows
for i=length(founder_bp):-1:1 
    a=cellfun('length',founder_data(i,:));
    d=max(a');
    if d>1
        founder_data(i,:)=[];
        raw_chromosome(i,:)=[];
        founder_pse_bp(i,:)=[];
        founder_bp(i,:)=[];
        founder_allele_total(i,:)=[];
        major_A_F(i,:)=[];
    end
end
% create usible vectors and matrices
i=1;
start=1;
a=strcmp(chromosome(i),indiv_chr(i));
while a==1
    %location_vector(i)=indiv_bp(i);
    j=start;
    match=0;
    while (match==0 && j<=length(raw_chromosome))
        if founder_pse_bp(j)==indiv_bp(i)
            location_vector(i)=indiv_bp;
            for n=1:no_founders
              founder_matrix(i,n)=strcmp(founder_data(j,n),indiv_SNP(i));
            end
            match=1;
            start=j+1;
        else
            j=j+1;
        end
    end
    i=i+1;
    a=strcmp(chromosome(i),indiv_chr(i)); 
end
edim=dim(founder_matrix);
%Emission Probability
delta=0.05;
error=0.1;
E=zeros(edim(1),no_founders+1);
for i=1:edim(1)
    for j=1:no_founders+1
        if j==no_founders+1
            E(i,j)=delta;
        elseif founder_matrix(i,j)==1
                E(i,j)=1-errorval-delta;
        else
            E(i,j)=errorval;
        end
    end
end

genome_length=10000;
no_cm=100;
filestart='Transion_matrix_'
for i=1:edim(1)
    if i==1
        last_SNP=0;
    else
        last_SNP=location_vector(i-1);
    end
    recombination_rate=(location_vector(i)-last_SNP)*no_cm/genome_length;
    %Transision matrix
    T=zeros(no_founders+1,no_founders+1);
    phi=exp(-recombination_rate);
    for j=1:no_founders+1
        for k=1:no_founders+1
            if j==k
                T(j,k)=phi+(1-phi)/(no_founders+1);
            else
                T(j,k)=(1-phi)/(no_founders+1);
            end
        end
        file2=num2str(i);
        new_filename = [filestart file2 ];
        save(new_filename,'T');
    end
    end