Analytische Berechnung von Derivaten auf C ++ - Vorlagen

    Hier haben sie neulich über das analytische Finden von Derivaten geschrieben, was mich an eine meiner kleinen C ++ - Bibliotheken erinnerte, die fast dasselbe tut, aber zur Kompilierungszeit.


    Was ist der Gewinn? Die Antwort ist einfach: Ich musste versuchen, das Minimum einer ziemlich komplizierten Funktion zu finden, die Ableitungen dieser Funktion anhand ihrer Parameter zu betrachten, da ein Stift auf einem Blatt Papier Faulheit war, und dann überprüfen, ob ich beim Schreiben des Codes nicht versiegelt war, und denselben Code beizubehalten, war doppelt faul, und so wurde beschlossen, zu schreiben das Ding, das es für mich tun wird. Nun, damit Sie so etwas in Ihren Code schreiben können:

    using Formula_t = decltype (k * (_1 - r0) / (_1 + r0) * (g0 / (alpha0 - logr0 / Num<300>) - _1));    // сама формула
    const auto residual = Formula_t::Eval (datapoint) - knownValue;    // регрессионный остаток
    // производные по параметрам:
    const auto dg0 = VarDerivative_t::Eval (datapoint);
    const auto dalpha0 = VarDerivative_t::Eval (datapoint);
    const auto dk = VarDerivative_t::Eval (datapoint);

    Anstelle von Krokodilen, die sich herausstellen werden, wenn wir zu Beginn teilweise abgeleitete Funktionen im Bild übernehmen (oder besser gesagt, eine vereinfachte Version davon, aber es sieht nicht so beängstigend aus).

    Es ist auch gut, sicher zu sein, dass der Compiler es so optimiert, als ob die entsprechenden Ableitungen und Funktionen von Hand geschrieben wurden. Und ich möchte sicher sein - man musste das Minimum oft finden (wirklich viel, irgendwo zwischen Hunderten von Millionen und einer Milliarde, das war die Essenz einer Art Computerexperiment), daher wäre die Berechnung von Derivaten ein Engpass, es geschah zur Laufzeit Jede Rekursion in der Baumstruktur. Wenn wir den Compiler zwingen, die Ableitung tatsächlich während der Kompilierung zu berechnen, besteht die Möglichkeit, dass sie vom Optimierer gemäß dem resultierenden Code übergeben wird, und wir verlieren nicht im Vergleich zum manuellen Ausschreiben aller Ableitungen. Die Chance wurde übrigens realisiert.

    Unter dem Schnitt - eine kurze Beschreibung, wie es dort funktioniert.

    Beginnen wir mit der Darstellung der Funktion im Programm. Aus irgendeinem Grund ist es so gekommen, dass jede Funktion ein Typ ist. Eine Funktion ist auch ein Ausdrucksbaum, und der Knoten dieses Baums wird durch einen Typ dargestellt Node:

    template
    struct Node;

    Hier NodeClassist der Knotentyp (Variable, Zahl, unäre Funktion, Binärfunktion), Argsdie Parameter dieses Knotens (Variablenindex, Wert der Zahl, untergeordnete Knoten).

    Knoten können sich für bestimmte Werte freier Variablen und Parameter differenzieren, drucken und berechnen. Wenn also ein Typ definiert ist, der einen Knoten mit einer regulären Nummer darstellt:

    using NumberType_t = long long;
    template
    struct Number {};

    dann ist die Spezialisierung des Knotens für Zahlen trivial:

    template
    struct Node>
    {
    	template
    	using Derivative_t = Node>;
    	static std::string Print ()
    	{
    		return std::to_string (N);
    	}
    	template
    	static typename Vec::value_type Eval (const Vec&)
    	{
    		return N;
    	}
    	constexpr Node () {}
    };

    Die Ableitung einer beliebigen Zahl in Bezug auf eine Variable ist Null (der Typ ist dafür verantwortlich Derivative_t, lassen wir seine Vorlagenparameter vorerst). Das Drucken einer Nummer ist ebenfalls einfach (siehe Print()). Berechnen Sie einen Knoten mit einer Nummer - geben Sie diese Nummer zurück (siehe Eval()den Vorlagenparameter wird später Vecerläutert).

    Die Variable wird auf ähnliche Weise angezeigt:

    template
    struct Variable {};

    Hier Familyund Indexist die "Familie" und der Index der Variablen. Also, denn sie werden gleich 'w'und 1und für - 'x'und 2dementsprechend.

    Der Knoten für die Variable ist etwas interessanter als für die Zahl bestimmt:

    template
    struct Node>
    {
    	template
    	using Derivative_t = std::conditional_t>,
    			Node>>;
    	static std::string Print ()
    	{
    		return std::string { Family, '_' } + std::to_string (Index);
    	}
    	template
    	static typename Vec::value_type Eval (const Vec& values)
    	{
    		return values (Node {});
    	}
    	constexpr Node () {}
    };

    Die Ableitung einer Variablen in Bezug auf sie ist also gleich Eins und zu jeder anderen gleich Null. Eigentlich sind die Parameter FPrimeund die IPrimeArt Derivative_t- eine Familie und eine Indexvariable , für die Sie die Ableitung nehmen wollen.

    Die Berechnung des Wertes einer aus einer Variablen bestehenden Funktion wird auf ihren Befund im Wertewörterbuch reduziert values, der an die Funktion übergeben wird Eval(). Das Wörterbuch selbst kann den Wert der gewünschten Variablen anhand ihres Typs ermitteln. Wir übergeben ihm also einfach den Typ unserer Variablen und geben den entsprechenden Wert zurück. Wie das Wörterbuch dies macht, werden wir später betrachten.

    Mit unären Funktionen werden die Dinge noch interessanter.

    enum class UnaryFunction
    {
    	Sin,
    	Cos,
    	Ln,
    	Neg
    };
    template
    struct UnaryFunctionWrapper;

    In der Spezialisierung werden UnaryFunctionWrapperwir die Logik vorantreiben, die Ableitungen jeder spezifischen unären Funktion zu übernehmen. Um die Codeduplizierung zu minimieren, nehmen wir die Ableitung der unären Funktion in Bezug auf ihr Argument. Der Aufrufer ist für die weitere Differenzierung des Arguments in Bezug auf die Zielvariable durch die Kettenregel verantwortlich:

    template<>
    struct UnaryFunctionWrapper
    {
    	template
    	using Derivative_t = Node;
    };
    template<>
    struct UnaryFunctionWrapper
    {
    	template
    	using Derivative_t = Node>;
    };
    template<>
    struct UnaryFunctionWrapper
    {
    	template
    	using Derivative_t = Node>, Child>;
    };
    template<>
    struct UnaryFunctionWrapper
    {
    	template
    	using Derivative_t = Node>;
    };

    Dann ist der Knoten selbst wie folgt:

    template
    struct Node, Node>
    {
    	using Child_t = Node;
    	template
    	using Derivative_t = Node::template Derivative_t,
    			typename Node::template Derivative_t>;
    	static std::string Print ()
    	{
    		return FunctionName (UF) + "(" + Node::Print () + ")";
    	}
    	template
    	static typename Vec::value_type Eval (const Vec& values)
    	{
    		const auto child = Child_t::Eval (values);
    		return EvalUnary (UnaryFunctionWrapper {}, child);
    	}
    };

    Wir betrachten die Ableitung durch die Kettenregel - es sieht beängstigend aus, die Idee ist einfach. Die Berechnung ist auch einfach: Wir berechnen den Wert des untergeordneten Knotens, dann berechnen wir den Wert unserer unären Funktion auf diesem Wert unter Verwendung der Funktion EvalUnary(). Vielmehr Funktionsfamilien: Das erste Argument für die Funktion ist der Typ, der unsere unäre Funktion definiert, um die Auswahl der gewünschten Überladung zur Kompilierungszeit zu gewährleisten. Ja, es wäre möglich, den Wert selbst zu vermitteln UF, und ein intelligenter Compiler würde mit ziemlicher Sicherheit alle erforderlichen konstanten Ausbreitungsdurchläufe ausführen, aber hier ist es einfacher, auf Nummer sicher zu gehen.

    Übrigens hätte keine separate unäre Negationsoperation eingeführt werden können, die durch Multiplikation mit minus eins ersetzt worden wäre.

    Bei binären Knoten ist alles ähnlich, nur Ableitungen sehen absolut beängstigend aus. Zur Teilung zum Beispiel:

    template<>
    struct BinaryFunctionWrapper
    {
    	template
    	using Derivative_t = Node,
    						V
    					>,
    					Node
    						>
    					>
    				>,
    				Node
    			>;
    };

    Dann wird die gewünschte Metafunktion VarDerivative_tganz einfach bestimmt, weil tatsächlich nur Derivative_tder an sie übergebene Knoten aufgerufen wird :

    template
    struct VarDerivative;
    template
    struct VarDerivative>>
    {
    	using Result_t = typename Expr::template Derivative_t;
    };
    template
    using VarDerivative_t = typename VarDerivative>::Result_t;

    Wenn Sie jetzt Hilfsvariablen und -typen definieren, zum Beispiel:

    // алиасы для типов унарных и бинарных функций:
    using Sin = UnaryFunctionWrapper;
    using Cos = UnaryFunctionWrapper;
    using Neg = UnaryFunctionWrapper;
    using Ln = UnaryFunctionWrapper;
    using Add = BinaryFunctionWrapper;
    using Mul = BinaryFunctionWrapper;
    using Div = BinaryFunctionWrapper;
    using Pow = BinaryFunctionWrapper;
    // variable template из C++14 для определения переменной в общем виде:
    template
    constexpr Node> Var {};
    // определим переменную x0 для удобства, авось, ей часто пользоваться будут:
    using X0 = Node>;
    constexpr X0 x0;
    // и так далее для других переменных
    // константа для единицы, единица часто встречается в формулах:
    constexpr Node> _1;
    // перегрузки операторов, им даже не нужно тело, достаточно типа:
    template
    Node, std::decay_t> operator+ (T1, T2);
    template
    Node, std::decay_t> operator* (T1, T2);
    template
    Node, std::decay_t> operator/ (T1, T2);
    template
    Node, Node>> operator- (T1, T2);
    // не совсем операторы, но тоже чтобы удобно писать было:
    template
    Node> Sin (T);
    template
    Node> Cos (T);
    template
    Node> Ln (T);

    Es wird möglich sein, Code ganz am Anfang des Beitrags zu schreiben.

    Was bleibt übrig?

    Behandeln Sie zunächst den Typ, der an die Funktion übergeben wird Eval(). Erwähnen Sie zweitens die Möglichkeit von Transformationen des gewünschten Ausdrucks durch Ersetzen eines Teilbaums durch einen anderen. Beginnen wir mit dem zweiten, es ist einfacher.

    Motivation (Sie können überspringen): Wenn Sie den Code, der mit der aktuellen Version funktioniert, leicht profilieren, fällt Ihnen auf, dass die Berechnung viel Zeit in Anspruch nimmt , was im Allgemeinen für jeden experimentellen Punkt gleich ist. Es spielt keine Rolle! Wir führen eine separate Variable ein, die wir einmal berechnen, bevor wir die Werte unserer Formel an jedem der experimentellen Punkte berechnen, und ersetzen alle Vorkommenzu dieser Variablen (tatsächlich wurde dies im Motivationscode ganz am Anfang bereits getan). Wenn wir jedoch die Ableitung in Bezug auf nehmen , müssen wir uns daran erinnern, dass im Allgemeinen kein freier Parameter, sondern eine Funktion von . Denken Sie daran , sehr einfach: ersetzen auf (die Meta-Funktionen ApplyDependency_t, obwohl es besser wäre , es zu nennen Rewrite_toder so ähnlich), wir unterscheiden, zurück zu zurück:

    using Unwrapped_t = ApplyDependency_t;
    using Derivative_t = VarDerivative_t;
    using CacheLog_t = ApplyDependency_t;

    Die Implementierung ist ausführlich, aber ideologisch einfach. Wir steigen rekursiv auf den Formelbaum ab und ersetzen das Baumelement, wenn es genau mit der Vorlage übereinstimmt. Andernfalls ändern wir nichts. Insgesamt gibt es drei Spezialisierungen: für den Abstieg entlang eines untergeordneten Knotens einer unären Funktion, für den Abstieg auf einen untergeordneten Knoten einer Binärfunktion und tatsächlich für den Ersatz, während Spezialisierungen für den Abstieg auf untergeordneten Knoten sicherstellen müssen, dass die Vorlage nicht mit dem Teilbaum übereinstimmt, der der betrachteten Unterfunktion entspricht:

    template
    struct ApplyDependency
    {
    	using Result_t = Formula;
    };
    template
    using ApplyDependency_t = typename ApplyDependency, std::decay_t, Formula>::Result_t;
    template
    struct ApplyDependency, Child>,
    		std::enable_if_t, Child>>::value>>
    {
    	using Result_t = Node<
    				UnaryFunctionWrapper,
    				ApplyDependency_t
    			>;
    };
    template
    struct ApplyDependency, FirstNode, SecondNode>,
    		std::enable_if_t, FirstNode, SecondNode>>::value>>
    {
    	using Result_t = Node<
    				BinaryFunctionWrapper,
    				ApplyDependency_t,
    				ApplyDependency_t
    			>;
    };
    template
    struct ApplyDependency
    {
    	using Result_t = Expr;
    };

    Puh. Es bleibt die Übertragung von Parameterwerten zu behandeln.

    Denken Sie daran, dass jeder Parameter seinen eigenen Typ hat. Wenn wir also eine Familie von Funktionen erstellen, die durch den Parametertyp überladen sind, von denen jeder den entsprechenden Wert zurückgibt, besteht wiederum (genau wie bei der Berechnung unärer Funktionen etwas früher) die Möglichkeit, dass der Compiler dies und minimiert optimiert (und er übrigens und optimiert so ein kluges Mädchen). Nun, so etwas wie:

    auto GetValue (Variable<'x', 0>)
    {
        return value_for_x0;
    }
    auto GetValue (Variable<'x', 1>)
    {
        return value_for_x1;
    }
    ...

    Nur wir wollen es schön machen, damit wir zum Beispiel schreiben können:

    BuildFunctor (g0, someValue,
            alpha0, anotherValue,
            k, yetOneMoreValue,
            r0, independentVariable,
            logr0, logOfTheIndependentVariable);

    Dabei g0sind alpha0und das Unternehmen Objekte mit den Typen der entsprechenden Variablen, gefolgt von den entsprechenden Werten.

    Wie können wir einen Igel und einen Igel kreuzen, nachdem wir im Allgemeinen eine Funktion erstellt haben, deren Parametertyp in der Comp-Zeit und der Wert in der Laufzeit festgelegt ist? Lambdas eilen zur Rettung!

    template
    auto BuildFunctor (NodeType, ValueType val)
    {
        return [val] (NodeType) { return val; };
    }

    Angenommen, wir haben zwei solche Funktionen. Wie können wir eine Funktionsfamilie in einem Namespace erhalten, sodass die gewünschte durch Überladung ausgewählt wird? Vererbung eilt zur Rettung!

    template
    struct Map : F, S
    {
    	using F::operator();
    	using S::operator();
    	Map (F f, S s)
    	: F { std::forward (f) }
    	, S { std::forward (s) }
    	{
    	}
    };

    Wir erben von beiden Lambdas (schließlich wird Lambda zu einer Struktur mit einem vom Compiler generierten Namen erweitert, was bedeutet, dass Sie davon erben können) und bringen deren Operator-Klammern in den Bereich.

    Darüber hinaus können Sie nicht nur von Lambdas erben, sondern auch von beliebigen Strukturen, die Klammeroperatoren haben. Hoppla, ich habe Algebra. Wenn es also N Lambdas gibt, können Sie das erste Mapvon den ersten beiden Lambdas erben , das nächste Mapvon den ersten Mapund nächsten Lambdas und so weiter. Lassen Sie es uns in Form von Code setzen:

    template
    auto Augment (F&& f)
    {
    	return f;
    }
    template
    auto Augment (F&& f, S&& s)
    {
    	return Map, std::decay_t> { f, s };
    }
    template
    auto BuildFunctor ()
    {
    	struct
    	{
    		ValueType operator() () const
    		{
    			return {};
    		}
    		using value_type = ValueType;
    	} dummy;
    	return dummy;
    }
    template
    auto BuildFunctor (NodeType, ValueType val, Tail&&... tail)
    {
    	return detail::Augment ([val] (NodeType) { return val; },
    			BuildFunctor (std::forward (tail)...));
    }

    Wir erhalten automatisch Vollständigkeit und Eindeutigkeit: Wenn einige Argumente nicht angegeben werden, ist dies ein Kompilierungsfehler, und wenn einige Argumente zweimal angegeben werden.

    Eigentlich ist das alles.

    Es sei denn, einer meiner Freunde, dem ich es damals gezeigt habe, schlug meiner Meinung nach eine elegantere Lösung für constexpr-Funktionen vor, aber ich habe sie seit 9 Monaten nicht mehr erreicht.

    Nun, der Link zur Bibliothek: I Am Mad . Nicht produktionsbereit, Pullrequests werden akzeptiert und so weiter.

    Nun, Sie können sich immer noch fragen, wie intelligent moderne Compiler sind, die diese Vorlagenebenen über Vorlagen über Lambdas über Vorlagen hinweg verarbeiten und ziemlich optimalen Code generieren können.

    Jetzt auch beliebt: