Adds a tool interface for Mathematica's AsymptoticDSolveValue command.

keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/btactics/ToolProvider.scala

@@ -68,6 +68,8 @@ object ToolProvider extends ToolProvider with Logging {
 
   def sosSolveTool(): Option[SOSsolveTool] = f.sosSolveTool()
 
+  def differentialSeriesApproxmationnTool() = f.differentialSeriesApproxmationnTool()
+
   def init(): Boolean = f.init()
 
   def shutdown(): Unit = f.shutdown()
@@ -147,6 +149,9 @@ trait ToolProvider {
   /** Returns a SOSsolve tool. */
   def sosSolveTool(): Option[SOSsolveTool]
 
+  /** Returns a series expansion tool. */
+  def differentialSeriesApproxmationnTool(): Option[DifferentialSolutionSeriesApproximationTool]
+
   /** Initializes the tools. */
   def init(): Boolean
 
@@ -174,6 +179,7 @@ class PreferredToolProvider[T <: Tool](val toolPreferences: List[T]) extends Too
   private[this] lazy val solver: Option[Tool with EquationSolverTool] = toolPreferences.find(_.isInstanceOf[EquationSolverTool]).map(_.asInstanceOf[Tool with EquationSolverTool])
   private[this] lazy val algebra: Option[Tool with AlgebraTool] = toolPreferences.find(_.isInstanceOf[AlgebraTool]).map(_.asInstanceOf[Tool with AlgebraTool])
   private[this] lazy val sossolve: Option[Tool with SOSsolveTool] = toolPreferences.find(_.isInstanceOf[SOSsolveTool]).map(_.asInstanceOf[Tool with SOSsolveTool])
+  private[this] lazy val diffSeriesApproximation: Option[Tool with DifferentialSolutionSeriesApproximationTool] = toolPreferences.find(_.isInstanceOf[DifferentialSolutionSeriesApproximationTool]).map(_.asInstanceOf[Tool with DifferentialSolutionSeriesApproximationTool])
 
   override def tools(): List[Tool] = toolPreferences
   override def defaultTool(): Option[Tool] = toolPreferences.headOption
@@ -193,6 +199,7 @@ class PreferredToolProvider[T <: Tool](val toolPreferences: List[T]) extends Too
   override def solverTool(): Option[EquationSolverTool] = ensureInitialized(solver)
   override def algebraTool(): Option[AlgebraTool] = ensureInitialized(algebra)
   override def sosSolveTool(): Option[SOSsolveTool] = ensureInitialized(sossolve)
+  override def differentialSeriesApproxmationnTool(): Option[DifferentialSolutionSeriesApproximationTool] = ensureInitialized(diffSeriesApproximation)
   override def init(): Boolean = false /* override to initialize tools in more specialized providers */
   override def shutdown(): Unit = toolPreferences.foreach(_.shutdown())
   override def isInitialized: Boolean = toolPreferences.forall(_.isInitialized)
@@ -221,6 +228,7 @@ class NoneToolProvider extends ToolProvider {
   override def solverTool(): Option[EquationSolverTool] = None
   override def algebraTool(): Option[AlgebraTool] = None
   override def sosSolveTool(): Option[SOSsolveTool] = None
+  override def differentialSeriesApproxmationnTool() = None
   override def init(): Boolean = true
   override def shutdown(): Unit = {}
   override def isInitialized: Boolean = true

keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/DifferentialSolutionSeriesApproximationTool.scala

@@ -0,0 +1,14 @@
+package edu.cmu.cs.ls.keymaerax.tools.ext
+
+import edu.cmu.cs.ls.keymaerax.core.{NamedSymbol, Number, ODESystem, Term, Variable}
+import edu.cmu.cs.ls.keymaerax.tools.ToolInterface
+
+trait DifferentialSolutionSeriesApproximationTool extends ToolInterface {
+  /**
+    *
+    * @param odes The ODEs
+    * @param ctx Context for any parameters in the ODEs.
+    * @return upper and lower bound series approximations of each primed variable in the ODEs.
+    */
+  def seriesApproximation(odes: ODESystem, ctx: Map[Term, Term]): Option[Map[Variable, Term]]
+}

keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/ExtMathematicaOpSpec.scala

@@ -105,6 +105,8 @@ object ExtMathematicaOpSpec {
 
   def d: BinaryMathOpSpec = BinaryMathOpSpec(symbol("D"))
 
+  def dsolveAsymptoticApproximation: NaryMathOpSpec = NaryMathOpSpec(symbol("AsymptoticDSolveValue"))
+
   //</editor-fold>
 
 }

keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/Mathematica.scala

@@ -21,6 +21,7 @@ import edu.cmu.cs.ls.keymaerax.tools.ext.SOSsolveTool.Result
 import edu.cmu.cs.ls.keymaerax.tools.qe.MathematicaConversion.{KExpr, MExpr}
 import edu.cmu.cs.ls.keymaerax.tools.qe.MathematicaOpSpec._
 import edu.cmu.cs.ls.keymaerax.tools.qe.MathematicaToKeYmaera
+import edu.cmu.cs.ls.keymaerax.tools.ext.MathematicaDifferentialSolutionSeriesApproximationTool
 
 import scala.annotation.tailrec
 import scala.collection.immutable.{Map, Seq}
@@ -37,7 +38,7 @@ import scala.collection.mutable.ListBuffer
 class Mathematica(private[tools] val link: MathematicaLink, override val name: String) extends Tool
     with QETacticTool with InvGenTool with ODESolverTool with CounterExampleTool
     with SimulationTool with DerivativeTool with EquationSolverTool with SimplificationTool with AlgebraTool
-    with PDESolverTool with SOSsolveTool with ToolOperationManagement {
+    with PDESolverTool with SOSsolveTool with ToolOperationManagement with DifferentialSolutionSeriesApproximationTool {
 
   /** Indicates whether the tool is initialized. */
   private var initialized = false
@@ -52,6 +53,7 @@ class Mathematica(private[tools] val link: MathematicaLink, override val name: S
   private val mSolve = new MathematicaEquationSolverTool(link)
   private val mAlgebra = new MathematicaAlgebraTool(link)
   private val mSimplify = new MathematicaSimplificationTool(link)
+  private val mExpand = new MathematicaDifferentialSolutionSeriesApproximationTool(link)
 
   private val qeInitialTimeout = Integer.parseInt(Configuration(Configuration.Keys.QE_TIMEOUT_INITIAL))
   private val qeCexTimeout = Integer.parseInt(Configuration(Configuration.Keys.QE_TIMEOUT_CEX))
@@ -70,6 +72,7 @@ class Mathematica(private[tools] val link: MathematicaLink, override val name: S
     mSolve.memoryLimit = memoryLimit
     mAlgebra.memoryLimit = memoryLimit
     mSimplify.memoryLimit = memoryLimit
+    mExpand.memoryLimit = memoryLimit
 
     // initialze tool thread pools
     mQE.init()
@@ -82,6 +85,8 @@ class Mathematica(private[tools] val link: MathematicaLink, override val name: S
     mAlgebra.init()
     mSimplify.init()
     mSOSsolve.init()
+    mExpand.init()
+
 
     initialized = link match {
       case l: JLinkMathematicaLink =>
@@ -113,6 +118,7 @@ class Mathematica(private[tools] val link: MathematicaLink, override val name: S
     mSolve.shutdown()
     mAlgebra.shutdown()
     mSimplify.shutdown()
+    mExpand.shutdown()
     //@note last, because we want to shut down all executors (tool threads) before shutting down the JLink interface
     link match {
       case l: JLinkMathematicaLink => l.shutdown()
@@ -334,4 +340,8 @@ class Mathematica(private[tools] val link: MathematicaLink, override val name: S
 
   /** @inheritdoc */
   override def getAvailableWorkers: Int = 1
+
+  /** @inheritdoc */
+  override def seriesApproximation(odes: ODESystem, ctx: Map[Term, Term]): Option[Map[Variable, Term]] =
+    mExpand.seriesApproximation(odes, ctx)
 }

keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaLink.scala

@@ -313,6 +313,14 @@ class JLinkMathematicaLink(val engineName: String) extends MathematicaLink with
     */
   override def run[T](cmd: () => T, executor: ToolExecutor): T = {
     if (ml == null) throw new IllegalStateException("No MathKernel set")
+    if (executor == null) throw new IllegalStateException(
+      """
+        |No Executor was set.
+        |
+        |Likely explanation: a tool was used before initialization. Remember to call .init() on all tools or subtools
+        |before use. E.g., in edu.cmu.cs.ls.keymaerax.tools.ext.Mathematica, a call to .init() should be added for
+        |every new subtool.
+      """.stripMargin)
     val taskId = executor.schedule(_ => { ml.synchronized { cmd() } })
 
     executor.wait(taskId) match {

keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaTools.scala

@@ -10,12 +10,13 @@ import edu.cmu.cs.ls.keymaerax.btactics.InvariantGenerator
 import edu.cmu.cs.ls.keymaerax.btactics.helpers.DifferentialHelper
 import edu.cmu.cs.ls.keymaerax.core.{Variable, _}
 import edu.cmu.cs.ls.keymaerax.infrastruct.ExpressionTraversal.{ExpressionTraversalFunction, StopTraversal}
-import edu.cmu.cs.ls.keymaerax.infrastruct.{ExpressionTraversal, FormulaTools, PosInExpr}
+import edu.cmu.cs.ls.keymaerax.infrastruct.{ExpressionTraversal, FormulaTools, PosInExpr, SubstitutionHelper}
 import edu.cmu.cs.ls.keymaerax.tools.qe.MathematicaConversion.{KExpr, _}
 import edu.cmu.cs.ls.keymaerax.tools.ext.SimulationTool.{SimRun, SimState, Simulation}
-import edu.cmu.cs.ls.keymaerax.tools.qe.{BinaryMathOpSpec, ExprFactory, K2MConverter, KeYmaeraToMathematica, M2KConverter, MathematicaNameConversion, MathematicaOpSpec, MathematicaToKeYmaera, NaryMathOpSpec, UnaryMathOpSpec}
+import edu.cmu.cs.ls.keymaerax.tools.qe.{DifferentialProgramToMathematica, BinaryMathOpSpec, ExprFactory, K2MConverter, KeYmaeraToMathematica, M2KConverter, MathematicaNameConversion, MathematicaOpSpec, MathematicaToKeYmaera, NaryMathOpSpec, UnaryMathOpSpec}
 import edu.cmu.cs.ls.keymaerax.tools.qe.MathematicaOpSpec._
 import edu.cmu.cs.ls.keymaerax.tools._
+import edu.cmu.cs.ls.keymaerax.infrastruct.Augmentors._
 
 import scala.collection.immutable
 import scala.math.BigDecimal
@@ -780,6 +781,76 @@ class MathematicaEquationSolverTool(override val link: MathematicaLink) extends
   }
 }
 
+/**
+  * Uses Mathematica's AsymptoticDSolveValue function to construct power series expansions of the solution to a system of ODEs.
+  *
+  * @author Nathan Fulton
+  */
+class MathematicaDifferentialSolutionSeriesApproximationTool(override val link: MathematicaLink) extends BaseKeYmaeraMathematicaBridge[KExpr](link, new UncheckedBaseK2MConverter, PegasusM2KConverter) with DifferentialSolutionSeriesApproximationTool {
+
+  private def pairToList(t: Term): List[Term] = t match {
+    case Pair(left, right) => {
+      assert(!left.isInstanceOf[Pair], "Expected canonical ordering so that tail recursion is possible.")
+      left :: pairToList(right)
+    }
+    case _ => t :: Nil
+  }
+
+  private def convertResult(odes: ODESystem, result: Expression) = {
+    result match {
+      case t: Term => {
+        // Convert nested pairs into a list of terms.
+        val approximations = pairToList(t)
+
+        // Replace all constant symbols in the terms with old(x_i) where x_i is the i^th primed variable in the ODEs.
+        val primedVariables = DifferentialHelper.getPrimedVariables(odes)
+        val approximationsWithInitialConditions = approximations.map(approximation => {
+          Range(1, primedVariables.length + 1).foldLeft(approximation)((approximation, i) => {
+            val constantSymbol = FuncOf(Function("C", None, Real, Real, true), Number(i)) // C[i]
+            val x_i = primedVariables(i - 1) //the i^th primed variable in the system.
+            val old_x_i = FuncOf(Function("old", None, Real, Real), x_i) // old(x_i)
+            approximation.replaceAll(constantSymbol, old_x_i)
+          })
+        })
+
+        Some(
+          DifferentialHelper.getPrimedVariables(odes)
+            .zip(approximationsWithInitialConditions)
+            .toMap
+        )
+      }
+      case _ => throw new ConversionException(s"Expected AsymptoticDSolveValue to give a term or a list of terms, but found ${result}")
+    }
+  }
+
+  /** @inheritdoc */
+  override def seriesApproximation(odes: ODESystem, ctx: Map[Term, Term]): Option[Map[Variable, Term]] = {
+    //See https://reference.wolfram.com/language/ref/AsymptoticDSolveValue.html for the command we're going to use.
+
+    // Apply the context to the ODEs.
+    val odesInCtx = ctx.foldLeft(odes)((currExpr, x) => {
+      SubstitutionHelper.replaceFree(currExpr)(x._1, x._2).asInstanceOf[ODESystem]
+    })
+
+    // Identify a unique time variable.
+    val TIME_VAR = Variable("t") //@todo do something here that's less stupid.
+
+    // Convert the ODEs into Mathematica expressions.
+    // @todo move this out into a utility class independent of the k2m converter interface.
+    val (mODEs, ivs, _, _) = DifferentialProgramToMathematica(k2m).apply(odesInCtx.ode, TIME_VAR, Map())
+
+    //construct the Mathematica expression that will execute AsymptoticDSolveValue.
+    val input = ExtMathematicaOpSpec.dsolveAsymptoticApproximation(
+      MathematicaOpSpec.list(mODEs:_*),
+      MathematicaOpSpec.list(DifferentialHelper.getPrimedVariables(odesInCtx).map(k2m):_*),
+      MathematicaOpSpec.list(k2m(TIME_VAR), k2m(Number(0)), k2m(Number(10))) //@todo I have no clue what these bounds should be. Make these arguments at least.
+    )
+
+    val (_, result) = run(input)
+    convertResult(odes, result)
+  }
+}
+
 /**
   * A link to Mathematica using the JLink interface.
   *

keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/DifferentialProgramToMathematica.scala

@@ -0,0 +1,81 @@
+package edu.cmu.cs.ls.keymaerax.tools.qe
+
+import edu.cmu.cs.ls.keymaerax.core.{AtomicODE, BaseVariable, DifferentialProduct, DifferentialProgram, DifferentialSymbol, Equal, FuncOf, Function, Number, ODESystem, Program, Term, Variable}
+import edu.cmu.cs.ls.keymaerax.infrastruct.ExpressionTraversal.{ExpressionTraversalFunction, StopTraversal}
+import edu.cmu.cs.ls.keymaerax.infrastruct.{ExpressionTraversal, PosInExpr}
+import edu.cmu.cs.ls.keymaerax.tools.ConversionException
+import edu.cmu.cs.ls.keymaerax.tools.ext.ExtMathematicaOpSpec
+import edu.cmu.cs.ls.keymaerax.tools.qe.MathematicaConversion.MExpr
+
+import scala.collection.immutable.{List, Map}
+import scala.math.BigDecimal
+
+case class DifferentialProgramToMathematica(k2m: K2MConverter[MathematicaConversion.KExpr]) {
+  /**
+    * Converts a DifferentialProgram to a Mathematica expression.
+    * This code is largely copy/pasted from MathematicaTools.MathematicalODESolvertool, so that we can KeYmaera
+    * translations of ODEs in other contexts.
+    * @return a 4-tuple containing:
+    *         a list of MExprs specifying the ODEs,
+    *         a list of MExprs specifying symbolic initial conditions,
+    *         a list of MExprs specifying functions (IDK what this is actually -- check the dsolve code and Wolfram documentaiton to figure out)
+    *         the variable used as the time variable (diffArg, as passed in initially)
+    * @author Nathan Fulton
+    */
+  def apply(dp: DifferentialProgram, diffArg: Variable, iv: Map[Variable, Variable]): (List[MExpr], List[MExpr], List[MExpr], Variable) = {
+    /**
+      * @note copied from MathematicaTools.MathematicalODESolvertool.
+      */
+    def toDiffSys(diffSys: DifferentialProgram, diffArg: Variable): List[(Variable, Term)] = {
+      var result = List[(Variable, Term)]()
+      ExpressionTraversal.traverse(new ExpressionTraversalFunction {
+        override def preP(p: PosInExpr, e: Program): Either[Option[StopTraversal], Program] = e match {
+          case AtomicODE(DifferentialSymbol(x), theta) if x != diffArg => result = result :+ (x, theta); Left(None)
+          case AtomicODE(DifferentialSymbol(x), _) if x == diffArg => Left(None)
+          case ODESystem(_, _) => Left(None)
+          case DifferentialProduct(_, _) => Left(None)
+        }
+      }, diffSys)
+      result
+    }
+
+    /** @note coped from dsolve tool. */
+    def functionalizeVars(t: Term, arg: Term, vars: Variable*) = ExpressionTraversal.traverse(
+      new ExpressionTraversalFunction {
+        override def postT(p: PosInExpr, e: Term): Either[Option[StopTraversal], Term] = e match {
+          case v@BaseVariable(name, idx, sort) if vars.isEmpty || vars.contains(v) =>
+            Right(FuncOf(Function(name, idx, arg.sort, sort), arg))
+          case _ => Left(None)
+        }
+      }, t) match {
+      case Some(resultTerm) => resultTerm
+      case None => throw ConversionException("Unable to functionalize " + t)
+    }
+
+    val diffSys = toDiffSys(dp, diffArg)
+
+    val primedVars = diffSys.map(_._1)
+    val functionalizedTerms = diffSys.map{ case (x, theta) => ( x, functionalizeVars(theta, diffArg, primedVars:_*)) }
+    val mathTerms = functionalizedTerms.map({case (x, theta) =>
+      (ExtMathematicaOpSpec.dx(ExtMathematicaOpSpec.primed(k2m(x)))(k2m(diffArg)), k2m(theta))})
+    val convertedDiffSys = mathTerms.map({case (x, theta) => MathematicaOpSpec.equal(x, theta)})
+
+    val functions = diffSys.map(t => k2m(functionalizeVars(t._1, diffArg)))
+
+    //@todo allows for partial initial conditions, but this will probably cause issues if IVs are used.
+    val initialValues = diffSys
+      .map(t =>
+        iv.get(t._1) match {
+          case Some(definedInitialValue) =>
+            Some(
+              Equal(functionalizeVars(t._1, Number(BigDecimal(0)), primedVars:_*), definedInitialValue)
+            )
+          case None => None //@todo perhaps add an error here?
+        }
+      )
+      .filterNot(x => x.isEmpty)
+      .map(x => k2m(x.get))
+
+    (convertedDiffSys, initialValues, functions, diffArg)
+  }
+}

keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaOpSpec.scala

@@ -129,6 +129,8 @@ object MathematicaOpSpec {
 
   def power: BinaryMathOpSpec = BinaryMathOpSpec(symbol("Power"))
 
+  def C: UnaryMathOpSpec = UnaryMathOpSpec(symbol("C")) //@todo document the meaning of this symbol.
+
   // implicit function application name[args]
   def func: NameMathOpSpec = NameMathOpSpec(
     (name: NamedSymbol, args: Array[Expr]) => {

keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaToKeYmaera.scala

@@ -41,6 +41,12 @@ class MathematicaToKeYmaera extends M2KConverter[KExpr] {
     //@note self-created MExpr with head RATIONAL are not rationalQ (type identifiers do not match)
     else if (MathematicaOpSpec.rational.applies(e)) convertBinary(e, Divide.apply)
 
+    // Constant symbols, typically as constants of integration.
+    // Should NOT be leaked out of tooling code. Marking as "interpreted" because that should prevent its unsound use in
+    // at least the most important soundness-critical contexts.
+    // @todo either enforce this invariant or give the function an obviously obnoxious name as a warning to the user.
+    else if (MathematicaOpSpec.C.applies(e))     convertUnary (e, (t: Term) => FuncOf(Function("C", None, Real, Real, true), t))
+
     // Arith expressions
     else if (MathematicaOpSpec.plus.applies(e))   convertNary  (e, Plus.apply)
     else if (MathematicaOpSpec.minus.applies(e))  convertBinary(e, Minus.apply)

keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala

@@ -67,7 +67,8 @@ object KeYmaeraX {
     val CONVERT: String = edu.cmu.cs.ls.keymaerax.cli.KeYmaeraX.Modes.CONVERT
     val SETUP: String = edu.cmu.cs.ls.keymaerax.cli.KeYmaeraX.Modes.SETUP
     val UI: String = "ui"
-    val modes: Set[String] = Set(CODEGEN, CONVERT, MODELPLEX, PROVE, REPL, UI, SETUP)
+    val TAYLORIZE: String = "taylorize"
+    val modes: Set[String] = Set(CODEGEN, CONVERT, MODELPLEX, PROVE, REPL, UI, SETUP, TAYLORIZE)
   }
 
   /** Usage -help information. */
@@ -95,6 +96,21 @@ object KeYmaeraX {
     try {
       //@todo allow multiple passes by filter architecture: -prove bla.key -tactic bla.scal -modelplex -codegen
       options.get('mode) match {
+        case Some(Modes.TAYLORIZE) => {
+          val toolConfig =
+            if (options.contains('quantitative)) {
+              configFromFile(Tools.MATHEMATICA) //@note quantitative ModelPlex uses Mathematica to simplify formulas
+            } else {
+              configFromFile("z3")
+            }
+          initializeProver(combineConfigs(options, toolConfig), usage)
+
+          val filename = options.get('file)
+          filename match {
+            case Some(s) => println(TaylorizeMain(s.asInstanceOf[String]))
+            case None => println("FAILED.")
+          }
+        }
         case Some(Modes.CODEGEN) =>
           val toolConfig =
             if (options.contains('quantitative)) {
@@ -161,6 +177,8 @@ object KeYmaeraX {
       case Nil => map
       case "-help" :: _ => println(usage); exit(1)
       // actions
+      case "-taylorize" :: value :: tail =>
+        nextOption(map ++ Map('mode -> Modes.TAYLORIZE, 'file -> value), tail)
       case "-sandbox" :: tail =>
         nextOption(map ++ Map('sandbox -> true), tail)
       case "-modelplex" :: value :: tail =>