1:- module(type_ndarray,
    2	[ ndarray/1,
    3	  op(500, yfx, '.+'),      % arithmetic operators on arrays
    4	  op(200, fy,  '.-'),      % same precedence and associativity as standard ops
    5	  op(500, yfx, '.-'),
    6	  op(400, yfx, '.*'),
    7	  op(400, yfx, './') , 
    8	  op(200, xfx, '.**')  
    9	]).   10
   11:- current_module(arithmetic_types) -> true ; use_module(library(arithmetic_types)).   12% for indexing and slicing support, arange/n content
   13:- current_module(type_list) 
   14	-> true 
   15	;  reexport(library(type_list)).  % reexport for operators
   16
   17:- arithmetic_function(new/2).        % new from shape, uninitialized
   18:- arithmetic_function(ndarray/1).    % new from ListOfLists
   19:- arithmetic_function(to_list/1).    % list from ndarray
   20:- arithmetic_function(ndim/1).       % number of dimensions
   21:- arithmetic_function(shape/1).      % shape (a list)
   22:- arithmetic_function(size/1).       % number of values
   23:- arithmetic_function([]/1).         % block indexing and slicing
   24:- arithmetic_function([]/2).   25:- arithmetic_function(init/2).       % initialize any variables
   26:- arithmetic_function(\\ /2).        % row concat
   27:- arithmetic_function(transpose/1).  % transpose array
   28% numeric functions
   29:- arithmetic_function(arange/2).     % new from range
   30:- arithmetic_function(arange/3).     % new from range
   31:- arithmetic_function(arange/4).     % new from range
   32:- arithmetic_function(identity/2).   % new NxN identity array
   33:- arithmetic_function(sum/1).        % sum of elements
   34:- arithmetic_function(min/1).        % minimum of elements
   35:- arithmetic_function(max/1).        % maximum of elements
   36:- arithmetic_function('.+'/2).       % numeric array addition
   37:- arithmetic_function('.-'/1).       % numeric array unary minus
   38:- arithmetic_function('.-'/2).       % numeric array subtraction
   39:- arithmetic_function('.*'/2).       % numeric array product
   40:- arithmetic_function('./'/2).       % numeric array division
   41:- arithmetic_function('.**'/2).      % numeric array power
   42:- arithmetic_function(apply/2).      % apply function/1 to array
   43:- arithmetic_function(cross/2).      % cross product of two 3D vectors
   44:- arithmetic_function(inner/2).      % inner product of two vectors
   45:- arithmetic_function(outer/2).      % outer product of two vectors
   46:- arithmetic_function(dot/2).        % dot product of two vectors/matrices
   47:- arithmetic_function(determinant/1).  % determinant of square matrix
   48:- arithmetic_function(inverse/1).    % inverse of square matrix
   49
   50% definition of a N dimensional array
   51ndarray(Arr) :- ndarray_(Arr,_).
   52
   53% for internal use, D is first dimension	
   54ndarray_(Arr,D) :- (compound(Arr) ; integer(D)), !, 
   55	functor(Arr,#,D).
   56
   57%
   58% helper functions for N dimensional array
   59%
   60% create uninitialized
   61new_ndarray([D|Dn],Arr) :- integer(D), D>0,
   62	ndarray_(Arr,D),
   63	(Dn=[] -> true ; new_fill(D,Arr,Dn)).
   64
   65new_fill(L,V,Dn) :-
   66	(L>0
   67	 ->	I is L-1,
   68		new_ndarray(Dn,Vs),
   69		[]([I],V,Vs),  %index_arr(V[I],Vs),  % V[I] is Vs[],  %%[]([I],V,Vs),
   70		new_fill(I,V,Dn)
   71	 ;	true
   72	).
   73
   74% array dimensions (shape)
   75ndarray_dim([D|Dn],Arr) :-
   76	ndarray_(Arr,D), !,
   77	[]([0],Arr,Arr0),  %index_arr(Arr[0],Arr0),
   78	ndarray_dim(Dn,Arr0).
   79ndarray_dim([],_Arr).
   80
   81% map to/from list form
   82map_ndarray(List,Arr) :- compound(Arr), Arr=..[#|Vals], !,
   83	map_elems_(List,Vals).
   84map_ndarray(List,Arr) :- is_list(List),
   85	map_elems_(List,Vals),
   86	Arr=..[#|Vals], !.
   87	
   88map_elems_([],[]).
   89map_elems_([E|Es],[E|Vs]) :- atomic(E), !,  % optimize
   90	map_elems_(Es,Vs).
   91map_elems_([E|Es],[V|Vs]) :- 
   92	(map_ndarray(E,V) -> true ; E=V),       % map sub-element of pass thru	
   93	map_elems_(Es,Vs).
   94
   95%
   96% Function: create new array
   97%
   98new(ndarray,Dim,Arr) :- var(Arr), is_list(Dim), !,
   99	new_ndarray(Dim,Arr).
  100	
  101%
  102% Function: create new array from nested list
  103%
  104ndarray(List,Arr) :- is_list(List), List \= [],
  105	map_ndarray(List,Arr).
  106		
  107%
  108% Function: create new array from nested list
  109%
  110to_list(Arr,List) :- ndarray(Arr),
  111	map_ndarray(List,Arr).
  112		
  113%
  114% Function: indexing and slicing - basically uses vector indexing and slicing
  115%
  116[](Arr, Arr) :- 
  117	ndarray(Arr).
  118[]([B:E],Arr,R) :-                   % slicing evaluates to array
  119	ndarray_(Arr,Len),
  120	slice_parameters(B:E,Len,Start,SLen), SLen>0, !,
  121	ndarray_(R,SLen),
  122	sub_term(0,SLen,Start,Arr,R).
  123[]([Ix],Arr,R) :-                      % indexing evaluates to element
  124	ndarray_(Arr,Len),	
  125	index_parameters(Ix,Len,I),
  126	AIx is I+1,  % index_parameters are 0 based, arg is 1 based
  127	arg(AIx,Arr,R).
  128
  129sub_term(Ix,Len,Start,From,To) :-
  130	(Ix < Len
  131	 -> ToIx is Ix+1,
  132	 	FromIx is Start+ToIx,
  133		arg(FromIx,From,El),
  134		arg(ToIx,To,El),
  135		sub_term(ToIx,Len,Start,From,To)
  136	 ;  true
  137	).
  138
  139%
  140% Function: shape
  141%
  142shape(Arr,S) :- ndarray(Arr),
  143	ndarray_dim(S,Arr).
  144
  145%
  146% Function: number of dimensions (length of shape)
  147%
  148ndim(Arr,N) :- 
  149	shape(Arr,S),
  150	length(S,N).
  151
  152%
  153% Function: number of elements (product of shape elements)
  154%
  155size(Arr,Sz) :- 
  156	shape(Arr,S),
  157	size_(S,1,Sz).
  158	
  159size_([],In,In).
  160size_([D|DN],In,Out) :-
  161	Nxt is D*In,
  162	size_(DN,Nxt,Out).
  163
  164%
  165% Function: binds any vars to a value
  166%
  167init(Arr,Val,RArr) :-
  168	ndarray_(Arr,D),	
  169	fill_each(D,Arr,Val),
  170	RArr=Arr.
  171
  172fill_each(N,Arr,Value) :-
  173	(N>0
  174	 ->	Ix is N-1,
  175	 	[]([Ix],Arr,Item),  %index_arr(Arr[Ix],Item),
  176		(ndarray_(Item,D)
  177		 -> fill_each(D,Item,Value)
  178		 ;	(var(Item) -> Item=Value ; true)
  179		),
  180		fill_each(Ix,Arr,Value)
  181	 ;	true
  182	).
  183
  184%
  185% Function: row concat, dimensions > 0 must be the same
  186%    Note: Concat of two vectors is a single vector.
  187%          Use transpose (twice) for column concat
  188%
  189\\(A1, A2, AR) :-  ndarray(A1), ndarray(A2),
  190	ndarray_dim([R1|S],A1),
  191	ndarray_dim([R2|S],A2),
  192	RR is R1+R2,
  193	new_ndarray([RR|S],AR),
  194	AR[0:R1] is A1[],
  195	AR[R1:_] is A2[].
  196	
  197%
  198% Function: transpose (1D=noop, 2D=matrix transpose, 3+D=2D matrix of array)
  199%
  200transpose(Arr,TArr)	:- 
  201	shape(Arr,S),
  202	transpose_(S,Arr,TArr).
  203
  204transpose_([C],Arr,TArr) :- !,       % 1D case result will be Cx1 array
  205	transpose1_([1,C],#(Arr),TArr).
  206transpose_([R,1|S],Arr,TArr) :- !,     % result will be a 1xR vector
  207	transpose1_([R,1|S],Arr,#(TArr)).
  208transpose_(Shp,Arr,TArr) :-          % general 2D case
  209	transpose1_(Shp,Arr,TArr).
  210	
  211transpose1_([R,C|S],Arr,TArr) :-
  212	new_ndarray([C,R|S],TArr),
  213	RIx is R-1,
  214	CIx is C-1,
  215	transpose_mtrx(CIx,RIx,RIx,Arr,TArr).
  216
  217transpose_mtrx(C,R,Rows,Arr,TArr) :-
  218	V is Arr[R,C],
  219	V is TArr[C,R],
  220	(matrix_iterate(C,R,Rows,NxtC,NxtR) -> transpose_mtrx(NxtC,NxtR,Rows,Arr,TArr) ; true).
  221
  222% Arithmetic functions on arrays, assumed to be numbers
  223%
  224% Function: arange
  225%
  226arange(ndarray,N,Arr) :- 
  227	Vs is arange(list,N),
  228	ndarray(Vs,Arr).
  229
  230arange(ndarray,B,E,Arr) :- number(B), number(E),
  231	Vs is arange(list,B,E),
  232	ndarray(Vs,Arr).
  233
  234arange(ndarray,B,E,S,Arr) :- number(B), number(E), number(S),
  235	Vs is arange(list,B,E,S),
  236	ndarray(Vs,Arr).
  237	
  238%
  239% Function: NxN identity matrix 
  240%
  241identity(ndarray,1,IdArr) :-  ndarray([1],IdArr).
  242identity(ndarray,N,IdArr) :-  integer(N), N>1,
  243	(var(IdArr) -> new(ndarray,[N,N],IdArr) ; true),
  244	Ix is N-1,
  245	diagonal_(Ix,IdArr),   % fill diagonal
  246	fill_each(N,IdArr,0).  % fill rest with 0
  247	
  248
  249diagonal_(N,Arr) :-
  250	(N>=0
  251	 ->	Arr[N,N] is 1,
  252		Nxt is N-1,
  253		diagonal_(Nxt,Arr)
  254	 ;  true
  255	).
  256
  257%
  258% Function: sum of array elements
  259%
  260sum(A,ASum) :- ndarray(A), !,
  261	ndarray_dim(Shp,A),
  262	sum_array(Shp,A,0,ASum).
  263	
  264sum_array([],N,Acc,R) :-
  265	catch(add(Acc,N,R),Err, constrain(Err,R is Acc+N)).  % final dimension, use accumulate number
  266sum_array([D|Ds],A,Acc,R) :-
  267	acc_subarray(D,sum_array,Ds,A,Acc,R).
  268
  269%
  270% Function: min of array elements
  271%
  272min(A,AMin) :- ndarray(A), !,
  273	ndarray_dim(Shp,A),
  274	min_array(Shp,A,inf,AMin).
  275	
  276min_array([],N,Acc,R) :-
  277	catch(min(Acc,N,R), Err, constrain(Err,R is min(Acc,N))).  % final dimension, use min
  278min_array([D|Ds],A,Acc,R) :-
  279	acc_subarray(D,min_array,Ds,A,Acc,R).
  280
  281%
  282% Function: max of array elements
  283%
  284max(A,AMax) :- ndarray(A), !,
  285	ndarray_dim(Shp,A),
  286	max_array(Shp,A,-inf,AMax).
  287	
  288max_array([],N,Acc,R) :-
  289	catch(max(Acc,N,R),Err, constrain(Err,R is max(Acc,N))).  % final dimension, use max
  290max_array([D|Ds],A,Acc,R) :-
  291	acc_subarray(D,max_array,Ds,A,Acc,R).
  292
  293%
  294% Function: array addition
  295%
  296'.+'(N,A,AR) :- number(N), ndarray(A), !,
  297	ndarray_dim(Shp,A),
  298	new_ndarray(Shp,AR),
  299	add_arrays(Shp,N,A,AR).
  300
  301'.+'(A,N,AR) :- number(N), ndarray(A), !,
  302	ndarray_dim(Shp,A),
  303	new_ndarray(Shp,AR),
  304	add_arrays(Shp,A,N,AR).
  305
  306'.+'(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  307	ndarray_dim(Shp,A1), ndarray_dim(Shp,A2),  % arrays have same shape
  308	new_ndarray(Shp,AR),
  309	add_arrays(Shp,A1,A2,AR).
  310	
  311add_arrays([],N1,N2,R) :- !,
  312	catch(add(N1,N2,R), Err, constrain(Err, R is N1 + N2)).  % final dimension, use numeric addition
  313add_arrays([D|Ds],A1,A2,AR) :-
  314	do_subarrays(D,add_arrays,Ds,A1,A2,AR).
  315
  316%
  317% Function: array unary minus
  318%
  319'.-'(A,R) :- ndarray(A),
  320	'.*'(-1,A,R).  % multiply by -1
  321
  322%
  323% Function: array subtraction 
  324%
  325'.-'(N,A,AR) :- number(N), ndarray(A), !,
  326	ndarray_dim(Shp,A),
  327	new_ndarray(Shp,AR),
  328	sub_arrays(Shp,N,A,AR).
  329
  330'.-'(A,N,AR) :- number(N), ndarray(A), !,
  331	ndarray_dim(Shp,A),
  332	new_ndarray(Shp,AR),
  333	sub_arrays(Shp,A,N,AR).
  334
  335'.-'(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  336	ndarray_dim(Shp,A1), ndarray_dim(Shp,A2),  % arrays have same shape
  337	new_ndarray(Shp,AR),
  338	sub_arrays(Shp,A1,A2,AR).
  339	
  340sub_arrays([],N1,N2,R) :- !,
  341	catch(sub(N1,N2,R), Err, constrain(Err, R is N1 - N2)).  % final dimension, use numeric subtraction
  342sub_arrays([D|Ds],A1,A2,AR) :-
  343	do_subarrays(D,sub_arrays,Ds,A1,A2,AR).
  344
  345%
  346% Function: array multiplication
  347%
  348'.*'(N,A,AR) :- number(N), ndarray(A), !,
  349	ndarray_dim(Shp,A),
  350	new_ndarray(Shp,AR),
  351	mul_arrays(Shp,N,A,AR).
  352
  353'.*'(A,N,AR) :- number(N), ndarray(A), !,
  354	ndarray_dim(Shp,A),
  355	new_ndarray(Shp,AR),
  356	mul_arrays(Shp,A,N,AR).
  357
  358'.*'(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  359	ndarray_dim(Shp,A1), ndarray_dim(Shp,A2),  % arrays have same shape
  360	new_ndarray(Shp,AR),
  361	mul_arrays(Shp,A1,A2,AR).
  362	
  363mul_arrays([],N1,N2,R) :- !,
  364	catch(mul(N1,N2,R), Err, constrain(Err, R is N1 * N2)).  % final dimension, use numeric multiplication
  365mul_arrays([D|Ds],A1,A2,AR) :-
  366	do_subarrays(D,mul_arrays,Ds,A1,A2,AR).
  367
  368%
  369% Function: array division
  370%
  371'./'(N,A,AR) :- number(N), ndarray(A), !,
  372	ndarray_dim(Shp,A),
  373	new_ndarray(Shp,AR),
  374	div_arrays(Shp,N,A,AR).
  375
  376'./'(A,N,AR) :- number(N), ndarray(A), !,
  377	ndarray_dim(Shp,A),
  378	new_ndarray(Shp,AR),
  379	div_arrays(Shp,A,N,AR).
  380
  381'./'(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  382	ndarray_dim(Shp,A1), ndarray_dim(Shp,A2),  % arrays have same shape
  383	new_ndarray(Shp,AR),
  384	div_arrays(Shp,A1,A2,AR).
  385	
  386div_arrays([],N1,N2,R) :- !,
  387	catch(div(N1,N2,R), Err, constrain(Err, R is N1 / N2)).  % final dimension, use numeric division
  388div_arrays([D|Ds],A1,A2,AR) :-
  389	do_subarrays(D,div_arrays,Ds,A1,A2,AR).
  390
  391%
  392% Function: array power
  393%
  394'.**'(N,A,AR) :- number(N), ndarray(A), !,
  395	ndarray_dim(Shp,A),
  396	new_ndarray(Shp,AR),
  397	pow_arrays(Shp,N,A,AR).
  398
  399'.**'(A,N,AR) :- number(N), ndarray(A), !,
  400	ndarray_dim(Shp,A),
  401	new_ndarray(Shp,AR),
  402	pow_arrays(Shp,A,N,AR).
  403
  404'.**'(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  405	ndarray_dim(Shp,A1), ndarray_dim(Shp,A2),  % arrays have same shape
  406	new_ndarray(Shp,AR),
  407	pow_arrays(Shp,A1,A2,AR).
  408	
  409pow_arrays([],N1,N2,R) :- !,
  410	catch(pow(N1,N2,R), Err, constrain(Err, R is N1 ** N2)).  % final dimension, use numeric power
  411pow_arrays([D|Ds],A1,A2,AR) :-
  412	do_subarrays(D,pow_arrays,Ds,A1,A2,AR).
  413	
  414%
  415% Function: apply arithmetic function	
  416%
  417apply(F,A,AR) :- ndarray(A),
  418	ndarray_dim(Shp,A),
  419	new_ndarray(Shp,AR),
  420	apply_arrays(Shp,A,F,AR).
  421
  422apply_arrays([],N,F,R) :- !,   % final dimension, apply function
  423	E=..[F,N], 
  424	catch(R is E, Err, constrain(Err, R is E)).
  425apply_arrays([D|Ds],A,F,AR) :-
  426	apply_subarrays(D,apply_arrays,Ds,A,F,AR).
  427
  428%
  429% Function: cross product of two 3 dimensional vectors	
  430%
  431cross(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  432	ndarray_dim([3],A1), ndarray_dim([3],A2), !,  % two 3D vectors
  433	new_ndarray([3],AR),
  434	cross_(A1,A2,AR).
  435	
  436cross(A1,A2,R) :- ndarray(A1), ndarray(A2),
  437	ndarray_dim([2],A1), ndarray_dim([2],A2),     % two 2D vectors, convert to 3D
  438	new_ndarray([3],D3A1), new_ndarray([3],D3A2),
  439	D3A1[0:2] is A1[], D3A1[2] is 0,  % copy with Z axis = 0
  440	D3A2[0:2] is A2[], D3A2[2] is 0,
  441	new_ndarray([3],AR),	
  442	cross_(D3A1,D3A2,AR),
  443	vector3d(AR,_,_,R).               % result is Z (a scaler, X and Y are 0)
  444
  445cross_(A1,A2,AR) :-                   % cross product of two 3D vectors
  446	vector3d(A1, A1_0,A1_1,A1_2),
  447	vector3d(A2, A2_0,A2_1,A2_2),
  448	vector3d(AR, AR_0,AR_1,AR_2),
  449	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)),
  450	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)),
  451	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)).
  452
  453vector3d(#(N0,N1,N2),N0,N1,N2).
  454
  455%
  456% Function: inner product of two vectors of same size	
  457%
  458inner(A1,A2,R) :- ndarray(A1), ndarray(A2),
  459	ndim(A1,1), ndim(A2,1),  % vectors only
  460	R is sum(A1 .* A2).      % fails if not same size
  461
  462%
  463% Function: outer product of two vectors
  464%
  465outer(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  466	ndarray_dim([S],A1),  % vector of size S
  467	ndarray_dim([_],A2),  % vector
  468	new_ndarray([S],R),
  469	Ix is S-1,
  470	outer_(Ix,A1,A2,R),
  471	((ndarray_dim([1,N],R), N>1) -> R = #(AR) ; R = AR).  % reduce if vector
  472
  473outer_(Ix,A1,A2,AR) :- Ix >= 0,
  474	R is A1[Ix] .* A2[],
  475	(R = #(V) ->  AR[Ix] is V ; AR[Ix] is R[]),
  476	Nxt is Ix-1,
  477	(Nxt >= 0 -> outer_(Nxt,A1,A2,AR) ; true).
  478
  479%
  480% Function: dot product of two arrays	
  481%
  482dot(A1,A2,AR) :- ndarray(A1), ndarray(A2),
  483	ndarray_dim(S1,A1),  % vector of same size S
  484	ndarray_dim(S2,A2),
  485	dot_(S1,S2,A1,A2,AR).
  486	
  487dot_([D], [D], A1, A2, AR) :-  !,        % vectors of same size, take inner product
  488	inner(A1,A2,AR).
  489dot_([D], [D,P], A1, A2, AR) :-  !,      % vector (A1), convert to matrix 
  490	dot_([1,D],[D,P],#(A1),A2,AR).
  491dot_([M,N], [N,P] , A1, A2, AR) :-  % matrix multiply
  492	new_ndarray([M,P],AR),
  493	transpose(A2,TA2),   % shape(TA2)  = [P,N]
  494	(P = 1 -> TTA2 = #(TA2) ; TTA2 = TA2),  % if vector(A2), convert to matrix for inner
  495	MIx is M-1, PIx is P-1,
  496	matrix_mul(MIx,PIx,PIx,A1,TTA2,AR).
  497
  498matrix_mul(M,P,Rows,A1,A2,AR) :-
  499	AR[M,P] is inner(A1[M],A2[P]),
  500	(matrix_iterate(M,P,Rows,NxtM,NxtP) -> matrix_mul(NxtM,NxtP,Rows,A1,A2,AR) ; true).
  501
  502%
  503% Function: determinant of a square matrix	
  504%
  505determinant(A,Det) :- ndarray(A),
  506	ndarray_dim(S,A),
  507	determinant_(S,A,Det).
  508
  509determinant_([1],A,Det) :- 
  510	[]([0],A,Det).  %index_arr(A[0],Det).
  511determinant_([2,2],A,Det) :-
  512	A = #( #(A_00,A_01) , #(A_10,A_11) ),
  513	catch(det2x2(A_00,A_01,A_10,A_11,Det), Err, constrain(Err, Det is A_00*A_11 - A_01*A_10)).
  514determinant_([N,N],A,Det) :-
  515	map_ndarray([Row|Rest],A),  % convert to list form, use top row
  516	MN is N-1,                  % minor dimension
  517	determinant_(0,N,MN,Row,Rest,0,Det).
  518	
  519determinant_(I,N,MN,Row,Rest,Acc,Det) :-
  520	(I<N
  521	 ->	El is Row[I],
  522	 	(number(El), El =:= 0
  523	 	 -> NxtA = Acc  % element is 0, no use calculating determinant
  524	 	 ;	minor_determinant_(Rest,I,MN,MDet),
  525 			I1 is 1 - 2*(I mod 2), % alternating 1 and -1 
  526	 	 	catch(acc_det(Acc,El,I1,MDet,NxtA), Err, constrain(Err, NxtA is Acc + El*I1*MDet))
  527	 	),
  528		NxtI is I+1,
  529		determinant_(NxtI,N,MN,Row,Rest,NxtA,Det)
  530	 ;	Det=Acc
  531	).
  532
  533minor_determinant_(List,I,2,Det) :- !,
  534	minor(List,I,[ [A_00,A_01], [A_10,A_11] ]),
  535	catch(det2x2(A_00,A_01,A_10,A_11,Det), Err, constrain(Err, Det is A_00*A_11 - A_01*A_10)).
  536minor_determinant_(List,I,MN,Det) :-
  537	minor(List,I,[Row|Rest]),      % list form of minor
  538	MMN is MN-1,
  539	determinant_(0,MN,MMN,Row,Rest,0,Det).	
  540
  541minor([],_I,[]).
  542minor([X|Xs],I,[M|Ms]) :-
  543	delete_item(I,X,M),
  544	minor(Xs,I,Ms).
  545	
  546delete_item(0,[_|Xs],Xs) :- !.
  547delete_item(I,[X|Xs],[X|Ms]) :-
  548	Nxt is I-1,	
  549	delete_item(Nxt,Xs,Ms).
  550
  551%
  552% Function: inverse of a square matrix	
  553%
  554inverse(A,Inv) :- ndarray(A),
  555	ndarray_dim(S,A),
  556	inverse_(S,A,Inv).
  557	
  558inverse_([1],A,Inv) :- !,
  559	A0 is A[0],
  560	catch(inv(A0,IN), Err, constrain(Err, IN is 1/A0)),  % fails if det=0
  561	ndarray([IN],Inv).
  562inverse_([2,2],A,Inv) :- !,     % special case 2x2 to avoid vector edge cases in NxN
  563	determinant_([2,2],A,ADet),
  564	catch(inv(ADet,IDet), Err, constrain(Err, IDet is 1/ADet)),  % fails if det=0
  565	A = #( #(A_00,A_01) , #(A_10,A_11) ),
  566	catch(mul(A_11,IDet,I_00), Err00, constrain(Err00,I_00 is  A_11*IDet)),
  567	catch(mml(A_01,IDet,I_01), Err01, constrain(Err01,I_01 is -A_01*IDet)),
  568	catch(mml(A_10,IDet,I_10), Err10, constrain(Err10,I_10 is -A_10*IDet)),
  569	catch(mul(A_00,IDet,I_11), Err11, constrain(Err11,I_11 is  A_00*IDet)),
  570	ndarray([[I_00, I_01], [I_10, I_11]],Inv).
  571inverse_([N,N],A,Inv) :-
  572	determinant_([N,N],A,ADet), 
  573	catch(inv(ADet,IDet), Err, constrain(Err, IDet is 1/ADet)),  % fails if det=0
  574	new_ndarray([N,N],Inv),
  575	map_ndarray(AList,A),       % use list form for minor/3
  576	D is N-1,  % for indexing
  577	inverse_iterate_(D,D,D,IDet,AList,Inv).
  578	
  579inverse_iterate_(R,C,D,IDet,AList,Inv) :-
  580	sub_array_(R,D,AList,Sub),
  581	minor_determinant_(Sub,C,D,SDet),
  582	X is Inv[C,R],             % inverted array element (transpose(A[R,C])
  583	Sgn is 1-2*((R+C)mod 2),   % alternating -1,1	
  584	catch(mul3(Sgn,SDet,IDet,X), Err, constrain(Err, X is Sgn*SDet*IDet)),
  585	(matrix_iterate(C,R,D,NxtC,NxtR)
  586	 -> inverse_iterate_(NxtR,NxtC,D,IDet,AList,Inv)
  587	 ;  true
  588	).
  589	
  590sub_array_(0,_,A,Sub) :- !,
  591	Sub is A[1:_].
  592sub_array_(D,D,A,Sub) :- !,
  593	Sub is A[0:D].
  594sub_array_(R,_,A,Sub) :-
  595	R1 is R+1,
  596	Sub is A[0:R]\\A[R1:_].
  597
  598%
  599% helpers for iterating through nested array levels
  600%
  601do_subarrays(N,Op,Ds,V1,V2,VR) :-
  602	(N>0
  603	 ->	I is N-1,
  604		(number(V1) -> V1i=V1 ; []([I],V1,V1i)),  %index_arr(V1[I],V1i)),  %  VS1 is V1[I]),  % simple "broadcast" 
  605		(number(V2) -> V2i=V2 ; []([I],V2,V2i)),  %index_arr(V2[I],V2i)),  %  VS2 is V2[I]),
  606		[]([I],VR,VRi),  % index_arr(VR[I],VRi),  % VSR is VR[I],
  607		SubOp =.. [Op,Ds,V1i,V2i,VRi], SubOp,   % SubOp(Ds,VS1,VS2,VSR)
  608		do_subarrays(I,Op,Ds,V1,V2,VR)
  609	 ;	true
  610	).
  611
  612apply_subarrays(N,Op,Ds,V,F,VR) :-
  613	(N>0
  614	 ->	I is N-1,
  615		(number(V) -> Vi=V ; []([I],V,Vi)),  % index_arr(V[I],Vi)),  %  VS is V[I]),  % simple "broadcast" 
  616		[]([I],VR,VRi),  % index_arr(VR[I],VRi),  % VSR is VR[I],
  617		SubOp =.. [Op,Ds,Vi,F,VRi], SubOp,          % SubOp(Ds,VS1,VS2,VSR)
  618		apply_subarrays(I,Op,Ds,V,F,VR)
  619	 ;	true
  620	).
  621
  622acc_subarray(N,Op,Ds,V,Acc,R) :-
  623	(N>0
  624	 ->	I is N-1,
  625		(number(V) -> Vi=V ; []([I],V,Vi)),        % simple "broadcst" 
  626		SubOp =.. [Op,Ds,Vi,Acc,AccR], SubOp,      % SubOp(Ds,VS1,VS2,VSR)
  627		acc_subarray(I,Op,Ds,V,AccR,R)
  628	 ;	R=Acc
  629	).
  630
  631%
  632% arithmetic --> clpBNR constraint
  633%
  634constrain(error(instantiation_error,_), Goal) :-
  635	current_module(clpBNR),            % clpBNR constraints available?
  636	Mod=clpBNR, Mod:constrain_(Goal),  % can't use module name directly as it invalidates current_module test
  637	!.
  638constrain(Err, _Goal) :-
  639	throw(Err).
  640
  641%
  642% Compiled arithmetic 
  643%
  644:- set_prolog_flag(optimise,true).  645
  646matrix_iterate(C,0,Rows,NxtC,Rows) :- !, C>0,
  647	NxtC is C-1.
  648matrix_iterate(C,R,_Rows,C,NxtR) :- R>0,
  649	NxtR is R-1.	
  650
  651add(N1,N2,R) :- R is  N1+N2.
  652sub(N1,N2,R) :- R is  N1-N2.
  653mul(N1,N2,R) :- R is  N1*N2.
  654mml(N1,N2,R) :- R is -N1*N2.
  655div(N1,N2,R) :- R is  N1/N2.
  656pow(N1,N2,R) :- R is  N1**N2.
  657max(N1,N2,R) :- R is  max(N1,N2).
  658min(N1,N2,R) :- R is  min(N1,N2).
  659mul3(N1,N2,N3,R) :- R is N1*N2*N3.
  660inv(N,R)     :- N=\=0, R is 1/N.
  661det2x2(A0_0,A0_1,A1_0,A1_1,R) :- R is A0_0*A1_1 - A0_1*A1_0.
  662acc_det(Acc,El,I1,MDet,NxtA)  :- NxtA is Acc + El*I1*MDet