/* Written by Marianne Rooman, Department of BioSystems, BioModeling & 
   BioProcesses, Université Libre de Bruxelles, Belgium 
*/

#include <R.h>
#include <Rmath.h>

#include <math.h>
#include <limits.h>
/*----------------------------*/
void LR_distance(dat1,dat2,nel1,nel2,dd)  
/*----------------------------*/
     double dat1[],dat2[],*dd;
     int *nel1,*nel2;
{
  int v;
  double m1,m2,cc,s1,s2,dd1,dd2;

  s1=0.0;
  s2=0.0;
  m1=0.0;
  m2=0.0;
  cc=0;
  *dd=(-1.0);
  if(*nel1!=(*nel2))
    {Rprintf("Error, the input vectors do not have the same length\n");
      return;
    } 
  if(*nel1<=0)
    {Rprintf("Error, length of the vectors = %d\n",*nel1);
      return;
    } 
  if(*nel1==1)
    {*dd=0.0;
      return;
    } 
  for(v=0;v<*nel1;v++)
    if(dat1[v]>=(double)INT_MAX)
      {Rprintf("Error, element %d of vector 1 is too large (= %f)\n",v+1,dat1[v]);
	return;
      } 
    else if(dat2[v]>=(double)INT_MAX)
      {Rprintf("Error, element %d of vector 2 is too large (= %f)\n",v+1,dat2[v]);
	return;
      } 

  for(v=0;v<*nel1;v++)
    {m1+=dat1[v];
      m2+=dat2[v];
      s1+=dat1[v]*dat1[v];
      s2+=dat2[v]*dat2[v];
      cc++;
    }
  m1/=cc;
  m2/=cc;
  s1=sqrt(s1/cc-m1*m1);
  s2=sqrt(s2/cc-m2*m2);

  if(fabs(s1)<=1.0/(double)INT_MAX||fabs(s2)<=1.0/(double)INT_MAX||s1*s2<0)
    {Rprintf("Error, one of the vectors contains almost constant values\n");
      return;
    }

  dd1=0.0;
  dd2=0.0;
  cc=0.0;
  for(v=0;v<*nel1;v++)
    {dd1+=pow((dat1[v]-m1)/s1-(dat2[v]-m2)/s2,2);
      dd2+=pow((dat1[v]-m1)/s1+(dat2[v]-m2)/s2,2);
      cc++;
    }
  dd1=sqrt(s1*s2*dd1/cc);
  dd2=sqrt(s1*s2*dd2/cc);
  if(dd1<dd2)*dd=dd1;
  else *dd=dd2;
  if(fabs(*dd)>=(double)INT_MAX)
    {Rprintf("Error, the computed distance is too large\n");
      *dd=(-1.0);
    }
}
/*----------------------------*/