Edit detail for SandBoxSurfaceComplex revision 1 of 1

changed:
-
\begin{axiom}
)abbrev domain CMAP CellMap
++
CellMap(R,n) : Exports == Implementation where

  R: Join(Ring,Comparable)
  n: NonNegativeInteger
  
  X   ==> Expression R
  DP  ==> DirectProduct
  OF  ==> OutputForm
  NNI ==> NonNegativeInteger
  MAP ==> List X -> List X
  DOM ==> List(Segment X)

  Exports == Join(CoercibleTo OF,SetCategory,Evalable X) with
  
    _= : (%,%) -> Boolean
      ++ f1=f2 checks if two given cell maps are equal, that is if they have
      ++ the same domain D and the same mapping from D into X^n.
    cellMap : (DOM,MAP) -> %
      ++ cellMap(D,f) is the constructor. Usually one has to specify the
      ++ dimension of the target space. For example, let Q=[a..b,c..d], then
      ++ cellMap(Q,Z+->[sin(Z.1),cos(Z.2),Z.1*Z.2])$CMAP(INT,3) defines a
      ++ 2-surface in X^3.
    getDom : % -> DOM
      ++ getDom(f) extracts the domain of f.
    getMap : % -> MAP
      ++ getMap(f) extracts the map of f.
    faces  : % -> List List(%)
      ++ faces(f) returns the faces of f, that means the images of the boundary
      ++ of the domain. Note: the returned list contains pairs of faces 
      ++ corresponding to the endpoints of intervals.
    coords : (Symbol,PositiveInteger) -> List X
      ++ coords(s,m) provides a sample of coordinates s[1],..,s[m] as a list.
    coordSymbols : (Symbol,PositiveInteger) -> List Symbol
      ++ coordSymbols(s,m) provides a sample of coordinates s[1],..,s[m] as a
      ++ list of symbols.
    jacobianMatrix : % -> (List X -> Matrix X)
      ++ jacobianMatrix(f) returns the Jacobian matrix as a marix valued
      ++ function defined on the same cell as the cellMap.
    tangentSpace : % -> (List(X) -> List(Vector X))
      ++ tangentSpace(f) returns a 
    coerce : % -> OutputForm
      ++ coerce(f) gives the output representation.



  Implementation ==  add 
    
    Rep := Record(d:DOM,f:MAP)


    (x:% = y:%):Boolean == 
      l:NNI:=min(#(x.d),#(y.d))
      v:List X
      for j in 1..l repeat
        s:X:=subscript('z,[j::OF])::X
        v:=concat(v,s::X)
      x.d =y.d and (x.f) v = (y.f) v => true
      false


    cellMap(dd:DOM,ff:MAP):% == 
      #dd > n => error concat("#DOM > ",string n)
      v:List X:=[1::X for j in 1..#dd]
      ~test(#ff(v)=n) => error concat("#Range ~= ", string n) 
      construct(dd,ff)


    faceLoHi(x:%,i:NNI,lo:Boolean):% ==
      l:NNI:=#(x.d)
      v:List X
      for j in 1..l repeat
        if j=i then
          if lo then
            s:X:=lo(x.d.i)
          else
            s:X:=hi(x.d.i)
        else
          if j>i then
            s:X:=subscript('%,[(j-1)::OF])::X
          else
            s:X:=subscript('%,[j::OF])::X
        v:=concat(v,s::X)
      vv:=delete(v,i..i)
      dd:List(Segment X):=delete(x.d,i..i)
      ff:MAP:=vv+->(x.f) v
      cellMap(dd,ff)
      
     
    faces(x:%):List List(%) ==
      l:NNI:=#(x.d)
      [[faceLoHi(x,j,true), faceLoHi(x,j,false)] for j in 1..l]
        
        
    getDom(x) == x.d
    getMap(x) == x.f


    coordSymbols(s:Symbol,m:PositiveInteger):List Symbol == 
      [subscript(s,[j::OF]) for j in 1..m]
      
      
    coords(s:Symbol,m:PositiveInteger):List X ==
      xs:=[subscript(s,[j::OF]) for j in 1..m]
      [coerce(xs.j)$X for j in 1..#xs]
    

    jacobianMatrix(S:%):List(X) -> Matrix(X) ==
      --xs:List Symbol:=v:=[subscript('x,[j::OF]) for j in 1..#(getDom S)]
      --x:List X:=[coerce(xs.j)$X for j in 1..#xs]
      xs:List Symbol:=coordSymbols('x,#(getDom S)::PositiveInteger)
      x:List X:=coords('x,#xs::PositiveInteger)
      F:List X:=(getMap S) x
      J:Matrix(X):=matrix [[D(ff,u) for u in xs] for ff in F]
      if Matrix(X) has Join(SetCategory,Evalable(X)) then
        (y:List X):Matrix(X)+-> eval(J,x,y)
      else
        (y:List X):Matrix(X)+-> J


    tangentSpace(S:%):List(X) -> List(Vector X) ==
      J:=jacobianMatrix(S)
      x:List X:=coords('x,#(getDom S)::PositiveInteger)
      if Vector(X) has Join(SetCategory,Evalable(X)) then
        if X has EuclideanDomain then 
          cs:List(Vector X):=columnSpace(J x)
          (y:List X):List Vector(X)+-> [eval(t,x,y) for t in cs]
      

    coerce(x) == 
      v:List X
      for j in 1..#(x.d) repeat
        s:X:=subscript('%,[j::OF])::X
        v:=concat(v,s::X)
      r:List X:=(x.f) v
      hconcat ["|",x.d::OF," -> ",r::OF,"|"]
  
    
      
)abbrev domain SCMPLX SurfaceComplex
++
SurfaceComplex(R,n) : Exports == Implementation where

  NNI ==> NonNegativeInteger
  INT ==> Integer
  
  n : NNI
  R : Join(Ring,Comparable)
  
  CMAP ==> CellMap(R,n)
  CTOF ==> CoercibleTo OutputForm
  X    ==> Expression R
  OF   ==> OutputForm
  MAP  ==> List X -> List X
  DOM  ==> List(Segment X)

  Exports == Join(AbelianGroup ,CTOF, RetractableTo CMAP) with
 
    bdry : % -> %
      ++ bdry(S) computes the boundary of the surface complex S.
    size : % -> NNI
      ++ size(S) returns the number of "pieces" of the surface complex S.
    nthCoef : (%,Integer) -> Integer
      ++ nthCoef(x, n) returns the coefficient of the n^th term of x.
    nthFactor : (%,Integer) -> CMAP
      ++ nthFactor(x, n) returns the factor of the n^th term of x.
    zero? : % -> Boolean 
      ++ zero?(S) returns true if S is the empty surface complex.
    _= : (%,%) -> Boolean
      ++ S=S' checks if the surface complexes S and S' are equal.
    terms : % -> List(Record(gen: CMAP,exp: Integer))
      ++ terms(S) returns all terms of S as a record.
    mapGen : ((CMAP -> CMAP),%) -> %
      ++ mapGen(f, e1 a1 +...+ en an) returns
      ++ \spad{e1 f(a1) +...+ en f(an)}.
    mapCoef : ((Integer -> Integer),%) -> %
      ++ mapCoef(f, e1 a1 +...+ en an) returns
      ++ \spad{f(e1) a1 +...+ f(en) an}.
    construct : (DOM,MAP) -> %
      ++ construct(d,f) constructs a term (piece) of a k-surface, where
      ++ d is the domain (a k-cell) and f is a mapping from d to a vector
      ++ space of dimension n.
    
    
    --coerce : % -> OutputForm
    
  Implementation == FreeAbelianGroup(CMAP) add

    Rep:=FreeAbelianGroup(CMAP)
      

    bdry(c:%):% == 
      if size(c) = 1 then
        s:=nthFactor(c,1)
        l:=faces(s)
        fs:=[(a.2::Rep-a.1::Rep) for a in l]
        sgn:=(j:INT):INT+->if even? (j-1) then 1 else -1
        nthCoef(c,1)*reduce("+",[sgn(j)*fs.j::Rep for j in 1..#fs])
      else
        ct:=[(nthCoef(c,j)*nthFactor(c,j))::Rep for j in 1..size(c)]
        reduce("+",map(bdry,ct))
        
 
    construct(d:DOM,f:MAP):% == cellMap(d,f)$CMAP::%

\end{axiom}

\begin{axiom}
)clear all
R ==> EXPR INT
OF ==> OutputForm

-- Cell map
R2 ==> CellMap(INT,2)
R3 ==> CellMap(INT,3)
R4 ==> CellMap(INT,4)

Q2 ==> [0..1,0..1::R]
Q3 ==> concat(Q2,[0..1::R])

--xs:List Symbol:=coordSymbols('x,4)$R4

----------------------------------------------------------------
-- https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant
----------------------------------------------------------------
-- Example 1
F1:=cellMap(Q2,X+->[X.1^2*X.2,5*X.1+sin(X.2)])$R2
J:=jacobianMatrix(F1)
x:=coords('x,2)$R2
J x
determinant(J x)

test(J x = matrix [[2*x.1*x.2,x.1^2],[5,cos(x.2)]])
test(determinant(J x) = 2*x.1*x.2*cos(x.2)-5*x.1^2)  


-- Example 2
F2:=cellMap(Q2,X+->[X.1*cos(X.2),X.1*sin(X.2)])$R2
J:=jacobianMatrix(F2)
x:=[r::R,phi::R]
(getMap F2) x
J x
determinant(J x)

test( J x = matrix [[cos(x.2),-x.1*sin(x.2)],[sin(x.2),x.1*cos(x.2)]])
test( normalize determinant(J x) = x.1)


-- Example 3
F3:=cellMap(Q3,Z+->[Z.1*sin(Z.2)*cos(Z.3),Z.1*sin(Z.2)*sin(Z.3),Z.1*cos(Z.2)])$R3
J:=jacobianMatrix(F3)
z:=[r::R,th::R,phi::R]
(getMap F3) z
J z
determinant(J z)

M:=[[sin(z.2)*cos(z.3),z.1*cos(z.2)*cos(z.3),-z.1*sin(z.2)*sin(z.3)],_
    [sin(z.2)*sin(z.3),z.1*cos(z.2)*sin(z.3),z.1*sin(z.2)*cos(z.3)],_
    [cos(z.2),-z.1*sin(z.2),0]]

test( J z = matrix M)
test( simplify determinant(J z) = z.1^2*sin(z.2) )


-- Example 4
F4:=cellMap(Q3,X+->[X.1,5*X.3,4*X.2^2-2*X.3,X.3*sin(X.1)])$R4
J:=jacobianMatrix(F4)
x:=coords('x,4)$R4
J x
nullSpace (J x)
rank (J x)
T:=tangentSpace(F4)$R4
T x


test(J x = matrix [[1,0,0],[0,0,5],[0,8*x.2,-2],[x.3*cos(x.1),0,sin(x.1)]])
test( rank (J x) = 3)
test( J x = transpose matrix (T x))

-- Example 5
F5:=cellMap(Q3,X+->[5*X.2,4*X.1^2-2*sin(X.2*X.3),X.2*X.3])$R3
J:=jacobianMatrix(F5)
x:=coords('x,3)$R3
J x
determinant (J x)

M:=[[0,5,0],[8*x.1,-2*x.3*cos(x.2*x.3),-2*x.2*cos(x.2*x.3)],[0,x.3,x.2]]
T:=tangentSpace(F5)$R3

test(J x = matrix M)
test(determinant (J x) = -40*x.1*x.2)
test( J x = transpose matrix (T x))
\end{axiom}

)abbrev domain CMAP CellMap
++
CellMap(R,n) : Exports == Implementation where

  R: Join(Ring,Comparable)
  n: NonNegativeInteger

  X   ==> Expression R
  DP  ==> DirectProduct
  OF  ==> OutputForm
  NNI ==> NonNegativeInteger
  MAP ==> List X -> List X
  DOM ==> List(Segment X)

  Exports == Join(CoercibleTo OF,SetCategory,Evalable X) with

    _= : (%,%) -> Boolean
      ++ f1=f2 checks if two given cell maps are equal, that is if they have
      ++ the same domain D and the same mapping from D into X^n.
    cellMap : (DOM,MAP) -> %
      ++ cellMap(D,f) is the constructor. Usually one has to specify the
      ++ dimension of the target space. For example, let Q=[a..b,c..d], then
      ++ cellMap(Q,Z+->[sin(Z.1),cos(Z.2),Z.1*Z.2])$CMAP(INT,3) defines a
      ++ 2-surface in X^3.
    getDom : % -> DOM
      ++ getDom(f) extracts the domain of f.
    getMap : % -> MAP
      ++ getMap(f) extracts the map of f.
    faces  : % -> List List(%)
      ++ faces(f) returns the faces of f, that means the images of the boundary
      ++ of the domain. Note: the returned list contains pairs of faces 
      ++ corresponding to the endpoints of intervals.
    coords : (Symbol,PositiveInteger) -> List X
      ++ coords(s,m) provides a sample of coordinates s[1],..,s[m] as a list.
    coordSymbols : (Symbol,PositiveInteger) -> List Symbol
      ++ coordSymbols(s,m) provides a sample of coordinates s[1],..,s[m] as a
      ++ list of symbols.
    jacobianMatrix : % -> (List X -> Matrix X)
      ++ jacobianMatrix(f) returns the Jacobian matrix as a marix valued
      ++ function defined on the same cell as the cellMap.
    tangentSpace : % -> (List(X) -> List(Vector X))
      ++ tangentSpace(f) returns a 
    coerce : % -> OutputForm
      ++ coerce(f) gives the output representation.

  Implementation ==  add 

    Rep := Record(d:DOM,f:MAP)

    (x:% = y:%):Boolean == 
      l:NNI:=min(#(x.d),#(y.d))
      v:List X
      for j in 1..l repeat
        s:X:=subscript('z,[j::OF])::X
        v:=concat(v,s::X)
      x.d =y.d and (x.f) v = (y.f) v => true
      false

    cellMap(dd:DOM,ff:MAP):% == 
      #dd > n => error concat("#DOM > ",string n)
      v:List X:=[1::X for j in 1..#dd]
      ~test(#ff(v)=n) => error concat("#Range ~= ", string n) 
      construct(dd,ff)

    faceLoHi(x:%,i:NNI,lo:Boolean):% ==
      l:NNI:=#(x.d)
      v:List X
      for j in 1..l repeat
        if j=i then
          if lo then
            s:X:=lo(x.d.i)
          else
            s:X:=hi(x.d.i)
        else
          if j>i then
            s:X:=subscript('%,[(j-1)::OF])::X
          else
            s:X:=subscript('%,[j::OF])::X
        v:=concat(v,s::X)
      vv:=delete(v,i..i)
      dd:List(Segment X):=delete(x.d,i..i)
      ff:MAP:=vv+->(x.f) v
      cellMap(dd,ff)

    faces(x:%):List List(%) ==
      l:NNI:=#(x.d)
      [[faceLoHi(x,j,true), faceLoHi(x,j,false)] for j in 1..l]

    getDom(x) == x.d
    getMap(x) == x.f

    coordSymbols(s:Symbol,m:PositiveInteger):List Symbol == 
      [subscript(s,[j::OF]) for j in 1..m]

    coords(s:Symbol,m:PositiveInteger):List X ==
      xs:=[subscript(s,[j::OF]) for j in 1..m]
      [coerce(xs.j)$X for j in 1..#xs]

    jacobianMatrix(S:%):List(X) -> Matrix(X) ==
      --xs:List Symbol:=v:=[subscript('x,[j::OF]) for j in 1..#(getDom S)]
      --x:List X:=[coerce(xs.j)$X for j in 1..#xs]
      xs:List Symbol:=coordSymbols('x,#(getDom S)::PositiveInteger)
      x:List X:=coords('x,#xs::PositiveInteger)
      F:List X:=(getMap S) x
      J:Matrix(X):=matrix [[D(ff,u) for u in xs] for ff in F]
      if Matrix(X) has Join(SetCategory,Evalable(X)) then
        (y:List X):Matrix(X)+-> eval(J,x,y)
      else
        (y:List X):Matrix(X)+-> J

    tangentSpace(S:%):List(X) -> List(Vector X) ==
      J:=jacobianMatrix(S)
      x:List X:=coords('x,#(getDom S)::PositiveInteger)
      if Vector(X) has Join(SetCategory,Evalable(X)) then
        if X has EuclideanDomain then 
          cs:List(Vector X):=columnSpace(J x)
          (y:List X):List Vector(X)+-> [eval(t,x,y) for t in cs]

    coerce(x) == 
      v:List X
      for j in 1..#(x.d) repeat
        s:X:=subscript('%,[j::OF])::X
        v:=concat(v,s::X)
      r:List X:=(x.f) v
      hconcat ["|",x.d::OF," -> ",r::OF,"|"]

      

)abbrev domain SCMPLX SurfaceComplex
++
SurfaceComplex(R,n) : Exports == Implementation where

  NNI ==> NonNegativeInteger
  INT ==> Integer

  n : NNI
  R : Join(Ring,Comparable)

  CMAP ==> CellMap(R,n)
  CTOF ==> CoercibleTo OutputForm
  X    ==> Expression R
  OF   ==> OutputForm
  MAP  ==> List X -> List X
  DOM  ==> List(Segment X)

  Exports == Join(AbelianGroup ,CTOF, RetractableTo CMAP) with

    bdry : % -> %
      ++ bdry(S) computes the boundary of the surface complex S.
    size : % -> NNI
      ++ size(S) returns the number of "pieces" of the surface complex S.
    nthCoef : (%,Integer) -> Integer
      ++ nthCoef(x, n) returns the coefficient of the n^th term of x.
    nthFactor : (%,Integer) -> CMAP
      ++ nthFactor(x, n) returns the factor of the n^th term of x.
    zero? : % -> Boolean 
      ++ zero?(S) returns true if S is the empty surface complex.
    _= : (%,%) -> Boolean
      ++ S=S' checks if the surface complexes S and S' are equal.
    terms : % -> List(Record(gen: CMAP,exp: Integer))
      ++ terms(S) returns all terms of S as a record.
    mapGen : ((CMAP -> CMAP),%) -> %
      ++ mapGen(f, e1 a1 +...+ en an) returns
      ++ \spad{e1 f(a1) +...+ en f(an)}.
    mapCoef : ((Integer -> Integer),%) -> %
      ++ mapCoef(f, e1 a1 +...+ en an) returns
      ++ \spad{f(e1) a1 +...+ f(en) an}.
    construct : (DOM,MAP) -> %
      ++ construct(d,f) constructs a term (piece) of a k-surface, where
      ++ d is the domain (a k-cell) and f is a mapping from d to a vector
      ++ space of dimension n.

    --coerce : % -> OutputForm

  Implementation == FreeAbelianGroup(CMAP) add

    Rep:=FreeAbelianGroup(CMAP)

    bdry(c:%):% == 
      if size(c) = 1 then
        s:=nthFactor(c,1)
        l:=faces(s)
        fs:=[(a.2::Rep-a.1::Rep) for a in l]
        sgn:=(j:INT):INT+->if even? (j-1) then 1 else -1
        nthCoef(c,1)*reduce("+",[sgn(j)*fs.j::Rep for j in 1..#fs])
      else
        ct:=[(nthCoef(c,j)*nthFactor(c,j))::Rep for j in 1..size(c)]
        reduce("+",map(bdry,ct))

    construct(d:DOM,f:MAP):% == cellMap(d,f)$CMAP::%

Compiling FriCAS source code from file 
      /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/9220523246559970704-25px.001.spad
      using old system compiler.
   CMAP abbreviates domain CellMap 
------------------------------------------------------------------------
   initializing NRLIB CMAP for CellMap 
   compiling into NRLIB CMAP 
****** Domain: R already in scope
   compiling exported = : ($,$) -> Boolean
Time: 0.09 SEC.

   compiling exported cellMap : (List Segment Expression R,List Expression R -> List Expression R) -> $
Time: 0.02 SEC.

   compiling local faceLoHi : ($,NonNegativeInteger,Boolean) -> $
****** comp fails at level 9 with expression: ******
error in function faceLoHi 

(SEQ (|:=| (|:| |l| (|NonNegativeInteger|)) (|#| (|x| |d|)))
     (|:| |v| (|List| (|Expression| R)))
     (REPEAT (IN |j| (SEGMENT 1 |l|))
             (SEQ
              (IF (= |j| |i|)
                  (IF |lo|
                      (|:=| (|:| |s| (|Expression| R)) (|lo| ((|x| |d|) |i|)))
                      (|:=| (|:| |s| (|Expression| R)) (|hi| ((|x| |d|) |i|))))
                  (IF (> |j| |i|)
                      (|:=| (|:| |s| (|Expression| R))
                       (|::|
                        (|subscript| '$
                         (|construct| (|::| (- |j| 1) (|OutputForm|))))
                        (|Expression| R)))
                      (|:=| (|:| |s| (|Expression| R))
                       (|::|
                        (|subscript| '$
                         (|construct| (|::| |j| (|OutputForm|))))
                        (|Expression| R)))))
              (|exit| 1
               (|:=| |v| (|concat| |v| (|::| |s| (|Expression| R)))))))
     (|:=| |vv| (|delete| |v| (SEGMENT |i| |i|)))
     (|:=| (|:| |dd| (|List| (|Segment| (|Expression| R))))
      (|delete| (|x| |d|) (SEGMENT |i| |i|)))
     (|:=|
      (|:| |ff|
           (|Mapping| (|List| (|Expression| R)) (|List| (|Expression| R))))
      (+-> |vv| ((|x| |f|) |v|)))
     (|exit| 1 (|cellMap| |dd| |ff|)))
****** level 9  ******
$x:= lo
$m:= $
$f:=
((((|s| #) (|j| # #) (|v| #) (|l| # #) ...)))

   >> Apparent user error:
   Cannot coerce lo 
      of mode (Boolean) 
      to mode (Record (: d (List (Segment (Expression R)))) (: f (Mapping (List (Expression R)) (List (Expression R)))))

)clear all

   All user variables and function definitions have been cleared.
R ==> EXPR INT

OF ==> OutputForm

-- Cell map
R2 ==> CellMap(INT,2)

R3 ==> CellMap(INT,3)

R4 ==> CellMap(INT,4)

Q2 ==> [0..1,0..1::R]

Q3 ==> concat(Q2,[0..1::R])

--xs:List Symbol:=coordSymbols('x,4)$R4

----------------------------------------------------------------
-- https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant
----------------------------------------------------------------
-- Example 1
F1:=cellMap(Q2,X+->[X.1^2*X.2,5*X.1+sin(X.2)])$R2

   CellMap is an unknown constructor and so is unavailable. Did you 
      mean to use -> but type something different instead?