function [dr,dirvec] = hdd1SM(data,trainp,trainn,imeth,methpar) ;
% HDD1SM, High Dimensional Discrimination
%   Steve Marron's matlab function
%     Computes various High Dimensional Discrimination methods 
%
% Inputs:
%     data   - d x nd column vector of data
%     trainp - d x n1 training data for the class "positive"
%     trainn - d x n2 training data for the class "negative"
%     imeth  - index for the method:
%                  1   - simple mean difference
%                               (methpar has no effect)
%                  2   - correlation with class means
%                               (methpar has no effect)
%                  3   - Nearest Neighbor method
%                               (methpar = number of nearest neighbors,
%                                    1 is default, recommendation:  odd)
%                  4   - Support Vector Machine
%                               (methpar = C, when >= 0,
%                                methpar = penalty factor, when < 0
%                                    (will adjust by median pairwise dist.)
%                                    1000 is default)
%                  11  - Distance Weighted Discrimination, 1
%                           based on Mike Todd's SDPT3 - sep Optimization
%                               (methpar = penalty factor,
%                                    (will adjust by median pairwise dist.)
%                                    100 is default)
%                  21  - Orthogonal Subspace Projection, 1
%                           uses SV representation,
%                           for sv's > thresh * max sv
%                               (methpar = threshfact,
%                                    10^(-12) is default)
%    methpar - method parameter (discussed above for each method)
%     
% Outputs:
%         dr - discrimination result, 1 x nd row vector
%                  of +1 or -1, with assignment of each
%                  column vector in data, to a class
%     dirvec - direction vector pointing towards positive class,
%                  unit vector (i.e. length 1)
%                  (this return is empty for imeth = 2)
%

%    Copyright (c) J. S. Marron 2002


if imeth == 1 ;    %  simple mean difference

  vmeanp = mean(trainp,2) ;
  vmeanm = mean(trainn,2) ;

  vdir = vmeanp - vmeanm ;
      %  direction vector from mean +, to mean -

  vcent = mean([vmeanp, vmeanm],2) ;
      %  centerpoint

  cutoff = vdir' * vcent ;
      %  inner product of direction vector with center

  dr = vdir' * data ;
      %  row vector of inner products with direction vector
  dr = dr >= cutoff ;
      %  1 where choose plus class, 0 otherwise
  dr = 2 * dr - 1 ;
      %  1 where choose plus class, -1 otherwise

  if nargout >= 2 ;
    dirvec = vdir / norm(vdir) ;
  end ;


elseif imeth == 2 ;    %  correlation with class means

  vmeanp = mean(trainp,2) ;
  vmeanm = mean(trainn,2) ;

  nd = size(data,2) ;
  mcorr = corrcoef([vmeanp, data]) ;
      %  (nd+1) x (nd+1) matrix of all correlations
  vcorrp = mcorr(1,2:(nd+1)) ;
      %  corrs of vmeanp with others

  mcorr = corrcoef([vmeanm, data]) ;
      %  (nd+1) x (nd+1) matrix of all correlations
  vcorrm = mcorr(1,2:(nd+1)) ;
      %  corrs of vmeanm with others

  dr = vcorrp >= vcorrm ;
      %  1 where choose plus class, 0 otherwise
  dr = 2 * dr - 1 ;
      %  1 where choose plus class, -1 otherwise

  if nargout >= 2 ;
    dirvec = [] ;
  end ;


elseif imeth == 3 ;    %  Nearest Neighbor method 

  if nargin > 4 ;    %  then have input a threshfact, so use it
    numnn = methpar ;
  else ;    %  then use default threshfact
    numnn = 1 ;
  end ;

  np = size(trainp,2) ;
  nn = size(trainn,2) ;

  mtrainall = [trainp trainn] ;
  vclassall = [ones(1,np) -ones(1,nn)] ;

  nd = size(data,2) ;
  dr = [] ;
  for idata = 1:nd ;

    dist2 = sum((vec2matSM(data(:,idata),np+nn) - mtrainall).^2,1) ;
        %  squared distances from this data point to training sample
    [temp,vind] = sort(dist2) ;
    vclassnear = vclassall(vind(1:numnn)) ;
        %  vector of classes of nearest neighbors
    dr = [dr (sum(vclassnear) >= 0)] ;
        %  discrimination result: 1 for positive, 0 for negative

  end ;
  dr = 2 * dr - 1 ;
      %  map from 0,1 to -1,1

  dirvec = [] ;
      %  no direction vector for NN discrimination


