RAI seems to be the closest to what I need to do. It has
a DSL with arrays and matrices, it generates C code,
and it even has automatic differentiation, according
to the docs. It is designed for DSP, but probably
can be extended to non-DSP programming. I should
look at it closer.
For the record: here is what I have got with RAI.
#lang racket/base
(require rai/tools
rai/ai-autodiff
rai/ai-array
rai/ai-array-c
rai/ai-symbolic
rai/stream-syntax
rai/stream-lib
rai/stream-meta)
;; Something defined in Racket with a RAI macro
(define poly
(ai-lambda-meta (x coefs)
((+ ai-add)
(* ai-mul))
(begin
(let ((rcoefs (reverse coefs)))
(for/fold
((acc (car rcoefs)))
((coef (cdr rcoefs)))
(+ coef (* x acc)))))))
(module test rai/stream
;; Here goes RAI code
(provide test-fn test-poly test-dotproduct)
;; function of two variables
(define (test-fn x y)
(/ (sin x) y))
;; polynomial of one variable (Horner's rule)
(define (test-poly x)
(poly x (vector 4 5 6 7)))
;; dot product of two 3-vectors
(define (test-dotproduct xs ys)
(loop (i (n 3))
((accu 0))
((x xs) (y ys))
(+ accu (* x y)))))
(require 'test)
interactions:
;; symbolic derivative of test-fn w.r.t. x
> ((ai-symbolic (ai-deriv test-fn 0)) 'x 'y)
'(/ (cos x) y)
;; symbolic second derivative of test-fn w.r.t. y
> ((ai-symbolic (ai-deriv (ai-deriv test-fn 1) 1)) 'x 'y)
'(- 0 (- 0 (/ (* (sin x) (+ y y)) (* (* y y) (* y y)))))
;; test-fn in dual numbers for automatic differentiation
> ((ai-symbolic (ai-dual test-fn)) (make-dual 'x0 'dx0) (make-dual 'y0
'dy0))
(dual '(p_div (sin x0) y0) '(p_sub (p_div (p_mul (cos x0) dx0) y0)
(p_div (p_mul (sin x0) dy0) (p_mul y0 y0))))
;; C code for the first derivative
> (display (ai-array-c (ai-deriv test-fn)))
[...]
static void proc_loop (
for (int t = 0; t < t_endx; t ++) {
struct proc_si * restrict si = (struct proc_si *)(&state[(t^0)&1]);
struct proc_so * restrict so = (struct proc_so *)(&state[(t^1)&1]);
float32_t r0 = p_copy(t);
float32_t r1 = p_copy(t_endx);
float32_t r3 = p_copy(1.0);
float32_t r4 = p_sub(r3,1.0);
so->r6 = p_sub(si->r5,1.0);
float32_t r7 = p_sin(in->x[t]);
float32_t r8 = p_cos(in->x[t]);
float32_t r9 = p_div(r7,in->y[t]);
float32_t r10 = p_div(r8,in->y[t]);
float32_t r11 = p_mul(in->y[t],in->y[t]);
out->r12[t] = p_copy(r10);
}
;; C code for dual numbers (does not quite work)
> (display (ai-array-c (ai-dual test-fn)))
[...]
struct proc_si * restrict si = (struct proc_si *)(&state[(0)&1]);
struct proc_so * restrict so = (struct proc_so *)(&state[(1)&1]);
for (int t = 0; t < t_endx; t ++) {
struct proc_si * restrict si = (struct proc_si *)(&state[(t^0)&1]);
struct proc_so * restrict so = (struct proc_so *)(&state[(t^1)&1]);
float32_t r0 = p_copy(t);
float32_t r1 = p_copy(t_endx);
float32_t r3 = p_copy(1.0);
float32_t r4 = p_sub(r3,1.0);
so->r6 = p_sub(si->r5,1.0);
float32_t r7 = p_sin(in->x[t]);
float32_t r8 = p_cos(in->x[t]);
float32_t r9 = p_div(r7,in->y[t]);
float32_t r10 = p_mul(in->y[t],in->y[t]);
out->r11[t] = p_copy(#(struct:dual r9 0));
}
;; Symbolic representation of the polynomial
> ((ai-symbolic test-poly) 't)
'(+ 4 (* t (+ 5 (* t (+ 6 (* t 7))))))
;; and its derivative
> ((ai-symbolic (ai-deriv test-poly)) 't)
'(+ (+ 5 (* t (+ 6 (* t 7)))) (* t (+ (+ 6 (* t 7)) (* t 7))))
;; and its derivative in C
> (display (ai-array-c (ai-deriv test-poly)))
[...]
for (int t = 0; t < t_endx; t ++) {
struct proc_si * restrict si = (struct proc_si *)(&state[(t^0)&1]);
struct proc_so * restrict so = (struct proc_so *)(&state[(t^1)&1]);
float32_t r0 = p_copy(t);
float32_t r1 = p_copy(t_endx);
float32_t r3 = p_copy(1.0);
float32_t r4 = p_sub(r3,1.0);
so->r6 = p_sub(si->r5,1.0);
float32_t r7 = p_mul(param->x,7.0);
float32_t r8 = p_add(6.0,r7);
float32_t r9 = p_mul(param->x,r8);
float32_t r10 = p_mul(param->x,7.0);
float32_t r11 = p_add(r8,r10);
float32_t r12 = p_add(5.0,r9);
float32_t r13 = p_mul(param->x,r12);
float32_t r14 = p_mul(param->x,r11);
float32_t r15 = p_add(r12,r14);
float32_t r16 = p_add(4.0,r13);
out->r17[t] = p_copy(r15);
}
;; dot product in C
> (display (ai-array-c test-dotproduct))
[...]
for (int t = 0; t < t_endx; t ++) {
struct proc_si * restrict si = (struct proc_si *)(&state[(t^0)&1]);
struct proc_so * restrict so = (struct proc_so *)(&state[(t^1)&1]);
float32_t r0 = p_copy(t);
float32_t r1 = p_copy(t_endx);
float32_t r3 = p_copy(1.0);
float32_t r4 = p_sub(r3,1.0);
so->r6 = p_sub(si->r5,1.0);
float32_t r8 = p_copy(0.0);
for (int r7 = 0; r7 < 3; r7 ++) {
float32_t r9 = p_copy(in->xs[r7][t]);
float32_t r10 = p_copy(in->ys[r7][t]);
float32_t r11 = p_mul(r9,r10);
float32_t r12 = p_add(r8,r11);
r8 = p_copy(r12);
}
out->r13[t] = p_copy(r8);
}
;; symbolic dot product?
> ((ai-symbolic test-dotproduct) 'x 'y)
ai-symbolic-not-inplemented: reduce/n
;; C dot product with dual numbers?
> (display (ai-array-c (ai-dual test-dotproduct)))
ai-not-defined: for/n
Obviously RAI is currently implemented as only a DSP tool
and is not quite complete, but still promising.
Regards,
Dmitry
--
You received this message because you are subscribed to the Google Groups "Racket
Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to racket-users+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
