add multidual number

docs/source/api.rst

@@ -26,13 +26,45 @@ Dual Numbers
 
       dual part of the dual number 
 
-.. function:: proc prim(a: DualNumber)
+.. record:: MultiDual
+
+   A multidual number is a number if the form :math:`a + \sum_{i=1}^nb_i\epsilon_i`,
+   for which it holds :math:`\epsilon_i^2=0` and :math:`\epsilon_i\epsilon_j=0` for any indices :math:`i, j`. 
+
+
+   .. attribute:: var dom: domain(1)
+
+      domain for the array of dual parts
+
+   .. attribute:: var prim: real
+
+      primal part 
+
+   .. attribute:: var dual: [dom] real
+
+      dual parts 
+
+   .. method:: proc init(val: real, grad: [?dom])
+
+
+      constructor to create the multidual number from the primal part and array of dual parts, automatically inferring the domain.
+
+
+.. function:: proc todual(val: real, der: real)
+
+   Converts a pair of real numbers to dual number 
+
+.. function:: proc todual(val: real, grad: [?D] real)
+
+   Converts a real number and array of reals to a multidual number. 
+
+.. function:: proc prim(a)
 
 
    For dual numbers, it returns the primal part. For real numbers, it returns the number itself.
 
 
-.. function:: proc dual(a: DualNumber)
+.. function:: proc dual(a)
 
 
    For dual numbers, it returns the dual part, for real numbers it returns zero.

src/ForwardModeAD.chpl

