[Apologies for cross-posting]
For the algorithmic differentiation: I have started an implementation
(attached) of the jet machinery using a domain rather than the domain+package
used in Smith et al.
This provides considerably faster differentiation than symbolic
differentiation. For example:
g(x) == (1+x)^10000
eval(D(g(x),x), x=0.1) -- slow
g(jet(0.1)) -- fast!
I was not certain how to get x=jet(1.0) or eval(x^x, x=jet(2.0)) to work: any
suggestions?
Sincerely, Mark,
On 02/19/2017 12:21 PM, Gabriel Dos Reis wrote:
It was available in a branch of OpenAxiom. Jacob wanted to rewrite for merging
in the main trunk. He graduated and got a job, so did not get to rewrite it. I
added the AST part of the work to the trunk, but the jet machinery remained in
Jacob's branch.
Jacob: do you still have that?
-- Gaby
On Feb 19, 2017 3:15 AM, "Mark Clements"
<mark.cleme...@ki.se<mailto:mark.cleme...@ki.se>> wrote:
Is the Axiom code for Smith's work on algorithmic differentiation
available?
http://www.axiomatics.org/~gdr/ad/issac07.pdf<http://www.axiomatics.org/%7Egdr/ad/issac07.pdf>
http://oaktrust.library.tamu.edu/bitstream/handle/1969.1/ETD-TAMU-2010-05-7823/SMITH-DISSERTATION.pdf
Sincerely, Mark Clements.
------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, SlashDot.org! <http://sdm.link/slashdot>
http://sdm.link/slashdot
_______________________________________________
open-axiom-devel mailing list
open-axiom-devel@lists.sourceforge.net<mailto:open-axiom-devel@lists.sourceforge.net>
https://lists.sourceforge.net/lists/listinfo/open-axiom-devel
)abbrev domain JET Jet
Jet(T: Ring): Public == Private where
Public ==> CoercibleTo(OutputForm) with
jet: (T, T) -> %
++ construct a jet value from a pair
jet: (T) -> %
++ construct a jet value with delta=1
Zero: () -> %
++ zero value
One: () -> %
++ one value
elt: (%, "value") -> T
++ retrieve the value field of a jet
elt: (%, "delta") -> T
++ retrieve the delta field of a jet
setelt: (%, "value", T) -> T
++ set the value field of a jet
setelt: (%, "delta", T) -> T
++ set the delta field of a jet
"+": (%, %) -> %
"+": (T, %) -> %
"+": (%, T) -> %
"-": (%, %) -> %
"-": (T, %) -> %
"-": (%, T) -> %
"*": (%, %) -> %
"*": (T, %) -> %
"*": (%, T) -> %
"-": % -> %
if T has OrderedRing then
abs: % -> %
if T has PartialOrder then
"<": (%, %) -> Boolean
">": (%, %) -> Boolean
"<=": (%, %) -> Boolean
">=": (%, %) -> Boolean
"=": (%, %) -> Boolean
if T has Field then
"/": (%, %) -> %
"/": (T, %) -> %
"/": (%, T) -> %
if T has TranscendentalFunctionCategory and T has Field then
log: % -> %
if T has TranscendentalFunctionCategory and T has Field and T has
RadicalCategory then
asin: % -> %
acos: % -> %
atan: % -> %
if T has TranscendentalFunctionCategory then
exp: % -> %
"^": (%, %) -> %
"^": (%, T) -> %
"^": (T, %) -> %
sin: % -> %
cos: % -> %
tan: % -> %
sinh: % -> %
cosh: % -> %
tanh: % -> %
if T has RadicalCategory and T has Field then
sqrt: % -> %
Private ==> add
x, y : %
z : T
Rep := Record(val: T, der: T)
jet(v,d) == [v, d]
jet(v) == [v, 1$T]
Zero() == jet(0$T, 0$T)
One() == jet(1$T, 1$T)
elt(jet: %, key: "value") == jet.val
elt(jet: %, key: "delta") == jet.der
setelt(jet: %, key: "value", jetValue: T) == jet.val := jetValue
setelt(jet: %, key: "delta", jetDelta: T) == jet.der := jetDelta
coerce(x:%):OutputForm ==
hconcat(('jet)::OutputForm, paren [ (x.val)::OutputForm,
(x.der)::OutputForm ])
x + y == jet(x.value + y.value, x.delta + y.delta)
x + z == jet(x.value + z, x.delta)
z + x == jet(z + x.value, x.delta)
x - y == jet(x.value - y.value, x.delta - y.delta)
x - z == jet(x.value - z, x.delta)
z - x == jet(z - x.value, -x.delta)
x * y == jet(x.value * y.value, x.delta*y.value + x.value*y.delta)
x * z == jet(x.value * z, x.delta*z)
z * x == jet(z * x.value, z*x.delta)
-x == jet(-x.value, -x.delta)
if T has OrderedRing then
abs x == jet(abs(x.value), sign(x.value)*x.delta)
if T has PartialOrder then
x < y == x.value < y.value
x > y == x.value > y.value
x <= y == x.value <= y.value
x >= y == x.value >= y.value
x = y == x.value = y.value
if T has Field then
x / y ==
r : T := x.value / y.value
jet(r, (x.delta - r * y.delta) / y.value)
x / z ==
jet(x.value / z, x.delta / z)
z / x ==
jet(z / x.value, -z*x.delta / x.value / x.value)
if T has TranscendentalFunctionCategory and T has Field then
log x == jet(log(x.value), x.delta/x.value)
if T has TranscendentalFunctionCategory and T has Field and T has
RadicalCategory then
asin x == jet(asin(x.value), x.delta/sqrt(1-x.value*x.value))
acos x == jet(acos(x.value), -x.delta/sqrt(1-x.value*x.value))
atan x == jet(atan(x.value), x.delta/(1+x.value*x.value))
if T has TranscendentalFunctionCategory then
exp x == jet(exp(x.value), x.delta*exp(x.value))
x ^ z == jet(x.value ^ z, x.delta*z*x.value^(z-1))
x ^ y == jet(x.value ^ y.value, log(x.value)*x.value^y.value*y.delta +
x.delta*y.value*x.value^(y.value-1))
z ^ x == jet(z ^ x.value, log(z)*z^x.value*x.delta)
sin x == jet(sin(x.value), x.delta*cos(x.value))
cos x == jet(cos(x.value), -x.delta*sin(x.value))
tan x == jet(tan(x.value), x.delta*(1+tan(x.value)*tan(x.value)))
sinh x == jet(sinh(x.value), x.delta*cosh(x.value))
cosh x == jet(cosh(x.value), x.delta*sinh(x.value))
tanh x == jet(tanh(x.value), x.delta*(1-tanh(x.value)*tanh(x.value)))
if T has RadicalCategory and T has Field then
sqrt x == jet(sqrt(x.value), x.delta/sqrt(x.value)/(1$T+1$T))
)if false
)cd /home/marcle/Home/axiom
)co Jet2.spad
JF ==> Jet Float
-- constructors
jet(1.0,2.0)
jet(1.0)
Zero()$Jet Float
One()$Jet Float
-- check some operators
g := operator 'g
h := operator 'h
Dgx(f) == (f,D(f(g(x)),x))
[Dgx(log), Dgx(exp), Dgx(-), Dgx(h), Dgx(sin), Dgx(cos), Dgx(tan), Dgx(asin),
Dgx(acos), Dgx(atan), Dgx(sinh), Dgx(cosh), Dgx(tanh), Dgx(sqrt)]
Dgx2(f) == (f,D(f(g(x),h(x)),x))
[Dgx2(+),Dgx2(-),Dgx2(*),Dgx2(/),Dgx2(^)]
[f::(Float->Float) for f in [+,-,*,/,^]]
-- testing
testDeriv(f,value) == eval(D(f(),X),x=value) - f(jet(value)).delta
testDeriv(log,2.0)
testDeriv(exp,2.0)
testDeriv(-,2.0)
testDeriv(x +-> exp(exp(x)),2.0)
testDeriv(sin,2.0)
testDeriv(cos,2.0)
testDeriv(tan,2.0)
testDeriv(sinh,2.0)
testDeriv(cosh,2.0)
testDeriv(tanh,2.0)
testDeriv(asin,0.1)
testDeriv(acos,0.1)
testDeriv(atan,0.1)
testDeriv(x +-> x^(2*x),2.0)
testDeriv(sqrt,2.0)
testDeriv(x +-> exp(x)/sqrt(x),2.0)
testDeriv(x +-> exp(x)*sqrt(x),2.0)
testDeriv(x +-> x*x,2.0)
testDeriv(x +-> x+x,2.0)
testDeriv(x +-> 10*x,2.0)
jet(2.0,1.0)<=jet(2.0,1.0)
jet(1,2) <= jet(1,3)
g(x) == (1+x)^10000
eval(D(g(x),x), x=0.1) -- slow symbolic differentiation
g(jet(0.1)) -- fast algorithmic differentiation!
g(x) == x^exp(x)^exp(x)
eval(D(g(x),x),x=2.0)
g(jet(2.0))
fx(x, n) ==
if n=0 then return 0
else
a := 0
b := x
(a,b) := (0, x)
for k in 1..(n-1) repeat
h := b
b := log(a+b)
a := h
b
fx(jet(2.0,1.0), 7)
eval(D(fx(x,7),x), x=2.0)
-- partial likelihood
-- Assume distinct ordered times
event := [true,false,false,true,true,false]
times := [i for i in 1..#event]
x := [1,0,1,1,0,1]
eta := [beta*xi for xi in x]
index := [i for i in 1..#event for e in event | e]
L := reduce(*,[exp(eta(i))/reduce(+,[exp(eta(j)) for j in i..#event]) for i in
index])
l := reduce(+,[eta(i) - log(reduce(+,[exp(eta(j)) for j in i..#event])) for i
in index])
-- eval(l,beta=jet(1.0)) -- fails. How do I get beta=jet(1.0) to work?
-- time-varying effect
f(beta,gamma) == reduce(+,[x.i*beta+x.i*times.i*gamma -
log(reduce(+,[exp(x.j*beta+x.j*times.j*gamma) for j in i..#event])) for i in
index])
f(x,y)
f(jet(1.0))
f(jet(1.0),3.0)
)endif------------------------------------------------------------------------------
Announcing the Oxford Dictionaries API! The API offers world-renowned
dictionary content that is easy and intuitive to access. Sign up for an
account today to start using our lexical data to power your apps and
projects. Get started today and enter our developer competition.
http://sdm.link/oxford
_______________________________________________
open-axiom-devel mailing list
open-axiom-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/open-axiom-devel
