/**
         
 * Warranty & Liability
         
 * To the extent permitted by applicable law and unless explicitly
         
 * otherwise agreed upon, XLOG Technologies AG makes no warranties
         
 * regarding the provided information. XLOG Technologies AG assumes
         
 * no liability that any problems might be solved with the information
         
 * provided by XLOG Technologies AG.
         
 *
         
 * Rights & License
         
 * All industrial property rights regarding the information - copyright
         
 * and patent rights in particular - are the sole property of XLOG
         
 * Technologies AG. If the company was not the originator of some
         
 * excerpts, XLOG Technologies AG has at least obtained the right to
         
 * reproduce, change and translate the information.
         
 *
         
 * Reproduction is restricted to the whole unaltered document. Reproduction
         
 * of the information is only allowed for non-commercial uses. Selling,
         
 * giving away or letting of the execution of the library is prohibited.
         
 * The library can be distributed as part of your applications and libraries
         
 * for execution provided this comment remains unchanged.
         
 *
         
 * Restrictions
         
 * Only to be distributed with programs that add significant and primary
         
 * functionality to the library. Not to be distributed with additional
         
 * software intended to replace any components of the library.
         
 *
         
 * Trademarks
         
 * Jekejeke is a registered trademark of XLOG Technologies AG.
         
 */
         

         
:- ensure_loaded(library(lists)).
         

         
/**
         
 * Vector is represented as [A1,..,An]:
         
 * - Where A1 .. An are numbers.
         
 * Matrix is represeted as [V1,..,Vm]:
         
 * - Where V1 .. Vm are vectors.
         
 */
         

         
/*******************************************************************/
         
/* Vector Operations                                               */
         
/*******************************************************************/
         

         
/**
         
 * vecdot(V, W, X):
         
 * The predicate succeeds in X with the inner product V . W.
         
 */
         
% vecdot(+Vector, +Vector, -Number)
         
vecdot(V, W, X) :-
         
   foldl(sys_muladd, V, W, 0, X).
         

         
% sys_muladd(+Number, +Number, +Number, -Number)
         
sys_muladd(X, Y, Z, T) :- T is X*Y+Z.
         

         
/**
         
 * veccross(V, W, M):
         
 * The predicate succeeds in M with the outer product V (x) W.
         
 */
         
% veccross(+Vector, +Vector, -Matrix):
         
veccross(V, W, M) :-
         
   maplist(sys_vecscale(W), V, M).
         

         
% sys_vecscale(+Vector, +Number, -Vector)
         
sys_vecscale(V, K, W) :-
         
   maplist(sys_mul(K), V, W).
         

         
% sys_mul(+Number, +Number, -Number)
         
sys_mul(X, Y, Z) :- Z is X*Y.
         

         
/**
         
 * vecrand(N, V):
         
 * The predicate succeeds in V with a random vector of size N.
         
 */
         
% vecrand(+Integer, -Vector)
         
vecrand(N, L) :-
         
   length(L, N),
         
   maplist(sys_rand, L).
         

         
% sys_rand(-Number)
         
sys_rand(Y) :- random(X), Y is 2*X-1.
         

         
/*******************************************************************/
         
/* Matrix Operations                                               */
         
/*******************************************************************/
         

         
/**
         
 * matapply(M, V, W):
         
 * The predicate succeeeds in W with the linear mapping M * V.
         
 */
         
% matapply(+Matrix, +Vector, -Vector)
         
matapply(M, V, W) :-
         
   maplist(vecdot(V), M, W).
         

         
/**
         
 * mattran(M, N):
         
 * The predicate succeeds in N with the transpose T(M).
         
 */
         
% mattran(+Matrix, -Matrix)
         
mattran([H|T], R) :- sys_mattran(T, H, R).
         

         
% sys_mattran(+Matrix, +Vector, -Matrix)
         
sys_mattran([], U, R) :- maplist(sys_rdot([]), U, R).
         
sys_mattran([J|T], H, R) :- sys_mattran(T, J, TC), maplist(sys_rdot, TC, H, R).
         

         
% sys_rdot(+Term, +Term, -Term)
         
sys_rdot(T, H, [H|T]).
         

         
/**
         
 * matmul(M, N, R):
         
 * The predicate succeeds in R with the matrix multiplication M * N.
         
 */
         
% matmul(+Matrix, +Matrix, -Matrix)
         
matmul(M, N, R) :-
         
   mattran(N, H),
         
   maplist(matapply(H), M, R).
         

         
/**
         
 * matrand(R, C, M):
         
 * The predicate succeeds in M with a random matrix of size R x C.
         
 */
         
% matrand(+Integer, +Integer, -List)
         
matrand(R, C, M) :-
         
   length(M, C),
         
   maplist(vecrand(R), M).