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

)abbrev package INTPACK AnnaNumericalIntegrationPackage
++ Author: Brian Dupee
++ Date Created: August 1994
++ Date Last Updated: December 1997
++ Description:
++ \axiomType{AnnaNumericalIntegrationPackage} is a \axiom{package}
++ of functions for the \axiom{category} 
++ \axiomType{NumericalIntegrationCategory} 
++ with \axiom{measure}, and \axiom{integrate}.

AnnaNumericalIntegrationPackage() : SIG == CODE where 

 EDF  ==> Expression DoubleFloat
 DF  ==> DoubleFloat
 EF  ==> Expression Float
 F  ==> Float
 INT  ==> Integer
 SOCDF  ==> Segment OrderedCompletion DoubleFloat
 OCDF  ==> OrderedCompletion DoubleFloat
 SBOCF  ==> SegmentBinding OrderedCompletion Float
 LSOCF  ==> List Segment OrderedCompletion Float
 SOCF  ==> Segment OrderedCompletion Float
 OCF  ==> OrderedCompletion Float
 LS  ==> List Symbol
 S  ==> Symbol
 LST  ==> List String
 ST  ==> String
 RT  ==> RoutinesTable
 NIA  ==> Record(var:S, fn:EDF, range:SOCDF, abserr:DF, relerr:DF)
 MDNIA  ==> Record(fn:EDF,range:List SOCDF,abserr:DF,relerr:DF)
 IFL  ==> List(Record(ifail:Integer,instruction:String))
 Entry  ==> Record(chapter:String, type:String, domainName: String, 
                   defaultMin:F, measure:F, failList:IFL, explList:List String)
 Measure  ==> Record(measure:F, name:ST, explanations:LST, extra:Result)

 SIG ==> with

  integrate : (EF,SOCF,F,F,RT) -> Result
    ++ integrate(exp, a..b, epsrel, routines) is a top level ANNA function 
    ++ to integrate an expression, {\tt exp}, over a given range {\tt a} 
    ++ to {\tt b} to the required absolute and relative accuracy using 
    ++ the routines available in the RoutinesTable provided.
    ++ 
    ++ It iterates over the \axiom{domains} of 
    ++ \axiomType{NumericalIntegrationCategory} 
    ++ to get the name and other 
    ++ relevant information of the the (domain of the) numerical 
    ++ routine likely to be the most appropriate, 
    ++ have the best \axiom{measure}.
    ++ 
    ++ It then performs the integration of the given expression 
    ++ on that \axiom{domain}.  

  integrate : NumericalIntegrationProblem -> Result
    ++ integrate(IntegrationProblem) is a top level ANNA function 
    ++ to integrate an expression over a given range or ranges 
    ++ to the required absolute and relative accuracy.
    ++ 
    ++ It iterates over the \axiom{domains} of 
    ++ \axiomType{NumericalIntegrationCategory} to get the name and other 
    ++ relevant information of the the (domain of the) numerical 
    ++ routine likely to be the most appropriate, 
    ++ have the best \axiom{measure}.
    ++ 
    ++ It then performs the integration of the given expression 
    ++ on that \axiom{domain}.  

  integrate : (EF,SOCF,F,F) -> Result
    ++ integrate(exp, a..b, epsabs, epsrel) is a top level ANNA function 
    ++ to integrate an expression, {\tt exp}, over a given range {\tt a} 
    ++ to {\tt b} to the required absolute and relative accuracy.
    ++ 
    ++ It iterates over the \axiom{domains} of 
    ++ \axiomType{NumericalIntegrationCategory} to get the name and other 
    ++ relevant information of the the (domain of the) numerical 
    ++ routine likely to be the most appropriate, 
    ++ have the best \axiom{measure}.
    ++ 
    ++ It then performs the integration of the given expression 
    ++ on that \axiom{domain}.  

  integrate : (EF,SOCF,F) -> Result
    ++ integrate(exp, a..b, epsrel) is a top level ANNA 
    ++ function to integrate an expression, {\tt exp}, over a given
    ++ range {\tt a} to {\tt b} to the required relative accuracy.
    ++ 
    ++ It iterates over the \axiom{domains} of 
    ++ \axiomType{NumericalIntegrationCategory} to get the name and other 
    ++ relevant information of the the (domain of the) numerical 
    ++ routine likely to be the most appropriate, 
    ++ have the best \axiom{measure}.
    ++ 
    ++ It then performs the integration of the given expression 
    ++ on that \axiom{domain}.  
    ++
    ++ If epsrel = 0, a default absolute accuracy is used.
 
  integrate : (EF,SOCF) -> Result
    ++ integrate(exp, a..b) is a top 
    ++ level ANNA function to integrate an expression, {\tt exp},
    ++ over a given range {\tt a} to {\tt b}.
    ++ 
    ++ It iterates over the \axiom{domains} of 
    ++ \axiomType{NumericalIntegrationCategory} to get the name and other 
    ++ relevant information of the the (domain of the) numerical 
    ++ routine likely to be the most appropriate, 
    ++ have the best \axiom{measure}.
    ++ 
    ++ It then performs the integration of the given expression 
    ++ on that \axiom{domain}.  
    ++
    ++ Default values for the absolute and relative error are used.

  integrate : (EF,LSOCF) -> Result
    ++ integrate(exp, [a..b,c..d,...]) is a top 
    ++ level ANNA function to integrate a multivariate expression, {\tt exp},
    ++ over a given set of ranges.
    ++ 
    ++ It iterates over the \axiom{domains} of 
    ++ \axiomType{NumericalIntegrationCategory} to get the name and other 
    ++ relevant information of the the (domain of the) numerical 
    ++ routine likely to be the most appropriate, 
    ++ have the best \axiom{measure}.
    ++ 
    ++ It then performs the integration of the given expression 
    ++ on that \axiom{domain}.  
    ++
    ++ Default values for the absolute and relative error are used.

  integrate : (EF,LSOCF,F) -> Result
    ++ integrate(exp, [a..b,c..d,...], epsrel) is a top 
    ++ level ANNA function to integrate a multivariate expression, {\tt exp},
    ++ over a given set of ranges to the required relative
    ++ accuracy.
    ++ 
    ++ It iterates over the \axiom{domains} of 
    ++ \axiomType{NumericalIntegrationCategory} to get the name and other 
    ++ relevant information of the the (domain of the) numerical 
    ++ routine likely to be the most appropriate, 
    ++ have the best \axiom{measure}.
    ++ 
    ++ It then performs the integration of the given expression 
    ++ on that \axiom{domain}.  
    ++
    ++ If epsrel = 0, a default absolute accuracy is used.

  integrate : (EF,LSOCF,F,F) -> Result
    ++ integrate(exp, [a..b,c..d,...], epsabs, epsrel) is a top 
    ++ level ANNA function to integrate a multivariate expression, {\tt exp},
    ++ over a given set of ranges to the required absolute and relative
    ++ accuracy.
    ++ 
    ++ It iterates over the \axiom{domains} of 
    ++ \axiomType{NumericalIntegrationCategory} to get the name and other 
    ++ relevant information of the the (domain of the) numerical 
    ++ routine likely to be the most appropriate, 
    ++ have the best \axiom{measure}.
    ++ 
    ++ It then performs the integration of the given expression 
    ++ on that \axiom{domain}.

  integrate : (EF,LSOCF,F,F,RT) -> Result
    ++ integrate(exp, [a..b,c..d,...], epsabs, epsrel, routines) is a top 
    ++ level ANNA function to integrate a multivariate expression, {\tt exp},
    ++ over a given set of ranges to the required absolute and relative
    ++ accuracy, using the routines available in the RoutinesTable provided.
    ++ 
    ++ It iterates over the \axiom{domains} of 
    ++ \axiomType{NumericalIntegrationCategory} to get the name and other 
    ++ relevant information of the the (domain of the) numerical 
    ++ routine likely to be the most appropriate, 
    ++ have the best \axiom{measure}.
    ++ 
    ++ It then performs the integration of the given expression 
    ++ on that \axiom{domain}.

  measure : NumericalIntegrationProblem -> Measure
    ++ measure(prob) is a top level ANNA function for identifying the most
    ++ appropriate numerical routine for solving the numerical integration
    ++ problem defined by \axiom{prob}.
    ++
    ++ It calls each \axiom{domain} of \axiom{category}
    ++ \axiomType{NumericalIntegrationCategory} in turn to calculate all measures
    ++ and returns the best 
    ++ the name of the most appropriate domain and any other relevant
    ++ information.

  measure : (NumericalIntegrationProblem,RT) -> Measure
    ++ measure(prob,R) is a top level ANNA function for identifying the most
    ++ appropriate numerical routine from those in the routines table
    ++ provided for solving the numerical integration
    ++ problem defined by \axiom{prob}.
    ++
    ++ It calls each \axiom{domain} listed in \axiom{R} of \axiom{category}
    ++ \axiomType{NumericalIntegrationCategory} in turn to calculate all measures
    ++ and returns the best 
    ++ the name of the most appropriate domain and any other relevant
    ++ information.

  integrate : (EF,SBOCF,ST) -> Union(Result,"failed")
    ++ integrate(exp, x = a..b, "numerical") is a top level ANNA function to 
    ++ integrate an expression, {\tt exp}, over a given range, {\tt a}
    ++ to {\tt b}.
    ++ 
    ++ It iterates over the \axiom{domains} of 
    ++ \axiomType{NumericalIntegrationCategory} to get the name and other 
    ++ relevant information of the the (domain of the) numerical 
    ++ routine likely to be the most appropriate, 
    ++ have the best \axiom{measure}.
    ++ 
    ++ It then performs the integration of the given expression 
    ++ on that \axiom{domain}.
    ++ 
    ++ Default values for the absolute and relative error are used.
    ++
    ++ It is an error of the last argument is not {\tt "numerical"}.

  integrate : (EF,SBOCF,S) -> Union(Result,"failed")
    ++ integrate(exp, x = a..b, numerical) is a top level ANNA function to 
    ++ integrate an expression, {\tt exp}, over a given range, {\tt a}
    ++ to {\tt b}.
    ++ 
    ++ It iterates over the \axiom{domains} of 
    ++ \axiomType{NumericalIntegrationCategory} to get the name and other 
    ++ relevant information of the the (domain of the) numerical 
    ++ routine likely to be the most appropriate, 
    ++ have the best \axiom{measure}.
    ++ 
    ++ It then performs the integration of the given expression 
    ++ on that \axiom{domain}.
    ++ 
    ++ Default values for the absolute and relative error are used.
    ++
    ++ It is an error if the last argument is not {\tt numerical}.

 CODE ==> add

  zeroMeasure: Measure -> Result
  scriptedVariables?: MDNIA -> Boolean
  preAnalysis:(Union(nia:NIA,mdnia:MDNIA),RT) -> RT
  measureSpecific:(ST,RT,Union(nia:NIA,mdnia:MDNIA)) -> _
      Record(measure:F,explanations:LST,extra:Result)
  changeName:(Result,ST) -> Result
  recoverAfterFail:(Union(nia:NIA,mdnia:MDNIA),RT,Measure,INT,Result) -> _
      Record(a:Result,b:Measure)
  better?:(Result,Result) -> Boolean
  integrateConstant:(EF,SOCF) -> Result
  integrateConstantList: (EF,LSOCF) -> Result
  integrateArgs:(NumericalIntegrationProblem,RT) -> Result
  integrateSpecific:(Union(nia:NIA,mdnia:MDNIA),ST,Result) -> Result

  import ExpertSystemToolsPackage

  integrateConstantList(exp:EF,ras:LSOCF):Result ==
    c:OCF := ((retract(exp)@F)$EF)::OCF
    b := [hi(j)-lo(j) for j in ras]
    c := c*reduce((x,y) +-> x*y,b)
    a := coerce(c)$AnyFunctions1(OCF)
    text := coerce("Constant Function")$AnyFunctions1(ST)
    construct([[result@S,a],[method@S,text]])$Result
    
  integrateConstant(exp:EF,ra:SOCF):Result ==
    c := (retract(exp)@F)$EF
    r:OCF := (c::OCF)*(hi(ra)-lo(ra))
    a := coerce(r)$AnyFunctions1(OCF)
    text := coerce("Constant Function")$AnyFunctions1(ST)
    construct([[result@S,a],[method@S,text]])$Result
    
  zeroMeasure(m:Measure):Result ==
    a := coerce(0$DF)$AnyFunctions1(DF)
    text := coerce("Constant Function")$AnyFunctions1(String)
    r := construct([[result@Symbol,a],[method@Symbol,text]])$Result
    concat(measure2Result m,r)$ExpertSystemToolsPackage

  scriptedVariables?(mdnia:MDNIA):Boolean ==
    vars:List Symbol := variables(mdnia.fn)$EDF
    var1 := first(vars)$(List Symbol)
    not scripted?(var1) => false
    name1 := name(var1)$Symbol
    for i in 2..# vars repeat
      not ((scripted?(vars.i)$Symbol) and (name1 = name(vars.i)$Symbol)) => 
         return false
    true

  preAnalysis(args:Union(nia:NIA,mdnia:MDNIA),t:RT):RT ==
    import RT
    r:RT := selectIntegrationRoutines t
    args case nia => 
      arg:NIA := args.nia
      rangeIsFinite(arg)$d01AgentsPackage case finite => 
        selectFiniteRoutines r
      selectNonFiniteRoutines r
    selectMultiDimensionalRoutines r
    
  changeName(ans:Result,name:ST):Result ==
    sy:S := coerce(name "Answer")$S
    anyAns:Any := coerce(ans)$AnyFunctions1(Result)
    construct([[sy,anyAns]])$Result

  measureSpecific(name:ST,R:RT,args:Union(nia:NIA,mdnia:MDNIA)):
      Record(measure:F,explanations:ST,extra:Result) ==
    args case nia => 
      arg:NIA := args.nia
      name = "d01ajfAnnaType" => measure(R,arg)$d01ajfAnnaType
      name = "d01akfAnnaType" => measure(R,arg)$d01akfAnnaType
      name = "d01alfAnnaType" => measure(R,arg)$d01alfAnnaType
      name = "d01amfAnnaType" => measure(R,arg)$d01amfAnnaType
      name = "d01anfAnnaType" => measure(R,arg)$d01anfAnnaType
      name = "d01apfAnnaType" => measure(R,arg)$d01apfAnnaType
      name = "d01aqfAnnaType" => measure(R,arg)$d01aqfAnnaType
      name = "d01asfAnnaType" => measure(R,arg)$d01asfAnnaType
      name = "d01TransformFunctionType" => 
                     measure(R,arg)$d01TransformFunctionType
      error("measureSpecific","invalid type name: " name)$ErrorFunctions
    args case mdnia => 
      arg2:MDNIA := args.mdnia
      name = "d01gbfAnnaType" => measure(R,arg2)$d01gbfAnnaType
      name = "d01fcfAnnaType" => measure(R,arg2)$d01fcfAnnaType
      error("measureSpecific","invalid type name: " name)$ErrorFunctions
    error("measureSpecific","invalid type name")$ErrorFunctions

  measure(a:NumericalIntegrationProblem,R:RT):Measure ==
    args:Union(nia:NIA,mdnia:MDNIA) := retract(a)$NumericalIntegrationProblem
    sofar := 0$F
    best := "none" :: ST
    routs := copy R
    routs := preAnalysis(args,routs)
    empty?(routs)$RT => 
      error("measure", "no routines found")$ErrorFunctions
    rout := inspect(routs)$RT
    e := retract(rout.entry)$AnyFunctions1(Entry)
    meth:LST := ["Trying " e.type " integration routines"]
    ext := empty()$Result
    for i in 1..# routs repeat
      rout := extract!(routs)$RT
      e := retract(rout.entry)$AnyFunctions1(Entry)
      n := e.domainName
      if e.defaultMin > sofar then
        m := measureSpecific(n,R,args)
        if m.measure > sofar then
          sofar := m.measure
          best := n
        ext := concat(m.extra,ext)$ExpertSystemToolsPackage
        str:LST := [string(rout.key)$S "measure: " outputMeasure(m.measure) 
                     " - " m.explanations]
      else
        str:LST :=  [string(rout.key)$S " is no better than other routines"]
      meth := append(meth,str)$LST
    [sofar,best,meth,ext]

  measure(a:NumericalIntegrationProblem):Measure ==
    measure(a,routines()$RT)

  integrateSpecific(args:Union(nia:NIA,mdnia:MDNIA),n:ST,ex:Result):Result ==
    args case nia => 
      arg:NIA := args.nia
      n = "d01ajfAnnaType" => numericalIntegration(arg,ex)$d01ajfAnnaType
      n = "d01TransformFunctionType" =>
        numericalIntegration(arg,ex)$d01TransformFunctionType
      n = "d01amfAnnaType" => numericalIntegration(arg,ex)$d01amfAnnaType
      n = "d01apfAnnaType" => numericalIntegration(arg,ex)$d01apfAnnaType
      n = "d01aqfAnnaType" => numericalIntegration(arg,ex)$d01aqfAnnaType
      n = "d01alfAnnaType" => numericalIntegration(arg,ex)$d01alfAnnaType
      n = "d01akfAnnaType" => numericalIntegration(arg,ex)$d01akfAnnaType
      n = "d01anfAnnaType" => numericalIntegration(arg,ex)$d01anfAnnaType
      n = "d01asfAnnaType" => numericalIntegration(arg,ex)$d01asfAnnaType
      error("integrateSpecific","invalid type name: " n)$ErrorFunctions
    args case mdnia => 
      arg2:MDNIA := args.mdnia
      n = "d01gbfAnnaType" => numericalIntegration(arg2,ex)$d01gbfAnnaType
      n = "d01fcfAnnaType" => numericalIntegration(arg2,ex)$d01fcfAnnaType
      error("integrateSpecific","invalid type name: " n)$ErrorFunctions
    error("integrateSpecific","invalid type name: " n)$ErrorFunctions

  better?(r:Result,s:Result):Boolean ==
    a1 := search("abserr"::S,r)$Result
    a1 case "failed" => false
    abserr1 := retract(a1)$AnyFunctions1(DF)
    negative?(abserr1) => false
    a2 := search("abserr"::S,s)$Result
    a2 case "failed" => true
    abserr2 := retract(a2)$AnyFunctions1(DF)
    negative?(abserr2) => true
    (abserr1 < abserr2) -- true if r.abserr better than s.abserr

  recoverAfterFail(n:Union(nia:NIA,mdnia:MDNIA),routs:RT,m:Measure,iint:INT,
                         r:Result):Record(a:Result,b:Measure) ==
    bestName := m.name
    while positive?(iint) repeat
      routineName := m.name
      s := recoverAfterFail(routs,routineName(1..6),iint)$RoutinesTable
      s case "failed" => iint := 0
      if s = "changeEps" then
        nn := n.nia
        zero?(nn.abserr) =>
          nn.abserr := 1.0e-8 :: DF
          m := measure(n::NumericalIntegrationProblem,routs)
          zero?(m.measure) => iint := 0
          r := integrateSpecific(n,m.name,m.extra)
          iint := 0
      rn := routineName(1..6)
      buttVal := getButtonValue(rn,"functionEvaluations")$AttributeButtons
      if (s = "incrFunEvals") and (buttVal < 0.8) then
        increase(rn,"functionEvaluations")$AttributeButtons
      if s = "increase tolerance" then
        (n.nia).relerr := (n.nia).relerr*(10.0::DF)
      if s = "decrease tolerance" then
        (n.nia).relerr := (n.nia).relerr/(10.0::DF)
      fl := coerce(s)$AnyFunctions1(ST)
      flrec:Record(key:S,entry:Any):=[failure@S,fl]
      m2 := measure(n::NumericalIntegrationProblem,routs)
      zero?(m2.measure) => iint := 0
      r2:Result := integrateSpecific(n,m2.name,m2.extra)
      better?(r,r2) => 
        m.name := m2.name
        insert!(flrec,r)$Result
      bestName := m2.name
      m := m2
      insert!(flrec,r2)$Result
      r := concat(r2,changeName(r,routineName))$ExpertSystemToolsPackage
      iany := search(ifail@S,r2)$Result
      iany case "failed" => iint := 0
      iint := retract(iany)$AnyFunctions1(INT)
    m.name := bestName
    [r,m]

  integrateArgs(prob:NumericalIntegrationProblem,t:RT):Result ==
    args:Union(nia:NIA,mdnia:MDNIA):= retract(prob)$NumericalIntegrationProblem
    routs := copy(t)$RT
    if args case mdnia then
      arg := args.mdnia
      v := (# variables(arg.fn))
      not scriptedVariables?(arg) => 
        error("MultiDimensionalNumericalIntegrationPackage",
                "invalid variable names")$ErrorFunctions
      (v ~= # arg.range)@Boolean =>
        error("MultiDimensionalNumericalIntegrationPackage",
          "number of variables do not match number of ranges")$ErrorFunctions
    m := measure(prob,routs)
    zero?(m.measure) => zeroMeasure m
    r := integrateSpecific(args,m.name,m.extra)
    iany := search(ifail@S,r)$Result
    iint := 0$INT
    if (iany case Any) then
      iint := retract(iany)$AnyFunctions1(INT)
    if positive?(iint) then
      tu:Record(a:Result,b:Measure) := recoverAfterFail(args,routs,m,iint,r)
      r := tu.a
      m := tu.b
    r := concat(measure2Result m,r)$ExpertSystemToolsPackage
    n := m.name
    nn:ST := 
      (# n > 14) => "d01transform"
      n(1..6)
    expl := getExplanations(routs,nn)$RoutinesTable
    expla := coerce(expl)$AnyFunctions1(LST)
    explaa:Record(key:Symbol,entry:Any) := ["explanations"::Symbol,expla]
    r := concat(construct([explaa]),r)
    args case nia =>
      att := showAttributes(args.nia)$IntegrationFunctionsTable
      att case "failed" => r
      concat(att2Result att,r)$ExpertSystemToolsPackage
    r

  integrate(args:NumericalIntegrationProblem):Result ==
    integrateArgs(args,routines()$RT)

  integrate(exp:EF,ra:SOCF,epsabs:F,epsrel:F,r:RT):Result ==
    Var:LS := variables(exp)$EF
    empty?(Var)$LS => integrateConstant(exp,ra)
    args:NIA:= [first(Var)$LS,ef2edf exp,socf2socdf ra,f2df epsabs,f2df epsrel]
    integrateArgs(args::NumericalIntegrationProblem,r)

  integrate(exp:EF,ra:SOCF,epsabs:F,epsrel:F):Result ==
    integrate(exp,ra,epsabs,epsrel,routines()$RT)

  integrate(exp:EF,ra:SOCF,err:F):Result ==
    positive?(err)$F => integrate(exp,ra,0$F,err)
    integrate(exp,ra,1.0E-5,err)

  integrate(exp:EF,ra:SOCF):Result == integrate(exp,ra,0$F,1.0E-5)

  integrate(exp:EF,sb:SBOCF, st:ST) ==
    st = "numerical" => integrate(exp,segment sb)
    "failed"

  integrate(exp:EF,sb:SBOCF, s:S) ==
    s = (numerical::Symbol) => integrate(exp,segment sb)
    "failed"

  integrate(exp:EF,ra:LSOCF,epsabs:F,epsrel:F,r:RT):Result ==
    vars := variables(exp)$EF
    empty?(vars)$LS => integrateConstantList(exp,ra)
    args:MDNIA := [ef2edf exp,convert ra,f2df epsabs,f2df epsrel]
    integrateArgs(args::NumericalIntegrationProblem,r)

  integrate(exp:EF,ra:LSOCF,epsabs:F,epsrel:F):Result ==
    integrate(exp,ra,epsabs,epsrel,routines()$RT)

  integrate(exp:EF,ra:LSOCF,epsrel:F):Result ==
    zero? epsrel => integrate(exp,ra,1.0e-6,epsrel)
    integrate(exp,ra,0$F,epsrel)

  integrate(exp:EF,ra:LSOCF):Result == integrate(exp,ra,1.0e-4)