source: trunk/sources/ALGLIB/lda.cs @ 2636

/*************************************************************************
Copyright (c) 2008, Sergey Bochkanov (ALGLIB project).

>>> SOURCE LICENSE >>>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation (www.fsf.org); either version 2 of the 
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

A copy of the GNU General Public License is available at
http://www.fsf.org/licensing/licenses

>>> END OF LICENSE >>>
*************************************************************************/

using System;

namespace alglib
{
    public class lda
    {
        /*************************************************************************
        Multiclass Fisher LDA

        Subroutine finds coefficients of linear combination which optimally separates
        training set on classes.

        INPUT PARAMETERS:
            XY          -   training set, array[0..NPoints-1,0..NVars].
                            First NVars columns store values of independent
                            variables, next column stores number of class (from 0
                            to NClasses-1) which dataset element belongs to. Fractional
                            values are rounded to nearest integer.
            NPoints     -   training set size, NPoints>=0
            NVars       -   number of independent variables, NVars>=1
            NClasses    -   number of classes, NClasses>=2


        OUTPUT PARAMETERS:
            Info        -   return code:
                            * -4, if internal EVD subroutine hasn't converged
                            * -2, if there is a point with class number
                                  outside of [0..NClasses-1].
                            * -1, if incorrect parameters was passed (NPoints<0,
                                  NVars<1, NClasses<2)
                            *  1, if task has been solved
                            *  2, if there was a multicollinearity in training set,
                                  but task has been solved.
            W           -   linear combination coefficients, array[0..NVars-1]

          -- ALGLIB --
             Copyright 31.05.2008 by Bochkanov Sergey
        *************************************************************************/
        public static void fisherlda(ref double[,] xy,
            int npoints,
            int nvars,
            int nclasses,
            ref int info,
            ref double[] w)
        {
            double[,] w2 = new double[0,0];
            int i_ = 0;

            fisherldan(ref xy, npoints, nvars, nclasses, ref info, ref w2);
            if( info>0 )
            {
                w = new double[nvars-1+1];
                for(i_=0; i_<=nvars-1;i_++)
                {
                    w[i_] = w2[i_,0];
                }
            }
        }