elseif imeth == 4 ;    %  Support Vector Machine

  if nargin > 4 ;    %  then have input a threshfact, so use it
    threshfact = methpar ;
  else ;    %  then use default threshfact
    threshfact = 1000 ;
  end ;

  d = size(trainp,1) ;
  np = size(trainp,2) ;
  nn = size(trainn,2) ;


  if threshfact < 0 ;    %  then have signalled should adjust using
                         %  median pairwise dist.

    %  Compute median of pairwise distances squared between classes
    %
    vpwdist2 = [] ;
    for ip = 1:np ;
      pwdist2 = sum((vec2matSM(trainp(:,ip),nn) - trainn).^2,1) ;
      vpwdist2 = [vpwdist2 pwdist2] ;
    end ;
    medianpwdist2 = median(vpwdist2) ;

    C = -threshfact / medianpwdist2 ;
        %  threshfact "makes this large", 
        %  and 1 / medianpwdist2 "puts on correct scale"
        %      [recall minus sign unencodes parameter]

  else ;

    C = threshfact ;

  end ;


  mtrainall = [trainp trainn]' ;
  vclassall = [ones(1,np) -ones(1,nn)]' ;
      %  since SVM wants data vectors as rows

  %  Run SVM algorithm
  %
  %  Using lines from program:  
  %      \Research\ComplexPopn\Old\PolyEmbed\PEegToyDisc1.m
  %
  addpath \Research\MatlabSoftware\SVM\SVM1 -end ;

    ker = 'linear' ;
  [nsv,alpha,bias] = SVC(mtrainall,vclassall,ker,C) ;
      %  classification version of SVM

  %  Find normal vector to separating hyperplane
  %
  w = sum(vec2matSM(alpha,d) .* ...
                  mtrainall .* ...
          vec2matSM(vclassall,d),1) ;
      %  sum_i alpha_i x_i y_i

  %  Find intercept value of separating hyperplane
  %
  [temp,r] = max(alpha .* (vclassall == 1)) ;
      %  r is index of biggest alpha from class 1
  [temp,s] = max(alpha .*(vclassall == -1 )) ;
      %  s is index of biggest alpha from class -1
  b = 0.5 * w * (mtrainall(r,:) + mtrainall(s,:))' ;

  dr = (data') * w' - b ;
      %  use dot product with SVM vector as discriminator
      %  need to transpose data, since SVM treats rows as data
  dr = (dr > 0)' ;
      %  one where think class seems positive, 0 otherwise
      %  tranpose to put back in "cols as data" form
  dr = 2 * dr - 1 ;
      %  map from 0,1 to -1,1

  dirvec = w' / norm(w) ;
      %  make this a unit vector, and a column vector


elseif imeth == 11 ;    %  Distance Weighted Discrimination, 1
                        %  Using lines from ToddEmail2.txt
                        %  and the version of sep.m in ToddEmaail12.txt

  if nargin > 4 ;    %  then have input a threshfact, so use it
    threshfact = methpar ;
  else ;    %  then use default threshfact
    threshfact = 100 ;
  end ;
  

  addpath '\Research\MatlabSoftware\HDD\SDPT3\Solver' -end ;   
  addpath '\Research\MatlabSoftware\HDD\SDPT3\Solver\mexexec' -end ;
  parameters ;
      %% specify default parameters for sqlp.m,
      %% they are specified in the structure called OPTIONS.
  global CACHE_SIZE   % cache size in kbytes
  global LOOP_LEVEL   % loop unrolling level
  CACHE_SIZE = 256;
  LOOP_LEVEL = 8;
      %% set global variables for functions imported from LIPSOL


  %  Compute median of pairwise distances squared between classes
  %
  np = size(trainp,2) ;
  nn = size(trainn,2) ;
  vpwdist2 = [] ;
  for ip = 1:np ;
    pwdist2 = sum((vec2matSM(trainp(:,ip),nn) - trainn).^2,1) ;
    vpwdist2 = [vpwdist2 pwdist2] ;
  end ;
  medianpwdist2 = median(vpwdist2) ;

  penalty = threshfact / medianpwdist2 ;
      %  threshfact "makes this large", 
      %  and 1 / medianpwdist2 "puts on correct scale"
  [w,beta,residp,residn,alp,totalviolation,dualgap,flag] = sep(trainp,trainn,penalty) ;
%  [w,beta,residp,residn,totalviolation,X,y,Z,err0,flag] = sepn(trainp,trainn,penalty) ;
%  [w,beta,residp,residn,totalviolation,X,y,Z,err0,flag] = sepOld(trainp,trainn,penalty) ;

  if flag == -1 ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Warning from HDD1SM:                           !!!') ;
    disp('!!!   sep optimization gave an inaccurate solution   !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  elseif flag == -2 ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Error from HDD1SM:                             !!!') ;
    disp('!!!   Infeasible or unbounded optimization problem   !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    dr = [] ;
    dirvec = [] ;
    return ;
  end ;

  dr = w'* data + beta >= 0 ;
      %  1 where choose plus class, 0 otherwise
  dr = 2 * dr - 1 ;
      %  1 where choose plus class, -1 otherwise

  if nargout >= 2 ;
    dirvec = w / norm(w) ;
  end ;


elseif imeth == 21 ;    %  Orthogonal Subspace Projection, 1

  if nargin > 4 ;    %  then have input a threshfact, so use it
    threshfact = methpar ;
  else ;    %  then use default threshfact
    threshfact = 10^(-12) ;
  end ;

  d = size(data,1) ;
  nd = size(data,2) ;

  np = size(trainp,2) ;
  meanp = mean(trainp,2) ;
      %  mean column vector
  trainpmm = trainp - vec2matSM(meanp,np) ;
  datammp = data - vec2matSM(meanp,nd) ;
      %  mean recentered versions
  [up,sp,vp] = svd(trainpmm,0) ;
      %  Singular value decomposition, to get projection matrix
      %  0 for "economy size" SVD, useful for HDLSS
  vsvp = diag(sp) ;
      %  vector of singular values
  flagp = (vsvp >= threshfact * max(vsvp)) ;
  if sum(flagp) >= d ;
    %  then training data are full rank, 
    %  so can't do OSP
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Error from HDD1SM:                       !!!') ;
    disp('!!!   OSP failed because trainp is full rank   !!!') ;
    disp('!!!   giving empty returns                     !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    dr = [] ;
    dirvec = [] ;
    return ;
  else ;    %  then can proceed with OSP calculation
    up = up(:,flagp) ;
        %  just keep cols with large enough sv's
        %  projection matrix is up * up', 
        %  but don't compute, since could be too big (dxd)
    projdatap = up' * datammp ;
    projdatap = up * projdatap ;
        %  projection of data onto plane generated by trainp
    projdataop = datammp - projdatap ;
        %  projection of data onto plane orthogonal to trainp
    vnop = sum(projdataop.^2,1) ;
        %  vector of squared norms of data 
        %  projected orthogonal to trainp
  end ;

  nn = size(trainn,2) ;
  meann = mean(trainn,2) ;
      %  mean column vector
  trainnmm = trainn - vec2matSM(meann,nn) ;
  datammn = data - vec2matSM(meann,nd) ;
      %  mean recentered versions
  [un,sn,vn] = svd(trainnmm,0) ;
      %  Singular value decomposition, to get projection matrix
      %  0 for "economy size" SVD, useful for HDLSS
  vsvn = diag(sn) ;
      %  vector of singular values
  flagn = (vsvn >= threshfact * max(vsvn)) ;
  if sum(flagn) >= d ;
    %  then training data are full rank, 
    %  so can't do OSP
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Error from HDD1SM:                       !!!') ;
    disp('!!!   OSP failed because trainn is full rank   !!!') ;
    disp('!!!   giving empty returns                     !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    dr = [] ;
    dirvec = [] ;
    return ;
  else ;    %  then can proceed with OSP calculation
    un = un(:,flagn) ;
        %  just keep cols with large enough sv's
        %  projection matrix is un * un', 
        %  but don't compute, since could be too big (dxd)
    projdatan = un' * datammn ;
    projdatan = un * projdatan ;
        %  projection of data onto plane generated by trainn
    projdataon = datammn - projdatan ;
        %  projection of data onto plane orthogonal to trainn
    vnon = sum(projdataon.^2,1) ;
        %  vector of squared norms of data 
        %  projected orthogonal to trainn
  end ;


  dr = (vnop <= vnon) ;
      %  1 where  data are "farther from negative class"
      %  i.e. choose plus class, 0 otherwise
  dr = 2 * dr - 1 ;
      %  1 where choose plus class, -1 otherwise


  dirvec = []  ;


else ;

  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Error from hdd1SM.m:         !!!') ;
  disp('!!!   invalid imeth, empty return  !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

  dr = [] ;
  dirvec = [] ;


end ;



