Metodi generati per la valutazione polinomiale

Question

Metodi generati per la valutazione polinomiale

algebra c# expression-trees linq math

Domanda

Sto cercando di trovare un modo elegante per gestire alcuni polinomi generati. Ecco la situazione su cui ci concentreremo (in esclusiva) per questa domanda:

ordine è un parametro nel generare un n ° ordine polinomio, dove n: = ordine + 1.
i è un parametro intero compreso nell'intervallo 0..n
Il polinomio ha zeri in x_j, dove j = 1..n e j â ‰ i (dovrebbe essere chiaro a questo punto che StackOverflow ha bisogno di una nuova funzione o è presente e io non lo so)
Il polinomio restituisce 1 a x_i.

Poiché questo particolare esempio di codice genera x_1 .. x_n, spiegherò come vengono trovati nel codice. I punti sono equamente distanziati x_j = j * elementSize / order apart, dove n = order + 1 .

Func<double, double> un Func<double, double> per valutare questo polinomio¹.

private static Func<double, double> GeneratePsi(double elementSize, int order, int i)
{
    if (order < 1)
        throw new ArgumentOutOfRangeException("order", "order must be greater than 0.");

    if (i < 0)
        throw new ArgumentOutOfRangeException("i", "i cannot be less than zero.");
    if (i > order)
        throw new ArgumentException("i", "i cannot be greater than order");

    ParameterExpression xp = Expression.Parameter(typeof(double), "x");

    // generate the terms of the factored polynomial in form (x_j - x)
    List<Expression> factors = new List<Expression>();
    for (int j = 0; j <= order; j++)
    {
        if (j == i)
            continue;

        double p = j * elementSize / order;
        factors.Add(Expression.Subtract(Expression.Constant(p), xp));
    }

    // evaluate the result at the point x_i to get scaleInv=1.0/scale.
    double xi = i * elementSize / order;
    double scaleInv = Enumerable.Range(0, order + 1).Aggregate(0.0, (product, j) => product * (j == i ? 1.0 : (j * elementSize / order - xi)));

    /* generate an expression to evaluate
     *   (x_0 - x) * (x_1 - x) .. (x_n - x) / (x_i - x)
     * obviously the term (x_i - x) is cancelled in this result, but included here to make the result clear
     */
    Expression expr = factors.Skip(1).Aggregate(factors[0], Expression.Multiply);
    // multiplying by scale forces the condition f(x_i)=1
    expr = Expression.Multiply(Expression.Constant(1.0 / scaleInv), expr);

    Expression<Func<double, double>> lambdaMethod = Expression.Lambda<Func<double, double>>(expr, xp);
    return lambdaMethod.Compile();
}

Il problema: ho anche bisogno di valutare Ï ² = dÏ / dx. Per fare ciò, posso riscrivere Ï = scaleÃ- (x_0 - x) (x_1 - x) Ã - .. Ã- (x_n - x) / (x_i - x) nella forma Ï = α ± _nÃ-x ^ n + Î ± _nÃ-x ^ (n-1) + .. + α ± _1Ã-x + α ± _0. Questo dà Ï â € œ nÃ-α ± _nÃ-x ^ (n-1) + (n-1) Ã-α ± _n-x ^ (n-2) + .. + 1Ã-α ± _1.

Per ragioni computazionali, possiamo riscrivere la risposta finale senza chiamate a Math.Pow scrivendo Ïâ² ² = xÃ- (xÃ- (xÃ- (..) - β_2) - β_1) - β_0.

Per fare tutto questo "inganno" (tutta algebra di base), ho bisogno di un modo pulito per:

Espandi Expression fattorizzata contenente le Expression ConstantExpression e ParameterExpression e le operazioni matematiche di base (termina BinaryExpression con il NodeType impostato sull'operazione) - il risultato qui può includere elementi InvocationExpression in MethodInfo per Math.Pow che gestiremo in modo speciale per tutto.
Quindi prendo la derivata rispetto ad alcuni ParameterExpression specificati. I termini nel risultato in cui il parametro del lato destro di un'invocazione di Math.Pow era la costante 2 sono sostituiti da ConstantExpression(2) moltiplicato per quello che era il lato sinistro (l'invocazione di Math.Pow(x,1) è rimosso). I termini nel risultato che diventano zero perché erano costanti rispetto a x vengono rimossi.
Quindi estrai le istanze di alcuni ParameterExpression specifici in cui si verificano come il parametro sul lato sinistro di un'invocazione di Math.Pow . Quando il lato destro dell'invocazione diventa un'espressione ConstantExpression con il valore 1 , sostituiamo l'invocazione con solo ParameterExpression .

