1/*	The MIT License (MIT)
    2 *
    3 *	Copyright (c) 2021-25 Rick Workman
    4 *
    5 *	Permission is hereby granted, free of charge, to any person obtaining a copy
    6 *	of this software and associated documentation files (the "Software"), to deal
    7 *	in the Software without restriction, including without limitation the rights
    8 *	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
    9 *	copies of the Software, and to permit persons to whom the Software is
   10 *	furnished to do so, subject to the following conditions:
   11 *
   12 *	The above copyright notice and this permission notice shall be included in all
   13 *	copies or substantial portions of the Software.
   14 *
   15 *	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
   16 *	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   17 *	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
   18 *	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   19 *	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
   20 *	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
   21 *	SOFTWARE.
   22 */
   23:- module(type_ndarray,
   24	[ ndarray/1,
   25	  op(500, yfx, '.+'),      % arithmetic operators on arrays
   26	  op(200, fy,  '.-'),      % same precedence and associativity as standard ops
   27	  op(500, yfx, '.-'),
   28	  op(400, yfx, '.*'),
   29	  op(400, yfx, './') , 
   30	  op(400, yfx, '.@'), 
   31	  op(200, xfx, '.**')  
   32	]).   33
   34%:- reexport(library(type_list),[op(_,_,_)]).  % requires fix circa 9.3.19 % for operators
   35:- reexport(library(type_list),except([slice_parameters/4,index_parameters/3])).  % for operators

arithmetic type ndarray (mulitidimensional array)

