Deeply embedded Automatic Differentiation in Lambda Prolog

diff.mod

module diff.

%% First attempt:
% Automatic differentiation, easy (i.e. non-composite) case.
% Lambda Prolog supports pattern-matching on lambdas,
% so we can use it to define "deep-embedded" autodiffs:

% We can exploit the higher-order expressiveness of Lambda Prolog
% to express that "F is a constant function".
ediff F (X : real\ 0.0) :-  pi Y\ pi X\ F X = F 0.0, !.
ediff (X\ X)     (X : real\ 1.0) :- !.
ediff (X\ ~ F X) (X\ ~ G X) :- ediff F G, !.
ediff sin cos         :- !.
ediff cos (X\~ sin X) :- !.
ediff (X\ F X + G X) (X\ F' X + G' X) :-
  !, ediff F F', ediff G G'.
ediff (X\ F X - G X) (X\ F' X - G' X) :-
  !, ediff F F', ediff G G'.
ediff (X\ F X * G X) (X\ F' X * G X + F X * G' X) :-
  !, ediff F F', ediff G G'.

ediff sqrt (X\ 1.0 / (2.0 * sqrt X)) :- !.
ediff (X\ F X / G X) (X\ (F' X / G X) - (F X * G' X / (G X * G X))) :-
  !, ediff F F', ediff G G'.

% If we add the following rule, computation explodes.
% ediff (X\F (G X)) (X \ F' (G X) * G' X) :-
%   ediff F F', ediff G G'.

%% Second attempt:
% Automatic differentiation, which supports composition.
% We declare each rule in a composite form to avoid explosion.
adiff (X\ X) (X\ 1.0) :- !.
adiff F (X : real\ 0.0) :-
  pi X\ F X = F 0.0, !.
