1Package $DesignOfExperiments
2
3-- D := [0..3, 0..2];
4-- D := [-1..1, -1..1];
5-- FacCl := [1, x, y, x^2, y^2];
6export FullDesignPoints;  -- FullDesignPoints(D);
7export FullDesignLevels;  -- FullDesignLevels(D);
8export FullDesignQB;      -- FullDesignQB(R, D);
9export FullDesignPolys;   -- FullDesignPolys(R, D);
10export corners;           -- corners(FD_QB, SmallerQB);
11export border;            -- border(FD_QB, SmallerQB);
12export GBRelations;       -- GBRelations(D, SmallerQB)
13export BBRelations;       -- BBRelations(D, SmallerQB)
14export AllGBases;         -- AllGBases(GBRelations(D,SmallerQB))
15export AllBBases;         -- AllBBases(BBRelations(D,SmallerQB))
16
17export FractionQB;        -- FractionQB(FD_QB, [pp1, pp2,...]);
18
19
20
21------------------------ [ Procedures ] ----------------------------------
22
23-- The number of parameters necessary for Corner, QB.
24
25Define MakeNumParams(Corners, QB)
26  Return sum([len([P in QB | P<T]) | T in Corners]);
27EndDefine; -- MakeNumParams
28
29
30-- For each term T_i in QB, gives the (parametric) polynomial G_i
31Define FillGB(Corners, QB)
32  P := RingOf(Corners[1]);
33  phi := CoeffEmbeddingHom(P);
34  params := apply(phi, indets(CoeffRing(P)));
35  Index := 0;
36  L := [];
37  For i := 1 To len(Corners) Do
38    K := [T In QB | T<Corners[i]];
39    append(ref L, Corners[i] + sum([params[Index+j]*K[j] | j in 1..len(K)]));
40    Index := Index + len(K);
41  EndFor;
42  Return L;
43EndDefine; -- FillGB
44
45
46Define FillBB(Border, QB)
47  P := RingOf(Border[1]);
48  phi := CoeffEmbeddingHom(P);
49  params := apply(phi, indets(CoeffRing(P)));
50  Index := 0;
51  L := [];
52  For i := 1 To len(Border) Do
53    append(ref L, record[head:=Border[i],
54			 tail:=sum([params[Index+j]*QB[j] | j in 1..len(QB)])]);
55    Index := Index + len(QB);
56  EndFor;
57  Return L;
58EndDefine; -- FillBB
59
60
61Define corners(BigQB, SmallQB)
62  return interreduced(diff(BigQB,SmallQB));
63EndDefine; -- corners
64
65
66Define border(BigQB, SmallQB)
67  IsInBorder :=
68  func(T)
69    foreach X in indets(RingOf(BigQB[1])) do
70      if IsDivisible(T,X) and (T/X isin SmallQB) then return true; endif;
71    endforeach;
72    return false;
73  endfunc;
74  return [ T in diff(BigQB, SmallQB) | IsInBorder(T) ];
75EndDefine; -- corners
76
77
78define FullDesignPoints(FD)
79  return CartesianProductList(FD);
80enddefine; -- FullDesignPoints
81
82
83define FullDesignLevels(FD)
84  return [ len(D) | D in FD];
85enddefine; -- FullDesignLevels
86
87
88Define FullDesignQB(P, FD)
89  if NumIndets(P) < len(FD) then error("Insufficient NumIndets in ring"); endif;
90  Return [MakeTerm(P, L) | L in CartesianProductList([0..len(D)-1 | D In FD])];
91EndDefine; -- FullDesign
92
93
94Define FullDesignPolys(P, FD)
95  if NumIndets(P) < len(FD) then error("Insufficient NumIndets in ring"); endif;
96  Return [product([ indet(P,i)-a | a in FD[i]]) | i in 1..len(FD)];
97EndDefine; -- FullDesignEq
98
99
100Define FractionQB(QB, Corners)
101  I := ideal(Corners);
102  Return [T in QB | not(T IsIn I)];
103EndDefine; -- CornerFullDesign
104
105
106-------------[ Point Solving ]
107
108-- Gives the list of fractions described by the specialization of Sys for
109-- every element of Sol.
110Define Points(equations, GenericEQ, ParameterSols)
111  Return [$.Point(equations, GenericEQ, Sol) | Sol in ParameterSols];
112EndDefine;
113
114
115Define Point(equations, GenericEQ, ParameterSol)
116  S := RingOf(GenericEQ[1]);
117  R := RingOf(equations[1]);
118  phiSR := PolyRingHom(S, R, PolyAlgebraHom(CoeffRing(S),RingQQ(),ParameterSol), indets(R));
119  ActualSys := concat(equations, apply(phiSR, GenericEQ));
120  Return RationalSolve(ActualSys);
121EndDefine;
122
123
124
125-----------------------[Buchberger procedures]-----------------------------
126
127
128------------[ Terms ]------------
129
130Define TermsAreCoprime(S,T)
131  L1 := exponents(S);
132  L2 := exponents(T);
133  For I := 1 To NumIndets(RingOf(T)) Do
134    if L1[I]*L2[I] <> 0 then return false; endif;
135  EndFor;
136  Return true;
137EndDefine;
138
139-----------------[ Pairs ]---------------
140
141Define BCriterion(ref Pairs, T) --> T is the LT of the new polynomial
142  Pairs := [P in Pairs | not(IsDivisible(P.lcm, T)
143                         and lcm(T,LT(P.f1)) <> P.lcm
144                         and lcm(T,LT(P.f2)) <> P.lcm) ];
145EndDefine;
146
147
148Define NewPair(F,G)
149  return Record[f1:=F,  f2:=G,  lcm:=lcm(LT(F),LT(G)),
150		IsCoprime := $.TermsAreCoprime(LT(F),LT(G))];
151EndDefine;
152
153
154Define InsertPairGM(ref Pairs, P)
155  LCM_P := P.lcm;
156  ToBeInserted := true;
157  for i := 1 to len(Pairs) do
158    LCM_i := Pairs[i].lcm;
159    if LCM_i = LCM_P Then
160      if P.IsCoprime then Pairs[i] := P; endif;
161      ToBeInserted := false;
162    elif IsDivisible(LCM_i, LCM_P) Then Pairs[i] := "remove";
163    elif IsDivisible(LCM_P, LCM_i) Then ToBeInserted := false;
164    endif;
165    if not(ToBeInserted) then break; endif;
166  endfor;
167  If ToBeInserted Then append(ref Pairs, P); EndIf;
168  Pairs := [ PPP In Pairs | type(PPP) <> STRING ];
169EndDefine; -- InsertPairGM
170
171
172-------------------[ Buchberger auxiliary ]-----------------
173
174Define SPoly(P)
175  F := P.f1;  G := P.f2; lcmP := P.lcm;
176  phi := CoeffEmbeddingHom(RingOf(F));
177  Return phi(LC(G))*(lcmP/LT(F))*F - phi(LC(F))*(lcmP/LT(G))*G;
178EndDefine;
179
180
181-------------------[ Buchberger Main ]-----------------
182
183Define BuildNewPairs(ref GB, F)
184  Pairs := [];
185  for i := 1 to len(GB) do
186    P := NewPair(GB[i], F);
187    if LT(GB[i]) = P.lcm then GB[i] := "remove"; endif;
188    $.InsertPairGM(ref Pairs, P);
189  endfor;
190  GB := [FF In GB | type(FF) <> STRING ];
191  Return [ PP In Pairs | not(PP.IsCoprime) ];
192EndDefine;
193
194
195Define UpdateBasisAndPairs(ref GB, ref Pairs, f)
196  $.BCriterion(ref Pairs, f);
197  Pairs_f := $.BuildNewPairs(ref GB, f);
198  append(ref GB, f);
199  Pairs := concat(Pairs, Pairs_f);
200EndDefine;
201
202
203Define relations_Buchberger(GB, GenericPolys)
204  Pairs := [];
205  Foreach f In GenericPolys Do
206    $.UpdateBasisAndPairs(ref GB, ref Pairs, f);
207  EndForeach;
208  PrintLn "len(Pairs) = ", len(Pairs);
209  PrintLn "GB = ", GB;
210  Out := [NR(SPoly(PP), GB) | PP in Pairs];
211  return record[GB := GB,
212	  rel := ConcatLists([coefficients(f) | f In Out and not(IsZero(f))])];
213EndDefine; -- relations_Buchberger
214
215
216Define GBRelations(FullDesign, FractionQB)
217  rings := diff(RingsOf(FractionQB), CurrentTypes());
218  if len(rings)>1 then error("more than 1 ring in FractionQB");endif;
219  P := rings[1];
220  FD_QB := FullDesignQB(P, FullDesign);
221  PrintLn "FD_QB = ", FD_QB;
222  foreach T in FractionQB do
223    if not(T IsIn FD_QB) then error(sprint(T)+" not contained in QB"); endif;
224  endforeach;
225  FractionCorners := corners(FD_QB, FractionQB);
226  PrintLn "FractionCorners = ", FractionCorners;
227  S := BuildParamRing(P, "a", MakeNumParams(FractionCorners, FractionQB));
228  phiPS := PolyRingHom(P,S, CanonicalHom(CoeffRing(P),CoeffRing(S)), indets(S));
229  GenericEqns := FillGB(apply(phiPS, FractionCorners), apply(phiPS, FractionQB));
230  PrintLn "GenericEqns = "; indent(GenericEqns);
231  FD_polys := FullDesignPolys(S, FullDesign);
232  return relations_Buchberger(FD_polys, GenericEqns);
233EndDefine; -- GBRelations
234
235
236define FindReducerBB(f, L)
237  for j:=1 to len(L) do
238    if IsDivisible(LT(f),L[j].head) then return j; endif;
239  endfor;
240  return 0;
241enddefine; -- FindReducerBB
242
243define NRBB(f, L)
244//  PrintLn "NRBB(f, L) ", f;
245  TmpRes := zero(RingOf(f));
246  if IsZero(f) then return TmpRes; endif;
247  i := FindReducerBB(f, L);
248  repeat
249    while i = 0 do
250      TmpRes := TmpRes + LM(f);      f := f - LM(f);
251      if IsZero(f) then return TmpRes; endif;
252      i := FindReducerBB(f, L);
253    endwhile;
254    f := f - (LM(f)/L[i].head)*(L[i].head + L[i].tail); -- NR
255//    PrintLn "f --> ", f;
256    if IsZero(f) then return TmpRes; endif;
257    i := FindReducerBB(f, L);
258  endrepeat;
259enddefine; -- NRBB
260
261
262Define MultiplicationMat(H, x, QB)
263//  PrintLn "QB = ", QB;
264  SortBy(ref H, func(a,b) return a.head>b.head; endfunc);
265  return transposed(mat([coefficients(NRBB(x*T, H), QB) | T in QB]));
266EndDefine; -- MultiplicationMat
267
268
269define eval(f,M)
270  v := ZeroMat(RingOf(M[1]), NumRows(M[1]), NumRows(M[1]));
271  foreach m in monomials(f) do
272    L := log(LT(m));
273    v := v+ RingElem(RingOf(v),LC(m)) * product([ M[i]^L[i] | i in 1..len(M)]);
274//    PrintLn "m = ", m;
275//    PrintLn "v = ", v;
276  endforeach;
277//  PrintLn "f = ", f;
278  return v;
279enddefine; -- eval
280
281
282define relations_MultMat(FD_polys, GenericEqnsBB, FractionQB, FractionBorder)
283  H := concat([ record[head:=LT(f), tail:=f-LT(f)] | f in FD_polys], GenericEqnsBB);
284  X := [indet(RingOf(f), UnivariateIndetIndex(f)) | f in FD_polys];
285  M := [MultiplicationMat(H, x, FractionQB) | x in X];
286//  PrintLn "M = ", M;
287  commutators := [M[i[1]]*M[i[2]]-M[i[2]]*M[i[1]] | i in subsets(1..len(M),2)];
288//  PrintLn "commutators = ", commutators;
289  ColMatrices := [ GetCol(eval(f,M),1) | f in FD_polys and IsZero(NR(LT(f), FractionBorder)) ];
290//  PrintLn "ColMatrices = ", ColMatrices;
291  AllMatEntries := flatten(concat([GetRows(N) | N in commutators],ColMatrices));
292  return record[BB := H,
293		rel := [ g in AllMatEntries | not(IsZero(g))]];
294enddefine; -- relations_MultMat
295
296
297Define BBRelations(FullDesign, FractionQB)
298  rings := diff(RingsOf(FractionQB), CurrentTypes());
299  if len(rings)>1 then error("more than 1 ring in FractionQB");endif;
300  P := rings[1];
301  FD_QB := FullDesignQB(P, FullDesign);
302//  PrintLn "FD_QB = ", FD_QB;
303  foreach T in FractionQB do
304    if not(T IsIn FD_QB) then error(sprint(T)+" not contained in QB"); endif;
305  endforeach;
306  FractionBorder := border(FD_QB, FractionQB);  -- C2
307//  PrintLn "FractionBorder = ", FractionBorder;
308  S := BuildParamRing(P, "a", len(FractionBorder)*len(FractionQB));
309  phiPS := PolyRingHom(P,S, CanonicalHom(CoeffRing(P),CoeffRing(S)), indets(S));
310  GenericEqnsBB := FillBB(apply(phiPS, FractionBorder), apply(phiPS, FractionQB));
311//  PrintLn "GenericEqnsBB = "; indent(GenericEqnsBB);
312  FD_polys := FullDesignPolys(S, FullDesign);
313//  PrintLn "FD_polys = "; indent(FD_polys);
314  return relations_MultMat(FD_polys, GenericEqnsBB,
315			   apply(phiPS,FractionQB), apply(phiPS,FractionBorder));
316EndDefine;
317
318
319define AllGBases(Rel)
320  TopLevel CurrentRing;
321  P := CurrentRing;
322  ParameterSol := RationalSolve(Rel.rel);
323  S := RingOf(Rel.GB[1]);
324  CoeffHoms := [PolyAlgebraHom(CoeffRing(S),CoeffRing(P),
325			 [AsRAT(x)|x in sol]) | sol in ParameterSol];
326  return [ apply(PolyRingHom(S,P,phi,indets(P)), Rel.GB) | phi in CoeffHoms];
327enddefine; -- AllGBases
328
329
330define AllBBases(Rel)
331  TopLevel CurrentRing;
332  P := CurrentRing;
333  ParameterSol := RationalSolve(Rel.rel);
334  BB := [ f.head + f.tail | f in Rel.BB];
335  S := RingOf(BB[1]);
336  CoeffHoms := [PolyAlgebraHom(CoeffRing(S),CoeffRing(P),
337			 [AsRAT(x)|x in sol]) | sol in ParameterSol];
338  return [ apply(PolyRingHom(S,P,phi,indets(P)), BB) | phi in CoeffHoms];
339enddefine; -- AllBBases
340
341
342-----------------------  [ Misc Procedures ] --------------
343
344Define Delta_Var(D_Var, L)
345  L:=concat([0],L);
346  S:=1;
347  For I:=2 To len(L) Do
348    S:=S*binomial(AsINT(D_Var-L[I-1]), AsINT(L[I]-L[I-1]));
349  EndFor;
350  Return S;
351EndDefine;
352
353
354-- L is a list of terms, D are the vars max degs
355Define Delta(L,D)
356  L:=GetCols(mat([exponents(T)|T In L]) );
357  L:=[sorted(R)|R In L];
358  Ris:=1;
359  For I:=1 To len(L) Do
360    Ris := Ris*$.Delta_Var(D[I], L[I]);
361  EndFor;
362  Return Ris;
363EndDefine;
364
365--------------------------- [ Imaging Procedures ] -------------------------
366
367Define BuildParamRing(P, ParamsName, NumParams)
368//  PrintLn NumParams;
369  R := NewPolyRing(RingQQ(), SymbolRange(ParamsName, 1, NumParams));
370//  K := NewFractionField(R);
371  return NewPolyRing(R, IndetSymbols(P));
372EndDefine;
373
374EndPackage; -- $DesignOfExperiments
375
376
377-- PrintLn "/**/  FD := [0..1, 0..1, 0..1];  QB := [1, x, y, z];";
378-- PrintLn "/**/  FD := [-1..1, -1..1];  QB := [1, x, y, x^2, y^2]; // ex (h)";
379-- PrintLn "/**/  FD := [-1..1, [-1,1]];  QB := [1, x, y]; // ex (g)";
380-- PrintLn;
381-- PrintLn "/**/  Rel := GBRelations(FD, QB);";
382-- PrintLn "/**/  multiplicity(RingOf(Rel.rel[1])/ideal(Rel.rel));";
383-- PrintLn "/**/  AllGBs := AllGBases(Rel);  indent(AllGBs);";
384-- PrintLn "/**/  indent([RationalSolve(B) | B in AllGBs]);";
385-- PrintLn;
386-- PrintLn "/**/  Rel := BBRelations(FD, QB);";
387-- PrintLn "/**/  multiplicity(RingOf(Rel.rel[1])/ideal(Rel.rel));";
388-- PrintLn "/**/  AllBBs := AllBBases(Rel);  indent(AllBBs);";
389-- PrintLn "/**/  indent([RationalSolve(B) | B in AllBBs]);";
390-- PrintLn "/**/  indent(diff(subsets(FullDesignPoints(FD),3), [RationalSolve(B) | B in AllBBs]));";
391