        /*************************************************************************
        N-dimensional multiclass Fisher LDA

        Subroutine finds coefficients of linear combinations which optimally separates
        training set on classes. It returns N-dimensional basis whose vector are sorted
        by quality of training set separation (in descending order).

        INPUT PARAMETERS:
            XY          -   training set, array[0..NPoints-1,0..NVars].
                            First NVars columns store values of independent
                            variables, next column stores number of class (from 0
                            to NClasses-1) which dataset element belongs to. Fractional
                            values are rounded to nearest integer.
            NPoints     -   training set size, NPoints>=0
            NVars       -   number of independent variables, NVars>=1
            NClasses    -   number of classes, NClasses>=2


        OUTPUT PARAMETERS:
            Info        -   return code:
                            * -4, if internal EVD subroutine hasn't converged
                            * -2, if there is a point with class number
                                  outside of [0..NClasses-1].
                            * -1, if incorrect parameters was passed (NPoints<0,
                                  NVars<1, NClasses<2)
                            *  1, if task has been solved
                            *  2, if there was a multicollinearity in training set,
                                  but task has been solved.
            W           -   basis, array[0..NVars-1,0..NVars-1]
                            columns of matrix stores basis vectors, sorted by
                            quality of training set separation (in descending order)

          -- ALGLIB --
             Copyright 31.05.2008 by Bochkanov Sergey
        *************************************************************************/
        public static void fisherldan(ref double[,] xy,
            int npoints,
            int nvars,
            int nclasses,
            ref int info,
            ref double[,] w)
        {
            int i = 0;
            int j = 0;
            int k = 0;
            int m = 0;
            double v = 0;
            int[] c = new int[0];
            double[] mu = new double[0];
            double[,] muc = new double[0,0];
            int[] nc = new int[0];
            double[,] sw = new double[0,0];
            double[,] st = new double[0,0];
            double[,] z = new double[0,0];
            double[,] z2 = new double[0,0];
            double[,] tm = new double[0,0];
            double[,] sbroot = new double[0,0];
            double[,] a = new double[0,0];
            double[,] xyproj = new double[0,0];
            double[,] wproj = new double[0,0];
            double[] tf = new double[0];
            double[] d = new double[0];
            double[] d2 = new double[0];
            double[] work = new double[0];
            int i_ = 0;

            
            //
            // Test data
            //
            if( npoints<0 | nvars<1 | nclasses<2 )
            {
                info = -1;
                return;
            }
            for(i=0; i<=npoints-1; i++)
            {
                if( (int)Math.Round(xy[i,nvars])<0 | (int)Math.Round(xy[i,nvars])>=nclasses )
                {
                    info = -2;
                    return;
                }
            }
            info = 1;
            
            //
            // Special case: NPoints<=1
            // Degenerate task.
            //
            if( npoints<=1 )
            {
                info = 2;
                w = new double[nvars-1+1, nvars-1+1];
                for(i=0; i<=nvars-1; i++)
                {
                    for(j=0; j<=nvars-1; j++)
                    {
                        if( i==j )
                        {
                            w[i,j] = 1;
                        }
                        else
                        {
                            w[i,j] = 0;
                        }
                    }
                }
                return;
            }
            
            //
            // Prepare temporaries
            //
            tf = new double[nvars-1+1];
            work = new double[Math.Max(nvars, npoints)+1];
            
            //
            // Convert class labels from reals to integers (just for convenience)
            //
            c = new int[npoints-1+1];
            for(i=0; i<=npoints-1; i++)
            {
                c[i] = (int)Math.Round(xy[i,nvars]);
            }
            
            //
            // Calculate class sizes and means
            //
            mu = new double[nvars-1+1];
            muc = new double[nclasses-1+1, nvars-1+1];
            nc = new int[nclasses-1+1];
            for(j=0; j<=nvars-1; j++)
            {
                mu[j] = 0;
            }
            for(i=0; i<=nclasses-1; i++)
            {
                nc[i] = 0;
                for(j=0; j<=nvars-1; j++)
                {
                    muc[i,j] = 0;
                }
            }
            for(i=0; i<=npoints-1; i++)
            {
                for(i_=0; i_<=nvars-1;i_++)
                {
                    mu[i_] = mu[i_] + xy[i,i_];
                }
                for(i_=0; i_<=nvars-1;i_++)
                {
                    muc[c[i],i_] = muc[c[i],i_] + xy[i,i_];
                }
                nc[c[i]] = nc[c[i]]+1;
            }
            for(i=0; i<=nclasses-1; i++)
            {
                v = (double)(1)/(double)(nc[i]);
                for(i_=0; i_<=nvars-1;i_++)
                {
                    muc[i,i_] = v*muc[i,i_];
                }
            }
            v = (double)(1)/(double)(npoints);
            for(i_=0; i_<=nvars-1;i_++)
            {
                mu[i_] = v*mu[i_];
            }
            
            //
            // Create ST matrix
            //
            st = new double[nvars-1+1, nvars-1+1];
            for(i=0; i<=nvars-1; i++)
            {
                for(j=0; j<=nvars-1; j++)
                {
                    st[i,j] = 0;
                }
            }
            for(k=0; k<=npoints-1; k++)
            {
                for(i_=0; i_<=nvars-1;i_++)
                {
                    tf[i_] = xy[k,i_];
                }
                for(i_=0; i_<=nvars-1;i_++)
                {
                    tf[i_] = tf[i_] - mu[i_];
                }
                for(i=0; i<=nvars-1; i++)
                {
                    v = tf[i];
                    for(i_=0; i_<=nvars-1;i_++)
                    {
                        st[i,i_] = st[i,i_] + v*tf[i_];
                    }
                }
            }
            
            //
            // Create SW matrix
            //
            sw = new double[nvars-1+1, nvars-1+1];
            for(i=0; i<=nvars-1; i++)
            {
                for(j=0; j<=nvars-1; j++)
                {
                    sw[i,j] = 0;
                }
            }
            for(k=0; k<=npoints-1; k++)
            {
                for(i_=0; i_<=nvars-1;i_++)
                {
                    tf[i_] = xy[k,i_];
                }
                for(i_=0; i_<=nvars-1;i_++)
                {
                    tf[i_] = tf[i_] - muc[c[k],i_];
                }
                for(i=0; i<=nvars-1; i++)
                {
                    v = tf[i];
                    for(i_=0; i_<=nvars-1;i_++)
                    {
                        sw[i,i_] = sw[i,i_] + v*tf[i_];
                    }
                }
            }
            
            //
            // Maximize ratio J=(w'*ST*w)/(w'*SW*w).
            //
            // First, make transition from w to v such that w'*ST*w becomes v'*v:
            //    v  = root(ST)*w = R*w
            //    R  = root(D)*Z'
            //    w  = (root(ST)^-1)*v = RI*v
            //    RI = Z*inv(root(D))
            //    J  = (v'*v)/(v'*(RI'*SW*RI)*v)
            //    ST = Z*D*Z'
            //
            //    so we have
            //
            //    J = (v'*v) / (v'*(inv(root(D))*Z'*SW*Z*inv(root(D)))*v)  =
            //      = (v'*v) / (v'*A*v)
            //
            if( !sevd.smatrixevd(st, nvars, 1, true, ref d, ref z) )
            {
                info = -4;
                return;
            }
            w = new double[nvars-1+1, nvars-1+1];
            if( (double)(d[nvars-1])<=(double)(0) | (double)(d[0])<=(double)(1000*AP.Math.MachineEpsilon*d[nvars-1]) )
            {
                
                //
                // Special case: D[NVars-1]<=0
                // Degenerate task (all variables takes the same value).
                //
                if( (double)(d[nvars-1])<=(double)(0) )
                {
                    info = 2;
                    for(i=0; i<=nvars-1; i++)
                    {
                        for(j=0; j<=nvars-1; j++)
                        {
                            if( i==j )
                            {
                                w[i,j] = 1;
                            }
                            else
                            {
                                w[i,j] = 0;
                            }
                        }
                    }
                    return;
                }
                
                //
                // Special case: degenerate ST matrix, multicollinearity found.
                // Since we know ST eigenvalues/vectors we can translate task to
                // non-degenerate form.
                //
                // Let WG is orthogonal basis of the non zero variance subspace
                // of the ST and let WZ is orthogonal basis of the zero variance
                // subspace.
                //
                // Projection on WG allows us to use LDA on reduced M-dimensional
                // subspace, N-M vectors of WZ allows us to update reduced LDA
                // factors to full N-dimensional subspace.
                //
                m = 0;
                for(k=0; k<=nvars-1; k++)
                {
                    if( (double)(d[k])<=(double)(1000*AP.Math.MachineEpsilon*d[nvars-1]) )
                    {
                        m = k+1;
                    }
                }
                System.Diagnostics.Debug.Assert(m!=0, "FisherLDAN: internal error #1");
                xyproj = new double[npoints-1+1, nvars-m+1];
                blas.matrixmatrixmultiply(ref xy, 0, npoints-1, 0, nvars-1, false, ref z, 0, nvars-1, m, nvars-1, false, 1.0, ref xyproj, 0, npoints-1, 0, nvars-m-1, 0.0, ref work);
                for(i=0; i<=npoints-1; i++)
                {
                    xyproj[i,nvars-m] = xy[i,nvars];
                }
                fisherldan(ref xyproj, npoints, nvars-m, nclasses, ref info, ref wproj);
                if( info<0 )
                {
                    return;
                }
                blas.matrixmatrixmultiply(ref z, 0, nvars-1, m, nvars-1, false, ref wproj, 0, nvars-m-1, 0, nvars-m-1, false, 1.0, ref w, 0, nvars-1, 0, nvars-m-1, 0.0, ref work);
                for(k=nvars-m; k<=nvars-1; k++)
                {
                    for(i_=0; i_<=nvars-1;i_++)
                    {
                        w[i_,k] = z[i_,k-(nvars-m)];
                    }
                }
                info = 2;
            }
            else
            {
                
                //
                // General case: no multicollinearity
                //
                tm = new double[nvars-1+1, nvars-1+1];
                a = new double[nvars-1+1, nvars-1+1];
                blas.matrixmatrixmultiply(ref sw, 0, nvars-1, 0, nvars-1, false, ref z, 0, nvars-1, 0, nvars-1, false, 1.0, ref tm, 0, nvars-1, 0, nvars-1, 0.0, ref work);
                blas.matrixmatrixmultiply(ref z, 0, nvars-1, 0, nvars-1, true, ref tm, 0, nvars-1, 0, nvars-1, false, 1.0, ref a, 0, nvars-1, 0, nvars-1, 0.0, ref work);
                for(i=0; i<=nvars-1; i++)
                {
                    for(j=0; j<=nvars-1; j++)
                    {
                        a[i,j] = a[i,j]/Math.Sqrt(d[i]*d[j]);
                    }
                }
                if( !sevd.smatrixevd(a, nvars, 1, true, ref d2, ref z2) )
                {
                    info = -4;
                    return;
                }
                for(k=0; k<=nvars-1; k++)
                {
                    for(i=0; i<=nvars-1; i++)
                    {
                        tf[i] = z2[i,k]/Math.Sqrt(d[i]);
                    }
                    for(i=0; i<=nvars-1; i++)
                    {
                        v = 0.0;
                        for(i_=0; i_<=nvars-1;i_++)
                        {
                            v += z[i,i_]*tf[i_];
                        }
                        w[i,k] = v;
                    }
                }
            }
            
            //
            // Post-processing:
            // * normalization
            // * converting to non-negative form, if possible
            //
            for(k=0; k<=nvars-1; k++)
            {
                v = 0.0;
                for(i_=0; i_<=nvars-1;i_++)
                {
                    v += w[i_,k]*w[i_,k];
                }
                v = 1/Math.Sqrt(v);
                for(i_=0; i_<=nvars-1;i_++)
                {
                    w[i_,k] = v*w[i_,k];
                }
                v = 0;
                for(i=0; i<=nvars-1; i++)
                {
                    v = v+w[i,k];
                }
                if( (double)(v)<(double)(0) )
                {
                    for(i_=0; i_<=nvars-1;i_++)
                    {
                        w[i_,k] = -1*w[i_,k];
                    }
                }
            }
        }
    }
}