/*
 *  IFD-R.c: solving the ideal-free distribution for a generalised functional
 *           response with exponential interference function
 *  
 *
 *  Created by Adriaan Dokter
 *  last modified: Apr-2015
 *
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <float.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

#define MAX(a,b) (((a)>(b))?(a):(b))

struct f_params
  {
	  int size;
	  float *areas;
	  float *ks;
	  float *FRs;
	  int Nbirds;
	  float p0;
  };


/********************************************************************/
/* Prototypes of local functions:                                   */
/********************************************************************/

double f_ifd (double x, void *params);
void IFDgsl(float *FRs,float *areas,int *Nbirds,float *kk,float *frandom,int *Ndata, float *distrib,float *maxw, double *logival, int *maxiter, int *type);

/********************************************************************/
/* local functions                                                  */
/********************************************************************/

double f_ifd (double x, void *params){
	struct f_params *p = (struct f_params *) params;
	double output = 0;
	int i;
	int size = p->size;
	float *areas = p->areas;
	float *ks = p->ks;
	float *FRs = p->FRs;
	int Nbirds = p->Nbirds;
	float p0 = p->p0;
	output=-Nbirds;
	for (i=0; i<size; i++) {
		output+=MAX(0,(areas[i]/ks[i])*(log(FRs[i])-x)-areas[i]*p0);
	}
	return output;
}

/* calculates an ideal free distribution for a generalised functional response with an exponential interference function.
          FRs: array of interference-free gain rates
        areas: array of patch areas
       Nbirds: number of birds to be distributed
           kk: array of interference-constants
      frandom: fraction of birds distributed randomly over patches 
        Ndata: array-length of the arrays FRs, areas, kk, distrib
      distrib: array to hold the calculated bird-densities
         maxw: logarithm of the maximised gain rate 
      logival: array of length 2. Exp(logival) gives the interval in which
               the maximised gain rates should be found
      maxiter: the maximum number of iterations of Brent`s root finding
               algorithm
	     type: determines whether the randomly distributed birds (given by the fraction frandom)
			   do (type=1) or do not (type=0) interfere with the ideal-free distributed birds.
*/
void IFDgsl(float *FRs,float *areas,int *Nbirds,float *kk,float *frandom,int *Ndata, float *distrib,float *maxw, double *logival, int *maxiter, int *type){
	int status;
	int iter = 0, max_iter = maxiter[0];

	int imax = Ndata[0];
	float fran=frandom[0];
	int Nbird=(1-frandom[0])*Nbirds[0];
	float k = kk[0];
    int i = 0;
	float areasum=0;
	for (i=0; i<imax; i++) areasum=areasum+areas[i];
    
    /* set all negative gain rates to zero */
    float FRs0[imax];
    int nFRs0=0;
    for (i=0; i<imax; i++){
        if (FRs[i]>0){
            FRs0[i]=FRs[i];
        }
        else {
            FRs0[i]=0;
            nFRs0++;
        }
    }
    
    /* if all gain rates zero or negative, distribute all birds randomly */
    if (nFRs0==imax) {
        for (i=0; i<imax; i++) distrib[i]=Nbirds[0]/areasum;
        return;
    }
    
	const gsl_root_fsolver_type *T;
	gsl_root_fsolver *s;
	double r = 0;
	double x_lo = logival[0], x_hi = logival[1];
	gsl_function F;
	struct f_params params;
	params.size=imax;
	params.areas=areas;
	params.ks=kk;
	params.FRs=FRs0;
	params.Nbirds=Nbird;
	if (type[0]==0) {
		params.p0=0;
	}
	else {
		params.p0=fran*Nbirds[0]/areasum;
	}
	F.function = &f_ifd;
	F.params = &params;
	
	T = gsl_root_fsolver_brent;
	s = gsl_root_fsolver_alloc (T);
	gsl_root_fsolver_set (s, &F, x_lo, x_hi);

	/*maximise the gain rate r */
	do
    {
		iter++;
		status = gsl_root_fsolver_iterate (s);
		r = gsl_root_fsolver_root (s);
		x_lo = gsl_root_fsolver_x_lower (s);
		x_hi = gsl_root_fsolver_x_upper (s);
		status = gsl_root_test_interval (x_lo, x_hi,
										 0, 0.001);
    }
	while (status == GSL_CONTINUE && iter < max_iter);
	
	gsl_root_fsolver_free (s);
	
	//calculate the ideal free density, given the maximised gain rate r
	if(r>logival[0]){
		for (i=0; i<imax; i++) {
			distrib[i]=MAX(0,(log(FRs0[i])-r)/kk[i]-params.p0);
		} 
	}
	else{
		//no food available, so distribute all birds evenly
		fran=1.0;
	}
	//add the non-foraging randomly distributed birds
	for (i=0; i<imax; i++) distrib[i]=distrib[i]+fran*Nbirds[0]/areasum;
	maxw[0]=r;
}