#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "lbl_init.h"

/*---------------------------------------------------------*/

int init_with_residual(double fix_residual)
{
  static int index=0,count=0;
  static double x[NB_RESIDUAL+1];
  double med;
  int flag_init;

  /*--- count residuals */
  count++;
  if (count > NB_RESIDUAL)
    count = NB_RESIDUAL;

  /*--- change index in vector x with modulus NB_RESIDUAL */
  index++;
  if (index > NB_RESIDUAL)
    index = 1;
  x[index] = fix_residual;

  /*--- compare current residual to median only if NB_RESIDUAL residuals 
    were already obtained */
  flag_init = 0;
  if (count == NB_RESIDUAL) {
    mediane(count,x,&med);
    if (fix_residual < med) {
      flag_init = 1;
    }
  }

  return(flag_init);
}

/*---------------------------------------------------------*/

int mediane(int n,double v[],double *med)
{
        int j;
        double arr[NB_RESIDUAL+1];

        /*--- copy v[] in arr[] */
        for (j=1; j<=n; j++)
                arr[j] = v[j];

        if (n & 1) { /* case odd */
                j = ((n+1) >> 1);
                (*med) = select(j,n,arr);
        }
        else { /* case even */
                j = n >> 1;
                (*med) = 0.5*(select(j,n,arr) + select(j+1,n,arr));
        }

        free(arr);

        return(1);
}

/*---------------------------------------------------------*/
/* select(): from numerical recipes in C                   */
/*---------------------------------------------------------*/

double select(int k,int n,double arr[])
{
        int i,ir,j,l,mid;
        double a,temp;

        l = 1;
        ir = n;
        while (1) {
                if (ir <= l+1) {
                        if (ir == l+1 && arr[ir] < arr[l]) {
                                swap(arr[l],arr[ir])
                        }
                        return arr[k];
                }
                else {
                        mid = (l+ir) >> 1;
                        swap(arr[mid],arr[l+1])
                        if (arr[l+1] > arr[ir]) {
                                swap(arr[l+1],arr[ir])
                        }
                        if (arr[l] > arr[ir]) {
                                swap(arr[l],arr[ir])
                        }
                        if (arr[l+1] > arr[l]) {
                                swap(arr[l+1],arr[l])
                        }
                        i = l+1;
                        j = ir;
                        a = arr[l];
                        while (1) {
                                do i++; while (arr[i] < a);
                                do j--; while (arr[j] > a);
                                if (j < i) break;
                                swap(arr[i],arr[j])
                        }
                        arr[l] = arr[j];
                        arr[j] = a;
                        if (j >= k) ir = j-1;
                        if (j <= k) l = i;
                }
        }
}