Â¹ In futuro, mi piacerebbe che il metodo eseguisse ParameterExpression e restituisse Expression che valuta in base a quel parametro. In questo modo posso aggregare le funzioni generate. Non sono ancora arrivato. Â² In futuro, spero di pubblicare una libreria generale per lavorare con LINQ Expressions come matematica simbolica.

Sam Harwell

Vedi di più su Stack Overflow

Sam Harwell · Accepted Answer · 10/28/2009 3:28:22 AM

Ho scritto le basi di diverse funzioni matematiche simboliche utilizzando il tipo ExpressionVisitor in .NET 4. Non è perfetto, ma sembra il fondamento di una soluzione praticabile.

Symbolic è una classe statica pubblica che espone metodi come Expand , Simplify e PartialDerivative
ExpandVisitor è un tipo di helper interno che espande le espressioni
SimplifyVisitor è un tipo di helper interno che semplifica le espressioni
DerivativeVisitor è un tipo di helper interno che prende la derivata di un'espressione
ListPrintVisitor è un tipo di helper interno che converte Expression in una notazione di prefisso con una sintassi Lisp

`Symbolic`

public static class Symbolic
{
    public static Expression Expand(Expression expression)
    {
        return new ExpandVisitor().Visit(expression);
    }

    public static Expression Simplify(Expression expression)
    {
        return new SimplifyVisitor().Visit(expression);
    }

    public static Expression PartialDerivative(Expression expression, ParameterExpression parameter)
    {
        bool totalDerivative = false;
        return new DerivativeVisitor(parameter, totalDerivative).Visit(expression);
    }

    public static string ToString(Expression expression)
    {
        ConstantExpression result = (ConstantExpression)new ListPrintVisitor().Visit(expression);
        return result.Value.ToString();
    }
}

Espressioni espandibili con `ExpandVisitor`

internal class ExpandVisitor : ExpressionVisitor
{
    protected override Expression VisitBinary(BinaryExpression node)
    {
        var left = Visit(node.Left);
        var right = Visit(node.Right);

        if (node.NodeType == ExpressionType.Multiply)
        {
            Expression[] leftNodes = GetAddedNodes(left).ToArray();
            Expression[] rightNodes = GetAddedNodes(right).ToArray();
            var result =
                leftNodes
                .SelectMany(x => rightNodes.Select(y => Expression.Multiply(x, y)))
                .Aggregate((sum, term) => Expression.Add(sum, term));

            return result;
        }

        if (node.Left == left && node.Right == right)
            return node;

        return Expression.MakeBinary(node.NodeType, left, right, node.IsLiftedToNull, node.Method, node.Conversion);
    }

    /// <summary>
    /// Treats the <paramref name="node"/> as the sum (or difference) of one or more child nodes and returns the
    /// the individual addends in the sum.
    /// </summary>
    private static IEnumerable<Expression> GetAddedNodes(Expression node)
    {
        BinaryExpression binary = node as BinaryExpression;
        if (binary != null)
        {
            switch (binary.NodeType)
            {
            case ExpressionType.Add:
                foreach (var n in GetAddedNodes(binary.Left))
                    yield return n;

                foreach (var n in GetAddedNodes(binary.Right))
                    yield return n;

                yield break;

            case ExpressionType.Subtract:
                foreach (var n in GetAddedNodes(binary.Left))
                    yield return n;

                foreach (var n in GetAddedNodes(binary.Right))
                    yield return Expression.Negate(n);

                yield break;

            default:
                break;
            }
        }

        yield return node;
    }
}

Prendendo un derivato con `DerivativeVisitor`

internal class DerivativeVisitor : ExpressionVisitor
{
    private ParameterExpression _parameter;
    private bool _totalDerivative;

    public DerivativeVisitor(ParameterExpression parameter, bool totalDerivative)
    {
        if (_totalDerivative)
            throw new NotImplementedException();

        _parameter = parameter;
        _totalDerivative = totalDerivative;
    }

    protected override Expression VisitBinary(BinaryExpression node)
    {
        switch (node.NodeType)
        {
        case ExpressionType.Add:
        case ExpressionType.Subtract:
            return Expression.MakeBinary(node.NodeType, Visit(node.Left), Visit(node.Right));

        case ExpressionType.Multiply:
            return Expression.Add(Expression.Multiply(node.Left, Visit(node.Right)), Expression.Multiply(Visit(node.Left), node.Right));

        case ExpressionType.Divide:
            return Expression.Divide(Expression.Subtract(Expression.Multiply(Visit(node.Left), node.Right), Expression.Multiply(node.Left, Visit(node.Right))), Expression.Power(node.Right, Expression.Constant(2)));

        case ExpressionType.Power:
            if (node.Right is ConstantExpression)
            {
                return Expression.Multiply(node.Right, Expression.Multiply(Visit(node.Left), Expression.Subtract(node.Right, Expression.Constant(1))));
            }
            else if (node.Left is ConstantExpression)
            {
                return Expression.Multiply(node, MathExpressions.Log(node.Left));
            }
            else
            {
                return Expression.Multiply(node, Expression.Add(
                    Expression.Multiply(Visit(node.Left), Expression.Divide(node.Right, node.Left)),
                    Expression.Multiply(Visit(node.Right), MathExpressions.Log(node.Left))
                    ));
            }

        default:
            throw new NotImplementedException();
        }
    }

    protected override Expression VisitConstant(ConstantExpression node)
    {
        return MathExpressions.Zero;
    }

    protected override Expression VisitInvocation(InvocationExpression node)
    {
        MemberExpression memberExpression = node.Expression as MemberExpression;
        if (memberExpression != null)
        {
            var member = memberExpression.Member;
            if (member.DeclaringType != typeof(Math))
                throw new NotImplementedException();

            switch (member.Name)
            {
            case "Log":
                return Expression.Divide(Visit(node.Expression), node.Expression);

            case "Log10":
                return Expression.Divide(Visit(node.Expression), Expression.Multiply(Expression.Constant(Math.Log(10)), node.Expression));

            case "Exp":
            case "Sin":
            case "Cos":
            default:
                throw new NotImplementedException();
            }
        }

        throw new NotImplementedException();
    }

    protected override Expression VisitParameter(ParameterExpression node)
    {
        if (node == _parameter)
            return MathExpressions.One;

        return MathExpressions.Zero;
    }
}

Semplificazione delle espressioni con `SimplifyVisitor`

internal class SimplifyVisitor : ExpressionVisitor
{
    protected override Expression VisitBinary(BinaryExpression node)
    {
        var left = Visit(node.Left);
        var right = Visit(node.Right);

        ConstantExpression leftConstant = left as ConstantExpression;
        ConstantExpression rightConstant = right as ConstantExpression;
        if (leftConstant != null && rightConstant != null
            && (leftConstant.Value is double) && (rightConstant.Value is double))
        {
            double leftValue = (double)leftConstant.Value;
            double rightValue = (double)rightConstant.Value;

            switch (node.NodeType)
            {
            case ExpressionType.Add:
                return Expression.Constant(leftValue + rightValue);
            case ExpressionType.Subtract:
                return Expression.Constant(leftValue - rightValue);
            case ExpressionType.Multiply:
                return Expression.Constant(leftValue * rightValue);
            case ExpressionType.Divide:
                return Expression.Constant(leftValue / rightValue);
            default:
                throw new NotImplementedException();
            }
        }

        switch (node.NodeType)
        {
        case ExpressionType.Add:
            if (IsZero(left))
                return right;
            if (IsZero(right))
                return left;
            break;

        case ExpressionType.Subtract:
            if (IsZero(left))
                return Expression.Negate(right);
            if (IsZero(right))
                return left;
            break;

        case ExpressionType.Multiply:
            if (IsZero(left) || IsZero(right))
                return MathExpressions.Zero;
            if (IsOne(left))
                return right;
            if (IsOne(right))
                return left;
            break;

        case ExpressionType.Divide:
            if (IsZero(right))
                throw new DivideByZeroException();
            if (IsZero(left))
                return MathExpressions.Zero;
            if (IsOne(right))
                return left;
            break;

        default:
            throw new NotImplementedException();
        }

        return Expression.MakeBinary(node.NodeType, left, right);
    }

    protected override Expression VisitUnary(UnaryExpression node)
    {
        var operand = Visit(node.Operand);

        ConstantExpression operandConstant = operand as ConstantExpression;
        if (operandConstant != null && (operandConstant.Value is double))
        {
            double operandValue = (double)operandConstant.Value;

            switch (node.NodeType)
            {
            case ExpressionType.Negate:
                if (operandValue == 0.0)
                    return MathExpressions.Zero;

                return Expression.Constant(-operandValue);

            default:
                throw new NotImplementedException();
            }
        }

        switch (node.NodeType)
        {
        case ExpressionType.Negate:
            if (operand.NodeType == ExpressionType.Negate)
            {
                return ((UnaryExpression)operand).Operand;
            }

            break;

        default:
            throw new NotImplementedException();
        }

        return Expression.MakeUnary(node.NodeType, operand, node.Type);
    }

