Edit detail for SandBoxLinearOperator revision 1 of 15

changed:
-
\begin{spad}
)abbrev domain LAZY LazyLinearOperator
)lib CARTEN
LazyLinearOperator(dim:NNI,gen:OrderedFinite,K:CommutativeRing): Exports == Implementation where
  NNI ==> NonNegativeInteger
  NAT ==> PositiveInteger
  T ==> CartesianTensor(1,dim,K)

  Exports ==> Join(Ring, BiModule(K,K), Monoidal NNI, RetractableTo K) with
    inp: List K -> %
      ++ incoming vector
    inp: List % -> %
    out: List K -> %
      ++ output vector
    out: List % -> %
    arity: % -> Prop %
    basisVectors: () -> List %
    basisForms: () -> List %
    tensor: % -> T
    map: (K->K,%) -> %
    if K has Evalable(K) then Evalable(K)
    ravel: % -> List K
    unravel: (Prop %,List K) -> %
    coerce:(x:List NAT) -> %
      ++ identity for composition and permutations of its products
    coerce:(x:List None) -> %
      ++ [] = 1
    elt: (%,%) -> %
    elt: (%,NAT) -> %
    elt: (%,NAT,NAT) -> %
    elt: (%,NAT,NAT,NAT) -> %
    _/: (Tuple %,Tuple %) -> %
    _/: (Tuple %,%) -> %
    _/: (%,Tuple %) -> %
      ++ yet another syntax for product
    ev: NAT -> %
      ++ (2,0)-tensor for evaluation
    co: NAT -> %
      ++ (0,2)-tensor for co-evaluation

  Implementation ==> add
    import List NNI
    import NAT
    L == Record(domain:NNI, codomain:NNI, data:T)
    -- FreeMonoid provides unevaluated products
    Rep == FreeMonoid L
    rep(x:%):Rep == x pretend Rep
    per(x:Rep):% == x pretend %

    -- Prop (arity)
    dom(f:%):NNI == reduce(_+,map((y:L):NNI +-> y.gen.domain*y.exp,factors(rep f)))
    cod(f:%):NNI == reduce(_+,map((y:L):NNI +-> y.gen.codomain*y.exp,factors(rep f)))
    dat(f:%):T == reduce(product,map((y:L):T +-> y.gen.data^y.exp,factors(rep f)))
    arity(f:%):Prop % == f::Prop %

    retractIfCan(f:%):Union(K,"failed") ==
      dom(f)=0 and cod(f)=0 => retract(dat f)$T
      return "failed"
    retract(f:%):K ==
      dom(f)=0 and cod(f)=0 => retract(dat f)$T
      error "failed"

    -- basis
    basisVectors():List % == [per [0,1,entries(row(1,i)$SquareMatrix(dim,K))::T] for i in 1..dim]
    basisForms():List % == [per [1,0,entries(row(1,i)$SquareMatrix(dim,K))::T] for i in 1..dim]
    ev(n:NAT):% ==
      dx:= basisForms()
      reduce(_+,[ (dx.i)^n * (dx.i)^n for i in 1..dim])
    co(n:NAT):% == 
      Dx:= basisVectors()
      reduce(_+,[ (Dx.i)^n * (Dx.i)^n for i in 1..dim])

    -- manipulation
    map(f:K->K, g:%):% == per [dom g,cod g,unravel(map(f,ravel dat g))$T]
    if K has Evalable(K) then
      eval(g:%,f:List Equation K):% == map((x:K):K+->eval(x,f),g)
    ravel(g:%):List K == ravel dat g
    unravel(p:Prop %,r:List K):% ==
      dim^(dom(p)+cod(p)) ~= #r => error "failed"
      per [dom(p),cod(p),unravel(r)$T]
    tensor(x:%):T == dat(x)

    -- sum
    (f:% + g:%):% ==
      dat(f)=0 => g
      dat(g)=0 => f
      dom(f) ~= dom(g) or cod(f) ~= cod(g) => error "arity"
      per [dom f,cod f,dat(f)+dat(g)]

    (f:% - g:%):% ==
      dat(f)=0 => g
      dat(g)=0 => f
      dom(f) ~= dom(f) or cod(g) ~= cod(g) => error "arity"
      per [dom f, cod f,dat(f)-dat(g)]

    _-(f:%):% == per [dom f, cod f,-dat(f)]

    -- repeated sum
    (p:NNI * f:%):% ==
      p=1 => f
      q:=subtractIfCan(p,1)
      q case NNI => q*f + f
      -- zero map (non-trivial)
      per [dom f,cod f,0*dat(f)]

    -- identity for sum (trivial zero map)
    0 == per [0,0,0]
    zero?(f:%):Boolean == dat f = 0 * dat f
    -- identity for product
    1:% == per [0,0,1]
    one?(f:%):Boolean == dat f = 1$T
    -- identity for composition
    I == per([1,1,kroneckerDelta()$T])
    -- inherited from Ring
    --(x:% = y:%):Boolean ==
    --  dom(x) ~= dom(y) or cod(x) ~= cod(y) => error "arity"
    --  dat(x) = dat(y)
    (x:% = y:%):Boolean == zero? (x - y)

    -- permutations and identities
    coerce(p:List NAT):% ==
      r:=I^#p
      #p = 1 and p.1 = 1 => return r
      p1:List Integer:=[i for i in 1..#p]
      p2:List Integer:=[#p+i for i in p]
      p3:=concat(p1,p2)
      per [#p,#p,reindex(dat r,p3)]
    coerce(p:List None):% == per [0,0,1]
    coerce(x:K):% == 1*x

    -- product
    elt(f:%,g:%):% == f * g
    elt(f:%,g:NAT):% == f * I^g
    elt(f:%,g1:NAT,g2:NAT):% == f * [g1 @ NAT,g2 @ NAT]::List NAT::%
    elt(f:%,g1:NAT,g2:NAT,g3:NAT):% == f * [g1 @ NAT,g2 @ NAT,g3 @ NAT]::List NAT::%
    apply(f:%,g:%):% == f * g
    (f:% * g:%):% ==
      r:T := product(dat f,dat g)
      -- dom(f) + cod(f) + dom(g) + cod(g)
      p:List Integer := concat                _
        [[i for i in 1..dom(f)],              _
        [dom(f)+cod(f)+i for i in 1..dom(g)], _
        [dom(f)+i for i in 1..cod(f)],        _
        [dom(f)+dom(g)+cod(f)+i for i in 1..cod(g)]]
        -- dom(f) + dom(g) + cod(f) + cod(g)
      per [dom(f)+dom(g),cod(f)+cod(g),reindex(r,p)]

    -- repeated product
    (f:% ^ p:NNI):% ==
      p=1 => f
      q:=subtractIfCan(p,1)
      q case NNI => f^q * f
      1

    -- composition:
    -- f/g : A^n -> A^p = f:A^n -> A^m / g:A^m -> A^p
    (ff:% / gg:%):% ==
      g:=gg; f:=ff
      -- partial application from the left
      n:=subtractIfCan(cod ff,dom gg)
      if n case NNI and n>0 then
      -- apply g on f from the left, pass extra f outputs on the right
        print(hconcat([message("arity warning: "), _
              over(arity(ff)::OutputForm, _
                   arity(gg)::OutputForm*(arity(I)::OutputForm)^n::OutputForm) ]))$OutputForm
        g:=gg*I^n
      m:=subtractIfCan(dom gg, cod ff)
      --  apply g on f from the left, add extra g inputs on the left
      if m case NNI and m>0 then
        print(hconcat([message("arity warning: "), _
              over((arity(I)::OutputForm)^m::OutputForm*arity(ff)::OutputForm, _
                   arity(gg)::OutputForm)]))$OutputForm
        f:=I^m*ff
      f1:Integer:=dom(f)+1
      r:T := contract(cod(f),dat f,f1, dat g,1)
      per [dom(f),cod(g),r]

    -- another notation for composition of products
    (t:Tuple % / x:%):% == t / construct([x])$PrimitiveArray(%)::Tuple(%)
    (x:% / t:Tuple %):% == construct([x])$PrimitiveArray(%)::Tuple(%) / t
    (f:Tuple % / g:Tuple %):% ==
      -- optimize leading and trailing identities ?
      f1:=0; f2:=length(f)-1
      --for i in 0..length(f)-1 repeat
      --  if select(f,i)~=I then
      --    if f1=0 then f1:=i
      --    f2:=i
      fs:List % := [select(f,i) for i in f1..f2]
      g1:=0; g2:=length(g)-1
      --for i in 0..length(g)-1 repeat
      --  if select(g,i)~=I then
      --    if g1=0 then g1:=i
      --    g2:=i
      gs:List % := [select(g,i) for i in g1..g2]
      fr:=reduce(elt@(%,%)->%,fs,1)
      gr:=reduce(elt@(%,%)->%,gs,1)
      fr / gr
      --f2:=length(f)-1-f2
      --g2:=length(g)-1-g2
      --if f1+cod(fr)+f2 ~= g1+dom(gr)+g2 then error "arity"
      --if cod(fr) < dom(gr) then  -- more inputs
      --  r:T := contract(cod(fr),dat fr,dom(fr)+1, dat gr,1+f1)
      --  -- move f1 inputs of gr before inputs of fr and f2 inputs of gr after inputs of fr
      --  -- r:=reindex(r,[])
      --  return per [f1+dom(fr)+f2,cod(gr),r]
      --else  -- more outputs?
      --  r:T := contract(dom(gr),dat fr,dom(fr)+1+g1, dat gr,1)
      --  -- move g1 outputs of fr before outputs of gr and g2 outputs of fr after outputs of gr
      --  -- r:=reindex(r,[])
      --  return per [dom(fr),g1+cod(gr)+g2,r]

    (x:K * y:%):% == per [dom y, cod y,x*dat(y)]
    (x:% * y:K):% == per [dom x,cod x,dat(x)*y]
    (x:Integer * y:%):% == per [dom y,cod y,x*dat(y)]

    -- constructors
    inp(x:List K):% == per [1,0,entries(x)::T]
    inp(x:List %):% ==
      #removeDuplicates([dom(y) for y in x]) ~= 1 or
        #removeDuplicates([cod(y) for y in x]) ~= 1 => error "arity"
      per [(dom(first x)+1),cod(first x),[dat(y) for y in x]::T]$Rep
    out(x:List K):% == per [0,1,entries(x)::T]
    out(x:List %):% ==
      #removeDuplicates([dom(y) for y in x])~=1 or
        #removeDuplicates([cod(y) for y in x])~=1 => error "arity"
      per [dom(first x),(cod(first x)+1),[dat(y) for y in x]::T]$Rep

    -- display operators using basis
    coerce(x:%):OutputForm ==
      dom(x)=0 and cod(x)=0 => return dat(x)::OutputForm
      if size()$gen > 0 then
        gens:List OutputForm:=[index(i::PositiveInteger)$gen::OutputForm for i in 1..dim]
      else
        -- default to numeric indices
        gens:List OutputForm:=[i::OutputForm for i in 1..dim]
      -- input basis
      inps:List OutputForm := []
      for i in 1..dom(x) repeat
        empty? inps => inps:=gens
        inps:=concat [[(inps.k * gens.j) for j in 1..dim] for k in 1..#inps]
      -- output basis
      outs:List OutputForm := []
      for i in 1..cod(x) repeat
        empty? outs => outs:=gens
        outs:=concat [[(outs.k * gens.j) for j in 1..dim] for k in 1..#outs]
      -- combine input (superscripts) and/or output(subscripts) to form basis symbols
      bases:List OutputForm
      if #inps > 0 and #outs > 0 then
        bases:=concat([[ scripts(message("|"),[i,j]) for i in outs] for j in inps])
      else if #inps > 0 then
        bases:=[super(message("|"),i) for i in inps]
      else if #outs > 0 then
        bases:=[sub(message("|"),j) for j in outs]
      else
        bases:= []
      -- merge bases with data to form term list
      terms:=[(k=1 => base;k::OutputForm*base)
        for base in bases for k in ravel dat(x) | k~=0]
      empty? terms => return 0::OutputForm
      -- combine the terms
      return reduce(_+,terms)
\end{spad}

)abbrev domain LAZY LazyLinearOperator
)lib CARTEN
LazyLinearOperator(dim:NNI,gen:OrderedFinite,K:CommutativeRing): Exports == Implementation where
  NNI ==> NonNegativeInteger
  NAT ==> PositiveInteger
  T ==> CartesianTensor(1,dim,K)

  Exports ==> Join(Ring, BiModule(K,K), Monoidal NNI, RetractableTo K) with
    inp: List K -> %
      ++ incoming vector
    inp: List % -> %
    out: List K -> %
      ++ output vector
    out: List % -> %
    arity: % -> Prop %
    basisVectors: () -> List %
    basisForms: () -> List %
    tensor: % -> T
    map: (K->K,%) -> %
    if K has Evalable(K) then Evalable(K)
    ravel: % -> List K
    unravel: (Prop %,List K) -> %
    coerce:(x:List NAT) -> %
      ++ identity for composition and permutations of its products
    coerce:(x:List None) -> %
      ++ [] = 1
    elt: (%,%) -> %
    elt: (%,NAT) -> %
    elt: (%,NAT,NAT) -> %
    elt: (%,NAT,NAT,NAT) -> %
    _/: (Tuple %,Tuple %) -> %
    _/: (Tuple %,%) -> %
    _/: (%,Tuple %) -> %
      ++ yet another syntax for product
    ev: NAT -> %
      ++ (2,0)-tensor for evaluation
    co: NAT -> %
      ++ (0,2)-tensor for co-evaluation

  Implementation ==> add
    import List NNI
    import NAT
    L == Record(domain:NNI, codomain:NNI, data:T)
    -- FreeMonoid provides unevaluated products
    Rep == FreeMonoid L
    rep(x:%):Rep == x pretend Rep
    per(x:Rep):% == x pretend %

    -- Prop (arity)
    dom(f:%):NNI == reduce(_+,map((y:L):NNI +-> y.gen.domain*y.exp,factors(rep f)))
    cod(f:%):NNI == reduce(_+,map((y:L):NNI +-> y.gen.codomain*y.exp,factors(rep f)))
    dat(f:%):T == reduce(product,map((y:L):T +-> y.gen.data^y.exp,factors(rep f)))
    arity(f:%):Prop % == f::Prop %

    retractIfCan(f:%):Union(K,"failed") ==
      dom(f)=0 and cod(f)=0 => retract(dat f)$T
      return "failed"
    retract(f:%):K ==
      dom(f)=0 and cod(f)=0 => retract(dat f)$T
      error "failed"

    -- basis
    basisVectors():List % == [per [0,1,entries(row(1,i)$SquareMatrix(dim,K))::T] for i in 1..dim]
    basisForms():List % == [per [1,0,entries(row(1,i)$SquareMatrix(dim,K))::T] for i in 1..dim]
    ev(n:NAT):% ==
      dx:= basisForms()
      reduce(_+,[ (dx.i)^n * (dx.i)^n for i in 1..dim])
    co(n:NAT):% == 
      Dx:= basisVectors()
      reduce(_+,[ (Dx.i)^n * (Dx.i)^n for i in 1..dim])

    -- manipulation
    map(f:K->K, g:%):% == per [dom g,cod g,unravel(map(f,ravel dat g))$T]
    if K has Evalable(K) then
      eval(g:%,f:List Equation K):% == map((x:K):K+->eval(x,f),g)
    ravel(g:%):List K == ravel dat g
    unravel(p:Prop %,r:List K):% ==
      dim^(dom(p)+cod(p)) ~= #r => error "failed"
      per [dom(p),cod(p),unravel(r)$T]
    tensor(x:%):T == dat(x)

    -- sum
    (f:% + g:%):% ==
      dat(f)=0 => g
      dat(g)=0 => f
      dom(f) ~= dom(g) or cod(f) ~= cod(g) => error "arity"
      per [dom f,cod f,dat(f)+dat(g)]

    (f:% - g:%):% ==
      dat(f)=0 => g
      dat(g)=0 => f
      dom(f) ~= dom(f) or cod(g) ~= cod(g) => error "arity"
      per [dom f, cod f,dat(f)-dat(g)]

    _-(f:%):% == per [dom f, cod f,-dat(f)]

    -- repeated sum
    (p:NNI * f:%):% ==
      p=1 => f
      q:=subtractIfCan(p,1)
      q case NNI => q*f + f
      -- zero map (non-trivial)
      per [dom f,cod f,0*dat(f)]

    -- identity for sum (trivial zero map)
    0 == per [0,0,0]
    zero?(f:%):Boolean == dat f = 0 * dat f
    -- identity for product
    1:% == per [0,0,1]
    one?(f:%):Boolean == dat f = 1$T
    -- identity for composition
    I == per([1,1,kroneckerDelta()$T])
    -- inherited from Ring
    --(x:% = y:%):Boolean ==
    --  dom(x) ~= dom(y) or cod(x) ~= cod(y) => error "arity"
    --  dat(x) = dat(y)
    (x:% = y:%):Boolean == zero? (x - y)

    -- permutations and identities
    coerce(p:List NAT):% ==
      r:=I^#p
      #p = 1 and p.1 = 1 => return r
      p1:List Integer:=[i for i in 1..#p]
      p2:List Integer:=[#p+i for i in p]
      p3:=concat(p1,p2)
      per [#p,#p,reindex(dat r,p3)]
    coerce(p:List None):% == per [0,0,1]
    coerce(x:K):% == 1*x

    -- product
    elt(f:%,g:%):% == f * g
    elt(f:%,g:NAT):% == f * I^g
    elt(f:%,g1:NAT,g2:NAT):% == f * [g1 @ NAT,g2 @ NAT]::List NAT::%
    elt(f:%,g1:NAT,g2:NAT,g3:NAT):% == f * [g1 @ NAT,g2 @ NAT,g3 @ NAT]::List NAT::%
    apply(f:%,g:%):% == f * g
    (f:% * g:%):% ==
      r:T := product(dat f,dat g)
      -- dom(f) + cod(f) + dom(g) + cod(g)
      p:List Integer := concat                _
        [[i for i in 1..dom(f)],              _
        [dom(f)+cod(f)+i for i in 1..dom(g)], _
        [dom(f)+i for i in 1..cod(f)],        _
        [dom(f)+dom(g)+cod(f)+i for i in 1..cod(g)]]
        -- dom(f) + dom(g) + cod(f) + cod(g)
      per [dom(f)+dom(g),cod(f)+cod(g),reindex(r,p)]

    -- repeated product
    (f:% ^ p:NNI):% ==
      p=1 => f
      q:=subtractIfCan(p,1)
      q case NNI => f^q * f
      1

    -- composition:
    -- f/g : A^n -> A^p = f:A^n -> A^m / g:A^m -> A^p
    (ff:% / gg:%):% ==
      g:=gg; f:=ff
      -- partial application from the left
      n:=subtractIfCan(cod ff,dom gg)
      if n case NNI and n>0 then
      -- apply g on f from the left, pass extra f outputs on the right
        print(hconcat([message("arity warning: "), _
              over(arity(ff)::OutputForm, _
                   arity(gg)::OutputForm*(arity(I)::OutputForm)^n::OutputForm) ]))$OutputForm
        g:=gg*I^n
      m:=subtractIfCan(dom gg, cod ff)
      --  apply g on f from the left, add extra g inputs on the left
      if m case NNI and m>0 then
        print(hconcat([message("arity warning: "), _
              over((arity(I)::OutputForm)^m::OutputForm*arity(ff)::OutputForm, _
                   arity(gg)::OutputForm)]))$OutputForm
        f:=I^m*ff
      f1:Integer:=dom(f)+1
      r:T := contract(cod(f),dat f,f1, dat g,1)
      per [dom(f),cod(g),r]

    -- another notation for composition of products
    (t:Tuple % / x:%):% == t / construct([x])$PrimitiveArray(%)::Tuple(%)
    (x:% / t:Tuple %):% == construct([x])$PrimitiveArray(%)::Tuple(%) / t
    (f:Tuple % / g:Tuple %):% ==
      -- optimize leading and trailing identities ?
      f1:=0; f2:=length(f)-1
      --for i in 0..length(f)-1 repeat
      --  if select(f,i)~=I then
      --    if f1=0 then f1:=i
      --    f2:=i
      fs:List % := [select(f,i) for i in f1..f2]
      g1:=0; g2:=length(g)-1
      --for i in 0..length(g)-1 repeat
      --  if select(g,i)~=I then
      --    if g1=0 then g1:=i
      --    g2:=i
      gs:List % := [select(g,i) for i in g1..g2]
      fr:=reduce(elt@(%,%)->%,fs,1)
      gr:=reduce(elt@(%,%)->%,gs,1)
      fr / gr
      --f2:=length(f)-1-f2
      --g2:=length(g)-1-g2
      --if f1+cod(fr)+f2 ~= g1+dom(gr)+g2 then error "arity"
      --if cod(fr) < dom(gr) then  -- more inputs
      --  r:T := contract(cod(fr),dat fr,dom(fr)+1, dat gr,1+f1)
      --  -- move f1 inputs of gr before inputs of fr and f2 inputs of gr after inputs of fr
      --  -- r:=reindex(r,[])
      --  return per [f1+dom(fr)+f2,cod(gr),r]
      --else  -- more outputs?
      --  r:T := contract(dom(gr),dat fr,dom(fr)+1+g1, dat gr,1)
      --  -- move g1 outputs of fr before outputs of gr and g2 outputs of fr after outputs of gr
      --  -- r:=reindex(r,[])
      --  return per [dom(fr),g1+cod(gr)+g2,r]

    (x:K * y:%):% == per [dom y, cod y,x*dat(y)]
    (x:% * y:K):% == per [dom x,cod x,dat(x)*y]
    (x:Integer * y:%):% == per [dom y,cod y,x*dat(y)]

    -- constructors
    inp(x:List K):% == per [1,0,entries(x)::T]
    inp(x:List %):% ==
      #removeDuplicates([dom(y) for y in x]) ~= 1 or
        #removeDuplicates([cod(y) for y in x]) ~= 1 => error "arity"
      per [(dom(first x)+1),cod(first x),[dat(y) for y in x]::T]$Rep
    out(x:List K):% == per [0,1,entries(x)::T]
    out(x:List %):% ==
      #removeDuplicates([dom(y) for y in x])~=1 or
        #removeDuplicates([cod(y) for y in x])~=1 => error "arity"
      per [dom(first x),(cod(first x)+1),[dat(y) for y in x]::T]$Rep

    -- display operators using basis
    coerce(x:%):OutputForm ==
      dom(x)=0 and cod(x)=0 => return dat(x)::OutputForm
      if size()$gen > 0 then
        gens:List OutputForm:=[index(i::PositiveInteger)$gen::OutputForm for i in 1..dim]
      else
        -- default to numeric indices
        gens:List OutputForm:=[i::OutputForm for i in 1..dim]
      -- input basis
      inps:List OutputForm := []
      for i in 1..dom(x) repeat
        empty? inps => inps:=gens
        inps:=concat [[(inps.k * gens.j) for j in 1..dim] for k in 1..#inps]
      -- output basis
      outs:List OutputForm := []
      for i in 1..cod(x) repeat
        empty? outs => outs:=gens
        outs:=concat [[(outs.k * gens.j) for j in 1..dim] for k in 1..#outs]
      -- combine input (superscripts) and/or output(subscripts) to form basis symbols
      bases:List OutputForm
      if #inps > 0 and #outs > 0 then
        bases:=concat([[ scripts(message("|"),[i,j]) for i in outs] for j in inps])
      else if #inps > 0 then
        bases:=[super(message("|"),i) for i in inps]
      else if #outs > 0 then
        bases:=[sub(message("|"),j) for j in outs]
      else
        bases:= []
      -- merge bases with data to form term list
      terms:=[(k=1 => base;k::OutputForm*base)
        for base in bases for k in ravel dat(x) | k~=0]
      empty? terms => return 0::OutputForm
      -- combine the terms
      return reduce(_+,terms)

   Compiling FriCAS source code from file 
      /var/zope2/var/LatexWiki/383248084152889850-25px001.spad using 
      old system compiler.
   LAZY abbreviates domain LazyLinearOperator 
   CartesianTensor is now explicitly exposed in frame initial 
   CartesianTensor will be automatically loaded when needed from 
      /var/zope2/var/LatexWiki/CARTEN.NRLIB/CARTEN
------------------------------------------------------------------------
   initializing NRLIB LAZY for LazyLinearOperator 
   compiling into NRLIB LAZY 
****** comp fails at level 1 with expression: ******
((|Monoidal| (|NonNegativeInteger|)))
****** level 1  ******
$x:= (Monoidal (NonNegativeInteger))
$m:= $EmptyMode
$f:=
((((K # #) (|gen| # #) (|dim| # #) (|LazyLinearOperator| #) ...)))

   >> Apparent user error:
   cannot compile (Monoidal (NonNegativeInteger))