--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

)abbrev package ESTOOLS ExpertSystemToolsPackage
++ Author: Brian Dupee
++ Date Created: May 1994
++ Date Last Updated: July 1996
++ Basic Operations:
++ Description:
++ \axiom{ExpertSystemToolsPackage} contains some useful functions for use
++ by the computational agents of numerical solvers.
ExpertSystemToolsPackage():E == I 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)

  E ==> 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 a 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.

  I ==> 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
        one?(numerator expr) => 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

)abbrev package ESTOOLS1 ExpertSystemToolsPackage1
++ Author: Brian Dupee
++ Date Created: February 1995
++ Date Last Updated: February 1995
++ Basic Operations: neglist
++ Description:
++ \axiom{ExpertSystemToolsPackage1} contains some useful functions for use
++ by the computational agents of Ordinary Differential Equation solvers.
ExpertSystemToolsPackage1(R1:OR): E == I where
  OR	==> OrderedRing
  E ==>  with
    neglist:List R1 -> List R1
      ++ neglist(l) returns only the negative elements of the list \spad{l}
  I ==> add
    neglist(l:List R1):List R1 == [u for u in l | negative?(u)$R1]

)abbrev package ESTOOLS2 ExpertSystemToolsPackage2
++ Author: Brian Dupee
++ Date Created: February 1995
++ Date Last Updated: July 1996
++ Basic Operations: map
++ Related Constructors: Matrix
++ Description:
++ \axiom{ExpertSystemToolsPackage2} contains some useful functions for use
++ by the computational agents of Ordinary Differential Equation solvers.
ExpertSystemToolsPackage2(R1:R,R2:R): E == I where
  R	==> Ring
  E ==>  with
    map:(R1->R2,Matrix R1) -> Matrix R2
      ++ map(f,m) applies a mapping f:R1 -> R2 onto a matrix
      ++ \spad{m} in R1 returning a matrix in R2
  I ==> add
    map(f:R1->R2,m:Matrix R1):Matrix R2 ==
      matrix([[f u for u in v] for v in listOfLists(m)$(Matrix R1)])$(Matrix R2)

