Picat machine learning

knn.pi

% Machine learning toolbox
module ml.
import util.

% calculates only half of the distances: symmetry
pdist(M) = M3 =>
    M3 = new_array(M.length, M.length),
    foreach(I in 1..M.length, J in I..M.length)
        M3[I,J] = sum([(M[I,K]-M[J,K])**2 : K in 1..M[1].length]),
		M3[J,I] = M3[I,J] 
    end.

% Calculates pairwise distances between each row in the first matrix vs. the 
% rows in the second matrix
pdist2(M1,M2) = M3 =>
    M3 = new_array(M1.length, M2.length),
    foreach(I in 1..M1.length, J in 1..M2.length)
        M3[I,J] = sum([(M1[I,K]-M2[J,K])**2 : K in 1..M1[1].length])
    end.

% Small utility printing a matrix
printmat(M) => foreach (X in 1 .. M.length) writeln(M[X]) end.

% Get indices for permutation turning the array into a sorted array
grade_up(L) = [ I : (E,I) in sort([(E,I) : {E,I} in zip(L,1..L.length)])].
% Same as grade_up, but also get the values
sortind(L)  = [ (I,E) : (E,I) in sort([(E,I) : {E,I} in zip(L,1..L.length)])].

% Center data around column averages
center(Dat) = Centered =>
	Dat_T = transpose(Dat),
	M = Dat[1].length,
	A = [avg(to_list(Dat_T[I])): I in 1 .. M],
	Centered = center(Dat,A).

% Center a data set around A.
center(Dat,A) = Centered =>
	[N,M] = [Dat.length, Dat[1].length],
	Centered = new_array(N, M),
	foreach (I in 1..N , J in 1..M)
		Centered[I,J] = Dat[I,J] - A[J]
	end.

% Sample Covariance
cov(Dat) = Cov =>
	Cen = center(Dat),
	[N,M] = [length(Dat),length(Dat[1])],
	Cov = matrix_multi(transpose(Cen),Cen),
	foreach(I in 1..M, J in 1..M)
		Cov[I,J]:=Cov[I,J]/(N-1)
	end.

% Identity matrix of size(N,N)
idmat(N) = I =>
	I = new_array(N,N),
	foreach ( K in 1..N, L in 1..N) 
		I[K,L]= cond( K==L, 1,0)
	end.

% Reduce [L | R] to row-echolon form
gaussjordan(L,R) = [LL,RR] =>
	[N,M1,M2] = [length(L),length(L[1]),length(R[1])],
	LL = new_array(N,M1),
	foreach ( I in 1..N, J in 1..M1) LL[I,J]=L[I,J] end,
	RR = new_array(N,M2),
	foreach ( I in 1..N, J in 1..M2) RR[I,J]=R[I,J] end,
	% Make Echelon form.
	foreach(J in 1..N, I in J+1..N) 
		C = LL[I,J]/LL[J,J],
		foreach(K in 1 .. M1)
			LL[I,K]:=LL[I,K]-C*LL[J,K]
		end,
		foreach(K in 1 .. M2)
			RR[I,K]:=RR[I,K]-C*RR[J,K]
		end
	end.

% Do back-substitution to obtain solutions for righ hand sides.
backsubst(L,R) = RR => 
	[N,M2] = [length(L),length(R[1])],
	RR = new_array(N,M2),
	foreach(I in 1..N, J in 1..M2) RR[I,J] = R[I,J] end,
	% Do backsubstitution, start from the bottom row, and work up.
	foreach (I in N .. -1 .. 1)
		foreach(K in 1..M2)
			Sum = 0,
			foreach(J in I+1 .. N)
				Sum := Sum + L[I,J]*RR[J,K]
			end,
			RR[I,K] := (RR[I,K]-Sum)/L[I,I]
		end
	end.	

% Matrix inverse via Gauss-Jordan.
inv(M) = Id =>
	N = length(M),
	MM = new_array(N,N), % Copy to avoid overwriting original
	foreach(I in 1..N,J in 1..N)
		MM[I,J]=M[I,J]
	end,
	[L,R] = gaussjordan(MM,idmat(N)),
	Id = backsubst(L,R).
	
