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

)abbrev package D02AGNT d02AgentsPackage
++ Author: Brian Dupee
++ Date Created: May 1994
++ Date Last Updated: January 1997
++ Description:

d02AgentsPackage() : SIG == CODE 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)

  SIG ==>  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].\indent{20}
      ++ 400 ``operation units'' -> 0.75 
      ++ 200 ``operation units'' -> 0.5 
      ++ 83 ``operation units'' -> 0.25 \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].

  CODE ==> 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]
      -- array === matrix
      realPartsMatrix:MDF := retract(realParts)$AnyFunctions1(MDF) 
      imagParts:Union(Any,"failed") := search(ri::Symbol,res)$Result
      imagParts case "failed" => [0$F,0$F]
      -- array === matrix
      imagPartsMatrix:MDF := retract(imagParts)$AnyFunctions1(MDF) 
      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