axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / ESTOOLS.spad

)abbrev package ESTOOLS ExpertSystemToolsPackage
++ Author: Brian Dupee
++ Date Created: May 1994
++ Date Last Updated: July 1996
++ Description:
++ \axiom{ExpertSystemToolsPackage} contains some useful functions for use
++ by the computational agents of numerical solvers.

ExpertSystemToolsPackage() : SIG == CODE where

  LEDF  ==> List Expression DoubleFloat
  KEDF  ==> Kernel Expression DoubleFloat
  LKEDF  ==> List Kernel Expression DoubleFloat
  VEDF  ==> Vector Expression DoubleFloat
  VEF  ==> Vector Expression Float
  VMF  ==> Vector MachineFloat
  EF2  ==> ExpressionFunctions2
  EFI  ==> Expression Fraction Integer
  MDF  ==> Matrix DoubleFloat
  LDF  ==> List DoubleFloat
  PDF  ==> Polynomial DoubleFloat
  EDF  ==> Expression DoubleFloat
  EF  ==> Expression Float
  SDF  ==> Stream DoubleFloat
  DF  ==> DoubleFloat
  F  ==> Float
  MF  ==> MachineFloat
  INT  ==> Integer
  NNI  ==> NonNegativeInteger
  LS  ==> List Symbol
  ST  ==> String
  LST  ==> List String
  SS  ==> Stream String
  FI  ==> Fraction Integer
  R  ==> Ring
  OR  ==> OrderedRing
  ON  ==> Record(additions:INT,multiplications:INT,exponentiations:INT,functionCalls:INT)
  RVE   ==> Record(val:EDF,exponent:INT)
  BO  ==> BasicOperator
  OCF  ==> OrderedCompletion Float
  OCDF  ==> OrderedCompletion DoubleFloat
  SOCF  ==> Segment OrderedCompletion Float
  SOCDF  ==> Segment OrderedCompletion DoubleFloat
  Measure  ==> Record(measure:F, name:String, explanations:List String)
  Measure2  ==> Record(measure:F, name:String, explanations:List String, extra:Result)
  CTYPE  ==> Union(continuous: "Continuous at the end points",
               lowerSingular: "There is a singularity at the lower end point",
                upperSingular: "There is a singularity at the upper end point",
                 bothSingular: "There are singularities at both end points",
                  notEvaluated: "End point continuity not yet evaluated")
  RTYPE  ==> Union(finite: "The range is finite",
                lowerInfinite: "The bottom of range is infinite",
                  upperInfinite: "The top of range is infinite",
                    bothInfinite: "Both top and bottom points are infinite",
                      notEvaluated: "Range not yet evaluated")
  STYPE  ==> Union(str:SDF,
                     notEvaluated:"Internal singularities not yet evaluated")
  ATT  ==> Record(endPointContinuity:CTYPE,singularitiesStream:STYPE,range:RTYPE)
  IFV  ==> Record(stiffness:F,stability:F,expense:F,accuracy:F,intermediateResults:F)

  SIG ==> with

    f2df : F -> DF
      ++ f2df(f) is a function to convert a \axiomType{Float} to a
      ++ \axiomType{DoubleFloat}

    ef2edf : EF -> EDF
      ++ ef2edf(f) is a function to convert an \axiomType{Expression Float} 
      ++ to an \axiomType{Expression DoubleFloat}

    ocf2ocdf : OCF -> OCDF
      ++ ocf2ocdf(a) is a function to convert an \axiomType{OrderedCompletion 
      ++ Float} to an \axiomType{OrderedCompletion DoubleFloat}

    socf2socdf : SOCF -> SOCDF
      ++ socf2socdf(a) is a function to convert a \axiomType{Segment OrderedCompletion Float}
      ++ to a \axiomType{Segment OrderedCompletion DoubleFloat}

    convert : List SOCF -> List SOCDF
      ++ convert(l) is a function to convert a \axiomType{Segment OrderedCompletion Float}
      ++ to a \axiomType{Segment OrderedCompletion DoubleFloat}

    df2fi : DF -> FI
      ++ df2fi(n) is a function to convert a \axiomType{DoubleFloat} to a
      ++ \axiomType{Fraction Integer}

    edf2fi : EDF -> FI
      ++ edf2fi(n) maps \axiomType{Expression DoubleFloat} to 
      ++ \axiomType{Fraction Integer}
      ++ It is an error if n is not coercible to Fraction Integer

    edf2df : EDF -> DF
      ++ edf2df(n) maps \axiomType{Expression DoubleFloat} to 
      ++ \axiomType{DoubleFloat}
      ++ It is an error if \spad{n} is not coercible to DoubleFloat

    isQuotient : EDF -> Union(EDF,"failed")
      ++ isQuotient(expr) returns the quotient part of the input
      ++ expression or \spad{"failed"} if the expression is not of that form.

    expenseOfEvaluation : VEDF -> F
      ++ expenseOfEvaluation(o) gives an approximation of the cost of
      ++ evaluating a list of expressions in terms of the number of basic
      ++ operations.
      ++ < 0.3 inexpensive ; 0.5 neutral ; > 0.7 very expensive
      ++ 400 `operation units' -> 0.75 
      ++ 200 `operation units' -> 0.5 
      ++ 83 `operation units' -> 0.25
      ++ ** = 4 units , function calls = 10 units.

    numberOfOperations : VEDF -> ON
      ++ numberOfOperations(ode) counts additions, multiplications, 
      ++ exponentiations and function calls in the input set of expressions.

    edf2efi : EDF -> EFI
      ++ edf2efi(e) coerces \axiomType{Expression DoubleFloat} into 
      ++ \axiomType{Expression Fraction Integer}

    dfRange : SOCDF -> SOCDF
      ++ dfRange(r) converts a range including 
      ++ \inputbitmap{\htbmdir{}/plusminus.bitmap} \infty 
      ++ to \axiomType{DoubleFloat} equavalents.

    dflist : List(Record(left:FI,right:FI)) -> LDF
      ++ dflist(l) returns a list of \axiomType{DoubleFloat} equivalents of list l

    df2mf : DF -> MF
      ++ df2mf(n) coerces a \axiomType{DoubleFloat} to \axiomType{MachineFloat}

    ldf2vmf : LDF -> VMF
      ++ ldf2vmf(l) coerces a \axiomType{List DoubleFloat} to
      ++ \axiomType{List MachineFloat}

    edf2ef : EDF -> EF
      ++ edf2ef(e) maps \axiomType{Expression DoubleFloat} to 
      ++ \axiomType{Expression Float}

    vedf2vef : VEDF -> VEF
      ++ vedf2vef(v) maps \axiomType{Vector Expression DoubleFloat} to 
      ++ \axiomType{Vector Expression Float}

    in? : (DF,SOCDF) -> Boolean
      ++ in?(p,range) tests whether point p is internal to the 
      ++ range range

    df2st : DF -> ST
      ++ df2st(n) coerces a \axiomType{DoubleFloat} to \axiomType{String}

    f2st : F -> ST
      ++ f2st(n) coerces a \axiomType{Float} to \axiomType{String}

    ldf2lst : LDF -> LST
      ++ ldf2lst(ln) coerces a \axiomType{List DoubleFloat} to \axiomType{List String}

    sdf2lst : SDF -> LST
      ++ sdf2lst(ln) coerces a \axiomType{Stream DoubleFloat} to \axiomType{String}

    getlo : SOCDF -> DF
      ++ getlo(u) gets the \axiomType{DoubleFloat} equivalent of
      ++ the first endpoint of the range \spad{u}

    gethi : SOCDF -> DF
      ++ gethi(u) gets the \axiomType{DoubleFloat} equivalent of
      ++ the second endpoint of the range \spad{u}

    concat : (Result,Result) -> Result
      ++ concat(a,b) adds two aggregates of type \axiomType{Result}.

    concat : (List Result) -> Result
      ++ concat(l) concatenates a list of aggregates of type \axiomType{Result}

    outputMeasure : F -> ST
      ++ outputMeasure(n) rounds \spad{n} to 3 decimal places and outputs
      ++ it as a string

    measure2Result : Measure -> Result
      ++ measure2Result(m) converts a measure record into a \axiomType{Result}

    measure2Result : Measure2 -> Result
      ++ measure2Result(m) converts a measure record into a \axiomType{Result}

    att2Result : ATT -> Result
      ++ att2Result(m) converts a attributes record into a \axiomType{Result}

    iflist2Result : IFV -> Result
      ++ iflist2Result(m) converts attributes record into a \axiomType{Result}

    pdf2ef : PDF -> EF
      ++ pdf2ef(p) coerces a \axiomType{Polynomial DoubleFloat} to 
      ++ \axiomType{Expression Float}

    pdf2df : PDF -> DF
      ++ pdf2df(p) coerces a \axiomType{Polynomial DoubleFloat} to 
      ++ \axiomType{DoubleFloat}. It is an error if \axiom{p} is not
      ++ retractable to DoubleFloat.

    df2ef : DF -> EF
      ++ df2ef(a) coerces a \axiomType{DoubleFloat} to 
      ++ \axiomType{Expression Float}

    fi2df : FI -> DF
      ++ fi2df(f) coerces a \axiomType{Fraction Integer} to 
      ++ \axiomType{DoubleFloat}

    mat : (LDF,NNI) -> MDF
      ++ mat(a,n) constructs a one-dimensional matrix of a.

  CODE ==> add

    mat(a:LDF,n:NNI):MDF ==
      empty?(a)$LDF => zero(1,n)$MDF
      matrix(list([i for i in a for j in 1..n])$(List LDF))$MDF

    f2df(f:F):DF == (convert(f)@DF)$F

    ef2edf(f:EF):EDF == map(f2df,f)$EF2(F,DF)

    fi2df(f:FI):DF == coerce(f)$DF

    ocf2ocdf(a:OCF):OCDF ==
      finite? a => (f2df(retract(a)@F))::OCDF
      a pretend OCDF

    socf2socdf(a:SOCF):SOCDF ==
      segment(ocf2ocdf(lo a),ocf2ocdf(hi a))

    convert(l:List SOCF):List SOCDF == [socf2socdf a for a in l]

    pdf2df(p:PDF):DF == retract(p)@DF

    df2ef(a:DF):EF ==
      b := convert(a)@Float
      coerce(b)$EF

    pdf2ef(p:PDF):EF == df2ef(pdf2df(p))

    edf2fi(m:EDF):FI == retract(retract(m)@DF)@FI

    edf2df(m:EDF):DF == retract(m)@DF

    df2fi(r:DF):FI == (retract(r)@FI)$DF

    dfRange(r:SOCDF):SOCDF ==
      if infinite?(lo(r))$OCDF then r := -(max()$DF :: OCDF)..hi(r)$SOCDF
      if infinite?(hi(r))$OCDF then r := lo(r)$SOCDF..(max()$DF :: OCDF)
      r

    dflist(l:List(Record(left:FI,right:FI))):LDF == [u.left :: DF for u in l]

    edf2efi(f:EDF):EFI == map(df2fi,f)$EF2(DF,FI)

    df2st(n:DF):String == (convert((convert(n)@Float)$DF)@ST)$Float

    f2st(n:F):String == (convert(n)@ST)$Float

    ldf2lst(ln:LDF):LST == [df2st f for f in ln]

    sdf2lst(ln:SDF):LST ==
      explicitlyFinite? ln => 
        m := map(df2st,ln)$StreamFunctions2(DF,ST)
        if index?(20,m)$SS then
          split!(m,20)
          m := concat(m,".......")
        m := complete(m)$SS 
        entries(m)$SS
      empty()$LST

    df2mf(n:DF):MF == (df2fi(n))::MF

    ldf2vmf(l:LDF):VMF ==
      m := [df2mf(n) for n in l]
      vector(m)$VMF

    edf2ef(e:EDF):EF == map(convert$DF,e)$EF2(DF,Float)

    vedf2vef(vedf:VEDF):VEF == vector([edf2ef e for e in members(vedf)])

    getlo(u:SOCDF):DF == retract(lo(u))@DF

    gethi(u:SOCDF):DF == retract(hi(u))@DF
  
    in?(p:DF,range:SOCDF):Boolean ==
      top := gethi(range)
      bottom := getlo(range)
      a:Boolean := (p < top)$DF
      b:Boolean := (p > bottom)$DF
      (a and b)@Boolean

    isQuotient(expr:EDF):Union(EDF,"failed") ==
      (k := mainKernel expr) case KEDF =>
        (expr = inv(f := k :: KEDF :: EDF)$EDF)$EDF => f
        (numerator expr) = 1 => denominator expr
        "failed"
      "failed"

    numberOfOperations1(fn:EDF,numbersSoFar:ON):ON ==
      (u := isQuotient(fn)) case EDF =>
        numbersSoFar := numberOfOperations1(u,numbersSoFar)
      (p := isPlus(fn)) case LEDF =>
        p := coerce(p)@LEDF
        np := #p
        numbersSoFar.additions := (numbersSoFar.additions)+np-1
        for i in 1..np repeat
          numbersSoFar := numberOfOperations1(p.i,numbersSoFar)
        numbersSoFar
      (t:=isTimes(fn)) case LEDF => 
        t := coerce(t)@LEDF
        nt := #t
        numbersSoFar.multiplications := (numbersSoFar.multiplications)+nt-1
        for i in 1..nt repeat
          numbersSoFar := numberOfOperations1(t.i,numbersSoFar)
        numbersSoFar
      if (e:=isPower(fn)) case RVE then
        e := coerce(e)@RVE
        e.exponent>1 =>  
          numbersSoFar.exponentiations := inc(numbersSoFar.exponentiations)
          numbersSoFar := numberOfOperations1(e.val,numbersSoFar)
      lk := kernels(fn)
      #lk = 1 =>        -- #lk = 0 => constant found (no further action)
        k := first(lk)$LKEDF
        n := name(operator(k)$KEDF)$BO
        entry?(n,variables(fn)$EDF)$LS => numbersSoFar  -- solo variable found
        a := first(argument(k)$KEDF)$LEDF
        numbersSoFar.functionCalls := inc(numbersSoFar.functionCalls)$INT
        numbersSoFar := numberOfOperations1(a,numbersSoFar)
      numbersSoFar
      
    numberOfOperations(ode:VEDF):ON ==
      n:ON := [0,0,0,0]
      for i in 1..#ode repeat
        n:ON := numberOfOperations1(ode.i,n)
      n

    expenseOfEvaluation(o:VEDF):F ==
      ln:ON := numberOfOperations(o)
      a := ln.additions
      m := ln.multiplications
      e := ln.exponentiations
      f := 10*ln.functionCalls
      n := (a + m + 4*e + 10*e)
      (1.0-exp((-n::F/288.0))$F)

    concat(a:Result,b:Result):Result ==
      membersOfa := (members(a)@List(Record(key:Symbol,entry:Any)))
      membersOfb := (members(b)@List(Record(key:Symbol,entry:Any)))
      allMembers:=
        concat(membersOfa,membersOfb)$List(Record(key:Symbol,entry:Any))
      construct(allMembers)

    concat(l:List Result):Result ==
      import List Result
      empty? l => empty()$Result
      f := first l
      if empty?(r := rest l) then
        f
      else
        concat(f,concat r)

    outputMeasure(m:F):ST ==
      fl:Float := round(m*(f:= 1000.0))/f
      convert(fl)@ST

    measure2Result(m:Measure):Result ==
      mm := coerce(m.measure)$AnyFunctions1(Float)
      mmr:Record(key:Symbol,entry:Any) := [bestMeasure@Symbol,mm]
      mn := coerce(m.name)$AnyFunctions1(ST)
      mnr:Record(key:Symbol,entry:Any) := [nameOfRoutine@Symbol,mn]
      me := coerce(m.explanations)$AnyFunctions1(List String)
      mer:Record(key:Symbol,entry:Any) := [allMeasures@Symbol,me]
      mr := construct([mmr,mnr,mer])$Result
      met := coerce(mr)$AnyFunctions1(Result)
      meth:Record(key:Symbol,entry:Any):=[method@Symbol,met]
      construct([meth])$Result

    measure2Result(m:Measure2):Result ==
      mm := coerce(m.measure)$AnyFunctions1(Float)
      mmr:Record(key:Symbol,entry:Any) := [bestMeasure@Symbol,mm]
      mn := coerce(m.name)$AnyFunctions1(ST)
      mnr:Record(key:Symbol,entry:Any) := [nameOfRoutine@Symbol,mn]
      me := coerce(m.explanations)$AnyFunctions1(List String)
      mer:Record(key:Symbol,entry:Any) := [allMeasures@Symbol,me]
      mx := coerce(m.extra)$AnyFunctions1(Result)
      mxr:Record(key:Symbol,entry:Any) := [other@Symbol,mx]
      mr := construct([mmr,mnr,mer,mxr])$Result
      met := coerce(mr)$AnyFunctions1(Result)
      meth:Record(key:Symbol,entry:Any):=[method@Symbol,met]
      construct([meth])$Result

    att2Result(att:ATT):Result ==
      aepc := coerce(att.endPointContinuity)$AnyFunctions1(CTYPE)
      ar := coerce(att.range)$AnyFunctions1(RTYPE)
      as := coerce(att.singularitiesStream)$AnyFunctions1(STYPE)
      aa:List Any := [aepc,ar,as]
      aaa := coerce(aa)$AnyFunctions1(List Any)
      aar:Record(key:Symbol,entry:Any) := [attributes@Symbol,aaa]
      construct([aar])$Result

    iflist2Result(ifv:IFV):Result ==
      ifvs:List String := 
        [concat(["stiffness: ",outputMeasure(ifv.stiffness)]),
          concat(["stability: ",outputMeasure(ifv.stability)]),
           concat(["expense: ",outputMeasure(ifv.expense)]),
            concat(["accuracy: ",outputMeasure(ifv.accuracy)]),
             concat(["intermediateResults: ",_
                     outputMeasure(ifv.intermediateResults)])]
      ifa:= coerce(ifvs)$AnyFunctions1(List String)
      ifr:Record(key:Symbol,entry:Any) := [intensityFunctions@Symbol,ifa]
      construct([ifr])$Result