Differentiation and canonical ordering visit each shared subexpression once

CHANGELOG.md

@@ -74,6 +74,9 @@ Bug fixes for 4.0.
 ### Fixed
 
 - A `HomotopySolver` solve's records ask now includes its start points (#498).  It named the homotopy, target, settings and seed but not the given start points, so with recording on (the default) a second solve of the same homotopy from different start points was the same ask, and was answered from the record of the first: it silently returned the first solve's endpoints.  The start system now reports the identity of start points it was given (`StartSystem::GivenStartIdentity`, a digest of their exact values), and the solve adds it to its ask.  A solve from given start points therefore has a different run id than in 4.0.0, so records written by 4.0.0 for such solves are asked again rather than recalled.
+- Differentiating an expression differentiates each distinct subexpression once, however many paths in the expression reach it (#505).  It used to differentiate a shared subexpression once per path, so a determinant built by cofactor expansion with its sub-minors shared took time proportional to k! -- 24 s for the three partial derivatives of a 7 x 7 determinant of quadratics, all of it paid on a `System`'s first evaluation, which builds its Jacobian.  `Node::Differentiate` now keeps each derivative for the length of the outermost call, and a Jacobian is built under one `node::DerivativeScope`, so subexpressions its functions share are differentiated once among all of them.  A node type's own rule is now `DifferentiateImpl`; a Python subclass of `Node` still overrides `Differentiate`.
+- Building a sum or product no longer walks its operands (#505).  Putting the operands in canonical order asks each one's variables, multidegree and, on a tie, canonical encoding, and each answer was a walk of the operand's DAG -- of its whole expansion, for `Degree(v)` under a negation -- so building an expression cost the square of its size, or worse.  Each node now keeps a degree record, made once from its operands' records: its variables, its degree and multidegree in each variable it mentions, its total degree, and its answer for every variable it does not mention.  `Degree`, `MultiDegree` and `GatherVariables` read it.  On the Barth sextic's 50 deflated critical-point systems, the first evaluations took 121 s and take 0.6 s.  A node type's degree rule is now `DegreeImpl`; a Python subclass of `Node` still overrides `Degree`.
+- **The canonicalization tie-break is a content key, and the system encoding version is `b2sysenc/3`** (#505, ADR-0073).  Operands that tie on multidegree were ordered by their whole canonical encoding, which costs a walk of the operand.  They are now ordered by a content key: the SHA-256 of the operand's own canonical encoding with each of its operands written as that operand's key -- a Merkle hash, made from the operands' keys once per node.  The first evaluation of a 7 x 7 determinant of quadratics took 24.6 s in 4.0.0 and takes 0.04 s.  Where operands tie, their order in a system's canonical encoding changes, so the encoding version moves to `b2sysenc/3` and every system digest changes: records written by 4.0.0 are asked again rather than recalled, and `System.from_canonical` refuses `b2sysenc/2` text.
 
 _______________________________________________________________________________
 

core/include/bertini2/function_tree/canonical_encoding.hpp

@@ -101,5 +101,17 @@ void EncodeCanonical(std::shared_ptr<Node> const& root, std::ostream& out, Encod
 */
 std::string CanonicalEncoding(std::shared_ptr<Node> const& root);
 
+/**
+\brief The root's own canonical encoding, with each operand written as `@` and its content key.
+
+The text a content key is the SHA-256 of (Node::ContentKey): the root's kind, payload and
+operands in the format above, but an operand is not expanded -- it is `@<its content key>`.  So a
+node's key is made from its operands' keys, once per node, and equal content gives equal keys.
+
+\param root The node.
+\return The shallow encoding text.
+*/
+std::string ShallowCanonicalEncoding(std::shared_ptr<Node> const& root);
+
 } // namespace node
 } // namespace bertini

core/include/bertini2/function_tree/node.hpp

@@ -237,11 +237,17 @@ class Node : public VisitableBase<>, public std::enable_shared_from_this<Node>
     /**
     \brief Compute the derivative with respect to a single variable.
 
-    Virtual method for differentiating the node.  If no variable is passed, produces a Jacobian tree when all is said and done, which is used for evaluating the Jacobian.
+    If no variable is passed, produces a Jacobian tree when all is said and done, which is used for evaluating the Jacobian.
 
+    Each distinct node is differentiated once per variable, however many paths reach it: the
+    derivative of a subexpression is kept for as long as the outermost call (or an enclosing
+    DerivativeScope) lasts, and every later request for it is answered from there.  A node type
+    supplies its own rule as DifferentiateImpl, which calls this on its operands.
+
+    \param v The variable to differentiate with respect to.
     \return The result of differentiation.  Jacobian or regular Node depending on what you passed in.
     */
-    virtual std::shared_ptr<Node> Differentiate(std::shared_ptr<Variable> const& v = nullptr) const = 0;
+    std::shared_ptr<Node> Differentiate(std::shared_ptr<Variable> const& v = nullptr) const;
 
     /**
     \brief Differentiate repeatedly with respect to a single variable.
@@ -270,43 +276,71 @@ class Node : public VisitableBase<>, public std::enable_shared_from_this<Node>
     /**
     Compute the degree, optionally with respect to a single variable.
 
+    Answered from the node's degree record (see DegreeRecord), which is made once per node from
+    its operands' records, so asking costs a lookup however large the expansion beneath the node.
+    The per-type rule lives in DegreeImpl.
+
     \param v Shared pointer to variable with respect to which you want to compute the degree of the Node.
     \return The degree.  Will be negative if the Node is non-polynomial.
     */
-    virtual int Degree(std::shared_ptr<Variable> const& v = nullptr) const = 0;
-
+    int Degree(std::shared_ptr<Variable> const& v = nullptr) const;
 
+    /// \brief Per-node-type degree in one variable.  Asks its operands through Degree(), so
+    /// their records answer; a node's own record is made by calling this, so it must never ask
+    /// Degree() of its own node.
+    /// \param v The variable, or null for the total degree.
+    /// \return The degree.  Negative if the node is non-polynomial.
+    virtual int DegreeImpl(std::shared_ptr<Variable> const& v) const = 0;
 
-    /**
-     Compute the multidegree with respect to a variable group.  This is for homogenization, and testing for homogeneity.
-
-     \return A vector containing the degrees.  Negative entries indicate non-polynomiality.
-    */
     /**
      Compute the multidegree with respect to a variable group.
 
-     DAG-AWARE: this is a non-virtual wrapper that memoizes by node identity FOR THE DURATION
-     OF ONE TOP-LEVEL CALL, so a shared subtree is visited once however many times it is
-     reached.  The per-type work lives in MultiDegreeImpl.
-
-     Why it exists (b2#417): the node graph is hash-consed, so a small DAG can stand for an
-     enormous expansion -- measured, a 46-node expression whose printed form is 196 KB.  A
-     plain recursive walk visits the TREE, so it costs O(expansion) where the structure is
-     O(DAG).  `CanonicalizeNaryOperands` calls this for every operand of every n-ary
-     construction, which is how a 46-node expression came to cost 13 ms to multiply by a
-     variable, and how a cofactor expansion inside a deflation stage came to run for hours
-     with a DAG that never exceeded 1102 nodes.
-
-     Caching the RESULT between calls does not fix this and was measured 29% SLOWER -- filling
-     such a cache costs the same walk it is trying to avoid, so it only pays from a node's
-     second mention.  The memo has to be inside the traversal, which is what this is.
+     Answered from the node's degree record, like Degree.  The per-type rule lives in
+     MultiDegreeImpl.
+
+     Why it matters (b2#417, #505): the node graph is hash-consed, so a small DAG can stand for
+     an enormous expansion -- measured, a 46-node expression whose printed form is 196 KB.
+     `CanonicalizeNaryOperands` asks this of every operand of every sum and product it builds.
+     A plain recursive walk visits the TREE, so it cost O(expansion), which is how a cofactor
+     expansion inside a deflation stage came to run for hours with a DAG that never exceeded
+     1102 nodes.  A walk memoized within one call (#417) costs O(DAG), and asked once per
+     construction that is still O(DAG^2) for building an expression: the derivative of a 7 x 7
+     determinant took seconds.  The record makes each construction's question O(operands).
+
+     #417 measured a cache of RESULTS between calls 29% slower than the memoized walk, because
+     filling it cost the same walk.  The record does not walk: it is made from the operands'
+     records, which already exist for any operand that has been built into a sum or product.
     */
     std::vector<int> MultiDegree(VariableGroup const& vars) const;
 
-    /// \brief Per-node-type multidegree.  Recurses through MultiDegree(), so children are
-    /// memoized; never call this directly unless you want the unmemoized walk.
+    /// \brief Per-node-type multidegree.  Asks its operands through MultiDegree(), and its own
+    /// node only through DegreeImpl, never Degree() (see DegreeImpl).
+    /// \param vars The variable group.
+    /// \return One degree per variable of the group.
     virtual std::vector<int> MultiDegreeImpl(VariableGroup const& vars) const = 0;
 
+    /**
+    \brief The variables in this expression, distinct and sorted by name (GatherVariables).
+
+    Part of the node's degree record: made once per node, from its operands' variables.
+
+    \return The variables.
+    */
+    VariableGroup const& Variables() const;
+
+    /**
+    \brief This node's content key: the SHA-256, in hex, of its shallow canonical encoding.
+
+    The shallow encoding (ShallowCanonicalEncoding) is the node's own canonical encoding with
+    each operand written as its content key, so the key is a Merkle hash of the tree's content:
+    made from the operands' keys, once per node, and equal for equal content.  Canonicalization
+    orders the operands of a sum or product by it when their multidegrees tie (ADR-0073).
+    Computed on the first request and kept: a node is immutable.
+
+    \return The content key, 64 hex digits.
+    */
+    std::string const& ContentKey() const;
+
     /**
      Compute the overall degree with respect to a variable group.
 
@@ -376,10 +410,50 @@ class Node : public VisitableBase<>, public std::enable_shared_from_this<Node>
     /// EvalProgram().  Lazily populated by EvalExpression; transient, not serialized.
     mutable std::shared_ptr<const void> eval_program_;
 
+    /**
+    \brief What Variables, Degree and MultiDegree report about this node, made once per node.
+
+    Its variables; for each variable it mentions (its variables and its differentials'
+    variables), its degree and multidegree entry in that variable; its total degree; and one
+    degree and one multidegree entry for every variable it does not mention.  The last are a
+    single pair because the degree rules compare a variable only with the variables a subtree
+    mentions, so every unmentioned variable is answered alike -- not always 0: a negative
+    constant power reports -1 in every variable.  Made from the operands' records, so it costs
+    O(operands x mentioned variables), never a walk.  Defined in node.cpp.
+    */
+    struct DegreeRecord;
+
+    /// The degree record, once made (see Record()).  Published once, atomically, and never
+    /// replaced, so a reference into it lives as long as the node.  Transient, not serialized.
+    mutable std::shared_ptr<const DegreeRecord> degree_record_;
+
+    /// The content key, once computed (see ContentKey()).  Published like degree_record_.
+    mutable std::shared_ptr<const std::string> content_key_;
+
+    /// \brief This node's degree record, made on first request from its operands' records.
+    /// \return The record.
+    DegreeRecord const& Record() const;
+
+    /// \brief A variable that is in no expression, for the record's answer about unmentioned variables.
+    /// \return The variable.
+    static std::shared_ptr<Variable> const& Unmentioned();
+
     /// Compute this node's structural hash.  Default: identity (the object address), so
     /// distinct nodes hash distinctly.  Value/operator nodes override to be structural.
     virtual std::size_t HashImpl() const;
 
+    /**
+    \brief This node type's differentiation rule, applied to this node alone.
+
+    Called by Differentiate only when the derivative is not already known.  Differentiate an
+    operand by calling its Differentiate, never its DifferentiateImpl, so that a subexpression
+    shared by many operands is differentiated once.
+
+    \param v The variable to differentiate with respect to.
+    \return The derivative of this node.
+    */
+    virtual std::shared_ptr<Node> DifferentiateImpl(std::shared_ptr<Variable> const& v) const = 0;
+
     /// Combine an extra value into a running hash (boost::hash_combine recipe).
     static void HashCombine(std::size_t& seed, std::size_t value)
     {
@@ -402,6 +476,34 @@ class Node : public VisitableBase<>, public std::enable_shared_from_this<Node>
 
 }; // re: class node
 
+    /**
+    \brief While one of these is alive, every derivative computed on this thread is kept and reused.
+
+    Node::Differentiate opens one for the length of each outermost call, so a subexpression
+    reached by many paths is differentiated once, not once per path: the sub-minors of a
+    determinant built by cofactor expansion used to be differentiated k! times at the top (#505).
+    Open one around a loop that differentiates many expressions, such as a whole Jacobian, so that
+    subexpressions they share are differentiated once among all of them.  Scopes nest: an inner one
+    uses the outer one's derivatives.  It keeps every node it has differentiated, and every
+    derivative, alive until it closes.
+    */
+    class DerivativeScope
+    {
+    public:
+        /// \brief Start keeping derivatives on this thread, unless an enclosing scope already is.
+        DerivativeScope();
+        /// \brief Stop keeping them and release them, if this scope started it.
+        ~DerivativeScope();
+        DerivativeScope(DerivativeScope const&) = delete;
+        DerivativeScope& operator=(DerivativeScope const&) = delete;
+
+        /// \cond DERIVATIVE_SCOPE_DETAIL
+        struct Derivatives;
+        /// \endcond
+    private:
+        std::unique_ptr<Derivatives> owned_;   ///< the derivatives, when this scope is the outermost
+    };
+
     /**
      a single layer of indirection, though which to call the overridden virtual print() method which must be defined for each non-abstract Node type.
      */