// Random Normal Deviates.
// Includes functions adapted from the file http://www.mas.ncl.ac.uk/~ndjw1/teaching/sim/transf/norm.c
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include<string.h>
#include "reconstruction.h"

typedef struct {
  int i;
  double x;
} PAIR;

PHENOTYPE_TABLE *PermutePhenotypes( PHENOTYPE_TABLE *pt, int in_place ) {

  PHENOTYPE_TABLE *ptable = pt;
  int *idx = (int*)calloc(pt->nobs,sizeof(int));
  double *y = (double*)calloc(pt->nobs,sizeof(double));
  int i;

  if ( !in_place ) { //create a copy, otherwise return the original struct with y permuted
    int i;
    ptable = (PHENOTYPE_TABLE*)calloc(1,sizeof(PHENOTYPE_TABLE));
    ptable->nobs = pt->nobs;
    ptable->id = (char**)calloc(pt->nobs, sizeof(char*));
    ptable->x = (double*)calloc(pt->nobs, sizeof(double));
    ptable->name = strdup(pt->name);
    for(i=0;i<ptable->nobs;i++) {
      ptable->x[i] = pt->x[i];
      ptable->id[i] = strdup(pt->id[i]);
    }
  }

  for(i=0;i<pt->nobs;i++) { 
    idx[i] = i;
    y[i] = ptable->x[i];
  }
  PermuteArray( idx, pt->nobs);

  for(i=0;i<ptable->nobs;i++) {
    ptable->x[i] = y[idx[i]];
    //    printf( "%d %e %e\n", idx[i], ptable->x[i], y[i] );
  }
  free(y);
  free(idx);
  return(ptable);
}

int pair_cmp( const void *A, const void *B) {
  const PAIR *a = (PAIR*)A;
  const PAIR *b = (PAIR*)B;
  double x = a->x-b->x;
  return ( x > 0 ? 1 : ( x<0  ? -1 : 0 ));
}

int double_cmp( const void *A, const void *B) {
  const double *a = (double*)A;
  const double *b = (double*)B;
  double x = (*a) - (*b);
  return ( x > 0 ? 1 : ( x<0  ? -1 : 0 ));
}

void PermuteArray( int *array, int n ) { // permute an array of integers in place
  PAIR *pair = (PAIR*)calloc(n,sizeof(PAIR));
  int i;
  for (i=0;i<n;i++) {
    pair[i].i = array[i];
    pair[i].x = surand();
    //    printf("%d %e\n", pair[i].i, pair[i].x);
  }
  qsort( pair, n, sizeof(PAIR), pair_cmp );
  for (i=0;i<n;i++) {
    array[i] = pair[i].i;
    //    printf("%d %d %e\n", pair[i].i, array[i], pair[i].x);
  }
  free(pair);
}

double surand()
{
  return( (double) rand()/RAND_MAX );
}

double urand(double low, double high)
{
  return(low+(high-low)*surand());
}

double genexp(double lambda)
{
  double u,x;
  u=surand();
  x=(-1/lambda)*log(u);
  return(x);
}

double gennor()
{
  double theta,rsq,x;
  theta=urand(0,2*M_PI);
  rsq=genexp(0.5);
  x=sqrt(rsq)*cos(theta);
  return(x);
}