gap> START_TEST("GBNP test26");
gap> ######################### BEGIN COPYRIGHT MESSAGE #########################
gap> # GBNP - computing Gröbner bases of noncommutative polynomials
gap> # Copyright 2001-2010 by Arjeh M. Cohen, Dié A.H. Gijsbers, Jan Willem
gap> # Knopper, Chris Krook. Address: Discrete Algebra and Geometry (DAM) group
gap> # at the Department of Mathematics and Computer Science of Eindhoven
gap> # University of Technology.
gap> #
gap> # For acknowledgements see the manual. The manual can be found in several
gap> # formats in the doc subdirectory of the GBNP distribution. The
gap> # acknowledgements formatted as text can be found in the file chap0.txt.
gap> #
gap> # GBNP is free software; you can redistribute it and/or modify it under
gap> # the terms of the Lesser GNU General Public License as published by the
gap> # Free Software Foundation (FSF); either version 2.1 of the License, or
gap> # (at your option) any later version. For details, see the file 'LGPL' in
gap> # the doc subdirectory of the GBNP distribution or see the FSF's own site:
gap> #
https://www.gnu.org/licenses/lgpl.html
gap> ########################## END COPYRIGHT MESSAGE ##########################
gap>
gap> ### filename = "exampleColagen.g"
gap> ### authors Cohen & Gijsbers
gap>
gap> ### THIS IS A GAP PACKAGE FOR COMPUTING NON-COMMUTATIVE GROBNER BASES
gap>
gap> # <#GAPDoc Label="ExampleColagen">
gap> # <Section Label="ExampleColagen"><Heading>The cola gene puzzle</Heading>
gap> # <P/>
gap> # A prize question appearing in the January 2005 issue of
gap> # the Dutch journal "Natuur en Techniek" asked for a DNA change
gap> # of cows so that they could produce Cola instead of milk.
gap> # A team of genetic manipulators
gap> # has tools to perform the following
gap> # five DNA string operations. (Here the strings before and
gap> # after the equality sign can be interchanged at will.)
gap> # <P/>
gap> # operation 1: TCAT = T;
gap> # <P/>
gap> # operation 2: GAG = AG;
gap> # <P/>
gap> # operation 3: CTC = TC;
gap> # <P/>
gap> # operation 4: AGTA = A;
gap> # <P/>
gap> # operation 5: TAT = CT.
gap> # <P/>
gap> # The first question is to show how they can transform the
gap> # milk gene TAGCTAGCTAGCT to the
gap> # cola gene CTGACTGACT.
gap> # <P/>
gap> # A second question is to show that mad cow disease related
gap> # retro virus CTGCTACTGACT can be avoided at all times.
gap> # Can this be guaranteed?
gap> # <P/>
gap> # We answer these questions using the trace functions of the
gap> # Gröbner basis package GBNP in Section <Ref Sect="tracefun"/>.
gap> # <P/>
gap> # First load the package and set the standard infolevel <Ref
gap> # InfoClass="InfoGBNP" Style="Text"/> to 0 and the time infolevel <Ref
gap> # Func="InfoGBNPTime" Style="Text"/> to 0 to minimize the printing of
gap> # data regarding the computation (for more information about the info
gap> # level, see Chapter <Ref Chap="Info"/>).
gap>
gap> # <L>
gap> LoadPackage("gbnp", false);
true
gap> SetInfoLevel(InfoGBNP,0);
gap> SetInfoLevel(InfoGBNPTime,0);
gap> # </L>
gap>
gap>
gap> # We introduce the free algebra <C>ALG</C>
gap> # on the generators corresponding to the four letters in the DNA code
gap> # and express the milk gene and cola gene as monomials in this algebra.
gap>
gap> # <L>
gap> ALG:=FreeAssociativeAlgebraWithOne(Rationals, "A", "C", "G", "T");;
gap> g:=GeneratorsOfAlgebra(ALG);;
gap> A:=g[2];;
gap> C:=g[3];;
gap> G:=g[4];;
gap> T:=g[5];;
gap>
gap> milk := T*A*G*C*T*A*G*C*T*A*G*C*T;;
gap> cola := C*T*G*A*C*T*G*A*C*T;;
gap> # </L>
gap>
gap>
gap> # We next enter the set <M>K</M> of binomials <M>x-y</M>
gap> # corresponding to the five operations <M>x = y</M> listed above.
gap> # We enter the binomials as members of <C>ALG</C> and
gap> # let <C>KNP</C> be the corresponding set of NP polynomials.
gap>
gap> # <L>
gap> rule1 := T*C*A*T - T;;
gap> rule2 := G*A*G - A*G;;
gap> rule3 := C*T*C - T*C;;
gap> rule4 := A*G*T*A - A;;
gap> rule5 := T*A*T - C*T;;
gap> K := [rule1,rule2,rule3,rule4,rule5];;
gap> KNP := List(K,x-> GP2NP(x));;
gap> # </L>
gap>
gap> # We stipulate how the variables will be printed and print <C>KNP</C>.
gap> # See <Ref Func="PrintNPList" Style="Text"/>.
gap>
gap> # <L>
gap> GBNP.ConfigPrint("A","C","G","T");
gap> PrintNPList(KNP);
TCAT - T
GAG - AG
CTC - TC
AGTA - A
TAT - CT
gap> # </L>
gap>
gap>
gap> # Now calculate the usual Gröbner basis with <Ref Func="SGrobner"
gap> # Style="Text"/>.
gap>
gap> # <L>
gap> GB := SGrobner(KNP);;
gap> # </L>
gap>
gap>
gap> # Compare milk and cola after taking their strong normal forms
gap> # with respect to <C>GB</C>
gap> # using
gap> # <Ref Func="StrongNormalFormNP" Style="Text"/>.
gap> # We observe that they have the same normal form and so there is
gap> # a way to transform the milk gene into the cola gene.
gap>
gap> # <L>
gap> milkNP := GP2NP(milk);;
gap> colaNP := GP2NP(cola);;
gap> milkRed := NP2GP(StrongNormalFormNP(milkNP,GB),ALG);
(1)*T
gap> colaRed := NP2GP(StrongNormalFormNP(colaNP,GB),ALG);
(1)*T
gap> # </L>
gap>
gap> # But this information does not yet show us how to perform the transformation.
gap> # To this end we calculate the Gröbner basis with trace information,
gap> # using the function
gap> # <Ref Func="SGrobnerTrace" Style="Text"/>.
gap>
gap> # <L>
gap> GTrace := SGrobnerTrace(KNP);;
gap> # </L>
gap>
gap> # The full trace can be printed with <Ref Func="PrintTraceList" Style="Text"/>,
gap> # but we only print the relations (and no trace) by invoking <Ref
gap> # Func="PrintNPListTrace" Style="Text"/>.
gap>
gap> # <L>
gap> PrintNPListTrace(GTrace);
CT - T
GA - A
AGT - AT
ATA - A
TAT - T
TCA - TA
gap> # </L>
gap>
gap>
gap> # In order to display a proof that <M>milk-cola</M>
gap> # belongs to the ideal
gap> # we use <Ref Func="StrongNormalFormTraceDiff" Style="Text"/>,
gap> # The result is a record, <C>p</C> say, containing
gap> # <C>milk-cola</C> in the field <C>p.pol</C>
gap> # (the normal form will be subtracted from the argument <C>milk-cola</C>
gap> # to obtain <C>p.pol</C>, but the normal form is zero
gap> # because the argument belongs to the ideal generated by <C>K</C>).
gap> # The other field of the record <C>p</C> is
gap> # <C>p.trace</C>, the traced polynomial, which is best displayed
gap> # by use of <Ref Func="PrintTracePol" Style="Text"/>.
gap>
gap> # <L>
gap> p := StrongNormalFormTraceDiff(CleanNP(GP2NP(milk-cola)),GTrace);;
gap> NP2GP(p.pol,ALG);
(1)*(T*A*G*C)^3*T+(-1)*(C*T*G*A)^2*C*T
gap> PrintTracePol(p);
- TGATGAG(1) + TAGG(1)ATAT - TATATAGG(1) - TGAG(1)GACT + TGATGACG(1) - G(
1)GACTGACT - TAGCG(1)ATAT - TATAGG(1)AGT + TATATAGCG(1) + TGACG(1)GACT + CG(
1)GACTGACT - TAGG(1)AGTAGT + TAGTAGTAGG(1) + TATAGCG(1)AGT + TAGCG(
1)AGTAGT + TAGTAGG(1)AGCT - TAGTAGTAGCG(1) + TAGG(1)AGCTAGCT - TAGTAGCG(
1)AGCT - TAGCG(1)AGCTAGCT - TATG(2)TAT - TG(2)TATGAT - TGATGAG(3)AT + TAGG(
3)ATATAT - TATATAGG(3)AT - TGAG(3)ATGACT - G(3)ATGACTGACT - TATAGG(
3)ATAGT - TAGG(3)ATAGTAGT + TAGTAGTAGG(3)AT + TAGTAGG(3)ATAGCT + TAGG(
3)ATAGCTAGCT - TATG(4)T + TG(4)TAT + TATGG(4)T - TG(4)TGAT + TATATG(4)T + TGG(
4)TGAT - TG(4)TATAT + TATG(4)TAGT + TG(4)TAGTAGT + TAGG(5)ATAT - TATATAGG(
5) - TATAGG(5)AGT - TAGG(5)AGTAGT
gap> # </L>
gap>
gap> # In order to give a precise answer to the first question
gap> # we need to work a little on <C>p.trace</C>. To do so, we
gap> # introduce the following function, which creates the NP polynomial
gap> # corresponding
gap> # to the <C>i</C>-th term in the expression <C>p.trace</C>
gap> # of <C>p.pol</C> as a linear combination of members of <C>KNP</C>.
gap> # It is used to obtain the list <C>EvalList</C> of polynomials
gap> # for all <C>i</C>.
gap>
gap> # <L>
gap> EvalTracePol := function(i,p,KNP)
> local x,pi;
> pi := p.trace[i];
> x := BimulNP(pi[1],KNP[pi[2]],pi[3]);
> return [x[1],pi[4]*x[2]];
> end;;
gap>
gap> lev := Length(p.trace);;
gap> EvalList := List([1..lev], y -> CleanNP(EvalTracePol(y,p,KNP)));;
gap> # </L>
gap>
gap> # In order to find the rewrite from the milk gene to the cola gene
gap> # as required for an answer to the first question,
gap> # we match leading terms recursively.
gap>
gap> # <L>
gap> UnusedIndices := Set([1..lev]);;
gap> RunningTerm := milkNP[1][1];;
gap> stepno := 0;;
gap> NP2GP(milkNP,ALG);
(1)*(T*A*G*C)^3*T
gap> while Length(UnusedIndices) > 0 do
> i := 0;
> notfnd := true;
> while i < lev and notfnd do
> i := i+1;
> if EvalList[i][1][1] = RunningTerm and i in UnusedIndices then
> notfnd := false;
> RemoveSet(UnusedIndices, i);
> RunningTerm := EvalList[i][1][2];
> stepno := stepno+1;
> elif EvalList[i][1][2] = RunningTerm and i in UnusedIndices then
> notfnd := false;
> RemoveSet(UnusedIndices, i);
> RunningTerm := EvalList[i][1][1];
> stepno := stepno+1;
> fi;
> od;
> if i = lev and notfnd = true then Print("error not fnd in"); fi;
> Print(" -(",stepno,")- ");
> PrintNP([[p.trace[i][1]],[1]]);
> Print(" K[",p.trace[i][2],"]\n ");
> PrintNP([[p.trace[i][3]],[1]]);
> Print(" --> ");
> PrintNP([[EvalList[i][1][2]],[1]]);
> od;;
-(1)- TAGC
K[1]
AGCTAGCT
--> TAGCTAGCTAGCT
-(2)- TAG
K[3]
ATAGCTAGCT
--> TAGTCATAGCTAGCT
-(3)- TAG
K[1]
AGCTAGCT
--> TAGTAGCTAGCT
-(4)- TAGTAGC
K[1]
AGCT
--> TAGTAGCTAGCT
-(5)- TAGTAG
K[3]
ATAGCT
--> TAGTAGTCATAGCT
-(6)- TAGTAG
K[1]
AGCT
--> TAGTAGTAGCT
-(7)- TAGTAGTAGC
K[1]
1
--> TAGTAGTAGCT
-(8)- TAGTAGTAG
K[3]
AT
--> TAGTAGTAGTCAT
-(9)- TAGTAGTAG
K[1]
1
--> TAGTAGTAGT
-(10)- TAG
K[1]
AGTAGT
--> TAGTAGTAGT
-(11)- TAG
K[3]
ATAGTAGT
--> TAGTCATAGTAGT
-(12)- TAGC
K[1]
AGTAGT
--> TAGCTAGTAGT
-(13)- TAG
K[5]
AGTAGT
--> TAGCTAGTAGT
-(14)- T
K[4]
TAGTAGT
--> TATAGTAGT
-(15)- TATAG
K[1]
AGT
--> TATAGTAGT
-(16)- TATAG
K[3]
ATAGT
--> TATAGTCATAGT
-(17)- TATAGC
K[1]
AGT
--> TATAGCTAGT
-(18)- TATAG
K[5]
AGT
--> TATAGCTAGT
-(19)- TAT
K[4]
TAGT
--> TATATAGT
-(20)- TATATAG
K[1]
1
--> TATATAGT
-(21)- TATATAG
K[3]
AT
--> TATATAGTCAT
-(22)- TATATAGC
K[1]
1
--> TATATAGCT
-(23)- TATATAG
K[5]
1
--> TATATAGCT
-(24)- TATAT
K[4]
T
--> TATATAT
-(25)- T
K[4]
TATAT
--> TATATAT
-(26)- TAG
K[5]
ATAT
--> TAGCTATAT
-(27)- TAGC
K[1]
ATAT
--> TAGCTATAT
-(28)- TAG
K[3]
ATATAT
--> TAGTCATATAT
-(29)- TAG
K[1]
ATAT
--> TAGTATAT
-(30)- T
K[4]
TAT
--> TATAT
-(31)- TAT
K[4]
T
--> TATAT
-(32)- TAT
K[2]
TAT
--> TATAGTAT
-(33)- TATG
K[4]
T
--> TATGAT
-(34)- T
K[4]
TGAT
--> TATGAT
-(35)- T
K[2]
TATGAT
--> TAGTATGAT
-(36)- TG
K[4]
TGAT
--> TGATGAT
-(37)- TGATGA
K[1]
1
--> TGATGAT
-(38)- TGATGA
K[3]
AT
--> TGATGATCAT
-(39)- TGATGAC
K[1]
1
--> TGATGACT
-(40)- TGA
K[1]
GACT
--> TGATGACT
-(41)- TGA
K[3]
ATGACT
--> TGATCATGACT
-(42)- TGAC
K[1]
GACT
--> TGACTGACT
-(43)- 1
K[1]
GACTGACT
--> TGACTGACT
-(44)- 1
K[3]
ATGACTGACT
--> TCATGACTGACT
-(45)- C
K[1]
GACTGACT
--> CTGACTGACT
gap> NP2GP(colaNP,ALG);
(1)*(C*T*G*A)^2*C*T
gap> # </L>
gap>
gap>
gap> # And now the second question regarding the retro virus.
gap>
gap> # <L>
gap> retro := C*T*G*C*T*A*C*T*G*A*C*T;;
gap> # </L>
gap>
gap>
gap> # We compute the Strong Normal Form
gap> # <Ref Func="StrongNormalFormNP" Style="Text"/>
gap> # of <C>retro</C> with respect to <C>GB</C>.
gap> # As it is <C>TGT</C>, distinct to <C>T</C>, the strong normal form of milk,
gap> # there is no transformation from milk to retro.
gap> # <L>
gap> NP2GP(StrongNormalFormNP(CleanNP(GP2NP(retro)),GB), ALG);
(1)*T*G*T
gap> # </L>
gap>
gap>
gap> # Of course, here too we can verify the reduction,
gap> # by computing <Ref Func="StrongNormalFormTraceDiff" Style="Text"/>
gap> # with input the NP polynomial corresponding to <C>retro</C>
gap> # and with respect to <C>K</C>;
gap> # it is called <C>retroTrace</C>.
gap> # The symbol <C>G</C> in expression like <C>G(2)</C>
gap> # are not to be confused with the single symbols <C>G</C>
gap> # representing the DNA element.
gap>
gap> # <L>
gap> retroTrace := StrongNormalFormTraceDiff(CleanNP(GP2NP(retro)),GTrace);;
gap> PrintTracePol(retroTrace);
TGG(1) - TGC^2G(1) - TGTAG(1) + TGTACG(1) + TGTAGG(1)AT + TGTATGAG(
1) + TGTAG(1)GACT - TGTAGCG(1)AT - TGTATGACG(1) + TGG(1)ACTGACT - TGTACG(
1)GACT + G(1)GCTACTGACT - TGCG(1)ACTGACT - CG(1)GCTACTGACT + TGTATG(
2)TAT + TGG(3)AT + TGCG(3)AT - TGTAG(3)AT + TGTAGG(3)ATAT + TGTATGAG(
3)AT + TGTAG(3)ATGACT + TGG(3)ATACTGACT + G(3)ATGCTACTGACT + TGTG(
4)T + TGTATG(4)T - TGTG(4)TAT - TGTATGG(4)T + TGCG(5) + TGG(5)AT - TGTAG(
5) + TGTAGG(5)AT
gap> # </L>
gap>
gap> # </Section>
gap> # <#/GAPDoc>
gap>
gap> STOP_TEST("test26.g",10000);
