function [egv,discrim,coeff] = fda(X,gnum,reg)

// Fisher discriminant analysis with regularization

// X ... data matrix arranged as columns - variables and rows - observations
// across individual groups
// gnum ... vector of group marks, e.g. gnum = [1 1 1 1 1 2 2 2 2 2 2 2 ...
// 3 3 3 3]' of the length corresponding to number of rows of X
// reg 0 ... without regularization (default)
//     any numerical value ... parameter used within the regularization

// coeff ... discriminant coefficients - eigenvectors
// discrim ... discriminant variables
// egv ... eigenvalues

// Literature: Krzanowski, W.J. (2000). Principles of multivariate analysis,
// Univerity Press: Oxford
// Golub, G.H. and Van Loan, C.F. (1996). Matrix computations, The Johns
// Hopkins University Press: Baltimore

// Copyright (C): Lukáš Malec, Lenka Netušilová and Hana Tvrdíková;
// University of Business in Prague, Czech Republic
// Date: $2015-10-23$

[lhs,rhs]=argn(0)

[n,m] = size(X), g = max(unique(gnum))

Xm = mean(X,'r')
X = X-Xm(ones(n,1),:)
T = X'*X

W = zeros(T)

for i = 1:g
    f = find(gnum==i)
    nf = length(f)
    G = X(f,:)
    Xm = mean(G,'r')
    G = G-Xm(ones(nf,1),:)

    W = W+G'*G

end

B=T-W

if rhs==2
    R = chol(W)
  else R = chol(W+reg*eye(W))

end

D = R'\B/R
D = (D+D')/2
[evec,ev] = spec(D)
v = R\evec

[egv,eo] = gsort(diag(ev),'g','d')

v = v(:,eo)
vscal = diag((v'*W*v))' ./(n-g)
vscal(vscal<=0) = 1

coeff = v./repmat(sqrt(vscal),size(v,1),1)

for k = 1:m
    p = find((coeff(:,k)<0))
  if length(p) > m/2 then
    coeff(:,k) = -coeff(:,k)

  end
end

discrim = X*coeff

endfunction
