restrict.red 5.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209
  1. module restrict;
  2. % Restrict to a subset of a coframing
  3. % Author: David Hartley
  4. Comment. Data structures:
  5. rsx ::= list of prefix (usually !*sq)
  6. rmap ::= {map,rsx}
  7. Restrictions are store in rmap's, where the second part gives the
  8. restrictions to the coframing.
  9. endcomment;
  10. fluid '(xvars!* kord!* subfg!*);
  11. global '(!*sqvar!*);
  12. fluid '(cfrmcob!* cfrmcrd!* cfrmdrv!* cfrmrsx!*);
  13. % Type coersions
  14. symbolic procedure !*a2rmap u;
  15. % u:list of equation/inequality -> !*a2rmap:rmap
  16. % should remove x=x entries
  17. begin scalar map,rsx;
  18. for each j in getrlist u do
  19. if eqexpr j then
  20. << map := (!*a2k lhs j . rhs j) . map;
  21. rsx := addrsx(denr simp!* cdar map,rsx) >>
  22. else if eqcar(j := reval j,'neq) then
  23. rsx := addrsx(numr subtrsq(simp!* cadr j,simp!* caddr j),rsx)
  24. else typerr(j,"either equality or inequality");
  25. map := unrollmap map;
  26. rsx := pullbackrsx(rsx,map);
  27. return {map,rsx};
  28. end;
  29. symbolic procedure !*rmap2a u;
  30. % u:rmap -> !*rmap2a:prefix
  31. makelist append(foreach p in car u collect {'equal,car p,cdr p},
  32. foreach p in cadr u collect {'neq,p,0});
  33. symbolic procedure !*map2rmap x;
  34. % x:map -> !*map2rmap:rmap
  35. % Pick up denominators in x
  36. begin scalar rsx;
  37. for each s in x do
  38. rsx := addrsx(reorder denr simp!* cdr s,rsx);
  39. return {x,reversip rsx};
  40. end;
  41. % Restrict
  42. if not operatorp 'neq then
  43. mkop 'neq; % make it an algebraic operator so it can be reval'd
  44. put('restrict,'rtypefn,'getrtypecar);
  45. put('restrict,'cfrmfn,'restrictcfrm);
  46. put('restrict,'edsfn,'restricteds);
  47. put('restrict,'listfn,'restrictlist);
  48. put('restrict,'simpfn,'simprestrict);
  49. symbolic procedure restrictcfrm(m,x);
  50. % m:cfrm, x:list of equation/inequality -> restrictcfrm:cfrm
  51. % restricts m using rmap x
  52. restrictcfrm1(m,!*a2rmap x);
  53. symbolic procedure restricteds(s,x);
  54. % s:eds, x:list of equation/inequality -> restrict:eds
  55. % restricts s using rmap x
  56. restrict(s,!*a2rmap x);
  57. symbolic procedure restrictlist(u,v);
  58. % u:{prefix list of prefix,prefix list of equation/inequality}, v:bool
  59. % -> restrictlist:prefix list of prefix
  60. begin scalar x;
  61. x := car !*a2rmap reval cadr u;
  62. u := reval car u;
  63. return
  64. makelist foreach f in cdr u join
  65. if (f := !*pf2a1(restrictpf(xpartitop f,x),v)) neq 0 then {f};
  66. end;
  67. symbolic procedure simprestrict u;
  68. % u:{prefix,prefix list of equation/inequality} -> simprestrict:sq
  69. % just ignores inequalities
  70. !*pf2sq repartit restrictpf(f,x)
  71. where f = xpartitop car u,
  72. x = car !*a2rmap reval cadr u;
  73. symbolic procedure restrictcfrm1(m,x);
  74. % m:cfrm, x:rmap -> restrictcfrm:cfrm
  75. begin scalar kl,rl;
  76. if null car x and null cadr x then return m;
  77. m := copycfrm m;
  78. kl := union(!*map2cob car x,!*rsx2cob cadr x);
  79. % Get rsx restrictions from denominators of map part
  80. rl := purge foreach p in car x join
  81. if not cfrmconstant(p := denr simp!* cdr p) then
  82. {mk!*sq !*f2q p};
  83. % Put all rsx together and restrict
  84. rl := append(cfrm_rsx m,append(cadr x,rl));
  85. cfrm_rsx m := pullbackrsx(rl,car x);
  86. if not subsetp(kl,cfrm_cob m) or 0 member cfrm_rsx m then
  87. rerror(eds,000,
  88. "Map image not within target coframing in restrict");
  89. % Restrict derivatives
  90. cfrm_drv m := restrictdrv(cfrm_drv m,car x);
  91. return purgecfrm m;
  92. end;
  93. symbolic procedure !*map2cob x;
  94. % x:map -> !*map2cob:cob
  95. % Collect all 1-form variables in map x.
  96. begin scalar f,kl;
  97. foreach p in x do
  98. << f := simp!* cdr p;
  99. f := foreach k in union(kernels denr f,kernels numr f) join
  100. if exformp k then xpows exdfk k;
  101. f := append(xpows exdfk car p,f);
  102. kl := union(f,kl) >>;
  103. return kl;
  104. end;
  105. symbolic procedure !*rsx2cob x;
  106. % x:rsx -> !*rsx2cob:cob
  107. % Collect all 1-form variables in restrictions x.
  108. begin scalar f,kl;
  109. foreach p in x do
  110. << f := simp!* p;
  111. f := foreach k in union(kernels denr f,kernels numr f) join
  112. if exformp k then xpows exdfk k;
  113. kl := union(f,kl) >>;
  114. return kl;
  115. end;
  116. symbolic procedure restrictdrv(d,x);
  117. % d:drv, x:map -> restrictdrv:drv
  118. (foreach r in d collect
  119. {car r,cadr r,mk!*sq restrictsq(simp!* caddr r,x)})
  120. ; %%% where subfg!*=nil; %%% Why?
  121. symbolic procedure restrictsq(q,x);
  122. % q:sq, x:map -> restrictsq:sq
  123. !*pf2sq restrictpf(xpartitsq q,x);
  124. symbolic procedure restrict(s,x);
  125. % s:eds, x:rmap -> restrict:eds
  126. % restricts s using rmap x
  127. begin
  128. if null car x and null cadr x then return s;
  129. % Do coframing first (spot errors faster)
  130. s := copyeds s;
  131. eds_cfrm s := restrictcfrm1(eds_cfrm s,x);
  132. % Fix flags
  133. s := purgeeds!* s;
  134. foreach f in {'solved,'pfaffian,'quasilinear,'closed} do
  135. remfalseeds(s,f);
  136. rempropeds(s,'involutive);
  137. remkrns s;
  138. % Form new eds
  139. eds_sys s := foreach f in eds_sys s join
  140. if f := restrictpf(f,car x) then {f};
  141. eds_ind s := foreach f in eds_ind s join
  142. if f := restrictpf(f,car x) then {f}
  143. else <<edsverbose(
  144. "Restriction inconsistent with independence condition",
  145. nil,nil);>>;
  146. return normaleds s;
  147. end;
  148. symbolic procedure restrictpf(f,x);
  149. % f:pf, x:map -> restrictpf:pf
  150. % restricts f using map x
  151. % should watch out for partdf's
  152. if null f then nil
  153. else if null x then f % doesn't check let rules
  154. else (if numr c then lpow f .* c .+ restrictpf(red f,x)
  155. else restrictpf(red f,x)) where c = subsq(lc f,x);
  156. symbolic procedure pullbackrsx(rsx,x);
  157. % rsx:rsx, x:map -> pullbackrsx:rsx
  158. foreach p in rsx collect mk!*sq subf(numr simp!* p,x);
  159. endmodule;
  160. end;