1module transfrm;
2
3% Cobasis transformations
4
5% Author: David Hartley
6
7% Redistribution and use in source and binary forms, with or without
8% modification, are permitted provided that the following conditions are met:
9%
10%    * Redistributions of source code must retain the relevant copyright
11%      notice, this list of conditions and the following disclaimer.
12%    * Redistributions in binary form must reproduce the above copyright
13%      notice, this list of conditions and the following disclaimer in the
14%      documentation and/or other materials provided with the distribution.
15%
16% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
18% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
19% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR
20% CONTRIBUTORS
21% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27% POSSIBILITY OF SUCH DAMAGE.
28%
29
30
31COMMENT. Data structure:
32
33        xform   ::= list of 1-form kernel . 1-form pf
34
35endcomment;
36
37
38fluid '(xvars!* kord!* subfg!*);
39global '(!*sqvar!*);
40
41fluid '(cfrmcob!* cfrmcrd!* cfrmdrv!* cfrmrsx!* !*edssloppy);
42
43
44% Type coersions
45
46
47symbolic procedure !*a2xform u;
48   % u:list of equation -> !*a2xform:xform
49   % Turn off subfg!* to stop let rules being applied.
50   % should remove x=x entries
51   for each j in getrlist u collect
52      if eqexpr j then !*a2k lhs j . xpartitop rhs j where subfg!* = nil
53      else rerror(eds,000,"Incorrectly formed transform");
54
55
56symbolic procedure !*xform2map x;
57   % x:xform -> !*xform2map:map
58   foreach p in x collect
59      car p . mk!*sq !*pf2sq cdr p;
60
61
62symbolic procedure !*map2xform u;
63   % u:map -> !*map2xform:xform
64   % Turn off subfg!* to stop let rules being applied.
65   (for each x in u collect
66      car u . xpartitop cadr u) where subfg!* = nil;
67
68
69symbolic procedure !*xform2drv x;
70   % x:xform -> !*xform2drv:drv
71   % Turn off subfg!* to stop let rules being applied.
72   (foreach p in x collect
73      {'replaceby,car p,!*pf2a cdr p}) where subfg!* = nil;
74
75
76symbolic procedure !*xform2sys x;
77   % x:xform -> !*xform2sys:sys
78   foreach p in x collect addpf(!*k2pf car p,negpf cdr p);
79
80
81% Transform
82
83
84put('transform,'rtypefn,'getrtypecar);
85put('transform,'cfrmfn,'transformcfrm);
86put('transform,'edsfn,'transformeds);
87
88
89symbolic procedure transformcfrm(m,x);
90   % m:cfrm, x:list of equation -> transformcfrm:cfrm
91   % Transform m using map x.
92   checkcfrm xformcfrm(m,!*a2xform x);
93
94
95symbolic procedure transformeds(s,x);
96   % s:eds, x:list of equation -> transform:eds
97   % pulls back s using map x
98   xformeds(s,!*a2xform x);
99
100
101symbolic procedure xformcfrm(m,x);
102   % m:cfrm, x:xform -> xformcfrm:cfrm
103   % Apply transform x to m, where x may be either old=f(new,old) or
104   % new=f(old).
105   xformcfrm1(m,car u,cadr u,caddr u)
106      where u = getxform(x,cfrm_cob m);
107
108
109symbolic procedure xformcfrm1(m,x,y,new);
110   % m:cfrm, x,y:xform, new:cob -> xformcfrm1:cfrm
111   % Apply transform x to m, where x is old=f(new,old), y is x inverse,
112   % and new gives the new cobasis elements.
113   begin scalar p,z;
114   m := copycfrm m;
115   z := pair(foreach p in x collect car p,new);
116   cfrm_cob m :=            % replace old by new in-place
117      foreach k in cfrm_cob m collect % sublis here destroys kernels
118                                if p := atsoc(k,z) then cdr p else k;
119   cfrm_crd m :=        % retain all old coordinates (may appear in eds)
120       reverse union(foreach k in new join if exact k then {cadr k},
121                          reverse cfrm_crd m);
122   cfrm_drv m :=      % add new differentials and structure equations
123      append(xformdrv(cfrm_drv m,x),
124      append(!*xform2drv foreach p in x join if exact car p then {p},
125             structeqns(y,x)));
126   if !*edssloppy then m := updatersx m; % invxform may have added new
127                                         % rsx
128   m := purgecfrm m;
129   return m;
130   end;
131
132
133symbolic procedure xformcfrm0(m,x,new);
134   % m:cfrm, x:xform, new:cob -> xformcfrm0:cfrm
135   % Cut down version of xformcfrm1 which doesn't update structure
136   % equations. Useful when following operations are purely algebraic.
137   begin scalar p,z;
138   m := copycfrm m;
139   z := pair(foreach p in x collect car p,new);
140   cfrm_cob m :=            % replace old by new in-place
141      foreach k in cfrm_cob m collect % sublis here destroys kernels
142                                if p := atsoc(k,z) then cdr p else k;
143   cfrm_crd m :=        % retain all old coordinates (may appear in eds)
144       reverse union(foreach k in new join if exact k then {cadr k},
145                          reverse cfrm_crd m);
146   if !*edssloppy then m := updatersx m; % invxform may have added new
147                                         % rsx
148   m := purgecfrm m;
149   return m;
150   end;
151
152
153symbolic procedure xformdrv(d,x);
154   % d:drv, x:xform -> xformdrv:drv
155   % Apply xform to drv. Must suppress active rules, for example if d a
156   % => d b is active and x = {d b => d a}, then after applying x, it
157   % will be undone immediately.
158   pullbackdrv(d,!*xform2map x) where subfg!* = nil;
159
160
161symbolic procedure updatersx m;
162   % m:cfrm -> updatersx:cfrm
163   % Reinstall restrictions in s from global variable, typically
164   % after solvepfsys when !*edssloppy is t.
165   begin
166   m := copycfrm m;
167   cfrm_rsx m := foreach f in purge cfrmrsx!* collect !*pf2a f;
168   return m;
169   end;
170
171
172symbolic procedure xformeds(s,x);
173   % s:eds, x:xform -> xformeds:eds
174   % Apply transform x to m, where x may be either old=f(new,old) or
175   % new=f(old).
176   % possibly changes kernel order
177   xformeds1(s,car u,cadr u,caddr u)
178      where u = getxform(x,edscob s);
179
180
181
182symbolic procedure xformeds1(s,x,y,new);
183   % s:eds, x,y:xform, new:cob -> xformeds1:eds
184   % Apply transform x to m, where x is old=f(new,old), y is x inverse,
185   % and new gives the new cobasis elements.  Changes background
186   % coframing.
187   begin scalar k;
188   s := copyeds s;
189   % Transform coframing
190   eds_cfrm s := xformcfrm1(eds_cfrm s,x,y,new);
191   % Make sure old get eliminated (and add new to kord!* for safety)
192   k := updkordl append(foreach p in x collect car p,new);
193   x := !*xform2sys x;
194   % Transform rest of eds
195   eds_sys s := foreach f in eds_sys s collect xreduce(xreorder f,x);
196   eds_ind s := foreach f in eds_ind s collect xreduce(xreorder f,x);
197   remkrns s;
198   s := purgeeds!* s;
199   rempropeds(s,'jet0);
200   foreach f in {'solved,'reduced} do rempropeds(s,f);
201   setkorder k;
202   s := normaleds s; % Refine this a bit?
203   setcfrm eds_cfrm!* s;
204   return s;
205   end;
206
207
208symbolic procedure xformeds0(s,x,new);
209   % s:eds, x:xform, new:cob -> xformeds0:eds
210   % Cut down version of xformeds1 which doesn't care about structure
211   % equations (some are lost). Useful when following operations are
212   % purely algebraic. Changes background coframing.
213   begin scalar k;
214   s := copyeds s;
215   % Transform coframing (ignore structure equations)
216   eds_cfrm s := xformcfrm0(eds_cfrm s,x,new);
217   % Make sure old get eliminated (and add new to kord!* for safety)
218   k := updkordl append(foreach p in x collect car p,new);
219   x := !*xform2sys x;
220   % Transform rest of eds
221   eds_sys s := foreach f in eds_sys s collect xreduce(xreorder f,x);
222   eds_ind s := foreach f in eds_ind s collect xreduce(xreorder f,x);
223   remkrns s;
224   s := purgeeds!* s;
225   rempropeds(s,'jet0);
226   foreach f in {'solved,'reduced} do rempropeds(s,f);
227   setkorder k;
228   s := normaleds s; % Refine this a bit?
229   setcfrm eds_cfrm!* s;
230   return s;
231   end;
232
233
234symbolic procedure getxform(x,cob);
235   % x:xform, cob:cob -> getxform:{xform,xform,cob}
236   % Analyse transform x, which may be either old=f(new,old) or
237   % new=f(old). The sense is established by cob, which contains the old
238   % cobasis. Return value is {x,y,new} where x is in the sense old =
239   % f(new,old), and y is the inverse of x (ie new = f(old)).  The
240   % inverse y is calculated only if x is old = f(new,old) and there
241   % are anholonomic forms in new.
242   begin scalar old,new,y;
243   foreach p in x do
244   << new := union(xpows cdr p,new);
245      old := car p . old >>;
246   if not xnp(old,cob) then  % x is new=f(old), must invert
247   << y := x; x := invxform x;
248      new := old; old := foreach p in x collect car p >>;
249   new := sort(setdiff(new,cob),'termordp);
250   edsdebug("New cobasis elements...",new,'cob);
251   edsdebug("... replacing old cobasis elements",old,'cob);
252   if length new neq length old or not subsetp(old,cob) then
253      rerror(eds,000,"Bad cobasis transformation");
254   if not allexact new and null y then
255      y := invxform x; % for structure equations
256   return {x,y,new};
257   end;
258
259
260% Structure equations
261
262
263symbolic procedure xformdrveval u;
264   % u:{rlist,rlist} or {rlist} -> xformdrveval:rlist
265   begin scalar x,y,k;
266   y := !*a2xform car u;
267   x := if cdr u then !*a2xform cadr u else invxform y;
268   k := updkordl foreach p in x collect car p;
269   y := structeqns(y,x);
270   setkorder k;
271   return makelist y;
272   end;
273
274symbolic procedure xformdrveval u;
275   % u:{rlist,rlist} or {rlist} -> xformdrveval:rlist
276   begin scalar x,y,xvars!*;
277   y := !*a2xform car u;
278   x := if cdr u then !*a2xform cadr u else invxform y;
279   y := structeqns(y,x);
280   return makelist y;
281   end;
282
283
284symbolic procedure structeqns(y,x);
285   % y,x:xform -> structeqns:list of rule
286   % y is the inverse of x, and d lhs x are known.
287   % Returns rules for d lhs y.
288   begin scalar ok;
289   ok := updkordl foreach p in x collect car p;
290   x := !*xform2sys x;
291   y := foreach p in y join
292                 if not exact car p then
293              {{'replaceby,
294                {'d,car p},
295                    !*pf2a xreduce(exdfpf cdr p,x)}};
296   setkorder ok;
297   return y;
298   end;
299
300symbolic procedure structeqns(y,x);
301   % y,x:xform -> structeqns:list of rule
302   % y is the inverse of x, and d lhs x are known.
303   % Returns rules for d lhs y.
304   begin scalar ok;
305   ok := updkordl sort(foreach p in x collect car p,function ordop);
306   x := !*xform2sys x;
307   y := foreach p in y join
308                 if not exact car p then
309              {{'replaceby,
310                {'d,car p},
311                    !*pf2a xreduce(exdfpf cdr p,x)}};
312   setkorder ok;
313   return y;
314   end;
315
316% Inverting tranformations
317
318
319put('invert,'rtypefn,'quotelist);
320put('invert,'listfn,'inverteval);
321
322symbolic procedure inverteval(u,v);
323   % u:{prefix list of eqn} -> inverteval:{prefix list of eqn}
324   % u is unevaluated.
325   makelist foreach p in invxform !*a2xform(u := reval car u) collect
326      {'equal,car p,!*pf2a1(cdr p,v)};
327
328
329symbolic procedure invxform x;
330   % x:xform -> invxform:xform
331   % Returns inverse transformation. Selects kernels to eliminate based
332   % on prevailing order
333   begin scalar old,y,k, subfg!*;
334   subfg!* := nil;
335   foreach p in x do old := union(xpows cdr p,old);
336   old := sort(old,'termordp);
337   % ensure old eliminated, and add new to kord!* for safety
338   k := updkordl append(old,foreach p in x collect car p);
339   edsdebug("Inverting transform",nil,nil);
340   y := solvepfsys1(!*xform2sys x,old);       % invert transformation
341   if cadr y then
342      rerror(eds,000,"Cobasis transform could not be inverted");
343   setkorder k;
344   return foreach f in car y collect
345      lpow f . negpf xreorder red f;
346   end;
347
348
349symbolic procedure tmpind s;
350   % s:eds -> tmpind:{eds,xform}
351   % Returns s with eds_ind s all kernels, transforming to a new
352   % cobasis if necessary.  Second return value is nil if no change
353   % made, or the list of transformation relations.  Structure
354   % equations are not transformed, so s should be closed first if
355   % necessary.  NB.  Background coframing changed.
356   begin scalar new,x;
357   if singleterms eds_ind s then return {s,nil};
358   new := foreach f in eds_ind s collect mkform!*(intern gensym(),1);
359   x := invxform pair(new,eds_ind s);
360   updkordl foreach p in x collect car p;
361   return {xformeds0(s,x,new),x};
362   end;
363
364endmodule;
365
366end;
367