path: root/src/algebra/d02agents.spad.pamphlet

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra d02agents.spad}
\author{Brian Dupee}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{domain ODEIFTBL ODEIntensityFunctionsTable}
<<domain ODEIFTBL ODEIntensityFunctionsTable>>=
)abbrev domain ODEIFTBL ODEIntensityFunctionsTable
++ Author: Brian Dupee
++ Date Created: May 1994
++ Date Last Updated: January 1996
++ Basic Operations: showTheIFTable, insert!
++ Description:
++ \axiom{ODEIntensityFunctionsTable()} provides a dynamic table and a set of
++ functions to store details found out about sets of ODE's.

ODEIntensityFunctionsTable(): E == I where
  LEDF	==> List Expression DoubleFloat
  LEEDF	==> List Equation Expression DoubleFloat
  EEDF	==> Equation Expression DoubleFloat
  VEDF	==> Vector Expression DoubleFloat
  MEDF	==> Matrix Expression DoubleFloat
  MDF	==> Matrix DoubleFloat
  EDF	==> Expression DoubleFloat
  DF	==> DoubleFloat
  F	==> Float
  INT	==> Integer
  CDF	==> Complex DoubleFloat
  LDF	==> List DoubleFloat
  LF	==> List Float
  S	==> Symbol
  LS	==> List Symbol
  MFI	==> Matrix Fraction Integer
  LFI	==> List Fraction Integer
  FI	==> Fraction Integer
  ODEA	==> Record(xinit:DF,xend:DF,fn:VEDF,yinit:LDF,intvals:LDF,g:EDF,abserr:DF,relerr:DF)
  ON	==> Record(additions:INT,multiplications:INT,exponentiations:INT,functionCalls:INT)
  RVE 	==> Record(val:EDF,exponent:INT)
  RSS	==> Record(stiffnessFactor:F,stabilityFactor:F)
  ATT	==> Record(stiffness:F,stability:F,expense:F,accuracy:F,intermediateResults:F)
  ROA	==> Record(key:ODEA,entry:ATT)

  E ==> with
    showTheIFTable:() -> $
      ++ showTheIFTable() returns the current table of intensity functions.
    clearTheIFTable : () -> Void
      ++ clearTheIFTable() clears the current table of intensity functions.
    keys : $ -> List(ODEA)
      ++ keys(tab) returns the list of keys of f
    iFTable: List Record(key:ODEA,entry:ATT) -> $
      ++ iFTable(l) creates an intensity-functions table from the elements 
      ++ of l.
    insert!:Record(key:ODEA,entry:ATT) -> $
      ++ insert!(r) inserts an entry r into theIFTable
    showIntensityFunctions:ODEA -> Union(ATT,"failed")
      ++ showIntensityFunctions(k) returns the entries in the 
      ++ table of intensity functions k.

  I ==> add
    Rep := Table(ODEA,ATT)
    import Rep

    theIFTable:$ := empty()$Rep

    showTheIFTable():$ ==
      theIFTable

    clearTheIFTable():Void ==
      theIFTable := empty()$Rep
      void()$Void

    iFTable(l:List Record(key:ODEA,entry:ATT)):$ ==
      theIFTable := table(l)$Rep

    insert!(r:Record(key:ODEA,entry:ATT)):$ ==
      insert!(r,theIFTable)$Rep

    keys(t:$):List ODEA ==
      keys(t)$Rep

    showIntensityFunctions(k:ODEA):Union(ATT,"failed") ==
      search(k,theIFTable)$Rep

