splitpointcomputation.h

Go to the documentation of this file.

/*

Copyright (c) 2003, Cornell University
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

   - Redistributions of source code must retain the above copyright notice,
       this list of conditions and the following disclaimer.
   - Redistributions in binary form must reproduce the above copyright
       notice, this list of conditions and the following disclaimer in the
       documentation and/or other materials provided with the distribution.
   - Neither the name of Cornell University nor the names of its
       contributors may be used to endorse or promote products derived from
       this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.

*/

// -*- C++ -*-

/* Implement functions for split point computation together with
   necessary auxiliary functions
*/

#if !defined _CLUS_SPLITPOINTCOMPUTATION_H_
#define _CLUS_SPLITPOINTCOMPUTATION_H_

#include "general.h"
#include <math.h>
#include <iostream>
#include "trislv.h"
#include "fmat.h"
#include "cholesky.h"

#include <gsl/gsl_sf_gamma.h> // for incomplete regularized beta function
#include <gsl/gsl_errno.h>

#include <stdlib.h>

using namespace TNT;
using namespace std;

namespace CLUS
{

/// Determines in log time if a point is in a set.
/// The set is a vector of sorted values.
template< class T>
bool IsPointInSet(T value, Vector<T> set)
{
    int l=0, r=set.dim()-1;

    assert(r>=0);

    while ( l<=r )
    {
        int m=(l+r)/2;
        T vm=set
             [m];

        if ( vm ==value )
            return true;

        if ( vm < value )
            l=m+1;
        else
            r=m-1;
    }
    return false;
}


/** Compute P[X>=val] for X~Binomial(N,p).
    This is exactly the normalized incomplete beta function (see Eric's encyclopedia) .
 
    We use gsl to get this function. The result is IB_p(val,N-val-1.0)
*/
double PValueBinomialDistribution(double N, double p, double val)
{

    double mean=N*p;
    double std=sqrt(N*p*(1-p)); // standard deviantion

    // cout << "mean=" << mean << " std=" << std << endl;

    // if val is not withing 10 standard deviations do not bother with the approximation
    if (val<mean-10*std)
        return 1.0;

    if (val>mean+10*std)
        return 0.0;

    gsl_sf_result result;

    int status=gsl_sf_beta_inc_e(val+1,N-val,p,&result);

    switch(status)
    {
    case 0:
        // everything is fine. Do nothing
        break;

    default:
        // deal with the error
        // use the normal approximation to the binomial
        {
            cout << "gsl_sf_beta_inc_e failed for N=" << N << " val=" << val << " p=" << p << " .Approximating with normal" << endl;
            return .5*(1.0+erf( (mean-val)/(std*sqrt(2.0)) ));
        }
        break;
    }

    return result.val;

}

/** Computes int_{x>=eta} N(mu,var) dx */
double PValueNormalDistribution(double mu, double sigma, double eta)
{
    if (sigma==0.0)
        if (mu>eta)
            return 1.0;
        else
            if (mu==eta)
                return .5;
            else
                return 0.0;
    return .5*(1.0+erf( (mu-eta)/sqrt(2*pow2(sigma)) ));
}

/** Computes int_{n'*(x-xc)>=0} N(mu,Sigma) dx */
double PValueNormalDistribution(const Vector<double> mu,
                                const Fortran_Matrix<double> cholSigma,
                                const Vector<double> n,
                                const Vector<double> xc)
{
    /* Crafted after pvalue.m and thing.m (Matlab prototyping)
       The final formula is: 
                               1                   /       / sigma v \ \
           P=------------------------------------- | 1+ Erf| ------- | |
      2 sigma det(cholSigma) Sqrt[ det(S) ] \       \ Sqrt[2] / /
      
       where:
           v=n'*(mu-xc)/norm(n),
           sigma=Sqrt[ s-w'S^-1w ],
                                          / s w'\           
    (M'*cholSigma*cholSigma'*M)^-1=|     |
                                          \ w S /
       adn M is picked such that  M'n=e1

    */

    int d=n.dim(); // dimentionality of the space

    // compute v
    double v=dot_prod(n,mu-xc)/sqrt(dot_prod(n,n));

    // build initial M
    Fortran_Matrix<double> M(d,d); // already initialized to 0
    /* pick robustly the other d-1 vectors, eliminate ei, i=max(abs(n))
       since is the most "paralel" with n */
    // find i first
    int i=0;
    double max=fabs(n[i]);
    for(int j=1; j<d; j++)
    {
        double aux=fabs(n[j]);
        if (aux>max)
        {
            max=aux;
            i=j;
        }
    }
    i++; // adsust for the Fortran indexing
    // put n in the first line and 1 on the diagonal
    for (int j=1; j<=d; j++)
    {
        M(j,j)=1.0;
        M(j,1)=n[j-1];
    }
    // put e1 in row i
    if (i!=1)
    {
        M(i,i)=0.0;
        M(1,i)=1.0;
    }

    // Gram-Schmidt ortogonalize M
    for (int j=1; j<=d; j++)
    {
        for (int k=1; k<j; k++)
        {
            // mult=M(:,j)*M(:,k)'
            double mult=0.0;
            for (int l=1; l<=d; l++)
                mult+=M(l,j)*M(l,k);
            // M(:,j)=M(:,j)-mult*M(:,k);
            for (int l=1; l<=d; l++)
                M(l,j)-=mult*M(l,k);
        }
        // compute norm(M(j,:))
        double norm=0.0;
        for (int k=1; k<=d; k++)
            norm+=pow2(M(k,j));
        norm=sqrt(norm);
        // normalizde M(j,:)
        for (int k=1; k<=d; k++)
            M(k,j)/=norm;
    }

    // compute cholSigma^-1 M into M
    for (int i=1; i<=d; i++)
    {
        // for each column of M solve cholSigma^-1 z in z
        for (int j=1; j<=d; j++)
        {
            double tmp=0.0;
            for (int k=1; k<j; k++)
                tmp+=cholSigma(j,k)*M(k,i);
            M(j,i)=(M(j,i)-tmp)/cholSigma(j,j);
        }
    }

    // compute s
    double s=0.0;
    for (int i=1; i<=d; i++)
        s+=pow2(M(i,1));

    // compute w
    Vector<double> w(d-1);
    for (int i=0; i<d-1; i++)
    {
        // w[i]=<z1,zi>
        double tmp=0.0;
        for (int j=1; j<=d; j++)
            tmp+=M(j,1)*M(j,i+2);
        w[i]=tmp;
    }

    // compute S
    Fortran_Matrix<double> S(d-1,d-1);
    for (int i=1; i<=d-1; i++)
        for (int j=i; j<=d-1; j++)
        {
            // multiply z(i+1) with z(j+1)
            double tmp=0.0;
            for (int k=1; k<=d; k++)
                tmp+=M(k,i+1)*M(k,j+1);
            S(i,j)=tmp;
        }

    Cholesky_upper_factorization(S,S);
    // compute cholS^-1 * w
    Vector< double > aux = Lower_triangular_solve(S,w);

    double sigma=s-dot_prod(aux,aux);

    sigma=sqrt(sigma);
    double det_cholS=1.0;
    for (int i=1; i<=d-1; i++)
        det_cholS*=S(i,i);
    double det_cholSigma=1.0;
    for (int i=1; i<=d; i++)
        det_cholSigma*=cholSigma(i,i);

    // we can compute now the PValue
    return 1/( 2*sigma*det_cholSigma*det_cholS )*( 1+erf(sigma*v/sqrt(2.0)) );
}

/** auxiliary type */
struct T_array
{
    int index;
    double crit;
};

/** auxiliary function to sort elements */
int compare_array_elements(const void* x, const void* y)
{
    const T_array ax = *( (T_array*)x );
    const T_array ay = *( (T_array*)y );
    if (ax.crit<ay.crit)
        return -1;
    else
        if (ax.crit == ay.crit)
            return 0;
        else
            return 1;
}

/** Computes /Delta g(T).
    C1 is cluster 1; S1 is split 1
    
    @param p11     P[x /in C1 ^ x /in S1]
    @param p_1     P[x /in C1]
    @param p1_     P[x /in S1]
*/
double BinaryGiniGain(double p11, double p_1, double p1_)
{
    return p1_>0.0 ? 2*pow2(p11-p_1*p1_)/( p1_*(1-p1_) ) : 0.0;
}

/** Computes the maximum gain in gini by splitting on a discrete variable and the actual split
    Split: return the best split here
    Return: the new gini
    
    If d_s_p1[i]=d_N[i]=0 we don't know anything about value i. To avoid biases we distribute
    this values among the two splits.
 
    Theorem 9.4 from Breiman et al. justifies the linear algorithm.
 */
double DiscreteGiniGain(Vector<double>& d_s_p1, Vector<double>& d_N, double N, double alpha_1,
                        Vector<int>& Split)
{

    int n=d_N.dim();

#if 0

    Vector<int> B_order(n);
    for(int i=0; i<n; i++)
        B_order[i]=i;
    // order values in vector B_order. We use bubble sort since is the simplest one to implement and acceptable
    bool flag;
    do
    {
        flag=false;
        for(int i=0; i<n-1; i++)
        {
            int i1=B_order[i];
            double v1=d_N[i1]>0 ? d_s_p1[i1]/d_N[i1] : -1.0;
            int i2=B_order[i+1];
            double v2=d_N[i2]>0 ? d_s_p1[i2]/d_N[i2] : -1.0;
            if (v1>v2)
            {
                B_order[i]=i2;
                B_order[i+1]=i1;
                flag=true;
            }
        }
    }
    while (flag);
#endif

    T_array* B_order=new T_array[n];

    for (int i=0; i<n; i++)
    {
        B_order[i].index=i;
        B_order[i].crit=d_N[i]>0 ? d_s_p1[i]/d_N[i] : -1.0;
        ;
    }

    qsort((void*)B_order, n, sizeof(T_array), compare_array_elements );

    //    cout << "B_order ";
    //     for (int i=0; i<n; i++)
    //       cout << " [" << B_order[i].index << "," << B_order[i].crit << "] ";
    //     cout << endl;


    // find the splitting point
    int pos=-1;  // the splitting point is between pos and pos+1
    double p11=0.0;
    double p1_=0.0;
    double maxgini=0.0;
    bool partonefirst=true; // if true partition one is at the left
    for (int i=0; i<n-1; i++)
    {
        p11+=d_s_p1[B_order[i].index]/N;
        p1_ +=1.0*d_N[B_order[i].index]/N;
        double gini=BinaryGiniGain(p11,alpha_1,p1_);

        //cout << "XXXXXX " << p11 << " " << p1_ << " " << gini << endl;

        if (p1_>=0.0 && p1_<=1.0 && gini>=maxgini)
        {
            maxgini=gini;
            pos=i;
            if ( p11<(p1_-p11) )
                partonefirst=false;
            //cout << "maxgini=" << maxgini << " i=" << i << endl;

        }
    }

    double boundary;
    if (pos==-1 || B_order[pos].crit==-1.0)
        return 0.0;
    else
        boundary=d_s_p1[B_order[pos].index]/d_N[B_order[pos].index];

    // cout << "boundary=" << boundary << " pos=" << pos << endl;

    int nr=0;
    int buff[MAX_CATEGORIES];
    for (int i=0; i<n; i++)
        if ( (d_N[i]==0.0 && rand()*1.0/RAND_MAX > .5) ||
                ((!partonefirst) && d_s_p1[i]/d_N[i]>boundary)
                && NonEqual(d_s_p1[i]/d_N[i],boundary) ||
                (partonefirst &&
                 (d_s_p1[i]/d_N[i]<boundary ||
                  !NonEqual(d_s_p1[i]/d_N[i],boundary) )) )
        {
            // cout << "Adding " << i << " " << nr << endl;
            buff[nr]=i;
            nr++;
        }

    //cout << "BUFFER ";
    //for (int i=0; i<nr; i++)
    //cout << buff[i] << " ";
    //cout << endl;

    // put the result in Split
    Split.newsize(nr);
    for (int i=0; i<nr; i++)
        Split[i]=buff[i];

    //cout << "Split ";
    // for (int i=0; i<Split.size(); i++)
    //cout << Split[i] << " ";
    //cout << " ----- " << nr << endl;

    return (nr>0 && nr<=n-1) ? maxgini : 0.0; // to avoid pseudo nonzero ginis
}

/** Sister function of DiscreteGiniGain. Instead of finding a split set it finds
    probabilities that a point belongs to the left set 
*/
double ProbabilisticDiscreteGiniGain(const Vector<double>& d_s_p1, const Vector<double>& d_N,
                                     double N, double alpha_1,
                                     Vector<double>& probSet)
{

    int n=d_N.dim();

    T_array* B_order=new T_array[n];

    for (int i=0; i<n; i++)
    {
        B_order[i].index=i;
        B_order[i].crit=d_N[i]>0.0 ? d_s_p1[i]/d_N[i] : -1.0;
        ;
    }

    qsort((void*)B_order, n, sizeof(T_array), compare_array_elements );

    // find the splitting point
    int pos=-1;  // the splitting point is between pos and pos+1
    double p11=0.0;
    double p1_=0.0;
    double maxgini=0.0;
    bool partonefirst=true; // if true partition one is at the left
    for (int i=0; i<n-1; i++)
    {
        p11+=d_s_p1[B_order[i].index]/N;
        p1_ +=1.0*d_N[B_order[i].index]/N;
        double gini=BinaryGiniGain(p11,alpha_1,p1_);

        if (p1_>=0.0 && p1_<=1.0 && gini>=maxgini)
        {
            maxgini=gini;
            pos=i;
            if ( p11<(p1_-p11) )
                partonefirst=false;
        }
    }

    // take the boundary to be the average of the ratios where maximum occurs
    double boundary;
    probSet.newsize(n);
    if (pos==-1 || B_order[pos].crit==-1.0)
    {
        for (int i=0; i<n; i++)
            probSet[i]=0.5;
        return 0.0;
    }
    else
        boundary=(d_s_p1[B_order[pos].index]/d_N[B_order[pos].index]+
                  d_s_p1[B_order[pos+1].index]/d_N[B_order[pos+1].index])/2.0;

    delete [] B_order;

    // for every element in the domain compute the probability that the
    // observed ratio is at the left of the boundary point assuming Binomial distribution
    for (int i=0; i<n; i++)
        if (d_N[i]==0.0)
            probSet[i]=.5;
        else
        {
            double p=d_s_p1[i]/d_N[i];
            probSet[i]=1.0-PValueBinomialDistribution(d_N[i], p, boundary*d_N[i]);
            //cout << "i:" << i << " s_p1=" << d_s_p1[i] << " d_N=" << d_N[i]
            //     << " probSet[]=" << probSet[i] << endl;

        }

    return maxgini;
}



/**
   Form equation 
   eta^2(1/var1 - 1/var2)-2eta(eta1/var1-eta2/var2)+eta1^2/var1-eta2^2/var2
   = 2ln(alpha1/alpha2)-ln(var1/var2)
   and solve it.

   @param alpha_1
   @param eta1
   @param var1
   @param alpha_2
   @param eta2
   @param var2
   @param whichSol           set to 0 if first order equation solved, 1 if first sol of second 
                             order equation and 2 for second solution. Is set to 3 if the default
   @return                   weighted mean of averages
*/
double UnidimensionalQDA(double alpha_1, double eta1, double var1,
                         double alpha_2, double eta2, double var2, int& whichSol)
{
    whichSol=4; // special solution: eta1 or eta2

    if (isnan(alpha_1+eta1+var1))
    {
        assert( !isnan(alpha_2+eta2+var2) );
        return eta2;
    }

    if (isnan(alpha_2+eta2+var2))
    {
        assert( !isnan(alpha_1+eta1+var1) );
        return eta1;
    }

    if (!NonEqual(eta1,eta2))
        return eta1;

    double a=1.0/var1-1.0/var2;
    double b=-2.0*(eta1/var1-eta2/var2);
    double c=pow2(eta1)/var1-pow2(eta2)/var2-2.0*log(alpha_1/alpha_2)/log(exp(1.0))+
             log(var1/var2)/log(exp(1.0));


    if (var1<=SMALL_POZ_VALUE || var2<=SMALL_POZ_VALUE)
    {
        whichSol=5;
        return alpha_1*eta1+alpha_2*eta2;
    }

    // Alin:fixthis
    goto default_solution;

    if ( fabs(a) < 1e-8*(1.0/var1+1.0/var2) )
    {
        // first order equation

        if (NonZero(b))
        {
            double eta_l=-c/b;
            if ( (eta1-eta_l)*(eta2-eta_l)<0.0 )
            {
                whichSol=0;
                return eta_l;
            }
        } // otherwise return default solution; fallover
    }
    else
    {
        // second order equation
        double disc=pow2(b)-4.0*a*c;
        if (disc<0.0) // no solution, return default
            goto default_solution;
        else
        {
            // the two solutions
            double eta_1=(-b+sqrt(disc))/(2.0*a);
            double eta_2=(-b-sqrt(disc))/(2.0*a);
            // find the one in between
            if ( (eta1-eta_1)*(eta2-eta_1) < 0.0)
            {
                whichSol=1;
                return eta_1;
            }
            else
                if ( (eta1-eta_2)*(eta2-eta_2) < 0.0)
                {
                    whichSol=2;
                    return eta_2;
                } // otherwise return default solution; fallover
        }
    }

default_solution:
    whichSol=3; // default solution
    return (alpha_1*sqrt(var1)*eta1+alpha_2*eta2*sqrt(var2))/
           (alpha_1*sqrt(var1)+alpha_2*sqrt(var2));
    //return alpha_1*eta1+alpha_2*eta2;
}

/** Computes the variance of the split point.
   The prototype is in the Mathematica file StatisticsSplitPoint2.nb
   which has the computations of the variance using the delta method.
*/
double UnidimensionalQDAVariance(double n1, double m1, double v1,
                                 double n2, double m2, double v2, int whichSol)
{

    whichSol=6;

#ifdef DEBUG_PRINT

    cout << "UnidimentionalQDAVariance " << n1 << " " << n2 << " " << v1 << " " << v2;
    cout << " " << m1 << " " << m2 << "\t";
#endif

    double variance;

    if (v1<SMALL_POZ_VALUE  || v2<SMALL_POZ_VALUE )
    {
        if (n1<3 || n2<3)
            variance=pow2(m1-m2);
        else
            variance=SMALL_POZ_VALUE;
        return variance;
    }

    if (isnan(m1) || isnan(m2) || isnan(v1) || isnan(v2))
    {
        variance=LARGE_POZ_VALUE;
        return variance;
    }


    if ( IsEqual(v1,v2) && IsEqual(m1,m2) )
        variance = SMALL_POZ_VALUE;
    else
    {
        double CLog=log(n1/n2*sqrt(v2/v1));
        double SR=sqrt(v1*v2*( pow2(m1-m2)+2*(-v1+v2)*CLog ));
        double Q1=v1*(pow2(m1-m2)-v1+v2-2*(v1-2*v2)*CLog);
        double Q2=v2*(pow2(m1-m2)+v1-v2-2*(2*v1-v2)*CLog);
        double Q3=2*(pow(m2,4.0)-2*pow2(m2)*n2*v2+n2*pow2(v2));
        double Q4=2*(pow(m1,4.0)-2*pow2(m1)*n1*v1+n1*pow2(v1));
        double Q5=-m2*v1+m1*v2;

        // we compute all solutions even if not needed. any decent
        // compiler should get rid of unnecessary code

        double Sol1=
            ( v1*pow2(v2) ) / ( n1*pow2(v1-v2) ) * pow2(-1+(m1-m2)*v1/SR)+
            ( v2*pow2(v1) ) / ( n2*pow2(v1-v2) ) * pow2( 1-(m1-m2)*v2/SR)-
            Q4/( n1*(n1-1)*pow(v1-v2,4.0) ) * pow2(  Q5-SR+(v1-v2)*( m2+Q2/SR) )-
            Q3/( n2*(n2-1)*pow(v1-v2,4.0) ) * pow2( -Q5+SR+(v1-v2)*(-m1+Q1/SR) );

        double Sol2=
            ( v1*pow2(v2) ) / ( n1*pow2(v1-v2) ) * pow2( 1+(m1-m2)*v1/SR)+
            ( v2*pow2(v1) ) / ( n2*pow2(v1-v2) ) * pow2(-1-(m1-m2)*v2/SR)-
            Q4/( n1*(n1-1)*pow(v1-v2,4.0) ) * pow2( Q5+SR-(v1-v2)*(-m2+Q2/SR) )-
            Q3/( n2*(n2-1)*pow(v1-v2,4.0) ) * pow2( Q5+SR+(v1-v2)*( m1+Q1/SR) );

        // solution for the first order equation. Here v1=v2
        double QL=v1*log(n1/n2)/pow2(m1-m2);
        double Sol0=2*pow2(m1-m2)*pow2(QL)*(1.0/n1+1.0/n2)+v1/(4*n1)*pow2(1-2*QL)+
                    +v1/(4*n2)*pow2(1+2*QL);

        double SolDefault=(n1*pow2(v1)+n2*pow2(v2))/pow2(n1*sqrt(v1)+n2*sqrt(v2));
        //double SolDefault=(n1*v1+n2*v2)/pow2(n1+n2);

        switch(whichSol)
        {
        case 0:
            variance=Sol0;
            break;

        case 1:
            variance=Sol1;
            break;

        case 2:
            variance=Sol2;
            break;

        case 3:
            variance=SolDefault;
            break;

        case 5:
            variance=(n1*v1+n2*v2)/pow2(n1+n2);
            break;

        case 6: // paper version
            {
                double s1=sqrt(v1);
                double s2=sqrt(v2);
                variance=s1*s2/(s1+s2)*(s1/n1+s2/n2);
            }
            break;

        default: // very degenerate solution
            variance=LARGE_POZ_VALUE;
        }
#ifdef DEBUG_PRINT
        cout << "Sol0=" << Sol0 << "\tSol1=" << Sol1 << "\tSol2=" << Sol2;
        cout << "\tSolDefault=" << SolDefault << "\twhichSol=" << whichSol;
#endif

        // assert (isfinite(variance)); //NaN?
    }

#ifdef DEBUG_PRINT
    cout << "\tvariance=" << variance << endl;
#endif

    return (finite(variance)==true) ? variance : LARGE_POZ_VALUE;
}


/**
   The hyperplane is orthogonal on one of the axis not oblique.
   the best axis to split on is determined by comparing the value of 
   the t-test for the directions, i.e.
   
       (eta1-eta2)^2
   t=-----------------
     sigma1^2+sigma2^2
 
     the split point is determined by solvind QDA on the unidimentional 
     space like for the general case
*/
double ComputeSeparatingHyperplane_Anova(double mass /* the total mass of the two distributions*/,
        double alpha_1, Vector<double>& mu1, Fortran_Matrix<double>& S1,
        double alpha_2, Vector<double>& mu2, Fortran_Matrix<double>& S2,
        Vector<double>& SeparatingHyperplane)
{
    double maxt=0.0;
    int maxi=-1;
    int N=mu1.dim();

    for (int i=0; i<N; i++)
    {
        double t=pow2(mu1[i]-mu2[i])/(S1(i+1,i+1)+S2(i+1,i+1));

        //cout << "i=" << i << "\tt=" << t << "\tmaxt=" << maxt << endl;
        //cout << "mu1=" << mu1[i] << " mu2=" << mu2[i]
        // << " s1=" << S1(i+1,i+1) << " s2=" << S2(i+1,i+1) << endl;

        if (t>maxt)
        {
            maxi=i;
            maxt=t;
        }
    }

    if (maxi==-1)
        return 0.0; // none of the directions is any good

    // solve QDA on the projection maxi
    int whichSol; //
    double eta=UnidimensionalQDA(alpha_1, mu1[maxi], S1(maxi+1, maxi+1),
                                 alpha_1, mu2[maxi], S2(maxi+1, maxi+1), whichSol);
    double variance=UnidimensionalQDAVariance(mass*alpha_1, mu1[maxi], S1(maxi+1, maxi+1),
                    mass*alpha_1, mu2[maxi], S2(maxi+1, maxi+1), whichSol);

    assert(variance>0.0);

    // form the separating hyperplane
    // by divideing all components by sqrt(variance) we get the
    // quantities computed with the Separating Hyperplane to be
    // centered on 0 and have the standard distribution of the
    // split point 1.0
    SeparatingHyperplane.newsize(N+1);
    SeparatingHyperplane=0.0;
    double sgn= (mu1[maxi]>eta) ? 1.0 : -1.0;
    SeparatingHyperplane[0]=-sgn*eta/sqrt(variance);
    SeparatingHyperplane[maxi+1]=sgn/sqrt(variance);

    double p11=alpha_1*PValueNormalDistribution(mu1[maxi], sqrt(S1(maxi+1, maxi+1)), eta);
    double p1_=p11+alpha_2*PValueNormalDistribution(mu2[maxi], sqrt(S2(maxi+1, maxi+1)), eta);
    double gini=BinaryGiniGain(p11,alpha_1,p1_);
    return gini;
}

/**
    Compute the quadratic that separates two distributions
    and take the separating hyperplane to be the tangent to it
    in the intersection point with the line between the centers.

    @param mass              the total mass of the two distributions
    @param alpha_1            alpha, mu and S (sigma) describe a multinormal distribution
    @param mu1
    @param S1
    @param alpha_2            alpha, mu and S (sigma) describe a multinormal distribution
    @param mu2
    @param S2
    @param SeparatingHyperplane
    
    @return                  SeparatingHiperplane and gini of the split
*/
double ComputeSeparatingHyperplane_QDA(double mass,
                                       double alpha_1, Vector<double>& mu1, Fortran_Matrix<double>& S1,
                                       double alpha_2, Vector<double>& mu2, Fortran_Matrix<double>& S2,
                                       Vector<double>& SeparatingHyperplane)
{

    // Sice we got here we clearly have 2 clusters. So we should better produce a split at least
    // into the regressor space.

    int N=mu1.dim();

    /* Check if one of the matrices S1 or S2 has a 0 on the diagonal. If this is the case
       then the variable corresponding to that position is either extremelly predictive (one
       value of the variable determines exactly one of the clusters) or the variable is useless
       (all the values are identical). In the first case we split on it in the seccond we
       substibute the 0 with same small value in both matrices. */

    for (int i=1; i<=N; i++)
    {
        if (S1(i,i)==0.0 && S2(i,i)==0.0)
        {
            S1(i,i)=S2(i,i)=SMALL_POZ_VALUE;
            continue;
        }
        if (S1(i,i)==0.0 || S2(i,i)==0.0)
        {
            /* split on variable i exclusivelly. Do the split between the projections of the
               centers of the clusters on this direction proportional with the sizes of the clusters.
            */

            SeparatingHyperplane.newsize(N+1);
            SeparatingHyperplane[0]=-(alpha_1*mu1[i-1]+alpha_2*mu2[i-1]);
            for (int j=1; j<=N; j++)
                SeparatingHyperplane[j]= (i==j) ? 1.0 : 0.0;

            return 1-pow2(alpha_1)-pow2(alpha_2);
        }
    }

    // compute the Cholesky of S1 in lower part of S1
    if ( Cholesky_upper_factorization(S1,S1) || Cholesky_upper_factorization(S2,S2) )
    {
        // one of the matrices S1 and S2 is badly conditioned
        cerr << "Badly conditioned matrices in ComputeSeparatingHyperplane" << endl;

    }

    // compute the square root of determinant of S1 and S2
    double detS1=1.0;
    double detS2=1.0;
    for (int i=1; i<=N; i++)
    {
        detS1*=S1(i,i);
        detS2*=S2(i,i);
    }

    /* compute the quadratic that separates the distributions
    and take the separating hyperplane to be the tangent to it
    in the intersection point with the line between the centers */

    // Compute C
    double C=2.0*log(detS2/detS1)+log(alpha_1/alpha_2);

    Vector<double> dif=mu2-mu1;
    Vector<double> G1inv_dif=
        Lower_triangular_solve< Fortran_Matrix<double>,Vector<double> > (S1,dif);
    double I1=dot_prod(G1inv_dif, G1inv_dif);

    Vector<double> G2inv_dif=
        Lower_triangular_solve< Fortran_Matrix<double>,Vector<double> > (S2,dif);
    double I2=dot_prod(G2inv_dif, G2inv_dif);

    // compute t as solution of (1-t)^2 I1-t^2 I2=C
    double t, delta;
    if (fabs(I1-I2)<SMALL_POZ_VALUE)
    {
        // Liner equation
        t=(I1-C)/(2*I1);
    }
    else
    {
        // quadratic equation
        delta=C*(I1-I2)+I1*I2;
        if (delta<0.0)
        {
            t=-1.0;
        }
        else
        {
            t=(I1-sqrt(delta))/(I1-I2);
            if (t<0.0 || t>1.0)
                t=(I1+sqrt(delta))/(I1-I2);
        }
    }

    Vector<double> n; // the normal to the Separating Hyperplane

    /* check if the solution makes any sense. If not take a point proportional
       with the sizes (alpha) on the line that goes through the two centers
    */
    if (t<0.0 || t>1.0)
    {
        // fix t
        t=alpha_1/(alpha_1+alpha_2); // in case alphas are not normalized
        n=dif;
    }
    else
    {
        // compute n as the tangent to the the quadratic in point mu
        n=
            Transposed_Upper_triangular_solve(S1,G1inv_dif)*(-1.0+t)-
            Transposed_Upper_triangular_solve(S2,G2inv_dif)*t;
    }
    n=n*(1/sqrt(dot_prod(n,n)));  //normalize n

    // compute the point which the  Separating Hyperplane crosses
    Vector<double> mu=t*mu1+(1-t)*mu2;

    // compute etai, vari
    double eta1=dot_prod(n,mu1);
    double eta2=dot_prod(n,mu2);

    // compute vari=n^T*S1*n directly into vari
    double var1=0.0, var2=0.0;
    for (int j=1; j<=N; j++)
    {
        double aux1=0.0, aux2=0.0;
        // compute n^T*Si(:,j) into auxi
        for (int i=1; i<=j; i++)
        {
            aux1+=n(i)*S1(i,j);
            aux2+=n(i)*S2(i,j);
        }
        // use only upper part of S1 and S2
        for (int i=j+1; i<=N; i++)
        {
            aux1+=n(i)*S1(j,i);
            aux2+=n(i)*S2(j,2);
        }
        var1+=aux1*n(j);
        var2+=aux2*n(j);
    }

    double variance=UnidimensionalQDAVariance(mass*alpha_1, eta1, var1,
                    mass*alpha_2, eta2, var2, 5);

    assert(variance>0.0);

    // compute the equation of the SeparatingHyperplane
    double eta=dot_prod(n,mu);
    double sgn= (eta1>eta) ? 1.0 : -1.0;
    SeparatingHyperplane.newsize(N+1);
    SeparatingHyperplane[0]=-sgn*eta/sqrt(variance);
    for(int i=1; i<=N; i++)
        SeparatingHyperplane[i]=sgn*n[i-1]/sqrt(variance);

    double p11=alpha_1*PValueNormalDistribution(mu1, S1, n, mu);
    double p1_=p11+alpha_2*PValueNormalDistribution(mu2, S2, n, mu);
    double gini=BinaryGiniGain(p11,alpha_1,p1_);
    return gini;
}


/**
    The normal of the hyperplane is the best separating direction, i.e.
    the vector on which the projection of the two gaussians is as separated 
    as possible as measured by Fisher's discriminant.
 
    The equations are:
        - normal: n=(alpha_1*S1+alpha_2*S2)^{-1}(mu1-mu2)
        - separating equations: n^T*x-eta=0
        - eta is solution of the equation

    eta^2(1/var1 - 1/var2)-2eta(eta1/var1-eta2/var2)+eta1^2/var1-eta2^2/var2
           = 2ln(alpha1/alpha2)-ln(var1/var2)
 
    with etai=n^T*mu1 and vari=n^T*Si*n
 
    Parameters: alpha,mu and S (sigma) describe a multinormal distribution.
    There are 2 of them.
    Output: SeparatingHiperplane and gini of the split
 
    S1 and S2 are modified and have their respective Cholesky factorization in them
    
*/
double ComputeSeparatingHyperplane_LDA(double mass /* the total mass of the two distributions*/,
                                       double alpha_1, Vector<double>& mu1, Fortran_Matrix<double>& S1,
                                       double alpha_2, Vector<double>& mu2, Fortran_Matrix<double>& S2,
                                       Vector<double>& SeparatingHyperplane)
{

    if (alpha_1<TNNearlyZero || alpha_2<TNNearlyZero)
        return 0.0; // there's nothing we can do. One of the clusters is nonexistent

    int N = S1.num_rows();
    Vector<double> n(N); // normal to the hyperplane
    double eta1, eta2;
    double var1, var2;

    Fortran_Matrix<double> Sw(N,N);

    // Compute upper part of Sw
    for (int j=1; j<=N; j++)
        for (int i=1; i<=j; i++)
            Sw(i,j)=alpha_1*S1(i,j)+alpha_2*S2(i,j);

    // Put mu1-mu2 in n
    for (int i=1; i<=N; i++)
        n(i)=mu1(i)-mu2(i);

    // Solve Sw^{-1}*n into n
    // First compute Cholesky of Sw into lower Sw
    if ( Cholesky_upper_factorization(Sw,Sw) )
    {
        cerr << "Badly conditioned matrices in ComputeSeparatingHyperplane" << endl;
        cerr << "We use normal anova analysis to find the Separating Hyperplane" << endl;
        return ComputeSeparatingHyperplane_Anova(mass,alpha_1, mu1, S1, alpha_2, mu2, S2, SeparatingHyperplane);
    }

    Lower_triangular_solve(Sw,n);
    Transposed_Upper_triangular_solve(Sw,n);
    // normalize n
    if (IsZero(dot_prod(n,n)))
        return 0.0; // nothing we can do. Clusters are identical

    n=n*(1.0/sqrt(dot_prod(n,n)));

    // compute etai, vari
    eta1=dot_prod(n,mu1);
    eta2=dot_prod(n,mu2);

    // compute vari=n^T*S1*n directly into vari
    var1=0.0;
    var2=0.0;
    for (int j=1; j<=N; j++)
    {
        double aux1=0.0, aux2=0.0;
        // compute n^T*Si(:,j) into auxi
        for (int i=1; i<=j; i++)
        {
            aux1+=n(i)*S1(i,j);
            aux2+=n(i)*S2(i,j);
        }
        // use only upper part of S1 and S2
        for (int i=j+1; i<=N; i++)
        {
            aux1+=n(i)*S1(j,i);
            aux2+=n(i)*S2(j,2);
        }
        var1+=aux1*n(j);
        var2+=aux2*n(j);
    }


    int whichSol;
    double eta=UnidimensionalQDA(alpha_1, eta1, var1, alpha_2, eta2, var2, whichSol);
    double variance=UnidimensionalQDAVariance(mass*alpha_1, eta1, var1,
                    mass*alpha_2, eta2, var2, whichSol);

    assert(variance>0.0);

    // form the separating hyperplane
    // by divideing all components by sqrt(variance) we get the
    // quantities computed with the Separating Hyperplane to be
    // centered on 0 and have the standard distribution of the
    // split point 1.0
    SeparatingHyperplane.newsize(N+1);
    double sgn= (eta1>eta) ? 1.0 : -1.0;
    SeparatingHyperplane[0]=-sgn*eta/sqrt(variance);
    for (int i=1; i<=N; i++)
        SeparatingHyperplane[i]=sgn*n[i-1]/sqrt(variance);

    Vector<double> mu(N);
    mu=eta*n;


    // compute the Cholesky of S1 in lower part of S1
    if ( Cholesky_upper_factorization(S1,S1) || Cholesky_upper_factorization(S2,S2) )
    {
        // one of the matrices S1 and S2 is badly conditioned
        cerr << "Badly conditioned matrices in ComputeSeparatingHyperplane" << endl;
        return 0.0; // worst gini
    }

    double p11=alpha_1*PValueNormalDistribution(mu1, S1, n, mu);
    double p1_=p11+alpha_2*PValueNormalDistribution(mu2, S2, n, mu);
    double gini=BinaryGiniGain(p11,alpha_1,p1_);
    return gini;
}


}

#endif //_CLUS_SPLITPOINTCOMPUTATION_H_

---