@@ -1,46 +1,6 @@
 module ForwardModeAD {
-  /*
-  A dual number is a number in the form :math:`a + b\epsilon`, for which :math:`\epsilon^2 = 0`.
-  Via the algebra induced by this definition one can show that :math:`f(a + b\epsilon) = f(a) + bf'(a)\epsilon`
-  and hence dual numbers can be used for forward mode automatic differentiation by operators overloading.
-  */
-  record DualNumber {
-    /* primal part of the dual number */
-    var prim : real;
-    /* dual part of the dual number */
-    var dual : real;
-  }
-
-  /*
-  For dual numbers, it returns the primal part. For real numbers, it returns the number itself.
-  */
-  proc prim(a : DualNumber) {return a.prim;}
-
-  /*
-  For dual numbers, it returns the dual part, for real numbers it returns zero.
-  */
-  proc dual(a : DualNumber) {return a.dual;}
-
-  pragma "no doc"
-  proc prim(a : real) {return a;}
-
-  pragma "no doc"
-  proc dual(a : real) {return 0.0;}
-
-  pragma "no doc"
-  proc isDualNumberType(type t : DualNumber) param {return true;}
-
-  pragma "no doc"
-  proc isDualNumberType(type t) param {return false;}
-
-  pragma "no doc"
-  proc isEitherDualNumberType(type t, type s) param {return isDualNumberType(t) || isDualNumberType(s);}
 
-
-  pragma "no doc"
-  proc isclose(a, b, rtol=1e-5, atol=0.0) where isEitherDualNumberType(a.type, b.type) {
-      return isclose(prim(a), prim(b), rtol=rtol, atol=atol) && isclose(dual(a), dual(b), rtol=rtol, atol=atol);
-  }
+  public use dual;
 
   public use arithmetic;
 

src/arithmetic.chpl

@@ -1,37 +1,53 @@
 module arithmetic {
     use ForwardModeAD;
 
-    operator +(a : DualNumber) { return a; }
+    operator +(a) where isDualType(a.type) { return a; }
 
-    operator -(a : DualNumber) { return new DualNumber(-prim(a), -dual(a)); }
+    operator -(a) where isDualType(a.type) { 
+        var f = -prim(a),
+            df = -dual(a);
+        return todual(f, df); 
+    }
 
     operator +(a, b) where isEitherDualNumberType(a.type, b.type) {
-        return new DualNumber(prim(a) + prim(b), dual(a) + dual(b));
+        var f = prim(a) + prim(b);
+        var df = dual(a) + dual(b);
+        return todual(f, df);
     }
 
     operator -(a, b) where isEitherDualNumberType(a.type, b.type) {
-        return new DualNumber(prim(a) - prim(b), dual(a) - dual(b));
+        var f = prim(a) - prim(b);
+        var df = dual(a) - dual(b);
+        return todual(f, df);
     }
 
     operator *(a, b) where isEitherDualNumberType(a.type, b.type) {
-        return new DualNumber(prim(a) * prim(b), dual(a) * prim(b) + prim(a) * dual(b));
+        var f = prim(a) * prim(b),
+            df = dual(a) * prim(b) + prim(a) * dual(b);
+        return todual(f, df);
     }
 
     operator /(a, b) where isEitherDualNumberType(a.type, b.type) {
-        var f = prim(a) / prim(b);
-        var df = (dual(a) * prim(b) - prim(a) * dual(b)) / prim(b) ** 2;
-        return new DualNumber(f, df); 
+        var f = prim(a) / prim(b),
+            df = (dual(a) * prim(b) - prim(a) * dual(b)) / prim(b) ** 2;
+        return todual(f, df); 
     }
 
-    operator **(a : DualNumber, b : real) {
-        return new DualNumber(prim(a) ** b, b * (prim(a) ** (b - 1)) * dual(a));
+    operator **(a, b : real) where isDualType(a.type) {
+        var f = prim(a) ** b,
+            df = b * (prim(a) ** (b - 1)) * dual(a);
+        return todual(f, df);
     }
 
-    proc sqrt(a : DualNumber) {
-        return new DualNumber(sqrt(prim(a)), 0.5 * dual(a) / sqrt(prim(a)));
+    proc sqrt(a) where isDualType(a.type) {
+        var f = sqrt(prim(a)),
+            df = 0.5 * dual(a) / sqrt(prim(a));
+        return todual(f, df);
     }
 
-    proc cbrt(a : DualNumber) {
-        return new DualNumber(cbrt(prim(a)), 1.0 / 3.0 * dual(a) / cbrt(prim(a) ** 2));
+    proc cbrt(a) where isDualType(a.type) {
+        var f = cbrt(prim(a)),
+            df =  1.0 / 3.0 * dual(a) / cbrt(prim(a) ** 2);
+        return todual(f, df);
     }
 }

src/dual.chpl

@@ -0,0 +1,71 @@
+  /*
+  A dual number is a number in the form :math:`a + b\epsilon`, for which :math:`\epsilon^2 = 0`.
+  Via the algebra induced by this definition one can show that :math:`f(a + b\epsilon) = f(a) + bf'(a)\epsilon`
+  and hence dual numbers can be used for forward mode automatic differentiation by operators overloading.
+  */
+  record DualNumber {
+
+    /* primal part of the dual number */
+    var prim : real;
+
+    /* dual part of the dual number */
+    var dual : real;
+  }
+
+  /* A multidual number is a number if the form :math:`a + \sum_{i=1}^nb_i\epsilon_i`,
+     for which it holds :math:`\epsilon_i^2=0` and :math:`\epsilon_i\epsilon_j=0` for any indices :math:`i, j`. */
+  record MultiDual {
+
+    /* domain for the array of dual parts*/
+    var dom: domain(1);
+
+    /* primal part */
+    var prim: real;
+
+    /* dual parts */
+    var dual: [dom] real;
+
+    /*
+    constructor to create the multidual number from the primal part and array of dual parts, automatically inferring the domain.
+    */
+    proc init(val : real, grad : [?dom]) {
+      this.dom = dom; 
+      this.prim = val;
+      this.dual = grad;
+    }
+  }
+
+  /* Converts a pair of real numbers to dual number */
+  proc todual(val : real, der : real) {return new DualNumber(val, der);}
+
+  /* Converts a real number and array of reals to a multidual number. */
+  proc todual(val : real, grad : [?D] real) {return new MultiDual(val, grad);}
+
+  proc isDualType(type t : DualNumber) param {return true;}
+
+  proc isDualType(type t : MultiDual) param {return true;}
+
+  proc isDualType(type t) param {return false;}
+
+  proc isEitherDualNumberType(type t, type s) param {return isDualType(t) || isDualType(s);}
+
+  /*
+  For dual numbers, it returns the primal part. For real numbers, it returns the number itself.
+  */
+  proc prim(a) where isDualType(a.type) {return a.prim;}
+
+  /*
+  For dual numbers, it returns the dual part, for real numbers it returns zero.
+  */
+  proc dual(a) where isDualType(a.type) {return a.dual;}
+
+  pragma "no doc"
+  proc prim(a : real) {return a;}
+
+  pragma "no doc"
+  proc dual(a : real) {return 0.0;}
+
+  proc isclose(a, b, rtol=1e-5, atol=0.0) where isEitherDualNumberType(a.type, b.type) {
+      return isclose(prim(a), prim(b), rtol=rtol, atol=atol) && isclose(dual(a), dual(b), rtol=rtol, atol=atol);
+  }
+

src/hyperbolic.chpl

@@ -1,15 +1,39 @@
 module hyperbolic {
     use ForwardModeAD;
 
-    proc sinh(a : DualNumber) {return new DualNumber(sinh(prim(a)), dual(a) * cosh(prim(a)));}
+    proc sinh(a) where isDualType(a.type) {
+        var f = sinh(prim(a)), 
+        df = dual(a) * cosh(prim(a));
+        return todual(f, df);
+    }
 
-    proc cosh(a : DualNumber) {return new DualNumber(cosh(prim(a)), dual(a) * sinh(prim(a)));}
+    proc cosh(a) where isDualType(a.type) {
+        var f = cosh(prim(a)), 
+            df = dual(a) * sinh(prim(a));
+        return todual(f, df);
+    }
 
-    proc tanh(a : DualNumber) {return new DualNumber(tanh(prim(a)), dual(a) * (1 - tanh(prim(a))**2));}
+    proc tanh(a) where isDualType(a.type) {
+        var f = tanh(prim(a)), 
+            df = dual(a) * (1 - tanh(prim(a))**2);
+        return todual(f, df);
+    }
 
-    proc asinh(a : DualNumber) {return new DualNumber(asinh(prim(a)), dual(a) / sqrt(prim(a)**2 + 1));}
+    proc asinh(a) where isDualType(a.type) {
+        var f = asinh(prim(a)), 
+            df = dual(a) / sqrt(prim(a)**2 + 1);
+        return todual(f, df);
+    }
 
-    proc acosh(a : DualNumber) {return new DualNumber(acosh(prim(a)), dual(a) / sqrt(prim(a)**2 - 1));}
+    proc acosh(a) where isDualType(a.type) {
+        var f = acosh(prim(a)), 
+            df = dual(a) / sqrt(prim(a)**2 - 1);
+        return todual(f, df);
+    }
 
-    proc atanh(a : DualNumber) {return new DualNumber(atanh(prim(a)), dual(a) / (1 - prim(a) ** 2));}
+    proc atanh(a) where isDualType(a.type) {
+        var f = atanh(prim(a)), 
+            df = dual(a) / (1 - prim(a) ** 2);
+        return todual(f, df);
+    }
 }

src/transcendental.chpl

@@ -1,28 +1,58 @@
 module transcendental {
     use ForwardModeAD;
 
-    operator **(a : real, b : DualNumber) {
-        return new DualNumber(a ** prim(b), log(a) * (a ** prim(b)) * dual(b));
+    operator **(a : real, b) where isDualType(b.type) {
+        var f = a ** prim(b),
+            df = log(a) * (a ** prim(b)) * dual(b);
+        return todual(f, df);
     }
 
 
-    operator **(a : DualNumber, b : DualNumber) {
+    operator **(a, b) where isDualType(a.type) && a.type == b.type {
         var f = prim(a) ** prim(b);
         var df = f * (dual(b) * log(prim(a)) + prim(b) * dual(a) / prim(a)); 
-        return new DualNumber(f, df);
+        return todual(f, df);
     }
 
-    proc exp(a : DualNumber) { return new DualNumber(exp(prim(a)), dual(a) * exp(prim(a))); }
+    proc exp(a) where isDualType(a.type) {
+        var f = exp(prim(a)),
+            df = dual(a) * exp(prim(a));
+        return todual(f, df);
+    }
 
-    proc exp2(a : DualNumber) {return new DualNumber(exp2(prim(a)), ln_2 * exp2(prim(a)) * dual(a));}
+    proc exp2(a) where isDualType(a.type) {
+        var f = exp2(prim(a)), 
+            df = ln_2 * exp2(prim(a)) * dual(a);
+        return todual(f, df);
+    }
 
-    proc expm1(a : DualNumber) {return new DualNumber(expm1(prim(a)), exp(prim(a)) * dual(a));}
+    proc expm1(a) where isDualType(a.type) {
+        var f = expm1(prim(a)), 
+            df = exp(prim(a)) * dual(a);
+        return todual(f, df);
+    }
 
-    proc log(a : DualNumber) { return new DualNumber(log(prim(a)), dual(a) / prim(a)); }
+    proc log(a) where isDualType(a.type) {
+        var f = log(prim(a)), 
+            df = dual(a) / prim(a);
+        return todual(f, df);
+    }
 
-    proc log2(a : DualNumber) {return new DualNumber(log2(prim(a)), dual(a)/(prim(a) * ln_2));}
+    proc log2(a) where isDualType(a.type) {
+        var f = log2(prim(a)), 
+            df = dual(a)/(prim(a) * ln_2);
+        return todual(f, df);
+    }
 
-    proc log10(a : DualNumber) {return new DualNumber(log10(prim(a)), dual(a) / (prim(a) * ln_10));}
+    proc log10(a) where isDualType(a.type) {
+        var f = log10(prim(a)), 
+            df = dual(a) / (prim(a) * ln_10);
+        return todual(f, df);
+    }
 
-    proc log1p(a : DualNumber) {return new DualNumber(log1p(prim(a)), dual(a) / (prim(a) + 1));}
+    proc log1p(a) where isDualType(a.type) {
+        var f = log1p(prim(a)), 
+            df = dual(a) / (prim(a) + 1);
+        return todual(f, df);
+    }
 }