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(200, xfx, '.**')  
   31	]).   32
   33%:- reexport(library(type_list),[op(_,_,_)]).  % requires fix circa 9.3.19 % for operators
   34:- 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 Pyhon'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 array
:- 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(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. */

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