This module implements arithmetic type (see module `arithmetic_types) ndarry, immutable, N-dimensional arrays largely patterned after Python's NumPy library. Support from slicing and indexing (0 based) are imported from module type_list. A one dimensional array is represented a compound term with functor # and arguments being the elements of the array. N-dimensional arrays are supported by allowing arrays as array elements. The elements of the array can be any Prolog term, although a subset of the defined arithmetic functions on ndarray's only succeed if the elements are numbers. (There is one exception described below.)

The only exported predicate is the type checking ndarray/1. Also exported are a set of array arithmetic operations: `.+, .-, .*,` etc.; refer ot source for a complete list. The set of arithmetic functions defined by this module include:

:- arithmetic_function(new/2).        % new from shape, uninitialized
:- arithmetic_function(ndarray/1).    % new from ListOfLists
:- arithmetic_function(to_list/1).    % list from ndarray
:- arithmetic_function(ndim/1).       % number of dimensions
:- arithmetic_function(shape/1).      % shape (a list)
:- arithmetic_function(size/1).       % number of values
:- arithmetic_function([]/1).         % block indexing and slicing
:- arithmetic_function([]/2).
:- arithmetic_function(init/2).       % initialize any variables
:- arithmetic_function(\\ /2).        % row concat
:- arithmetic_function(flatten/1).    % flattened array
:- arithmetic_function(reshape/2).    % reshaped array
:- arithmetic_function(transpose/1).  % transpose array
% numeric functions
:- arithmetic_function(arange/2).     % new from range
:- arithmetic_function(arange/3).     % new from range
:- arithmetic_function(arange/4).     % new from range
:- arithmetic_function(identity/2).   % new NxN identity matrix
:- arithmetic_function(sum/1).        % sum of elements
:- arithmetic_function(min/1).        % minimum of elements
:- arithmetic_function(max/1).        % maximum of elements
:- arithmetic_function('.+'/2).       % numeric array addition
:- arithmetic_function('.-'/1).       % numeric array unary minus
:- arithmetic_function('.-'/2).       % numeric array subtraction
:- arithmetic_function('.*'/2).       % numeric array product
:- arithmetic_function('./'/2).       % numeric array division
:- arithmetic_function('.**'/2).      % numeric array power
:- arithmetic_function(apply/2).      % apply function/1 to array
:- arithmetic_function(cross/2).      % cross product of two 3D vectors
:- arithmetic_function(inner/2).      % inner product of two vectors
:- arithmetic_function(outer/2).      % outer product of two vectors
:- arithmetic_function('.@'/2).       % dot product of two vectors/matrices
:- arithmetic_function(dot/2).        % dot product of two vectors/matrices
:- arithmetic_function(determinant/1).  % determinant of square matrix
:- arithmetic_function(inverse/1).    % inverse of square matrix

type_ndarray supports the use of constrained numeric values as defined by library(clpBNR) If current_module(clpBNR) succeeds (i.e., the module has been loaded), instantiation errors from standard functional arithmetic will be interpreted as numeric constraints on the variable(s) in the sub-expression being evaluated. If clpBNR is not accessible, the original instantiation error will be propagated.

See the ReadMe for this pack for more documentation and examples. */

   84:- use_module(library(arithmetic_types)).   85% for indexing and slicing support, arange/n content
   86:- use_module(library(type_list)).   87
   88:- arithmetic_function(new/2).        % new from shape, uninitialized
   89:- arithmetic_function(ndarray/1).    % new from ListOfLists
   90:- arithmetic_function(to_list/1).    % list from ndarray
   91:- arithmetic_function(ndim/1).       % number of dimensions
   92:- arithmetic_function(shape/1).      % shape (a list)
   93:- arithmetic_function(size/1).       % number of values
   94:- arithmetic_function([]/1).         % block indexing and slicing
   95:- arithmetic_function([]/2).   96:- arithmetic_function(init/2).       % initialize any variables
   97:- arithmetic_function(\\ /2).        % row concat
   98%- arithmetic_function(flatten/1).    % flattened array (directive below)
   99:- arithmetic_function(reshape/2).    % reshaped array
  100:- arithmetic_function(transpose/1).  % transpose array
  101% numeric functions
  102:- arithmetic_function(arange/2).     % new from range
  103:- arithmetic_function(arange/3).     % new from range
  104:- arithmetic_function(arange/4).     % new from range
  105:- arithmetic_function(identity/2).   % new NxN identity matrix
  106:- arithmetic_function(sum/1).        % sum of elements
  107:- arithmetic_function(min/1).        % minimum of elements
  108:- arithmetic_function(max/1).        % maximum of elements
  109:- arithmetic_function('.+'/2).       % numeric array addition
  110:- arithmetic_function('.-'/1).       % numeric array unary minus
  111:- arithmetic_function('.-'/2).       % numeric array subtraction
  112:- arithmetic_function('.*'/2).       % numeric array product
  113:- arithmetic_function('./'/2).       % numeric array division
  114:- arithmetic_function('.**'/2).      % numeric array power
  115:- arithmetic_function(apply/2).      % apply function/1 to array
  116:- arithmetic_function(cross/2).      % cross product of two 3D vectors
  117:- arithmetic_function(inner/2).      % inner product of two vectors
  118:- arithmetic_function(outer/2).      % outer product of two vectors
  119:- arithmetic_function('.@'/2).       % dot product of two vectors/matrices
  120:- arithmetic_function(dot/2).        % dot product of two vectors/matrices
  121:- arithmetic_function(determinant/1).  % determinant of square matrix
  122:- arithmetic_function(inverse/1).    % inverse of square matrix
 ndarray(?X:array) is semidet
Succeeds if X is an array; otherwise fails. */
  129% definition of a N dimensional array
  130ndarray(Arr) :- ndarray_(Arr,_).
  131
  132% for internal use, D is first dimension
  133ndarray_(Arr,D) :- (compound(Arr) ; integer(D)), !, 
  134	functor(Arr,#,D).
  135
  136%
  137% helper functions for N dimensional array
  138%
  139% create uninitialized
  140new_ndarray([D|Dn],Arr) :- integer(D), D>0,
  141	ndarray_(Arr,D),
  142	(Dn=[] -> true ; new_fill(D,Arr,Dn)).
  143
  144new_fill(L,V,Dn) :-
  145	(L>0
  146	 ->	I is L-1,
  147	    new_ndarray(Dn,Vs),
  148	    []([I],V,Vs),  %index_arr(V[I],Vs),  % V[I] is Vs[],  %%[]([I],V,Vs),
  149	    new_fill(I,V,Dn)
  150	 ;  true
  151	).
  152
  153% array dimensions (shape)
  154ndarray_dim([D|Dn],Arr) :-
  155	ndarray_(Arr,D), !,
  156	[]([0],Arr,Arr0),  %index_arr(Arr[0],Arr0),
  157	ndarray_dim(Dn,Arr0).
  158ndarray_dim([],_Arr).
  159
  160% map to/from list form
  161map_ndarray(List,Arr) :- compound(Arr), Arr=..[#|Vals], !,
  162	map_elems_(List,Vals).
  163map_ndarray(List,Arr) :- is_list(List),
  164	map_elems_(List,Vals),
  165	Arr=..[#|Vals], !.
  166
  167map_elems_([],[]).
  168map_elems_([E|Es],[E|Vs]) :- atomic(E), !,  % optimize
  169	map_elems_(Es,Vs).
  170map_elems_([E|Es],[V|Vs]) :- 
  171	(map_ndarray(E,V) -> true ; E=V),       % map sub-element of pass thru	
  172	map_elems_(Es,Vs).
  173
  174%
  175% Function: create new array
  176%
  177new(ndarray,Dim,Arr) :- var(Arr), is_list(Dim), !,
  178	new_ndarray(Dim,Arr).
  179
  180%
  181% Function: create new array from nested list
  182%
  183ndarray(List,Arr) :- is_list(List), List \= [],
  184	map_ndarray(List,Arr).
  185
  186%
  187% Function: create new array from nested list
  188%
  189to_list(Arr,List) :- ndarray(Arr),
  190	map_ndarray(List,Arr).
  191
  192%
  193% Function: indexing and slicing - basically uses vector indexing and slicing
  194%
  195[](Arr, Arr) :- 
  196	ndarray(Arr).
  197[]([B:E],Arr,R) :-                   % slicing evaluates to array
  198	ndarray_(Arr,Len),
  199	slice_parameters(B:E,Len,Start,SLen), SLen>0, !,
  200	ndarray_(R,SLen),
  201	sub_term(0,SLen,Start,Arr,R).
  202[]([Ix],Arr,R) :-                      % indexing evaluates to element
  203	ndarray_(Arr,Len),
  204	index_parameters(Ix,Len,I),
  205	AIx is I+1,  % index_parameters are 0 based, arg is 1 based
  206	arg(AIx,Arr,R).
  207
  208sub_term(Ix,Len,Start,From,To) :-
  209	(Ix < Len
  210	 -> ToIx is Ix+1,
  211	    FromIx is Start+ToIx,
  212	    arg(FromIx,From,El),
  213	    arg(ToIx,To,El),
  214	    sub_term(ToIx,Len,Start,From,To)
  215	 ;  true
  216	).
  217
  218%
  219% Function: shape
  220%
  221shape(Arr,S) :- ndarray(Arr),
  222	ndarray_dim(S,Arr).
  223
  224%
  225% Function: number of dimensions (length of shape)
  226%
  227ndim(Arr,N) :- 
  228	shape(Arr,S),
  229	length(S,N).
  230
  231%
  232% Function: number of elements (product of shape elements)
  233%
  234size(Arr,Sz) :- 
  235	shape(Arr,S),
  236	size_(S,1,Sz).
  237
  238size_([],In,In).
  239size_([D|DN],In,Out) :-
  240	Nxt is D*In,
  241	size_(DN,Nxt,Out).
  242
  243%
  244% Function: binds any vars to a value
  245%
  246init(Arr,Val,RArr) :-
  247	ndarray_(Arr,D),	
  248	fill_each(D,Arr,Val),
  249	RArr=Arr.
  250
  251fill_each(N,Arr,Value) :-
  252	(N>0
  253	 -> Ix is N-1,
  254	    []([Ix],Arr,Item),  %index_arr(Arr[Ix],Item),
  255	    (ndarray_(Item,D)
  256	     -> fill_each(D,Item,Value)
  257	     ;  (var(Item) -> Item=Value ; true)
  258	    ),
  259	    fill_each(Ix,Arr,Value)
  260	 ;  true
  261	).
  262
  263%
  264% Function: row concat, dimensions > 0 must be the same
  265%    Note: Concat of two vectors is a single vector.
  266%          Use transpose (twice) for column concat
  267%
  268\\(A1, A2, AR) :-  ndarray(A1), ndarray(A2),
  269	ndarray_dim([R1|S],A1),
  270	ndarray_dim([R2|S],A2),
  271	RR is R1+R2,
  272	new_ndarray([RR|S],AR),
  273	AR[0:R1] is A1[],
  274	AR[R1:_] is A2[].
  275
  276%
  277% Function: flatten
  278%
  279flatten(Arr,FArr) :-
  280	% equivalent to:  size(Arr,S), reshape(Arr,[S],FArr)
  281	FList is flatten(to_list(Arr)),               % flattened input as list
  282	FArr =.. [#|FList].                           % new array from flat data
  283
  284:- arithmetic_function(flatten/1).    % flattened array from local flatten/2
  285
  286%
  287% Function: reshape
  288%
  289reshape(Arr,NShp,NArr) :- 
  290	new_ndarray(NShp,NArr),                       % target array of new shape
  291	length(NShp,N), length(Ixs,N), init_Ix(Ixs),  % starting indices for copy  
  292	FList is flatten(to_list(Arr)),               % flatten input array to list
  293	copy_list(Ixs,NShp,FList,[],NArr).            % copy data in list to target
  294
  295init_Ix([]).
  296init_Ix([0|Ix]) :- init_Ix(Ix). 
  297
  298copy_list([Ix],[D],[Val|Vals],NVals,NArr) :- !,
  299	[]([Ix],NArr,Val),
  300	NIx is Ix+1,
  301	(NIx < D
  302	 -> copy_list([NIx],[D],Vals,NVals,NArr)
  303	 ;  NVals = Vals
  304	).
  305copy_list([Ix|Ixs],[D|Ds],Vals,NVals,NArr) :-
  306	[]([Ix],NArr,NxtArr),
  307	copy_list(Ixs,Ds,Vals,NxtVals,NxtArr),
  308	NIx is Ix+1,
  309	(NIx < D
  310	 -> copy_list([NIx|Ixs],[D|Ds],NxtVals,NVals,NArr)
  311	 ; NVals = NxtVals
  312	).
  313
  314%
  315% Function: transpose
  316%
  317transpose(Arr,TArr)	:- 
  318	shape(Arr,S),
  319	(length(S,1)
  320	 -> TArr = Arr               % 1D case, numpy semantics (T.shape = reverse(A.shape)
  321	 ;  lists:reverse(S,TS),                                 % transposed shape is reverse
  322	    (shape(TArr,TS) -> true ; TArr is new(ndarray,TS)),  % target for transposed 
  323	    transpose_elements(0,[],Arr,TArr)                    % copy values to target
  324	).
  325
  326transpose_elements(Ix,Path,Arr,TArr) :-
  327	[]([Ix],Arr,Val), !,         % fail if Ix too big
  328	NPath = [Ix|Path],           % construct target path in reverse order
  329	(ndarray(Val)
  330	 -> transpose_elements(0,NPath,Val,TArr)  % copy sub array
  331	 ;  set_array_(NPath,TArr,Val)            % copy element
  332	),
  333	Ix1 is Ix+1,
  334	transpose_elements(Ix1,Path,Arr,TArr).
  335transpose_elements(_Ix,_,_,_).   % exit when Ix exceeds dimension.
  336
  337set_array_([Ix],TArr,Val):- !,
  338	[]([Ix],TArr,Val).
  339set_array_([Ix|Path],TArr,Val) :-
  340	[]([Ix],TArr,SubArr),
  341	set_array_(Path,SubArr,Val).
  342
  343% Arithmetic functions on arrays, assumed to be numbers
  344%
  345% Function: arange
  346%
  347arange(ndarray,N,Arr) :- 
  348	Vs is arange(list,N),      % Vs is flat
  349	Arr =.. [#|Vs].
  350
  351arange(ndarray,B,E,Arr) :- number(B), number(E),
  352	Vs is arange(list,B,E),    % Vs is flat
  353	Arr =.. [#|Vs].
  354
  355arange(ndarray,B,E,S,Arr) :- number(B), number(E), number(S),
  356	Vs is arange(list,B,E,S),  % Vs is flat
  357	Arr =.. [#|Vs].
  358
  359%
  360% Function: NxN identity matrix 
  361%
  362identity(ndarray,N,IdArr) :-  integer(N), N>=1,
  363	(var(IdArr) -> new(ndarray,[N,N],IdArr) ; true),
  364	Ix is N-1,
  365	diagonal_(Ix,IdArr),   % fill diagonal
  366	fill_each(N,IdArr,0).  % fill rest with 0
  367
  368diagonal_(N,Arr) :-
  369	(N>=0
  370	 -> Arr[N,N] is 1,
  371	    Nxt is N-1,
  372	    diagonal_(Nxt,Arr)
  373	 ;  true
  374	).
  375
  376%
  377% Function: sum of array elements
  378%
  379sum(A,ASum) :- ndarray(A), !,
  380	ndarray_dim(Shp,A),
  381	sum_array(Shp,A,0,ASum).
  382
  383sum_array([],N,Acc,R) :-
  384	catch(add(Acc,N,R),Err, constrain(Err,R is Acc+N)).  % final dimension, use accumulate number
  385sum_array([D|Ds],A,Acc,R) :-
  386	acc_subarray(D,sum_array,Ds,A,Acc,R).
  387
  388%
  389% Function: min of array elements
  390%
  391min(A,AMin) :- ndarray(A), !,
  392	ndarray_dim(Shp,A),
  393	min_array(Shp,A,inf,AMin).
  394
  395min_array([],N,Acc,R) :-
  396	catch(min(Acc,N,R), Err, constrain(Err,R is min(Acc,N))).  % final dimension, use min
  397min_array([D|Ds],A,Acc,R) :-
  398	acc_subarray(D,min_array,Ds,A,Acc,R).
  399
  400%
  401% Function: max of array elements
  402%
  403max(A,AMax) :- ndarray(A), !,
  404	ndarray_dim(Shp,A),
  405	max_array(Shp,A,-inf,AMax).
  406	
  407max_array([],N,Acc,R) :-
  408	catch(max(Acc,N,R),Err, constrain(Err,R is max(Acc,N))).  % final dimension, use max
  409max_array([D|Ds],A,Acc,R) :-
  410	acc_subarray(D,max_array,Ds,A,Acc,R).
  411
  412%
  413% Function: array addition
  414%
  415'.+'(N,A,AR) :- number(N), ndarray(A), !,
  416	ndarray_dim(Shp,A),
  417	new_ndarray(Shp,AR),
  418	add_arrays(Shp,N,A,AR).
  419
  420'.+'(A,N,AR) :- number(N), ndarray(A), !,
  421	ndarray_dim(Shp,A),
  422	new_ndarray(Shp,AR),
  423	add_arrays(Shp,A,N,AR).
  424
  425'.+'(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  426	ndarray_dim(Shp,A1), ndarray_dim(Shp,A2),  % arrays have same shape
  427	new_ndarray(Shp,AR),
  428	add_arrays(Shp,A1,A2,AR).
  429	
  430add_arrays([],N1,N2,R) :- !,
  431	catch(add(N1,N2,R), Err, constrain(Err, R is N1 + N2)).  % final dimension, use numeric addition
  432add_arrays([D|Ds],A1,A2,AR) :-
  433	do_subarrays(D,add_arrays,Ds,A1,A2,AR).
  434
  435%
  436% Function: array unary minus
  437%
  438'.-'(A,R) :- ndarray(A),
  439	'.*'(-1,A,R).  % multiply by -1
  440
  441%
  442% Function: array subtraction 
  443%
  444'.-'(N,A,AR) :- number(N), ndarray(A), !,
  445	ndarray_dim(Shp,A),
  446	new_ndarray(Shp,AR),
  447	sub_arrays(Shp,N,A,AR).
  448
  449'.-'(A,N,AR) :- number(N), ndarray(A), !,
  450	ndarray_dim(Shp,A),
  451	new_ndarray(Shp,AR),
  452	sub_arrays(Shp,A,N,AR).
  453
  454'.-'(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  455	ndarray_dim(Shp,A1), ndarray_dim(Shp,A2),  % arrays have same shape
  456	new_ndarray(Shp,AR),
  457	sub_arrays(Shp,A1,A2,AR).
  458
  459sub_arrays([],N1,N2,R) :- !,
  460	catch(sub(N1,N2,R), Err, constrain(Err, R is N1 - N2)).  % final dimension, use numeric subtraction
  461sub_arrays([D|Ds],A1,A2,AR) :-
  462	do_subarrays(D,sub_arrays,Ds,A1,A2,AR).
  463
  464%
  465% Function: array multiplication
  466%
  467'.*'(N,A,AR) :- number(N), ndarray(A), !,
  468	ndarray_dim(Shp,A),
  469	new_ndarray(Shp,AR),
  470	mul_arrays(Shp,N,A,AR).
  471
  472'.*'(A,N,AR) :- number(N), ndarray(A), !,
  473	ndarray_dim(Shp,A),
  474	new_ndarray(Shp,AR),
  475	mul_arrays(Shp,A,N,AR).
  476
  477'.*'(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  478	ndarray_dim(Shp,A1), ndarray_dim(Shp,A2),  % arrays have same shape
  479	new_ndarray(Shp,AR),
  480	mul_arrays(Shp,A1,A2,AR).
  481
  482mul_arrays([],N1,N2,R) :- !,
  483	catch(mul(N1,N2,R), Err, constrain(Err, R is N1 * N2)).  % final dimension, use numeric multiplication
  484mul_arrays([D|Ds],A1,A2,AR) :-
  485	do_subarrays(D,mul_arrays,Ds,A1,A2,AR).
  486
  487%
  488% Function: array division
  489%
  490'./'(N,A,AR) :- number(N), ndarray(A), !,
  491	ndarray_dim(Shp,A),
  492	new_ndarray(Shp,AR),
  493	div_arrays(Shp,N,A,AR).
  494
  495'./'(A,N,AR) :- number(N), ndarray(A), !,
  496	ndarray_dim(Shp,A),
  497	new_ndarray(Shp,AR),
  498	div_arrays(Shp,A,N,AR).
  499
  500'./'(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  501	ndarray_dim(Shp,A1), ndarray_dim(Shp,A2),  % arrays have same shape
  502	new_ndarray(Shp,AR),
  503	div_arrays(Shp,A1,A2,AR).
  504
  505div_arrays([],N1,N2,R) :- !,
  506	catch(div(N1,N2,R), Err, constrain(Err, R is N1 / N2)).  % final dimension, use numeric division
  507div_arrays([D|Ds],A1,A2,AR) :-
  508	do_subarrays(D,div_arrays,Ds,A1,A2,AR).
  509
  510%
  511% Function: array power
  512%
  513'.**'(N,A,AR) :- number(N), ndarray(A), !,
  514	ndarray_dim(Shp,A),
  515	new_ndarray(Shp,AR),
  516	pow_arrays(Shp,N,A,AR).
  517
  518'.**'(A,N,AR) :- number(N), ndarray(A), !,
  519	ndarray_dim(Shp,A),
  520	new_ndarray(Shp,AR),
  521	pow_arrays(Shp,A,N,AR).
  522
  523'.**'(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  524	ndarray_dim(Shp,A1), ndarray_dim(Shp,A2),  % arrays have same shape
  525	new_ndarray(Shp,AR),
  526	pow_arrays(Shp,A1,A2,AR).
  527
  528pow_arrays([],N1,N2,R) :- !,
  529	catch(pow(N1,N2,R), Err, constrain(Err, R is N1 ** N2)).  % final dimension, use numeric power
  530pow_arrays([D|Ds],A1,A2,AR) :-
  531	do_subarrays(D,pow_arrays,Ds,A1,A2,AR).
  532
  533%
  534% Function: apply arithmetic function
  535%
  536apply(F,A,AR) :- ndarray(A),
  537	ndarray_dim(Shp,A),
  538	new_ndarray(Shp,AR),
  539	apply_arrays(Shp,A,F,AR).
  540
  541apply_arrays([],N,F,R) :- !,   % final dimension, apply function
  542	E=..[F,N], 
  543	catch(R is E, Err, constrain(Err, R is E)).
  544apply_arrays([D|Ds],A,F,AR) :-
  545	apply_subarrays(D,apply_arrays,Ds,A,F,AR).
  546
  547%
  548% Function: cross product of two 3 dimensional vectors
  549%
  550cross(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  551	ndarray_dim([3],A1), ndarray_dim([3],A2), !,  % two 3D vectors
  552	new_ndarray([3],AR),
  553	cross_(A1,A2,AR).
  554	
  555cross(A1,A2,R) :- ndarray(A1), ndarray(A2),
  556	ndarray_dim([2],A1), ndarray_dim([2],A2),     % two 2D vectors, convert to 3D
  557	new_ndarray([3],D3A1), new_ndarray([3],D3A2),
  558	D3A1[0:2] is A1[], D3A1[2] is 0,  % copy with Z axis = 0
  559	D3A2[0:2] is A2[], D3A2[2] is 0,
  560	new_ndarray([3],AR),
  561	cross_(D3A1,D3A2,AR),
  562	vector3d(AR,_,_,R).               % result is Z (a scaler, X and Y are 0)
  563
  564cross_(A1,A2,AR) :-                   % cross product of two 3D vectors
  565	vector3d(A1, A1_0,A1_1,A1_2),
  566	vector3d(A2, A2_0,A2_1,A2_2),
  567	vector3d(AR, AR_0,AR_1,AR_2),
  568	catch(det2x2(A1_1,A1_2,A2_1,A2_2,AR_0), Err0, constrain(Err0, AR_0 is A1_1*A2_2 - A2_1*A1_2)),
  569	catch(det2x2(A2_0,A2_2,A1_0,A1_2,AR_1), Err1, constrain(Err1, AR_1 is A2_0*A1_2 - A1_0*A2_2)),
  570	catch(det2x2(A1_0,A1_1,A2_0,A2_1,AR_2), Err2, constrain(Err2, AR_2 is A1_0*A2_1 - A2_0*A1_1)).
  571
  572vector3d(#(N0,N1,N2),N0,N1,N2).
  573
  574%
  575% Function: inner product of two vectors of same size
  576%
  577inner(A1,A2,R) :- ndarray(A1), ndarray(A2),
  578	ndim(A1,1), ndim(A2,1),  % vectors only
  579	R is sum(A1 .* A2).      % fails if not same size
  580
  581%
  582% Function: outer product of two vectors
  583%
  584outer(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  585	ndarray_dim([S],A1),  % vector of size S
  586	ndarray_dim([_],A2),  % vector
  587	new_ndarray([S],R),
  588	Ix is S-1,
  589	outer_(Ix,A1,A2,R),
  590	((ndarray_dim([1,N],R), N>1) -> R = #(AR) ; R = AR).  % reduce if vector
  591
  592outer_(Ix,A1,A2,AR) :- Ix >= 0,
  593	R is A1[Ix] .* A2[],
  594	(R = #(V) ->  AR[Ix] is V ; AR[Ix] is R[]),
  595	Nxt is Ix-1,
  596	(Nxt >= 0 -> outer_(Nxt,A1,A2,AR) ; true).
  597
  598%
  599% Function: dot product of two arrays
  600%
  601'.@'(A1,A2,AR) :- dot(A1,A2,AR).
  602
  603dot(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  604	ndarray_dim(S1,A1),  % vector of same size S
  605	ndarray_dim(S2,A2),
  606	dot_(S1,S2,A1,A2,AR).
  607
  608dot_([D], [D], A1, A2, AR) :-  !,        % vectors of same size, take inner product
  609	inner(A1,A2,AR).
  610dot_([D], [D,P], A1, A2, AR) :-  !,      % vector (A1), convert to matrix 
  611	dot_([1,D],[D,P],#(A1),A2,#(AR)).
  612dot_([D,P], [P], A1, A2, AR) :-  !,      % vector (A2), convert to matrix 
  613	reshape(A2,[P,1],A21),
  614	dot_([D,P],[P,1],A1,A21,AR0),
  615	reshape(AR0,[P],AR).
  616dot_([M,N], [N,P] , A1, A2, AR) :-       % matrix multiply
  617	new_ndarray([M,P],AR),
  618	transpose(A2,TA2),   % shape(TA2)  = [P,N]
  619	MIx is M-1, PIx is P-1,
  620	matrix_mul(MIx,PIx,PIx,A1,TA2,AR).
  621
  622matrix_mul(M,P,Rows,A1,A2,AR) :-
  623	AR[M,P] is inner(A1[M],A2[P]),
  624	(matrix_iterate(M,P,Rows,NxtM,NxtP) -> matrix_mul(NxtM,NxtP,Rows,A1,A2,AR) ; true).
  625
  626%
  627% Function: determinant of a square matrix
  628%
  629determinant(A,Det) :- ndarray(A),
  630	ndarray_dim(S,A),
  631	determinant_(S,A,Det).
  632
  633determinant_([1,1],A, A0) :- 
  634	A = #( #(A0) ).
  635determinant_([2,2],A,Det) :-
  636	A = #( #(A_00,A_01) , #(A_10,A_11) ),
  637	catch(det2x2(A_00,A_01,A_10,A_11,Det), Err, constrain(Err, Det is A_00*A_11 - A_01*A_10)).
  638determinant_([N,N],A,Det) :-
  639	map_ndarray([Row|Rest],A),  % convert to list form, use top row
  640	MN is N-1,                  % minor dimension
  641	determinant_(0,N,MN,Row,Rest,0,Det).
  642
  643determinant_(I,N,MN,Row,Rest,Acc,Det) :-
  644	(I<N
  645	 ->	El is Row[I],
  646	 	(number(El), El =:= 0
  647	 	 -> NxtA = Acc  % element is 0, no use calculating determinant
  648	 	 ;	minor_determinant_(Rest,I,MN,MDet),
  649 			I1 is 1 - 2*(I mod 2), % alternating 1 and -1 
  650	 	 	catch(acc_det(Acc,El,I1,MDet,NxtA), Err, constrain(Err, NxtA is Acc + El*I1*MDet))
  651	 	),
  652		NxtI is I+1,
  653		determinant_(NxtI,N,MN,Row,Rest,NxtA,Det)
  654	 ;	Det=Acc
  655	).
  656
  657minor_determinant_(List,I,2,Det) :- !,
  658	minor(List,I,[ [A_00,A_01], [A_10,A_11] ]),
  659	catch(det2x2(A_00,A_01,A_10,A_11,Det), Err, constrain(Err, Det is A_00*A_11 - A_01*A_10)).
  660minor_determinant_(List,I,MN,Det) :-
  661	minor(List,I,[Row|Rest]),      % list form of minor
  662	MMN is MN-1,
  663	determinant_(0,MN,MMN,Row,Rest,0,Det).	
  664
  665minor([],_I,[]).
  666minor([X|Xs],I,[M|Ms]) :-
  667	delete_item(I,X,M),
  668	minor(Xs,I,Ms).
  669
  670delete_item(0,[_|Xs],Xs) :- !.
  671delete_item(I,[X|Xs],[X|Ms]) :-
  672	Nxt is I-1,	
  673	delete_item(Nxt,Xs,Ms).
  674
  675%
  676% Function: inverse of a square matrix
  677%
  678inverse(A,Inv) :- ndarray(A),
  679	ndarray_dim(S,A),
  680	inverse_(S,A,Inv).
  681
  682inverse_([1,1],A,Inv) :- !,
  683	determinant_([1,1],A,ADet),
  684	catch(inv(ADet,IN), Err, constrain(Err, IN is 1/ADet)),  % fails if det=0
  685	ndarray([[IN]],Inv).
  686inverse_([2,2],A,Inv) :- !,     % special case 2x2 to avoid vector edge cases in NxN
  687	determinant_([2,2],A,ADet),
  688	catch(inv(ADet,IDet), Err, constrain(Err, IDet is 1/ADet)),  % fails if det=0
  689	A = #( #(A_00,A_01) , #(A_10,A_11) ),
  690	catch(mul(A_11,IDet,I_00), Err00, constrain(Err00,I_00 is  A_11*IDet)),
  691	catch(mml(A_01,IDet,I_01), Err01, constrain(Err01,I_01 is -A_01*IDet)),
  692	catch(mml(A_10,IDet,I_10), Err10, constrain(Err10,I_10 is -A_10*IDet)),
  693	catch(mul(A_00,IDet,I_11), Err11, constrain(Err11,I_11 is  A_00*IDet)),
  694	ndarray([[I_00, I_01], [I_10, I_11]],Inv).
  695inverse_([N,N],A,Inv) :-
  696	determinant_([N,N],A,ADet), 
  697	catch(inv(ADet,IDet), Err, constrain(Err, IDet is 1/ADet)),  % fails if det=0
  698	new_ndarray([N,N],Inv),
  699	map_ndarray(AList,A),       % use list form for minor/3
  700	D is N-1,  % for indexing
  701	inverse_iterate_(D,D,D,IDet,AList,Inv).
  702
  703inverse_iterate_(R,C,D,IDet,AList,Inv) :-
  704	sub_array_(R,D,AList,Sub),
  705	minor_determinant_(Sub,C,D,SDet),
  706	X is Inv[C,R],             % inverted array element (transpose(A[R,C])
  707	Sgn is 1-2*((R+C)mod 2),   % alternating -1,1	
  708	catch(mul3(Sgn,SDet,IDet,X), Err, constrain(Err, X is Sgn*SDet*IDet)),
  709	(matrix_iterate(C,R,D,NxtC,NxtR)
  710	 -> inverse_iterate_(NxtR,NxtC,D,IDet,AList,Inv)
  711	 ;  true
  712	).
  713
  714sub_array_(0,_,A,Sub) :- !,
  715	Sub is A[1:_].
  716sub_array_(D,D,A,Sub) :- !,
  717	Sub is A[0:D].
  718sub_array_(R,_,A,Sub) :-
  719	R1 is R+1,
  720	Sub is A[0:R]\\A[R1:_].
  721
  722%
  723% helpers for iterating through nested array levels
  724%
  725do_subarrays(N,Op,Ds,V1,V2,VR) :-
  726	(N>0
  727	 ->	I is N-1,
  728		(number(V1) -> V1i=V1 ; []([I],V1,V1i)),  %index_arr(V1[I],V1i)),  %  VS1 is V1[I]),  % simple "broadcast" 
  729		(number(V2) -> V2i=V2 ; []([I],V2,V2i)),  %index_arr(V2[I],V2i)),  %  VS2 is V2[I]),
  730		[]([I],VR,VRi),  % index_arr(VR[I],VRi),  % VSR is VR[I],
  731		SubOp =.. [Op,Ds,V1i,V2i,VRi], SubOp,   % SubOp(Ds,VS1,VS2,VSR)
  732		do_subarrays(I,Op,Ds,V1,V2,VR)
  733	 ;	true
  734	).
  735
  736apply_subarrays(N,Op,Ds,V,F,VR) :-
  737	(N>0
  738	 ->	I is N-1,
  739		(number(V) -> Vi=V ; []([I],V,Vi)),  % index_arr(V[I],Vi)),  %  VS is V[I]),  % simple "broadcast" 
  740		[]([I],VR,VRi),  % index_arr(VR[I],VRi),  % VSR is VR[I],
  741		SubOp =.. [Op,Ds,Vi,F,VRi], SubOp,          % SubOp(Ds,VS1,VS2,VSR)
  742		apply_subarrays(I,Op,Ds,V,F,VR)
  743	 ;	true
  744	).
  745
  746acc_subarray(N,Op,Ds,V,Acc,R) :-
  747	(N>0
  748	 ->	I is N-1,
  749		(number(V) -> Vi=V ; []([I],V,Vi)),        % simple "broadcst" 
  750		SubOp =.. [Op,Ds,Vi,Acc,AccR], SubOp,      % SubOp(Ds,VS1,VS2,VSR)
  751		acc_subarray(I,Op,Ds,V,AccR,R)
  752	 ;	R=Acc
  753	).
  754
  755%
  756% arithmetic --> clpBNR constraint
  757%
  758constrain(error(instantiation_error,_), Goal) :-
  759	current_module(clpBNR),            % clpBNR constraints available?
  760	Mod=clpBNR, Mod:constrain_(Goal),  % can't use module name directly as it invalidates current_module test
  761	!.
  762constrain(Err, _Goal) :-
  763	throw(Err).
  764
  765%
  766% Compiled arithmetic (Note: adding/subtracting a constant already optimal)
  767%
  768:- set_prolog_flag(optimise,true).  769
  770matrix_iterate(C,0,Rows,NxtC,Rows) :- !, C>0,
  771	NxtC is C-1.
  772matrix_iterate(C,R,_Rows,C,NxtR) :- R>0,
  773	NxtR is R-1.
  774
  775add(N1,N2,R) :- R is  N1+N2.
  776sub(N1,N2,R) :- R is  N1-N2.
  777mul(N1,N2,R) :- R is  N1*N2.
  778mml(N1,N2,R) :- R is -N1*N2.
  779div(N1,N2,R) :- R is  N1/N2.
  780pow(N1,N2,R) :- R is  N1**N2.
  781max(N1,N2,R) :- R is  max(N1,N2).
  782min(N1,N2,R) :- R is  min(N1,N2).
  783mul3(N1,N2,N3,R) :- R is N1*N2*N3.
  784inv(N,R)     :- N=\=0, R is 1/N.
  785det2x2(A0_0,A0_1,A1_0,A1_1,R) :- R is A0_0*A1_1 - A0_1*A1_0.
  786acc_det(Acc,El,I1,MDet,NxtA)  :- NxtA is Acc + El*I1*MDet