diff --git a/Classes/DerivativeNumber.sc b/Classes/DerivativeNumber.sc index 73dabba1b7230a0bd754b4eddd24637596fa64cc..3f80baf98b4b2edca68aea5c3d2a9b6183f9ebaa 100644 --- a/Classes/DerivativeNumber.sc +++ b/Classes/DerivativeNumber.sc @@ -1,6 +1,6 @@ /* -Automatic differentiation, following +Automatic differentiation (forward mode), following http://conal.net/blog/posts/beautiful-differentiation http://conal.net/blog/posts/what-is-automatic-differentiation-and-why-does-it-work @@ -9,12 +9,14 @@ https://dl.acm.org/doi/10.1145/289423.289442 Barak A. Pearlmutter and Jeffrey Mark Siskind 2007 – Lazy Multivariate Higher-Order Forward-Mode AD - The differentiationStep is held for numerical differentiation, when no derivative is known. +This is based on the concept of dual numbers, due to William Clifford (1873) + terminology: -real, derivative -(same as in in Python, xad-py) +real, dual (the real part, the dual part) + +a subclass DerivativeNumber will implement higher derivatives not all operators are implemented. @@ -25,31 +27,34 @@ only those need to be lifted. But in sclang, we have many primitives that need t */ -DerivativeNumber : Number { +Dual : Number { - var <>real, <>derivative; + var <>real, <>dual; var <>differentiationStep = 1e-6; // for numerical approximation. This may need adjustment - *new { |real, derivative = 0| - ^super.newCopyArgs(real, derivative) + *new { |real, dual = 0| + ^super.newCopyArgs(real, dual) } + derivative { ^this.dual } + + // double dispatch performBinaryOpOnSimpleNumber { | aSelector, aNumber, adverb | - ^aNumber.asDerivative.perform(aSelector, this, adverb) + ^aNumber.asDual.perform(aSelector, this, adverb) } performBinaryOpOnSignal { | aSelector, aSignal, adverb | - ^aSignal.asDerivative.perform(aSelector, this) + ^aSignal.asDual.perform(aSelector, this) } performBinaryOpOnComplex { | aSelector, aComplex, adverb | - ^aComplex.asDerivative.perform(aSelector, this) + ^aComplex.asDual.perform(aSelector, this) } performBinaryOpOnUGen { | aSelector, aUGen, adverb | - ^aUGen.asDerivative.perform(aSelector, this, adverb) + ^aUGen.asDual.perform(aSelector, this, adverb) } performBinaryOpOnSomething { |aSelector, thing, adverb| - ^thing.asDerivative.perform(aSelector, this, adverb) + ^thing.asDual.perform(aSelector, this, adverb) } performBinaryOpOnSeqColl { | aSelector, thing, adverb | ^this.performBinaryOpOnSomething(aSelector, thing, adverb) @@ -58,15 +63,15 @@ DerivativeNumber : Number { // equality < { |other| - ^this.real < other.asDerivative.real + ^this.real < other.real } == { |other| - ^this.real == other.real and: { this.derivative == other.derivative } + ^this.real == other.real and: { this.dual == other.dual } } hash { - ^this.instVarHash(['real', 'derivative']) + ^this.instVarHash(['real', 'dual']) } // unary ops @@ -97,7 +102,7 @@ DerivativeNumber : Number { cubed { ^this.chainRule(cubed(_), { |x| x.squared * 3 }) } abs { ^this * this.sign } - // unary ops with derivative zero + // unary ops with dual zero ceil { ^this.class.new(this.real.ceil, 0.0) } floor { ^this.class.new(this.real.floor, 0.0) } @@ -106,7 +111,7 @@ DerivativeNumber : Number { sqrt { var r = this.real.sqrt; - ^this.class.new(r, this.derivative / (2 * r)) + ^this.class.new(r, this.dual / (2 * r)) } @@ -115,36 +120,31 @@ DerivativeNumber : Number { // using the simplified formula directly instead of the chain rule + { |other| - other = other.asDerivative; - ^this.class.new(this.real + other.real, this.derivative + other.derivative) + ^this.class.new(this.real + other.real, this.dual + other.dual) } - { |other| - other = other.asDerivative; - ^this.class.new(this.real - other.real, this.derivative - other.derivative) + ^this.class.new(this.real - other.real, this.dual - other.dual) } // Leibniz rule (ab)' = (a'b + b'a) * { |other| - other = other.asDerivative; - ^this.class.new(this.real * other.real, (other.derivative * this.real) + (this.derivative * other.real)) + ^this.class.new(this.real * other.real, (other.dual * this.real) + (this.dual * other.real)) } / { |other| - other = other.asDerivative; - ^this.class.new(this.real / other.real, ((other.real * this.derivative) - (this.real * other.derivative)) / other.real.squared) + ^this.class.new(this.real / other.real, ((other.real * this.dual) - (this.real * other.dual)) / other.real.squared) } pow { |other| var d1, d2; - other = other.asDerivative; // this uses the n-argument chain rule with partial derivatives over pow(_, y) and pow(x, _) - d1 = other.real * pow(this.real, other.real - 1) * this.derivative; - d2 = pow(this.real, other.real) * log(this.real) * other.derivative; + d1 = other.real * pow(this.real, other.real - 1) * this.dual; + d2 = pow(this.real, other.real) * log(this.real) * other.dual; ^this.class.new(pow(this.real, other.real), d1 + d2) } - mod { |other| ^this.class.new(mod(this.real, other.asDerivative.real), this.derivative) } + mod { |other| ^this.class.new(mod(this.real, other.real), this.dual) } max { |other| ^0.5 * (this + other + abs(this - other)) } min { |other| ^0.5 * (this + other - abs(this - other)) } amclip { |other| ^this * max(other, 0) } @@ -154,49 +154,47 @@ DerivativeNumber : Number { // binary ops with derivative zero - div { |other| ^this.class.new(div(this.real, other.asDerivative.real), 0) } - lcm { |other| ^this.class.new(lcm(this.real, other.asDerivative.real), 0) } - gcd { |other| ^this.class.new(gcd(this.real, other.asDerivative.real), 0) } - round { |other=1| ^this.class.new(round(this.real, other.asDerivative.real), 0) } - roundUp { |other=1| ^this.class.new(roundUp(this.real, other.asDerivative.real), 0) } - trunc { |other=1| ^this.class.new(trunc(this.real, other.asDerivative.real), 0) } + div { |other| ^this.class.new(div(this.real, other.real), 0) } + lcm { |other| ^this.class.new(lcm(this.real, other.real), 0) } + gcd { |other| ^this.class.new(gcd(this.real, other.real), 0) } + round { |other=1| ^this.class.new(round(this.real, other.real), 0) } + roundUp { |other=1| ^this.class.new(roundUp(this.real, other.real), 0) } + trunc { |other=1| ^this.class.new(trunc(this.real, other.real), 0) } hypot { |other| // (a * a - b * b).sqrt ^sqrt(this.squared + other.squared) } ring1 { |other| // (a * b) + a - ^(this * other.asDerivative) + this + ^(this * other) + this } ring2 { |other| // ((a*b) + a + b) - other = other.asDerivative; ^(this * other) + this + other } ring3 { |other| // (a * a * b) - ^this.squared * other.asDerivative + ^this.squared * other } ring4 { |other| // ((a*a *b) - (a*b*b)) - other = other.asDerivative; ^(this.squared * other) - (other.squared * this) } difsqr { |other| // (a*a) - (b*b) - ^this.squared - other.asDerivative.squared + ^this.squared - other.squared } sumsqr { |other| // (a*a) + (b*b) - ^this.squared + other.asDerivative.squared + ^this.squared + other.squared } sqrdif { |other| // (a - b) ** 2 - ^(this - other.asDerivative).squared + ^(this - other.asDual).squared } sqrsum { |other| // (a + b) ** 2 - ^(this + other.asDerivative).squared + ^(this + other).squared } absdif { |other| // (a - b).abs - ^abs(this - other.asDerivative) + ^abs(this - other) } moddif { |other| - other = other.asDerivative; - ^this.class.new(moddif(this.real - other.real), this.derivative - other.derivative) + other = other; + ^this.class.new(moddif(this.real - other.real), this.dual - other.dual) } excess { |other=1| ^this - clip2(this, other) @@ -278,16 +276,15 @@ DerivativeNumber : Number { chainRule { |func, derivativeFunc| ^this.class.new( func.(this.real), - derivativeFunc.(this.real) * this.derivative + derivativeFunc.(this.real) * this.dual ) } chainRuleN { |func, partialDerivativeFuncs, operands| - var reals, derivatives, r, d; + var reals, duals, r, d; - operands = operands.collect(_.asDerivative); reals = operands.collect(_.real); - derivatives = operands.collect(_.derivative); + duals = operands.collect(_.dual); if(operands.size != partialDerivativeFuncs.size) { Error("the number of partial derivative functions must match the number of operands.\n" @@ -296,7 +293,7 @@ DerivativeNumber : Number { r = func.valueArray(reals); d = operands.sum { |x, i| - derivatives[i] * partialDerivativeFuncs[i].valueArray(reals) + duals[i] * partialDerivativeFuncs[i].valueArray(reals) }; ^this.class.new(r, d) @@ -306,9 +303,8 @@ DerivativeNumber : Number { *approximate { |func, args| var y0, y1, delta, reals; - args = args.collect(_.asDerivative); reals = args.collect(_.real); - delta = args[0].differentiationStep; + delta = args[0].asDual.differentiationStep; // this is very rough still: we just move by one delta in all arguments y0 = func.valueArray(reals); y1 = func.valueArray(reals + delta); @@ -318,7 +314,7 @@ DerivativeNumber : Number { // conversions isValidUGenInput { ^this.real.isValidUGenInput } - isNumber { ^false } // not really a number, even though we subclass Number + isNumber { ^false } // like Complex, not a normal number, even though we subclass Number asInteger { ^this.real.asInteger } asFloat { ^this.real.asFloat } @@ -326,27 +322,72 @@ DerivativeNumber : Number { asPolar { ^Polar.new(this, this.class.new(0, 0)) } asPoint { ^Point.new(this, this.class.new(0, 0)) } - asDerivative { ^this } + asDual { ^this } + asDerivative { ^DerivativeNumber(this.real, this.dual) } + asUGenInput { ^this.real.asUGenInput } printOn { |stream| - stream << this.class.name << "(" <<* [real, derivative] << ")" + stream << this.class.name << "(" <<* [real, dual] << ")" + } + + +} + + +// DerivativeNumber will serve as a class that adds higher derivatives + +DerivativeNumber : Dual { + + // overriding the optimizations of Dual with the chain rule, + // so that higher derivatives can be implemented here + + * { |other| ^this.chainRuleN(_ * _, [other.real, this.real], [this, other]) } + - { |other| ^this.chainRuleN(_ - _, [1, -1], [this, other]) } + + { |other| ^this.chainRuleN(_ + _, [1, 1], [this, other]) } + + mod { |other| ^this.chainRuleN(mod(_, _), [1, 0], [this, other]) } + div { |other| ^this.chainRuleN(div(_, _), [0, 0], [this, other]) } + + // unary ops + sqrt { ^this.chainRule(sqrt(_), { |x| reciprocal(sqrt(x) * 2) }) } + + + asDual { + ^Dual(this.real, this.dual) } + asDerivative { ^this } + } + Number { + dual { ^0 } + derivative { ^0 } + + asDual { + ^Dual(this, 0) + } + asDerivative { ^DerivativeNumber(this, 0) } - } + SequenceableCollection { + + dual { ^0 } + derivative { ^0 } + + + asDual { + ^Dual(this, 0) + } + asDerivative { ^DerivativeNumber(this, 0) } @@ -355,26 +396,25 @@ DerivativeNumber : Number { + UGen { + dual { ^Slope.perform(this.methodSelectorForRate, this) } + derivative { ^this.dual } + + + asDual { + ^Dual(this, this.dual) + } + asDerivative { - ^DerivativeNumber(this, Slope.perform(this.methodSelectorForRate, this)) + ^DerivativeNumber(this, this.dual) } } -/* -for a future nth deriviative support, we conveniently put everything through the chainRule and then just iterate it. - - sqrtChain { ^this.chainRule(sqrt(_), { |x| reciprocal(sqrt(x) * 2) }) } - mulChain { |other| ^this.chainRuleN(_ * _, [other.real, this.real], [this, other]) } - diffChain { |other| ^this.chainRuleN(_ - _, [1, -1], [this, other]) } - addChain { |other| ^this.chainRuleN(_ + _, [1, 1], [this, other]) } - - modChain { |other| ^this.chainRuleN(mod(_, _), [1, 0], [this, other]) } - divChain { |other| ^this.chainRuleN(div(_, _), [0, 0], [this, other]) } +/* +// possible fallback when we have no definition: - // fallback when we have no definition. doesNotUnderstand { |selector ... args| if(this.real.respondsTo(\selector).not) { DoesNotUnderstandError(this.real, selector, args) diff --git a/Classes/extFunction.sc b/Classes/extFunction.sc index 325d8a49f0d5838979cdcf1d98c8093e3f350978..7723a1745baf7db8cbe3c25721b4ca81eb9fd56c 100644 --- a/Classes/extFunction.sc +++ b/Classes/extFunction.sc @@ -26,8 +26,8 @@ var defaults = this.defaultArgs; var argNames = this.argNames; args = argNames.collect { |name, i| args[i] ?? { defaults[i] } }; - args[argIndex] = DerivativeNumber(args[argIndex], 1); - this.valueArray(args).asDerivative.derivative + args[argIndex] = Dual(args[argIndex], 1); + this.valueArray(args).asDual.dual } } diff --git a/HelpSource/DerivativeNumber.schelp b/HelpSource/DerivativeNumber.schelp index baef36aae0ac950ebae46b28f84952b9133307b7..a779ef381c43443a47cefd8b9f0c7e4404fc2851 100644 --- a/HelpSource/DerivativeNumber.schelp +++ b/HelpSource/DerivativeNumber.schelp @@ -1,12 +1,13 @@ -TITLE:: DerivativeNumber -SUMMARY::Automatic Differentiation +TITLE:: Dual +SUMMARY::Automatic Differentiation (forward mode) CATEGORIES::Math +RELATED::Complex DESCRIPTION:: -DerivativeNumber implements forward automatic differentiation. +Because automatic differentiation coincides with the algebra of the so called dual numbers (found by William Clifford in 1873), we can make use of this general concept here. Dual numbers in the form math::a + b\epsilon:: resemble complex numbers math::a + bi::, but while math::i^2 = -1::, the second component in dual numbers vanishes: math::\epsilon^2 = 0::. This "magically" leads to an algebra of derivatives. -It follows Jerzy Karczmarczuk 1998, "Functional differentiation of computer programs". +The implementation idea follows Jerzy Karczmarczuk 1998, "Functional differentiation of computer programs". In order to find the derivative of a function, there are three approaches: @@ -19,20 +20,21 @@ LIST:: note::This is a reasonably well-tested, but preliminary version. Names (like instance variables or argument names) and functionality may change.:: -Calling operations on derivatives apply the chain rule and keep track of resulting reals and derivatives: +Calling operations on derivatives apply the chain rule and keep track of resulting reals and duals: CODE:: -a = DerivativeNumber(4, 1); +a = Dual(4, 1); b = a.squared.cos; b.real; // cos(4 * 4) -b.derivative // neg(sin(4 * 4)) * (2 * 4) +b.dual // neg(sin(4 * 4)) * (2 * 4) +b.derivative; // same, this is a synonym :: This is valid for Arrays: CODE:: -a = DerivativeNumber((-100..100)/100 * pi, 1); +a = Dual((-100..100)/100 * pi, 1); b = a.squared.cos; -[b.real, b.derivative].plot(separately:true); +[b.real, b.dual].plot(separately:true); :: It is also valid for UGens: @@ -40,9 +42,9 @@ CODE:: ( { var a, b; - a = DerivativeNumber(Line.ar(-pi, pi, 0.01), 1); + a = Dual(Line.ar(-pi, pi, 0.01), 1); b = a.squared.cos; - [b.real, b.derivative] + [b.real, b.dual] }.plot ) :: @@ -57,14 +59,14 @@ ARGUMENT:: real The real part of the number, which may be any number class or a UGen. -ARGUMENT:: derivative +ARGUMENT:: dual -The derivative part of the number, which may be any number class or a UGen. +The dual (derivative) part of the number, which may be any number class or a UGen. In algebra, this is called math::\epsilon:: (epsilon). CODE:: -a = DerivativeNumber(2.5, -3.652); +a = Dual(2.5, -3.652); a.real; -a.derivative; +a.dual; :: @@ -75,9 +77,9 @@ Numerical approximation may fail to converge well in some cases, so you can't re CODE:: d = (0..100).normalize(-1, 1); // approximate a range of numbers -a = DerivativeNumber(d, 1); -b = DerivativeNumber.approximate({ |x, y, z| x.cylBesselJ(y, z) }, [a, 2, 2]); -[b.real, b.derivative].plot +a = Dual(d, 1); +b = Dual.approximate({ |x, y, z| x.cylBesselJ(y, z) }, [a, 2, 2]); +[b.real, b.dual].plot :: @@ -87,15 +89,15 @@ METHOD::real Return the real part of the number, which carries the result of any past operations. CODE:: -a = DerivativeNumber(2, 0.6).squared.real; // 4 +a = Dual(2, 0.6).squared.real; // 4 :: -METHOD::derivative -Return the derivative part of the number, which carries the derivative that follows from all past operations. +METHOD::dual +Return the dual part of the number, which carries the derivative that follows from all past operations. CODE:: -a = DerivativeNumber(2, 0.6).squared.derivative; // 2.4, which is (2 * 2 * 0.6) +a = Dual(2, 0.6).squared.dual; // 2.4, which is (2 * 2 * 0.6) :: METHOD::< @@ -114,49 +116,51 @@ CODE:: ( var steps = [0.0001, 0.001, 0.01, 0.1, 0.5]; d = (0..100).normalize(-1, 1); // approximate a range of numbers -a = steps.collect { |each| DerivativeNumber(d, 1).differentiationStep_(each) }; -c = a.collect { |each| DerivativeNumber.approximate({ |x, y, z| x.cylBesselJ(y, z) }, [each, 0.1, 2]) }; -c.collect { |each| each.derivative }.plot.superpose_(true); +a = steps.collect { |each| Dual(d, 1).differentiationStep_(each) }; +c = a.collect { |each| Dual.approximate({ |x, y, z| x.cylBesselJ(y, z) }, [each, 0.1, 2]) }; +c.collect { |each| each.dual }.plot.superpose_(true); ) :: SUBSECTION::Conversions +note::Mixed operations are not yet fully implemented and may therefore not be commutative.:: + METHOD::asDerivative -Convert receiver to a code::DerivativeNumber::. In this case, it is already one, so return code::this::. -For other classes, it returns code::DerivativeNumber(receiver, 0):: +Convert receiver to a code::Dual::. In this case, it is already one, so return code::this::. +For other classes, it returns code::Dual(receiver, 0):: CODE:: // applied in mixed expressions -5 * DerivativeNumber(2, 0.3); // DerivativeNumber(10, 1.5) -[1, 3, 5] * DerivativeNumber(2, 0.3); // DerivativeNumber([2, 6, 10], [0.3, 0.9, 1.5]) +5 * Dual(2, 0.3); // Dual(10, 1.5) +[1, 3, 5] * Dual(2, 0.3); // Dual([2, 6, 10], [0.3, 0.9, 1.5]) :: METHOD::asComplex -Returns a link::Classes/Complex:: with the real part as the receiver, and the imaginary part as code::DerivativeNumber(0, 0)::. +Returns a link::Classes/Complex:: with the real part as the receiver, and the imaginary part as code::Dual(0, 0)::. CODE:: // applied in mixed expressions -DerivativeNumber(2, 0.3) * Complex(5, 2); // DerivativeNumber(Complex(10.0, 4.0), Complex(1.5, 0.6)) -DerivativeNumber(2, 0.3) * Complex([1, 3, 5], 2); // DerivativeNumber(Complex([2.0, 6.0, 10.0], [4.0, 4.0, 4.0]), Complex([0.3, 0.9, 1.5], [0.6, 0.6, 0.6])) +Dual(2, 0.3) * Complex(5, 2); // Dual(Complex(10.0, 4.0), Complex(1.5, 0.6)) +Dual(2, 0.3) * Complex([1, 3, 5], 2); // Dual(Complex([2.0, 6.0, 10.0], [4.0, 4.0, 4.0]), Complex([0.3, 0.9, 1.5], [0.6, 0.6, 0.6])) :: METHOD::asPolar -Returns a link::Classes/Complex:: with code::rho:: (radius) as the receiver and code::theta:: as code::DerivativeNumber(0, 0)::. +Returns a link::Classes/Polar:: with code::rho:: (radius) as the receiver and code::theta:: as code::Dual(0, 0)::. CODE:: -DerivativeNumber(2, 0.3).asPolar; // Polar(DerivativeNumber(2, 0.3), DerivativeNumber(0, 0)) +Dual(2, 0.3).asPolar; // Polar(Dual(2, 0.3), Dual(0, 0)) // applied in mixed expressions, it returns a Complex (default in sclang) -Polar(5, 0) * DerivativeNumber(2, 0.3); // DerivativeNumber(Complex(10.0, 0.0), Complex(1.5, 0.0)) +Polar(5, 0) * Dual(2, 0.3); // Dual(Complex(10.0, 0.0), Complex(1.5, 0.0)) :: METHOD::asPoint -Returns a link::Classes/Point:: with the x coordinate as the receiver and the y coordinate as code::DerivativeNumber(0, 0)::. +Returns a link::Classes/Point:: with the x coordinate as the receiver and the y coordinate as code::Dual(0, 0)::. CODE:: // applied in mixed expressions -Point(5, 9) * DerivativeNumber(2, 0.3); // Point(DerivativeNumber(10, 1.5), DerivativeNumber(0, 0)) +Point(5, 9) * Dual(2, 0.3); // Point(Dual(10, 1.5), Dual(0, 0)) :: @@ -174,37 +178,37 @@ f.(d).plot; // now, we want to automatically find the derivative // for probing, we use a seed with derivative 1.0 -a = DerivativeNumber(2.0, 1.0); +a = Dual(2.0, 1.0); // at 2.0, the real part is 2^3 = 8.0 and the derivative is 3 * (2.squared) = 12: b = cubed(a); // test the difference: zero -b.derivative - (3 * a.real.squared) +b.dual - (3 * a.real.squared) -// the real and derivatives can be arrays: +// the real and duals can be arrays: d = (-100..100)/100 * 6pi; -a = DerivativeNumber(d, 1); +a = Dual(d, 1); // evaluate the function above for each value of the domain d -// this results in a DerivativeNumber of arrays +// this results in a Dual of arrays b = f.(a); // plot it -[b.real, b.derivative].plot +[b.real, b.dual].plot // the same with UGens ... ( { - a = DerivativeNumber(Line.ar(-6pi, 6pi, 0.01), 1); + a = Dual(Line.ar(-6pi, 6pi, 0.01), 1); b = f.(a); - [b.real, b.derivative] + [b.real, b.dual] }.plot; ) // ... and arrays of UGens ( { - a = DerivativeNumber(Line.ar(-6pi, 6pi, 0.01) * (1..3), 1); + a = Dual(Line.ar(-6pi, 6pi, 0.01) * (1..3), 1); b = f.(a); - [b.real, b.derivative].lace + [b.real, b.dual].lace }.plot; ) @@ -216,18 +220,18 @@ You need the link::Classes/UnitTest:: class to run these. CODE:: ( -UnitTest.runTest("TestDerivativeNumber:test_unaryOps_invariants"); -UnitTest.runTest("TestDerivativeNumber:test_binaryOps_invariants"); -UnitTest.runTest("TestDerivativeNumber:test_mappingOps_invariants"); - -UnitTest.runTest("TestDerivativeNumber:test_unaryOps_sanityCheckWithNumerical"); -UnitTest.runTest("TestDerivativeNumber:test_binaryOps_sanityCheckWithNumerical"); -UnitTest.runTest("TestDerivativeNumber:test_mappingOps_sanityCheckWithNumerical"); - -UnitTest.runTest("TestDerivativeNumber:test_chainRuleN_withPow"); -UnitTest.runTest("TestDerivativeNumber:test_parametricCurve_1"); -UnitTest.runTest("TestDerivativeNumber:test_parametricSurface_1"); -UnitTest.runTest("TestDerivativeNumber:test_partialDerivatives"); -UnitTest.runTest("TestDerivativeNumber:test_specialTrigonometric_invariants"); +UnitTest.runTest("TestDual:test_unaryOps_invariants"); +UnitTest.runTest("TestDual:test_binaryOps_invariants"); +UnitTest.runTest("TestDual:test_mappingOps_invariants"); + +UnitTest.runTest("TestDual:test_unaryOps_sanityCheckWithNumerical"); +UnitTest.runTest("TestDual:test_binaryOps_sanityCheckWithNumerical"); +UnitTest.runTest("TestDual:test_mappingOps_sanityCheckWithNumerical"); + +UnitTest.runTest("TestDual:test_chainRuleN_withPow"); +UnitTest.runTest("TestDual:test_parametricCurve_1"); +UnitTest.runTest("TestDual:test_parametricSurface_1"); +UnitTest.runTest("TestDual:test_partialDerivatives"); +UnitTest.runTest("TestDual:test_specialTrigonometric_invariants"); ) :: \ No newline at end of file diff --git a/UnitTests/TestAutomaticDifferentiation.sc b/UnitTests/TestAutomaticDifferentiation.sc index cc16906ad599129af66a5d0aa0ae66db76f4311a..b80fad0cb3c4a8a770fff6d3e67f5c954b6d890f 100644 --- a/UnitTests/TestAutomaticDifferentiation.sc +++ b/UnitTests/TestAutomaticDifferentiation.sc @@ -1,7 +1,7 @@ -TestDerivativeNumber : UnitTest { +TestDual : UnitTest { classvar unaryOps = #[ \sin, @@ -104,7 +104,7 @@ TestDerivativeNumber : UnitTest { var d = number.asDerivative; var a1 = perform(d, op).real; var a2 = perform(d.real, op); - this.assertFloatEquals(a1, a2, "Unary op '%' should have the same result on a derivative number as on a number".format(op)); + this.assertFloatEquals(a1, a2, "Unary op '%' should have the same result on a dual number as on a number".format(op)); } } @@ -115,7 +115,7 @@ TestDerivativeNumber : UnitTest { var d2 = n2.asDerivative; var a1 = perform(d1, op, d2).real; var a2 = perform(n1, op, n2); - this.assertFloatEquals(a1, a2, "Binary op '%' should have the same result on a derivative number as on a number".format(op)); + this.assertFloatEquals(a1, a2, "Binary op '%' should have the same result on a dual number as on a number".format(op)); } } @@ -125,7 +125,7 @@ TestDerivativeNumber : UnitTest { mappingOps.do { |op| var a1 = ds[0].perform(op, *ds[1..]).real; var a2 = ns[0].perform(op, *ns[1..]); - this.assertFloatEquals(a1, a2, "Mapping operation '%' should have the same result on a derivative number as on a number.\nArgs:%".format(op, ns)); + this.assertFloatEquals(a1, a2, "Mapping operation '%' should have the same result on a dual number as on a number.\nArgs:%".format(op, ns)); } } @@ -135,37 +135,37 @@ TestDerivativeNumber : UnitTest { var d = rrand(range[0], range[1]).asDerivative; var a1 = perform(d, op).real; var a2 = perform(d.real, op); - this.assertFloatEquals(a1, a2, "Unary op '%' should have the same result on a derivative number as on a number".format(op)); + this.assertFloatEquals(a1, a2, "Unary op '%' should have the same result on a dual number as on a number".format(op)); } } test_unaryOps_sanityCheckWithNumerical { - var d = DerivativeNumber(0.3, 1); + var d = Dual(0.3, 1); unaryOps.do { |op| - var a1 = DerivativeNumber.approximate({ |num| perform(num, op) }, [d]).derivative; - var a2 = perform(d, op).derivative; + var a1 = Dual.approximate({ |num| perform(num, op) }, [d]).dual; + var a2 = perform(d, op).dual; this.assertFloatEquals(a1, a2, "Numerical differentiation should approximate automatic differentiation:\n%(%)".format(op, d), within: 0.001); // low precision is ok, this is a sanity check } } test_binaryOps_sanityCheckWithNumerical { - var d1 = DerivativeNumber(0.3, 1); - var d2 = DerivativeNumber(0.8, 1); + var d1 = Dual(0.3, 1); + var d2 = Dual(0.8, 1); binaryOps.do { |op| - var a1 = DerivativeNumber.approximate({ |n1, n2| perform(n1, op, n2) }, [d1, d2]).derivative; - var a2 = perform(d1, op, d2).derivative; + var a1 = Dual.approximate({ |n1, n2| perform(n1, op, n2) }, [d1, d2]).dual; + var a2 = perform(d1, op, d2).dual; this.assertFloatEquals(a1, a2, "Numerical differentiation should approximate automatic differentiation:\n%(%, %)".format(op, d1, d2), within: 0.001); // low precision is ok, this is a sanity check } } test_mappingOps_sanityCheckWithNumerical { var ns = [rrand(-2.0, 2.0), 0.2, 0.8, 0.1, 2.3]; - var ds = ns.collect { |real| DerivativeNumber(real, 1) }; + var ds = ns.collect { |real| Dual(real, 1) }; mappingOps.do { |op| var call = { |...args| args[0].performList(op, args[1..]) }; - var a1 = DerivativeNumber.approximate(call, ns).derivative; - var a2 = call.valueArray(ds).derivative; + var a1 = Dual.approximate(call, ns).dual; + var a2 = call.valueArray(ds).dual; this.assertFloatEquals(a1, a2, "Numerical differentiation should approximate automatic differentiation:\n%(%)".format(op, ds), within: 0.001); // low precision is ok, this is a sanity check @@ -195,9 +195,9 @@ TestDerivativeNumber : UnitTest { var domain = (0..100).normalize(-3, 3); - var derivativeDomain = DerivativeNumber(domain, 1); + var dualDomain = Dual(domain, 1); var codomain = curve.(domain); - var automaticDerivative = curve.(derivativeDomain).collect { |x| x.derivative }; + var automaticDerivative = curve.(dualDomain).collect { |x| x.dual }; var symbolicDerivative = derivative.(domain); var error = absdif(automaticDerivative, symbolicDerivative); var fine = error.flat.every(_ < 0.000001); @@ -238,15 +238,15 @@ TestDerivativeNumber : UnitTest { var codomain = surface.(domain1, domain2); var symbolicDerivative = derivatives.(domain1, domain2); - var derivativeDomain1 = DerivativeNumber(domain1, 1); - var derivativeDomain2 = DerivativeNumber(domain2, 1); + var dualDomain1 = Dual(domain1, 1); + var dualDomain2 = Dual(domain2, 1); var automaticDerivativeNumbers = [ - surface.(derivativeDomain1, domain2).collect(_.asDerivative), - surface.(domain1, derivativeDomain2).collect(_.asDerivative), + surface.(dualDomain1, domain2).collect(_.asDerivative), + surface.(domain1, dualDomain2).collect(_.asDerivative), ]; - var automaticDerivatives = automaticDerivativeNumbers.collect { |d| d.collect(_.derivative) }; - var automaticReals = automaticDerivativeNumbers[0].collect(_.real); // use the first partial derivative (u) + var automaticDerivatives = automaticDerivativeNumbers.collect { |d| d.collect(_.dual) }; + var automaticReals = automaticDerivativeNumbers[0].collect(_.real); // use the first partial dual (u) var error1 = absdif(automaticDerivatives, symbolicDerivative); @@ -255,7 +255,7 @@ TestDerivativeNumber : UnitTest { var error2 = absdif(automaticReals, codomain); var fine2 = error2.flat.every(_ < 0.000001); - this.assert(fine1, "parametric surface 1: automatic and symbolic differentiation should match in derivative part" ); + this.assert(fine1, "parametric surface 1: automatic and symbolic differentiation should match in dual part" ); if(fine2.not) { automaticReals.flop.plot; @@ -266,7 +266,7 @@ TestDerivativeNumber : UnitTest { } test_chainRuleN_withPow { - var ds = [DerivativeNumber(1.0, 3.2), DerivativeNumber(3.0, 2.2), DerivativeNumber(0.3, -0.2)]; // no negative real, because -1.2 pow: 0.3 is nan + var ds = [Dual(1.0, 3.2), Dual(3.0, 2.2), Dual(0.3, -0.2)]; // no negative real, because -1.2 pow: 0.3 is nan var pairs = ds.dup.allTuples; pairs.do { |pair| var chained = pair[0].chainRuleN( @@ -278,7 +278,7 @@ TestDerivativeNumber : UnitTest { pair ); var direct = pow(pair[0], pair[1]); - this.assertFloatEquals(direct.derivative, chained.derivative, "n-arg chain rule derivatives should match to formula for pow:\n%\n%".format(*pair)); + this.assertFloatEquals(direct.dual, chained.dual, "n-arg chain rule derivatives should match to formula for pow:\n%\n%".format(*pair)); this.assertFloatEquals(direct.real, chained.real, "n-arg chain rule reals should match to formula for pow:\n%\n%".format(*pair)); }; diff --git a/examples/DerivativeNumber.scd b/examples/DerivativeNumber.scd deleted file mode 100644 index 676a26312fe650a6e6799240448107370c25a51b..0000000000000000000000000000000000000000 --- a/examples/DerivativeNumber.scd +++ /dev/null @@ -1,357 +0,0 @@ - - - -( -a = DerivativeNumber(1, 1); -b = DerivativeNumber(0, -2); -c = DerivativeNumber(-2, -1); -) - -// commutativity -(a * b) == (b * a) -(a + b) == (b + a) - -// basic ops -a * b -b / a -a / a -a / a -a - b - -// multichannel expansion -a * [1, 2, 3] + [10, 20] - -// mixing classes -a * Complex(8, 1) -k = Complex(8, 1) * a + b -k.squared -Complex(8, 1).squared - - -k = a + [1, 2]; -k.max(0.5) - - -// Not yet there: these should be all zero -( -a = DerivativeNumber(1.2, 1); -b = DerivativeNumber(2, 1); -[ - a.sqrt - a.sqrtChain, - a * b - a.mulChain(b), - a - b - a.diffChain(b), - a + b - a.addChain(b), - a mod: b - a.modChain(b), - a div: b - a.divChain(b), - a.pow(b) - a.powChain(b), - (a / b) - (a * b.reciprocal) -].selectIndices { |x| x != DerivativeNumber(0.0, 0) } // should be empty -) - -a/b -b/a - -// currently, no fallback for methods that are not implemented. -// one could do it easily. - -DerivativeNumber(1, 0.2).degreeToKey([0, 2, 4, 5, 7, 9], 12) -DerivativeNumber(1, 0.2).partition(2, 1) // this won't do it. -a = DerivativeNumber(1, 1).seriesIter(1.5, 4); -a.domain.nextN(8) -a.derivative.next // nil, well.. - - -(0..DerivativeNumber(1, 0.2)) // this fails. - - -// UGens -{ (DerivativeNumber(1, 15) + Line.ar(0, 1, 0.01)).derivative }.plot -{ (DerivativeNumber(1, SinOsc.ar(100) * 20) + Line.ar(0, 1, 0.01)).derivative }.plot - - - -// mappings -DerivativeNumber(0.5).linlin(0, 0.3, 0, 1) -DerivativeNumber(0.5).linexp(0.01, 1, 0.01, 1) - - -expexp(-0.7942271232605, 0.2, 0.8, 0.1, 2.3) -expexp(-0.7942271232605.asDerivative, 0.2, 0.8, 0.1, 2.3) -curvelin(1.9459619522095, 0.2, 0.8, 0.1, 2.3) -curvelin(1.9459619522095.asDerivative, 0.2, 0.8, 0.1, 2.3) - - -// we could add: -// D : DerivativeNumber {} - - -/* - -sanity checks - -*/ - -UnitTest.gui; - -DerivativeNumber.methods.do { |x| if(x.argNames.size == 1) { "%: %(_),".format(x.name, x.name).postln } }; - - -( -UnitTest.runTest("TestDerivativeNumber:test_unaryOps_invariants"); -UnitTest.runTest("TestDerivativeNumber:test_binaryOps_invariants"); -UnitTest.runTest("TestDerivativeNumber:test_mappingOps_invariants"); - -UnitTest.runTest("TestDerivativeNumber:test_unaryOps_sanityCheckWithNumerical"); -UnitTest.runTest("TestDerivativeNumber:test_binaryOps_sanityCheckWithNumerical"); -UnitTest.runTest("TestDerivativeNumber:test_mappingOps_sanityCheckWithNumerical"); - -UnitTest.runTest("TestDerivativeNumber:test_chainRuleN_withPow"); -UnitTest.runTest("TestDerivativeNumber:test_parametricCurve_1"); -UnitTest.runTest("TestDerivativeNumber:test_parametricSurface_1"); -UnitTest.runTest("TestDerivativeNumber:test_partialDerivatives"); -UnitTest.runTest("TestDerivativeNumber:test_specialTrigonometric_invariants"); -) - -/* - -Various checks - -Many of what follows exists as a unit test, too. - -*/ - -( -var r1 = 0.4; -var r2 = 0.2; -var m = 3; -var n = 5; -~curve = { |u| - [ - (r1 + (r2 * cos(n * u))) * cos(m * u), - (r1 + (r2 * cos(n * u))) * sin(m * u), - r2 * sin(n * u) - ] -}; - -~derivative = { |u| - [ - (r2 * n * sin(n * u) * cos(m * u).neg) + (r2 * m * cos(n * u) * sin(m * u).neg) - (r1 * m * sin(m * u)), - (r2 * n * sin(n * u) * sin(m * u).neg) + (r2 * m * cos(n * u) * cos(m * u)) + (r1 * m * cos(m * u)), - r2 * n * cos(n * u) - ] -}; -) - - -~domain = (0..100).normalize(-3, 3); -~derivativeDomain = DerivativeNumber(~domain, 1); -~codomain = ~curve.(~domain); -~automaticDerivative =~curve.(~derivativeDomain).collect { |x| x.derivative }; -~symbolicDerivative = ~derivative.(~domain); - -~symbolicDerivative.plot -~automaticDerivative.plot - - -~automaticDerivative - ~symbolicDerivative - - -( -var r1 = 0.4; -var r2 = 0.2; -~surface = { |u, v| - [ - (r1 + (r2 * cos(u))) * cos(v), - (r1 + (r2 * cos(u))) * sin(v), - r2 * sin(u) - ] - }; - -~derivatives = { |u, v| - [ - [ - r2 * sin(u) * cos(v).neg, - r2 * sin(u) * sin(v).neg, - r2 * cos(u) - ], - - [ - (r1 + (r2 * cos(u))) * sin(v).neg, - (r1 + (r2 * cos(u))) * cos(v), - 0 - ], - ] - }; -) - - -( -~domain1 = (0..100).normalize(-3, 3); -~domain2 = (0..100).normalize(2, 4); -~codomain = ~surface.(~domain1, ~domain2); -~symbolicDerivative = ~derivatives.(~domain1, ~domain2); - -~derivativeDomain1 = DerivativeNumber(~domain1, 1); -~derivativeDomain2 = DerivativeNumber(~domain2, 1); - -~automaticDerivative = [ - ~surface.(~derivativeDomain1, ~domain2).collect { |x| x.asDerivative.derivative }, - ~surface.(~domain1, ~derivativeDomain2).collect { |x| x.asDerivative.derivative }, -]; -) - - -~automaticDerivative[1] - -~automaticDerivative.shape -~symbolicDerivative.shape - - -~symbolicDerivative[0].plot -~symbolicDerivative[1].plot -~automaticDerivative[0].plot -~automaticDerivative[1].plot // the last one is simply zero, so it won't plot - - -~automaticDerivative[1].flop.flop.plot // this will do - -~symbolicDerivative[1].collect { |x| x.size } - -[[1, 2, 3], [6, 2, 1], 0].plot -[0, [1, 2, 3], [6, 2, 1]].plot -[[0, 4, 3], [1, 2, 3], [[6, 2, 1]]].plot(separately: false) -[[0, 4, 3], [1, 2, 3], [[6, 2, 1]]].plot(separately: true) - -[[1, 2, 3], [4, 5, 6, 7]].plot - -[0, [1, 2, 3], [[6, 2], 1]].plot - - - - -a = { |x| sin(x) }.partialDerivativesAt([0]); -a.value((0..100)/4).plot - -a = { |x, y| sin(x) * sin(y) }.partialDerivativesAt([0]); -a.value((0..100)/4, (0..100)/4).plot - -a = { |x, y| sin(x) * sin(y) }.partialDerivativesAt([0, 1]); -a.value((0..100)/4, (0..100)/2).plot - - - -a = { |x, y| sin(x) * sin(y) }.partialDerivatives; -a.value((0..100)/4, (0..100)/2).plot - -a = { |x, y| sin(x) * sin(y) }.partialDerivatives(withRespectTo: [\y]); -a.value((0..100)/4, (0..100)/2).plot - -a = { |x, y| sin(x) * sin(y) }.partialDerivatives(withRespectTo: [\x, \y]); -a.value((0..100)/4, (0..100)/2).plot - -f = pow(_,_); -a = f.partialDerivatives; -a.((1..100)/30000, (0..100)/300).plot; - - - -a.(100/30000, 100/300) - -g = { |x| x.scurve }.partialDerivativesAt([0]).value((-10..10)); -h = { |x| (6 * x) - (6 * x.squared) }.value((-10..10)); - - -~max = { |a, other| 0.5 * (a + other + abs(a - other)) }; -~max.(-20, -30) -~min = { |a, other| 0.5 * (a + other - abs(a - other)) }; -~min.(-20, -30) -~min.(-2.0, -1.0) - -DerivativeNumber(1.23, 0) min: DerivativeNumber(0, 0) -DerivativeNumber(1.23, 0) min: DerivativeNumber(0, 4) -DerivativeNumber(1.23, 0) max: DerivativeNumber(0, 4) - -~clip2 = { |a, other| max(min(a, other.neg), other) } - -~clip2.(DerivativeNumber(1.23, 0), 0.5) - -DerivativeNumber(1.23, 0) < 1 - -{ a = DerivativeNumber(SinOsc.ar, 0) > 0; a.real }.plot -{ a = DerivativeNumber(SinOsc.ar, 0); (a < 0).real }.plot - -UnitTest.gui; - - - -absdif(DerivativeNumber(0.3, 1), DerivativeNumber(0.8, 1)) - -abs(DerivativeNumber(0.3, 1) - DerivativeNumber(0.8, 1)) -absdif(0.3, 0.8) - absdif(0.3 + 0.001, 0.8 + 0.001) / 0.001 - - -DerivativeNumber(-1.0, 2.2) pow: DerivativeNumber(0.3, -0.2) --1.1 pow: 1.2 - -log(0.2) -lcm(0.4, 7.1) - -DerivativeNumber(0.3, 1).bernouliB2n -0.3.bernouliB2n - -[1, 2, 3].asDerivative.max(2.5) - -f = { |x, y| sin(x) * sin(y) }.partialDerivatives(withRespectTo: [\x]); -f.((0..100).normalize(-1, 1), (0..100).normalize(-1, 2)).plot; - - -// chain rule: -f.g' = f'g * g' -neg.g' = -1 * g' - -DerivativeNumber.approximate({ |x| x.neg }, [DerivativeNumber(0.2, 1)], 0.001) -DerivativeNumber(0.2, 1).chainRule({ |x| neg(x) }, { |x| neg(x) }) - -DerivativeNumber.approximate({ |x| x.cubed }, [DerivativeNumber(0.2, 1)], 0.001) -DerivativeNumber(0.2, 1).chainRule({ |x| cubed(x) }, { |x| 3 * x }) - - -DerivativeNumber.approximate({ |x| x.squared }, [DerivativeNumber(0.2, 1)], 0.001) -DerivativeNumber(0.2, 1).chainRule({ |x| squared(x) }, { |x| 2 * x }) - -DerivativeNumber.approximate({ |x, y| x.roundUp(y) }, [DerivativeNumber(0.2, 1), DerivativeNumber(0.5, 1)], 0.001) -DerivativeNumber(0.2, 1).chainRule({ |x| x.roundUp(0.5) }, 1) - -ring2(DerivativeNumber(0.3, 1), DerivativeNumber(0.8, 1)) - -DerivativeNumber(0.3, 1) * DerivativeNumber(0.8, 1) - -DerivativeNumber.approximate({ |x| x.round(0.5) }, [DerivativeNumber(2.2, 1)], 0.001) -DerivativeNumber(2.2, 1).chainRule({ |x| x.round(0.5) }, 0) -DerivativeNumber(2.2, 1).round(0.5) -DerivativeNumber(2.2, 1).roundUp(0.5) - -DerivativeNumber(0.3, 0.5).linlin(0, 1, 0, 2); -DerivativeNumber(0.3, 0.5).linexp(0, 1, 0.02, 2) - - -// UGens can also be directly converted, a Slope ugen measures the approximate derivative. -// but mind the numerical inaccuracies. - -{ var a = SinOsc.ar(440).asDerivative; [a.real, a.derivative * SampleDur.ir] }.plot(separately: true); -{ var a = SinOsc.kr(40).asDerivative; [a.real, a.derivative * SampleDur.ir] }.plot(0.3, separately: true); - -{ var a = SinOsc.ar(440).asDerivative; a = (a + 2).squared; [a.real, a.derivative * SampleDur.ir] }.plot(separately: true); - -( -{ - var a = Line.ar(0, 8pi, 0.01).asDerivative; - var b = SinOsc.ar(440).asDerivative; - a = (a.sin + 2).squared; - b = (b + 2).squared; - [a.real, a.derivative * SampleDur.ir, a.real, b.derivative] -}.plot(separately: true); -) - - - diff --git a/examples/Dual.scd b/examples/Dual.scd new file mode 100644 index 0000000000000000000000000000000000000000..c7c94de464a2e310e0d7d5b7fcd89ef55985e06b --- /dev/null +++ b/examples/Dual.scd @@ -0,0 +1,338 @@ + + + +( +a = Dual(1, 1); +b = Dual(0, -2); +c = Dual(-2, -1); +) + +// commutativity +(a * b) == (b * a) +(a + b) == (b + a) + +// basic ops +a * b +b / a +a / a +a / a +a - b + +// multichannel expansion +a * [1, 2, 3] + [10, 20] + +// mixing classes +a * Complex(8, 1) +k = Complex(8, 1) * a + b +k.squared +Complex(8, 1).squared + + +k = a + [1, 2]; +k.max(0.5) + + +// currently, no fallback for methods that are not implemented. +// one could do it easily. + +Dual(1, 0.2).degreeToKey([0, 2, 4, 5, 7, 9], 12) +Dual(1, 0.2).partition(2, 1) // this won't do it. +a = Dual(1, 1).seriesIter(1.5, 4); +a.domain.nextN(8) +a.dual.next // nil, well.. + + +(0..Dual(1, 0.2)) // this fails. + + +// UGens +{ (Dual(1, 15) + Line.ar(0, 1, 0.01)).dual }.plot +{ (Dual(1, SinOsc.ar(100) * 20) + Line.ar(0, 1, 0.01)).dual }.plot + + + +// mappings +Dual(0.5).linlin(0, 0.3, 0, 1) +Dual(0.5).linexp(0.01, 1, 0.01, 1) + + +expexp(-0.7942271232605, 0.2, 0.8, 0.1, 2.3) +expexp(-0.7942271232605.asDual, 0.2, 0.8, 0.1, 2.3) +curvelin(1.9459619522095, 0.2, 0.8, 0.1, 2.3) +curvelin(1.9459619522095.asDual, 0.2, 0.8, 0.1, 2.3) + + +// we could add: +// D : Dual {} + + +/* + +sanity checks + +*/ + +UnitTest.gui; + +Dual.methods.do { |x| if(x.argNames.size == 1) { "%: %(_),".format(x.name, x.name).postln } }; + + +( +UnitTest.runTest("TestDual:test_unaryOps_invariants"); +UnitTest.runTest("TestDual:test_binaryOps_invariants"); +UnitTest.runTest("TestDual:test_mappingOps_invariants"); + +UnitTest.runTest("TestDual:test_unaryOps_sanityCheckWithNumerical"); +UnitTest.runTest("TestDual:test_binaryOps_sanityCheckWithNumerical"); +UnitTest.runTest("TestDual:test_mappingOps_sanityCheckWithNumerical"); + +UnitTest.runTest("TestDual:test_chainRuleN_withPow"); +UnitTest.runTest("TestDual:test_parametricCurve_1"); +UnitTest.runTest("TestDual:test_parametricSurface_1"); +UnitTest.runTest("TestDual:test_partialDerivatives"); +UnitTest.runTest("TestDual:test_specialTrigonometric_invariants"); +) + +/* + +Various checks + +Many of what follows exists as a unit test, too. + +*/ + +( +var r1 = 0.4; +var r2 = 0.2; +var m = 3; +var n = 5; +~curve = { |u| + [ + (r1 + (r2 * cos(n * u))) * cos(m * u), + (r1 + (r2 * cos(n * u))) * sin(m * u), + r2 * sin(n * u) + ] +}; + +~derivative = { |u| + [ + (r2 * n * sin(n * u) * cos(m * u).neg) + (r2 * m * cos(n * u) * sin(m * u).neg) - (r1 * m * sin(m * u)), + (r2 * n * sin(n * u) * sin(m * u).neg) + (r2 * m * cos(n * u) * cos(m * u)) + (r1 * m * cos(m * u)), + r2 * n * cos(n * u) + ] +}; +) + + +~domain = (0..100).normalize(-3, 3); +~derivativeDomain = Dual(~domain, 1); +~codomain = ~curve.(~domain); +~automaticDerivative =~curve.(~derivativeDomain).collect { |x| x.dual }; +~symbolicDerivative = ~derivative.(~domain); + +~symbolicDerivative.plot +~automaticDerivative.plot + + +~automaticDerivative - ~symbolicDerivative + + +( +var r1 = 0.4; +var r2 = 0.2; +~surface = { |u, v| + [ + (r1 + (r2 * cos(u))) * cos(v), + (r1 + (r2 * cos(u))) * sin(v), + r2 * sin(u) + ] + }; + +~derivatives = { |u, v| + [ + [ + r2 * sin(u) * cos(v).neg, + r2 * sin(u) * sin(v).neg, + r2 * cos(u) + ], + + [ + (r1 + (r2 * cos(u))) * sin(v).neg, + (r1 + (r2 * cos(u))) * cos(v), + 0 + ], + ] + }; +) + + +( +~domain1 = (0..100).normalize(-3, 3); +~domain2 = (0..100).normalize(2, 4); +~codomain = ~surface.(~domain1, ~domain2); +~symbolicDerivative = ~derivatives.(~domain1, ~domain2); + +~derivativeDomain1 = Dual(~domain1, 1); +~derivativeDomain2 = Dual(~domain2, 1); + +~automaticDerivative = [ + ~surface.(~derivativeDomain1, ~domain2).collect { |x| x.asDual.dual }, + ~surface.(~domain1, ~derivativeDomain2).collect { |x| x.asDual.dual }, +]; +) + + +~automaticDerivative[1] + +~automaticDerivative.shape +~symbolicDerivative.shape + + +~symbolicDerivative[0].plot +~symbolicDerivative[1].plot +~automaticDerivative[0].plot +~automaticDerivative[1].plot // the last one is simply zero, so it won't plot + + +~automaticDerivative[1].flop.flop.plot // this will do + +~symbolicDerivative[1].collect { |x| x.size } + +[[1, 2, 3], [6, 2, 1], 0].plot +[0, [1, 2, 3], [6, 2, 1]].plot +[[0, 4, 3], [1, 2, 3], [[6, 2, 1]]].plot(separately: false) +[[0, 4, 3], [1, 2, 3], [[6, 2, 1]]].plot(separately: true) + +[[1, 2, 3], [4, 5, 6, 7]].plot + +[0, [1, 2, 3], [[6, 2], 1]].plot + + + + +a = { |x| sin(x) }.partialDerivativesAt([0]); +a.value((0..100)/4).plot + +a = { |x, y| sin(x) * sin(y) }.partialDerivativesAt([0]); +a.value((0..100)/4, (0..100)/4).plot + +a = { |x, y| sin(x) * sin(y) }.partialDerivativesAt([0, 1]); +a.value((0..100)/4, (0..100)/2).plot + + + +a = { |x, y| sin(x) * sin(y) }.partialDerivatives; +a.value((0..100)/4, (0..100)/2).plot + +a = { |x, y| sin(x) * sin(y) }.partialDerivatives(withRespectTo: [\y]); +a.value((0..100)/4, (0..100)/2).plot + +a = { |x, y| sin(x) * sin(y) }.partialDerivatives(withRespectTo: [\x, \y]); +a.value((0..100)/4, (0..100)/2).plot + +f = pow(_,_); +a = f.partialDerivatives; +a.((1..100)/30000, (0..100)/300).plot; + + + +a.(100/30000, 100/300) + +g = { |x| x.scurve }.partialDerivativesAt([0]).value((-10..10)); +h = { |x| (6 * x) - (6 * x.squared) }.value((-10..10)); + + +~max = { |a, other| 0.5 * (a + other + abs(a - other)) }; +~max.(-20, -30) +~min = { |a, other| 0.5 * (a + other - abs(a - other)) }; +~min.(-20, -30) +~min.(-2.0, -1.0) + +Dual(1.23, 0) min: Dual(0, 0) +Dual(1.23, 0) min: Dual(0, 4) +Dual(1.23, 0) max: Dual(0, 4) + +~clip2 = { |a, other| max(min(a, other.neg), other) } + +~clip2.(Dual(1.23, 0), 0.5) + +Dual(1.23, 0) < 1 + +{ a = Dual(SinOsc.ar, 0) > 0; a.real }.plot +{ a = Dual(SinOsc.ar, 0); (a < 0).real }.plot + +UnitTest.gui; + + + +absdif(Dual(0.3, 1), Dual(0.8, 1)) + +abs(Dual(0.3, 1) - Dual(0.8, 1)) +absdif(0.3, 0.8) - absdif(0.3 + 0.001, 0.8 + 0.001) / 0.001 + + +Dual(-1.0, 2.2) pow: Dual(0.3, -0.2) +-1.1 pow: 1.2 + +log(0.2) +lcm(0.4, 7.1) + +Dual(0.3, 1).bernouliB2n +0.3.bernouliB2n + +[1, 2, 3].asDual.max(2.5) + +f = { |x, y| sin(x) * sin(y) }.partialDerivatives(withRespectTo: [\x]); +f.((0..100).normalize(-1, 1), (0..100).normalize(-1, 2)).plot; + + +// chain rule: +f.g' = f'g * g' +neg.g' = -1 * g' + +Dual.approximate({ |x| x.neg }, [Dual(0.2, 1)], 0.001) +Dual(0.2, 1).chainRule({ |x| neg(x) }, { |x| neg(x) }) + +Dual.approximate({ |x| x.cubed }, [Dual(0.2, 1)], 0.001) +Dual(0.2, 1).chainRule({ |x| cubed(x) }, { |x| 3 * x }) + + +Dual.approximate({ |x| x.squared }, [Dual(0.2, 1)], 0.001) +Dual(0.2, 1).chainRule({ |x| squared(x) }, { |x| 2 * x }) + +Dual.approximate({ |x, y| x.roundUp(y) }, [Dual(0.2, 1), Dual(0.5, 1)], 0.001) +Dual(0.2, 1).chainRule({ |x| x.roundUp(0.5) }, 1) + +ring2(Dual(0.3, 1), Dual(0.8, 1)) + +Dual(0.3, 1) * Dual(0.8, 1) + +Dual.approximate({ |x| x.round(0.5) }, [Dual(2.2, 1)], 0.001) +Dual(2.2, 1).chainRule({ |x| x.round(0.5) }, 0) +Dual(2.2, 1).round(0.5) +Dual(2.2, 1).roundUp(0.5) + +Dual(0.3, 0.5).linlin(0, 1, 0, 2); +Dual(0.3, 0.5).linexp(0, 1, 0.02, 2) + + +// UGens can also be directly converted, a Slope ugen measures the approximate derivative. +// but mind the numerical inaccuracies. + +{ var a = SinOsc.ar(440).asDual; [a.real, a.dual * SampleDur.ir] }.plot(separately: true); +{ var a = SinOsc.kr(40).asDual; [a.real, a.dual * SampleDur.ir] }.plot(0.3, separately: true); + +{ var a = SinOsc.ar(440).asDual; a = (a + 2).squared; [a.real, a.dual * SampleDur.ir] }.plot(separately: true); + +( +{ + var a = Line.ar(0, 8pi, 0.01).asDual; + var b = SinOsc.ar(440).asDual; + a = (a.sin + 2).squared; + b = (b + 2).squared; + [a.real, a.dual * SampleDur.ir, a.real, b.dual] +}.plot(separately: true); +) + + +