adiff (X\ A * F X) (X\ A * F' X) :- !, adiff F F'.
adiff (X\ ~ (F X)) (X\ ~ (F' X)) :- !, adiff F F'.
adiff (X\ sin (F X)) (X\ cos (F X) * F' X) :- !, adiff F F'.
adiff (X\ cos (F X)) (X\ ~ sin (F X) * F' X) :- !, adiff F F'.
adiff (X\ F X + G X) (X\ F' X + G' X) :-
  !, adiff F F', adiff G G'.
adiff (X\ F X - G X) (X\ F' X - G' X) :-
  !, adiff F F', adiff G G'.
adiff (X\ F X * G X) (X\ F' X * G X + F X * G' X) :-
  !, adiff F F', adiff G G'.
adiff (X\ F X / G X) (X\ (F' X / G X) - (F X * G' X / (G X * G X))) :-
  !, adiff F F', adiff G G'.
adiff (X\ sqrt (F X)) (X\ F' X / (2.0 * sqrt (F X))) :-
  !, adiff F F'.
  
ndiff 0 F F  :- !.
ndiff N F Fn :-
  N > 0, M is N - 1, ndiff M F F', adiff F' Fn.


%% Taylor series.
taylor X0 F S :- maclaurin (X\ F (X - X0)) S.

%% Maclaurin series.
maclaurin F S :-
  machelper F S0,
  facts Fs,
  zipped macTerm Fs S0 S.

% local macTerm  int -> real -> real -> o.
macTerm N X Y :- Z is int_to_real N, Y is (X / Z).

% local machelper (real -> real) -> stream real -> o.
machelper F S :-
  adiff F F', Y is F 0.0, scons Y (machelper F') S.

% local facts  stream A -> o.
facts (1 # value S) :- facHelper 1 1 S.

facHelper N K Ns :-
  M is N * K, N' is N + 1, scons M (facHelper N' M) Ns.

%% Helper functions.
scons A P (A # thunk P).

force (value X) X.
force (thunk P) X :- P X.

take 0 _ [] :- !.
take N (X # P) (X :: Xs) :-
  N > 0, force P T, M is N - 1, take M T Xs.

zipped P (X # S) (Y # T) U :-
  P X Y Z,
  force S S', force T T',
  scons Z (zipped P S' T') U.

diff.sig

sig diff.

%% Automatic differentiation
% Deeply embedded automatic differentiation, without support for composition.
exportdef ediff ((real -> real) -> (real -> real) -> o).

% Deeply embedded automatic differentiation, with composition.
exportdef adiff ((real -> real) -> (real -> real) -> o).

exportdef lim-0     (real -> real) -> real -> o.

%% Taylor expansion around the given point.
exportdef taylor     real -> (real -> real) -> stream real -> o.
exportdef maclaurin  (real -> real) -> stream real -> o.

%% Inifnite stream type.
kind stream type -> type.
type empty  (stream A).
type #      A -> cell (stream A) -> stream A.
infixr # 90.

% Smart constructor
type scons A -> (stream A -> o) -> stream A -> o.

% Returns the list cosisting of initial n values of the given stream
exportdef take     int -> stream A -> list A -> o.

% Drops initial n elements of the given stream.
exportdef drop     int -> stream A -> stream A -> o.

% Maps stream.
exportdef zipped (A -> B -> C -> o) -> stream A -> stream B -> stream C -> o.

% Cell containing a (possibly delayed) value.
kind cell  type -> type.
type thunk (A -> o) -> cell A.
type value A        -> cell A.

% Force cell to the value.
exportdef force  (cell A -> A -> o).

exportdef facts  stream int -> o.
exportdef machelper (real -> real) -> stream real -> o.
exportdef macTerm  int -> real -> real -> o.
exportdef facHelper int -> int -> stream int -> o.

log.sh

$ tjcc diff && tjlink diff
$ tjsim diff
Compiling: diff
Linking...
Welcome to Teyjus
Copyright (C) 2008 A. Gacek, S. Holte, G. Nadathur, X. Qi, Z. Snow
Teyjus comes with ABSOLUTELY NO WARRANTY
This is free software, and you are welcome to redistribute it
under certain conditions. Please view the accompanying file
COPYING for more information
[diff] ?- adiff (X\ X * sin X) E.

The answer substitution:
E = W1\ 1.000000 * sin W1 + W1 * (cos W1 * 1.000000)

More solutions (y/n)? y

no (more) solutions

[diff] ?- ndiff 2 (X\ X * sin X) E.

The answer substitution:
E = W1\ 1.000000 * (cos W1 * 1.000000) + (1.000000 * (cos W1 * 1.000000) + W1 * (- sin W1 * 1.000000 * 1.000000 + cos W1 * 0.000000))

More solutions (y/n)? n

yes

[diff] ?- maclaurin (X\ X * sin X) S, take 5 S M.

The answer substitution:
M = 0.000000 :: 0.000000 :: 1.000000 :: 0.000000 :: -0.166667 :: nil
S = 0.000000 # thunk (zipped macTerm (1 # thunk (facHelper 2 1)) (0.000000 # thunk (machelper (W1\ 1.000000 * (cos W1 * 1.000000) + (1.000000 * (cos W1 * 1.000000) + W1 * (- sin W1 * 1.000000 * 1.000000 + cos W1 * 0.000000))))))

More solutions (y/n)? y

no (more) solutions

[diff] ?- taylor 1.0 (X\ X * sin (1.0 + X)) S, take 5 S M.

The answer substitution:
M = -0.000000 :: -1.000000 :: 1.000000 :: 0.166667 :: -0.166667 :: nil
S = -0.000000 # thunk (zipped macTerm (1 # thunk (facHelper 2 1)) (-1.000000 # thunk (machelper (W1\ (1.000000 - 0.000000) * (cos (1.000000 + (W1 - 1.000000)) * (0.000000 + (1.000000 - 0.000000))) + ((1.000000 - 0.000000) * (cos (1.000000 + (W1 - 1.000000)) * (0.000000 + (1.000000 - 0.000000))) + (W1 - 1.000000) * (- sin (1.000000 + (W1 - 1.000000)) * (0.000000 + (1.000000 - 0.000000)) * (0.000000 + (1.000000 - 0.000000)) + cos (1.000000 + (W1 - 1.000000)) * 0.000000))))))

More solutions (y/n)? y

no (more) solutions

[diff] ?- taylor 0.5 sin E, take 5 E S.

The answer substitution:
S = -0.479426 :: 0.877583 :: 0.239713 :: -0.146264 :: -0.019976 :: nil
E = -0.479426 # thunk (zipped macTerm (1 # thunk (facHelper 2 1)) (0.877583 # thunk (machelper (W1\ - sin (W1 - 0.500000) * (1.000000 - 0.000000) * (1.000000 - 0.000000) + cos (W1 - 0.500000) * 0.000000))))

More solutions (y/n)? y

no (more) solutions