#Andrew Baxter
#Cook.txt version 1.0, February 2nd, 2006
#Cook makes examples for the PROFILE program to handle.

Help:=proc(): print(`Cook(t,L,n,mistakes), CosmicRay(Lmer), CosmicRays(Lmer, n)`); end:

#cook(t,L,n) makes t sequences of DNA of length n with the same random L-mer in each one
Cook:=proc(t,L,n,mistakes) local ra, ra2, Lmer, T, i, j,k, position, DNA:

ra:=rand(1..4):

Lmer:=[seq(ra(), i=1..L)]: #the L-mer that will be inserted


DNA:=[seq([seq(ra(), i=1..n-L)],j=1..t)]; # the DNA before inserting the Lmers

DNA:=DNA:
ra2:=rand(0..n - L): 

for i from 1 to t do
 position:=ra2():
 DNA[i]:=[seq(DNA[i][j], j=1..position), op(CosmicRays(Lmer,mistakes)), seq(DNA[i][k], k=position+1..n-L)]: #inserts the Lmers mutated by CosmicRays
od:

DNA:
end:

#CosmicRay(Lmer) takes in list Lmer and causes a random mutation to it
CosmicRay:=proc(Lmer) local Lmer2, ra, ra2, site, mutateto:

Lmer2:=Lmer:

ra:=rand(1..nops(Lmer)):
site:=ra():

ra2:=rand(1..4):
mutateto:=ra2():

Lmer2[site]:=mutateto:

Lmer2:
end:

#CosmicRays(Lmer, n) takes in list Lmer and applies n applications of CosmicRay
CosmicRays:=proc(Lmer, n) local i, Lmer2:

Lmer2:=Lmer:

for i from 1 to n do
 Lmer2:=CosmicRay(Lmer2):
od:

Lmer2:
end: