Wygenerowane metody oceny wielomianowej

Question

Próbuję wymyślić elegancki sposób obsługi niektórych generowanych wielomianów. Tutaj jest sytuacja skupimy się na (wyłącznie) na to pytanie:

1. zamówienia jest parametrem generowania n tego rzędu wielomianu, gdzie n: = zlecenia + 1.
2. ja jest parametrem całkowitym z zakresu 0..n
3. Wielomian ma zera przy x_j, gdzie j = 1..n i j ≠ i (w tym miejscu powinno być jasne, że StackOverflow potrzebuje nowej funkcji lub jej jest i nie wiem)
4. Ocena wielomianowa es do 1 w x_i.

Ponieważ ten konkretny przykład kodu generuje x_1 .. x_n, wyjaśnię, jak zostały znalezione w kodzie. Punkty są równomiernie rozłożone na odległość x_j = j * elementSize/order, gdzie n = order + 1.

Generuję numer Func<double, double>, aby ocenić wielomian¹.

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(); 
} 

Problem: również trzeba ocenić * F '= dψ/DX. Aby to zrobić, mogę przepisać ψ = skala × (x_0 - x) (x_1 - x) × .. × (x_n - x)/(x_i - x) w postaci ψ = α_n × x^n + α_n × x^(n-1) + .. + α_1 × x + α_0. Daje to ψ '= n × α_n × x^(n-1) + (n-1) × α_n × x^(n-2) + .. + 1 × α_1.

Ze względów obliczeniowych możemy przepisać ostateczną odpowiedź bez wywoływania numeru Math.Pow, pisząc ψ '= x × (x × (x × (..) - β_2) - β_1) - β_0.

zrobić wszystko to „oszustwo” (bardzo podstawowy algebra), muszę czystą drogę do:

1. rozwinąć uwzględnione Expression zawierający ConstantExpression i ParameterExpression liści i podstawowe operacje matematyczne (skończyć obu BinaryExpression z NodeType zestaw do operacji) - wynik tutaj może zawierać InvocationExpression elementów do MethodInfo dla Math.Pow, które będziemy obsługiwać w sposób specjalny w całym tekście.
2. Następnie przyjmuję pochodną w odniesieniu do niektórych określonych ParameterExpression. Warunki w wynikach, w których parametr po prawej stronie wywołania Math.Pow był stały 2, zostały zastąpione przez ConstantExpression(2) pomnożone przez to, co było po lewej stronie (wywołanie Math.Pow(x,1) zostało usunięte). Warunki w wyniku, które stają się zerowe, ponieważ były stałe względem x, zostały usunięte.
3. Następnie należy odliczyć instancje niektórych specyficznych ParameterExpression, gdzie występują one jako parametr po lewej stronie do wywołania Math.Pow. Po prawej stronie wywołania staje się ConstantExpression o wartości 1, zastępujemy wywołanie tylko samym ParameterExpression.

đ W przyszłości ja jak metoda wziąć ParameterExpression i zwrócić Expression że ocenia na podstawie tego parametru. W ten sposób mogę agregować wygenerowane funkcje. Jeszcze mnie nie ma. ² W przyszłości mam nadzieję udostępnić ogólną bibliotekę do pracy z wyrażeniami LINQ jako symboliczną matematyką.

Answer 1

Napisałem podstawy kilku symbolicznych funkcji matematycznych przy użyciu typu ExpressionVisitor w .NET 4. Nie jest to doskonałe, ale wygląda na to, że jest podstawą realnego rozwiązania.

• Symbolic jest public static class odsłaniając metody jak Expand, Simplify i PartialDerivative
• ExpandVisitor jest wewnętrznym typ pomocnika, który rozszerza wyrażenia
• SimplifyVisitor jest wewnętrznym typ pomocnika, który upraszcza form ekspresji
• DerivativeVisitor IS wewnętrzny typ helpera, który pobiera pochodną wyrażenia, jest jednym z następujących: wewnętrznym pomocnikiem, który przekształca Expression do zapisu prefiksu z składni 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(); 
    } 
} 

rozwijanymi wyrażenia z 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; 
    } 
} 

Biorąc pochodną z 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; 
    } 
} 

Uproszczenie wyrażenia z 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; 
    } 
} 

formatowania wyrażeń dla wyświetlacza z 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); 
    } 
} 

Testowanie wyniki

[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)); 
}