    private static bool IsZero(Expression expression)
    {
        ConstantExpression constant = expression as ConstantExpression;
        if (constant != null)
        {
            if (constant.Value.Equals(0.0))
                return true;
        }

        return false;
    }

    private static bool IsOne(Expression expression)
    {
        ConstantExpression constant = expression as ConstantExpression;
        if (constant != null)
        {
            if (constant.Value.Equals(1.0))
                return true;
        }

        return false;
    }
}

Formattazione di espressioni da visualizzare con `ListPrintVisitor`

internal class ListPrintVisitor : ExpressionVisitor
{
    protected override Expression VisitBinary(BinaryExpression node)
    {
        string op = null;

        switch (node.NodeType)
        {
        case ExpressionType.Add:
            op = "+";
            break;
        case ExpressionType.Subtract:
            op = "-";
            break;
        case ExpressionType.Multiply:
            op = "*";
            break;
        case ExpressionType.Divide:
            op = "/";
            break;
        default:
            throw new NotImplementedException();
        }

        var left = Visit(node.Left);
        var right = Visit(node.Right);
        string result = string.Format("({0} {1} {2})", op, ((ConstantExpression)left).Value, ((ConstantExpression)right).Value);
        return Expression.Constant(result);
    }

    protected override Expression VisitConstant(ConstantExpression node)
    {
        if (node.Value is string)
            return node;

        return Expression.Constant(node.Value.ToString());
    }

    protected override Expression VisitParameter(ParameterExpression node)
    {
        return Expression.Constant(node.Name);
    }
}

Test dei risultati

[TestMethod]
public void BasicSymbolicTest()
{
    ParameterExpression x = Expression.Parameter(typeof(double), "x");
    Expression linear = Expression.Add(Expression.Constant(3.0), x);
    Assert.AreEqual("(+ 3 x)", Symbolic.ToString(linear));

    Expression quadratic = Expression.Multiply(linear, Expression.Add(Expression.Constant(2.0), x));
    Assert.AreEqual("(* (+ 3 x) (+ 2 x))", Symbolic.ToString(quadratic));

    Expression expanded = Symbolic.Expand(quadratic);
    Assert.AreEqual("(+ (+ (+ (* 3 2) (* 3 x)) (* x 2)) (* x x))", Symbolic.ToString(expanded));
    Assert.AreEqual("(+ (+ (+ 6 (* 3 x)) (* x 2)) (* x x))", Symbolic.ToString(Symbolic.Simplify(expanded)));

    Expression derivative = Symbolic.PartialDerivative(expanded, x);
    Assert.AreEqual("(+ (+ (+ (+ (* 3 0) (* 0 2)) (+ (* 3 1) (* 0 x))) (+ (* x 0) (* 1 2))) (+ (* x 1) (* 1 x)))", Symbolic.ToString(derivative));

    Expression simplified = Symbolic.Simplify(derivative);
    Assert.AreEqual("(+ 5 (+ x x))", Symbolic.ToString(simplified));
}

Metodi generati per la valutazione polinomiale

Domanda

Fastest Entity Framework Extensions

Risposta accettata

`Symbolic`

Espressioni espandibili con `ExpandVisitor`

Prendendo un derivato con `DerivativeVisitor`

Semplificazione delle espressioni con `SimplifyVisitor`

Formattazione di espressioni da visualizzare con `ListPrintVisitor`

Test dei risultati

Performance

Expression Evaluator

Metodi generati per la valutazione polinomiale

Domanda

Fastest Entity Framework Extensions

Risposta accettata

Symbolic

Espressioni espandibili con ExpandVisitor

Prendendo un derivato con DerivativeVisitor

Semplificazione delle espressioni con SimplifyVisitor

Formattazione di espressioni da visualizzare con ListPrintVisitor

Test dei risultati

Performance

Expression Evaluator

`Symbolic`

Espressioni espandibili con `ExpandVisitor`

Prendendo un derivato con `DerivativeVisitor`

Semplificazione delle espressioni con `SimplifyVisitor`

Formattazione di espressioni da visualizzare con `ListPrintVisitor`