@
\section{package D02AGNT d02AgentsPackage}
<<package D02AGNT d02AgentsPackage>>=
)abbrev package D02AGNT d02AgentsPackage
++ Author: Brian Dupee
++ Date Created: May 1994
++ Date Last Updated: January 1997
++ Basic Operations: stiffnessFactor, jacobian
++ Description:
++ \axiom{d02AgentsPackage} contains a set of computational agents 
++ for use with Ordinary Differential Equation solvers.
d02AgentsPackage(): E == I where
  LEDF	==> List Expression DoubleFloat
  LEEDF	==> List Equation Expression DoubleFloat
  EEDF	==> Equation Expression DoubleFloat
  VEDF	==> Vector Expression DoubleFloat
  MEDF	==> Matrix Expression DoubleFloat
  MDF	==> Matrix DoubleFloat
  EDF	==> Expression DoubleFloat
  DF	==> DoubleFloat
  F	==> Float
  INT	==> Integer
  CDF	==> Complex DoubleFloat
  LDF	==> List DoubleFloat
  LF	==> List Float
  S	==> Symbol
  LS	==> List Symbol
  MFI	==> Matrix Fraction Integer
  LFI	==> List Fraction Integer
  FI	==> Fraction Integer
  ODEA	==> Record(xinit:DF,xend:DF,fn:VEDF,yinit:LDF,intvals:LDF,g:EDF,abserr:DF,relerr:DF)
  ON	==> Record(additions:INT,multiplications:INT,exponentiations:INT,functionCalls:INT)
  RVE 	==> Record(val:EDF,exponent:INT)
  RSS	==> Record(stiffnessFactor:F,stabilityFactor:F)
  ATT	==> Record(stiffness:F,stability:F,expense:F,accuracy:F,intermediateResults:F)
  ROA	==> Record(key:ODEA,entry:ATT)

  E ==>  with
    combineFeatureCompatibility: (F,F) -> F
      ++ combineFeatureCompatibility(C1,C2) is for interacting attributes
    combineFeatureCompatibility: (F,LF) -> F
      ++ combineFeatureCompatibility(C1,L) is for interacting attributes
    sparsityIF: MEDF -> F
      ++ sparsityIF(m) calculates the sparsity of a jacobian matrix
    jacobian: (VEDF,LS) -> MEDF
      ++ jacobian(v,w) is a local function to make a jacobian matrix
    eval: (MEDF,LS,VEDF) -> MEDF
      ++ eval(mat,symbols,values) evaluates a multivariable matrix at given values
      ++ for each of a list of variables
    stiffnessAndStabilityFactor: MEDF -> RSS
      ++ stiffnessAndStabilityFactor(me) calculates the stability and
      ++ stiffness factor of a system of first-order differential equations
      ++ (by evaluating the maximum difference in the real parts of the
      ++ negative eigenvalues of the jacobian of the system for which O(10)
      ++ equates to mildly stiff wheras stiffness ratios of O(10^6) are not
      ++ uncommon) and whether the system is likely to show any oscillations
      ++ (identified by the closeness to the imaginary axis of the complex
      ++ eigenvalues of the jacobian).
    stiffnessAndStabilityOfODEIF:ODEA -> RSS
      ++ stiffnessAndStabilityOfODEIF(ode) calculates the intensity values
      ++ of stiffness of a system of first-order differential equations
      ++ (by evaluating the maximum difference in the real parts of the
      ++ negative eigenvalues of the jacobian of the system for which O(10)
      ++ equates to mildly stiff wheras stiffness ratios of O(10^6) are not
      ++ uncommon) and whether the system is likely to show any oscillations
      ++ (identified by the closeness to the imaginary axis of the complex
      ++ eigenvalues of the jacobian).
      ++
      ++ It returns two values in the range [0,1].
    systemSizeIF:ODEA -> F
      ++ systemSizeIF(ode) returns the intensity value of the size of
      ++ the system of ODEs.  20 equations corresponds to the neutral
      ++ value.  It returns a value in the range [0,1].
    expenseOfEvaluationIF:ODEA -> F
      ++ expenseOfEvaluationIF(o) returns the intensity value of the 
      ++ cost of evaluating the input ODE.  This is in terms of the number
      ++ of ``operational units''.  It returns a value in the range
      ++ [0,1].\newline\indent{20}
      ++ 400 ``operation units'' -> 0.75 \newline
      ++ 200 ``operation units'' -> 0.5 \newline
      ++ 83 ``operation units'' -> 0.25 \newline\indent{15}
      ++ exponentiation = 4 units , function calls = 10 units.
    accuracyIF:ODEA -> F
      ++ accuracyIF(o) returns the intensity value of the accuracy
      ++ requirements of the input ODE.  A request of accuracy of 10^-6 
      ++ corresponds to the neutral intensity.  It returns a value
      ++ in the range [0,1].
    intermediateResultsIF:ODEA -> F
      ++ intermediateResultsIF(o) returns a value corresponding to the 
      ++ required number of intermediate results required and, therefore, 
      ++ an indication of how much this would affect the step-length of the 
      ++ calculation.  It returns a value in the range [0,1].

  I ==> add

    import ExpertSystemToolsPackage

    accuracyFactor:ODEA -> F
    expenseOfEvaluation:ODEA -> F
    eval1:(LEDF,LEEDF) -> LEDF
    stiffnessAndStabilityOfODE:ODEA -> RSS
    intermediateResultsFactor:ODEA -> F
    leastStabilityAngle:(LDF,LDF) -> F

    intermediateResultsFactor(ode:ODEA):F ==
      resultsRequirement := #(ode.intvals)
      (1.0-exp(-(resultsRequirement::F)/50.0)$F)

    intermediateResultsIF(o:ODEA):F ==
      ode := copy o
      (t := showIntensityFunctions(ode)$ODEIntensityFunctionsTable) case ATT =>
        s := coerce(t)@ATT
        negative?(s.intermediateResults)$F => 
          s.intermediateResults := intermediateResultsFactor(ode)
          r:ROA := [ode,s]
          insert!(r)$ODEIntensityFunctionsTable
          s.intermediateResults
        s.intermediateResults
      a:ATT := [-1.0,-1.0,-1.0,-1.0,e:=intermediateResultsFactor(ode)]
      r:ROA := [ode,a]
      insert!(r)$ODEIntensityFunctionsTable
      e

    accuracyFactor(ode:ODEA):F ==
      accuracyRequirements := convert(ode.abserr)@F
      if zero?(accuracyRequirements) then
        accuracyRequirements := convert(ode.relerr)@F
      val := inv(accuracyRequirements)$F
      n := log10(val)$F
      (1.0-exp(-(n/(2.0))**2/(15.0))$F)

    accuracyIF(o:ODEA):F ==
      ode := copy o
      (t := showIntensityFunctions(ode)$ODEIntensityFunctionsTable) case ATT =>
        s := coerce(t)@ATT
        negative?(s.accuracy)$F => 
          s.accuracy := accuracyFactor(ode)
          r:ROA := [ode,s]
          insert!(r)$ODEIntensityFunctionsTable
          s.accuracy
        s.accuracy
      a:ATT := [-1.0,-1.0,-1.0,e:=accuracyFactor(ode),-1.0]
      r:ROA := [ode,a]
      insert!(r)$ODEIntensityFunctionsTable
      e

    systemSizeIF(ode:ODEA):F ==
      n := #(ode.fn)
      (1.0-exp((-n::F/75.0))$F)

    expenseOfEvaluation(o:ODEA):F ==
      -- expense of evaluation of an ODE -- <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.
      ode := copy o.fn
      expenseOfEvaluation(ode)

    expenseOfEvaluationIF(o:ODEA):F ==
      ode := copy o
      (t := showIntensityFunctions(ode)$ODEIntensityFunctionsTable) case ATT =>
        s := coerce(t)@ATT
        negative?(s.expense)$F => 
          s.expense := expenseOfEvaluation(ode)
          r:ROA := [ode,s]
          insert!(r)$ODEIntensityFunctionsTable
          s.expense
        s.expense
      a:ATT := [-1.0,-1.0,e:=expenseOfEvaluation(ode),-1.0,-1.0]
      r:ROA := [ode,a]
      insert!(r)$ODEIntensityFunctionsTable
      e

    leastStabilityAngle(realPartsList:LDF,imagPartsList:LDF):F ==
      complexList := [complex(u,v)$CDF for u in realPartsList for v in imagPartsList]
      argumentList := [abs((abs(argument(u)$CDF)$DF)-(pi()$DF)/2)$DF for u in complexList]
      sortedArgumentList := sort(argumentList)$LDF
      list := [u for u in sortedArgumentList | not zero?(u) ]
      empty?(list)$LDF => 0$F
      convert(first(list)$LDF)@F

    stiffnessAndStabilityFactor(me:MEDF):RSS ==

      -- search first for real eigenvalues of the jacobian (symbolically)
      -- if the system isn't too big
      r:INT := ncols(me)$MEDF  
      b:Boolean := ((# me) < 150)
      if b then
        mc:MFI := map(edf2fi,me)$ExpertSystemToolsPackage2(EDF,FI)
        e:LFI := realEigenvalues(mc,1/100)$NumericRealEigenPackage(FI)
        b := ((# e) >= r-1)@Boolean       
      b =>
        -- if all the eigenvalues are real, find negative ones
        e := sort(neglist(e)$ExpertSystemToolsPackage1(FI))
        -- if there are two or more, calculate stiffness ratio
        ((n:=#e)>1)@Boolean => [coerce(e.1/e.n)@F,0$F] 
        -- otherwise stiffness not present
        [0$F,0$F]

      md:MDF := map(edf2df,me)$ExpertSystemToolsPackage2(EDF,DF)

      -- otherwise calculate numerically the complex eigenvalues
      -- using NAG routine f02aff.

      res:Result := f02aff(r,r,md,-1)$NagEigenPackage
      realParts:Union(Any,"failed") := search(rr::Symbol,res)$Result
      realParts case "failed" => [0$F,0$F]
      realPartsMatrix:MDF := retract(realParts)$AnyFunctions1(MDF) -- array === matrix
      imagParts:Union(Any,"failed") := search(ri::Symbol,res)$Result
      imagParts case "failed" => [0$F,0$F]
      imagPartsMatrix:MDF := retract(imagParts)$AnyFunctions1(MDF) -- array === matrix
      imagPartsList:LDF := members(imagPartsMatrix)$MDF
      realPartsList:LDF := members(realPartsMatrix)$MDF
      stabilityAngle := leastStabilityAngle(realPartsList,imagPartsList)
      negRealPartsList := sort(neglist(realPartsList)$ExpertSystemToolsPackage1(DF))
      empty?(negRealPartsList)$LDF => [0$F,stabilityAngle]
      ((n:=#negRealPartsList)>1)@Boolean => 
        out := convert(negRealPartsList.1/negRealPartsList.n)@F
        [out,stabilityAngle]    -- calculate stiffness ratio
      [-convert(negRealPartsList.1)@F,stabilityAngle]
      
    eval1(l:LEDF,e:LEEDF):LEDF ==
      [eval(u,e)$EDF for u in l]

    eval(mat:MEDF,symbols:LS,values:VEDF):MEDF ==
      l := listOfLists(mat)
      ledf := entries(values)$VEDF
      e := [equation(u::EDF,v)$EEDF for u in symbols for v in ledf]
      l := [eval1(w,e) for w in l]
      matrix l

    combineFeatureCompatibility(C1:F,C2:F):F ==

      --                        C1 C2
      --    s(C1,C2) = -----------------------
      --               C1 C2 + (1 - C1)(1 - C2)

      C1*C2/((C1*C2)+(1$F-C1)*(1$F-C2))

    combineFeatureCompatibility(C1:F,L:LF):F ==

      empty?(L)$LF => C1
      C2 := combineFeatureCompatibility(C1,first(L)$LF)
      combineFeatureCompatibility(C2,rest(L)$LF)

    jacobian(v:VEDF,w:LS):Matrix EDF ==
      jacobian(v,w)$MultiVariableCalculusFunctions(S,EDF,VEDF,LS)

    sparsityIF(m:Matrix EDF):F ==
      l:LEDF :=parts m
      z:LEDF := [u for u in l | zero?(u)$EDF]
      ((#z)::F/(#l)::F)

    sum(a:EDF,b:EDF):EDF == a+b

    stiffnessAndStabilityOfODE(ode:ODEA):RSS ==
      odefns := copy ode.fn
      ls:LS := [subscript(Y,[coerce(n)])$Symbol for n in 1..# odefns]
      yvals := copy ode.yinit
      for i in 1..#yvals repeat
        zero?(yvals.i) => yvals.i := 0.1::DF
      yexpr := [coerce(v)@EDF for v in yvals]
      yv:VEDF := vector(yexpr)
      j1:MEDF := jacobian(odefns,ls)
      ej1:MEDF := eval(j1,ls,yv)
      ej1 := eval(ej1,variables(reduce(sum,members(ej1)$MEDF)),vector([(ode.xinit)::EDF]))
      ssf := stiffnessAndStabilityFactor(ej1)
      stability := 1.0-sqrt((ssf.stabilityFactor)*(2.0)/(pi()$F))
      stiffness := (1.0)-exp(-(ssf.stiffnessFactor)/(500.0))
      [stiffness,stability]

    stiffnessAndStabilityOfODEIF(ode:ODEA):RSS ==
      odefn := copy ode
      (t := showIntensityFunctions(odefn)$ODEIntensityFunctionsTable) case ATT =>
        s:ATT := coerce(t)@ATT
        negative?(s.stiffness)$F => 
          ssf:RSS := stiffnessAndStabilityOfODE(odefn)
          s := [ssf.stiffnessFactor,ssf.stabilityFactor,s.expense,
                  s.accuracy,s.intermediateResults]
          r:ROA := [odefn,s]
          insert!(r)$ODEIntensityFunctionsTable
          ssf
        [s.stiffness,s.stability]
      ssf:RSS := stiffnessAndStabilityOfODE(odefn)
      s:ATT := [ssf.stiffnessFactor,ssf.stabilityFactor,-1.0,-1.0,-1.0]
      r:ROA := [odefn,s]
      insert!(r)$ODEIntensityFunctionsTable
      ssf

@
\section{License}
<<license>>=
--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.
@
<<*>>=
<<license>>

<<domain ODEIFTBL ODEIntensityFunctionsTable>>
<<package D02AGNT d02AgentsPackage>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}