table
mvnpdf(X,M,S) = P => % Note, det(Sigma) not included
	% X should be 2D, M 1D, S 2D.
	K = length(X[1]),
	Mul = 1/(sqrt((2 * math.pi ) **K * det(S))),
	Xc  = center(X,M),
	E  = matrix_multi(matrix_multi(Xc,inv(S)),transpose(Xc)),
	P  = Mul * exp(-1/2* E[1,1]).

normpdf(X,M,S) = P =>
	P = 1/(S*sqrt(2*math.pi)) * exp( -(X -M)**2/(2*S**2)).

% Calculate the Levi-Civita permutation symbol aka. permutation parity.
table(+,-)
permpar(P,S),not permutation(P,1..length(P)) => S=0.	% not a permutation
permpar([],S) => S=1.					% Empty set is always sorted.
permpar([1|T],S) =>					% 1 is already in front, continue
	TT = [T[K]-1 : K in 1..T.len],
	permpar(TT,S).
permpar([H|T],S) =>					% All other cases:
	X=find_first_of(zip(T,1..T.len),{1,X}),		%  Find 1
	TT = [ T[K]-1: K in 1..T.len],			%  subtract 1 (for next iteration)
	TT[X] := H-1,					%  swap 1 and H (-1 because after prev)
	permpar(TT,S2),					%  Continue on the rest.
	S = -S2.					%  Swap happened: flip sign.

% Calculate the determinant of a matrix
det(A) = D,length(A)=length(A[1]) =>
	N = A.length,
	Per = permutations(1..N),
	NP = Per.len,
	Par = [S: K in 1 .. NP, permpar(Per[K],S)],
	Contrib = new_list(NP),
	foreach (K in 1..NP)
		Sig = Per[K],
		Comp = new_list(N),
		foreach (I in 1..N)
			Id = Sig[I],
			Comp[I] = A[I,Id]
		end,
		Contrib[K] = Par[K]*prod(Comp)
	end,
	D = sum(Contrib)

ml.pi

% Machine learning toolbox
module ml.
import util.

% calculates only half of the distances: symmetry
pdist(M) = M3 =>
    M3 = new_array(M.length, M.length),
    foreach(I in 1..M.length, J in I..M.length)
        M3[I,J] = sum([(M[I,K]-M[J,K])**2 : K in 1..M[1].length]),
		M3[J,I] = M3[I,J] 
    end.

% Calculates pairwise distances between each row in the first matrix vs. the 
% rows in the second matrix
pdist2(M1,M2) = M3 =>
    M3 = new_array(M1.length, M2.length),
    foreach(I in 1..M1.length, J in 1..M2.length)
        M3[I,J] = sum([(M1[I,K]-M2[J,K])**2 : K in 1..M1[1].length])
    end.

% Small utility printing a matrix
printmat(M) => foreach (X in 1 .. M.length) writeln(M[X]) end.

% Get indices for permutation turning the array into a sorted array
grade_up(L) = [ I : (E,I) in sort([(E,I) : {E,I} in zip(L,1..L.length)])].
% Same as grade_up, but also get the values
sortind(L)  = [ (I,E) : (E,I) in sort([(E,I) : {E,I} in zip(L,1..L.length)])].

% Center data around column averages
center(Dat) = Centered =>
	Dat_T = transpose(Dat),
	M = Dat[1].length,
	A = [avg(to_list(Dat_T[I])): I in 1 .. M],
	Centered = center(Dat,A).

% Center a data set around A.
center(Dat,A) = Centered =>
	[N,M] = [Dat.length, Dat[1].length],
	Centered = new_array(N, M),
	foreach (I in 1..N , J in 1..M)
		Centered[I,J] = Dat[I,J] - A[J]
	end.

% Sample Covariance
cov(Dat) = Cov =>
	Cen = center(Dat),
	[N,M] = [length(Dat),length(Dat[1])],
	Cov = matrix_multi(transpose(Cen),Cen),
	foreach(I in 1..M, J in 1..M)
		Cov[I,J]:=Cov[I,J]/(N-1)
	end.

% Identity matrix of size(N,N)
idmat(N) = I =>
	I = new_array(N,N),
	foreach ( K in 1..N, L in 1..N) 
		I[K,L]= cond( K==L, 1,0)
	end.

% Reduce [L | R] to row-echolon form
gaussjordan(L,R) = [LL,RR] =>
	[N,M1,M2] = [length(L),length(L[1]),length(R[1])],
	LL = new_array(N,M1),
	RR = new_array(N,M2),
	% Make Echelon form.
	foreach(J in 1..N, I in J+1..N) 
		C = LL[I,J]/LL[J,J],
		foreach(K in 1 .. M1)
			LL[I,K]:=LL[I,K]-C*LL[J,K]
		end,
		foreach(K in 1 .. M2)
			RR[I,K]:=RR[I,K]-C*RR[J,K]
		end
	end.

% Do back-substitution to obtain solutions for righ hand sides.
backsubst(L,R) = [RR] => 
	[N,M2] = [length(L),length(R[1])],
	RR = new_array(N,M2),
	% Do backsubstitution, start from the bottom row, and work up.
	foreach (I in N .. -1 .. 1)
		foreach(K in 1..M2)
			Sum = 0,
			foreach(J in I+1 .. N)
				Sum := Sum + L[I,J]*RR[J,K]
			end,
			RR[I,K] := (RR[I,K]-Sum)/L[I,I]
		end
	end.	

% Matrix inverse via Gauss-Jordan.
inv(M) = I =>
	N = length(M),
	MM = new_array(N,N), % Copy to avoid overwriting original
	foreach(I in 1..N,J in 1..N)
		MM[I,J]=M[I,J]
	end,
	[L,R] = gaussjordan(MM,idmat(N)),
	[_,I] = backsubst(L,R).
	
table
mvnpdf(X,M,S) = P => % Note, det(Sigma) not included
	% X should be 2D, M 1D, S 2D.
	K = length(X[1]),
	Mul = 1/(sqrt(2*math.pi)**K),
	Xc  = center(X,M),
	E  = matrix_multi(matrix_multi(Xc,inv(S)),transpose(Xc)),
	P  = Mul * exp(-1/2* E[1,1]).

normpdf(X,M,S) = P =>
	P = 1/(S*sqrt(2*math.pi)) * exp( -(X -M)**2/(2*S**2)).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Incomplete stuff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%permsign(P,S),P=1..length(P) => S=1.
%permsign(P,S),not permutation(P,1..length(P)) => S=0.
%permsign(P,S) =>
%	S = prod(findall(length(C),cycles(P,C))).
%
%repext(Pat,Len) = Res =>
%	N = length(Pat),
%	I = ceiling(Len/N),	% "Integer" part (in terms of N)
%	Res = slice(flatten([Pat: _ in 1..I]),1,Len).
%
%det(A,D) =>
%	throw($error(nyi)),
%	N = A.length,
%	% permutations(L) yields 1 0 1 0, 0 1 0 1, ... as 
%	Per = findall(P,permutation(1..N,P)),
%
%	% Per = permutations(1..N) % Problem: does not always follow the same
%	% structure (concerning permutation parity)
%	% Note, permutation does neither...
%
%	% Todo: implement permutation parity check.
%
%	Par = repext([1,-1,-1,1],factorial(N)),
%	Contrib = new_list(length(Par)),
%	foreach (K in 1..length(Par))
%		Sig = Per[K],
%		Comp = new_list(N),
%		foreach (I in 1..N)
%			Id = Sig[I],
%	    	Comp[I] = A[I,Id]
%		end,
%		Contrib[K] = Par[K]*prod(Comp)
%	end,
%	D = [sum(Contrib),Per,Par].
%