#include <math.h>
#define NRANSI
#include "nrutil.h"
#include <stdio.h>

extern long idum;
extern double tt;

double damotsa(double **p, double y[], double psum[], int ndim, double pb[],
	double *yb, double (*funk)(double []), int ihi, double *yhi, double fac)
{
	float ran1(long *idum);
	int j;
	double fac1,fac2,yflu,ytry,*ptry;
	double total=0.0;

	ptry=dvector(1,ndim);
	fac1=(1.0-fac)/ndim;
	fac2=fac1-fac;
	for (j=1;j<=ndim;j++)
		ptry[j]=psum[j]*fac1-p[ihi][j]*fac2;
	for(j=1;j<=ndim;j++) {
		if(ptry[j]<0.0) { /*printf("exist at least one element less than 0\n"); fflush(stdout);*/ return 1.0e10; } 
		else total+=ptry[j];
	}
	if(total>1.0) { /*printf("sum of all elements greater than 1.0\n"); fflush(stdout);*/ return 1.0e10; }
	ytry=(*funk)(ptry);
	if (ytry <= *yb) {
		for (j=1;j<=ndim;j++) pb[j]=ptry[j];
		*yb=ytry;
//		for(int i=1;i<=20;i++) printf("%8.5f",pb[i]);
//		printf("   LOG likelihood: %12.5f\n",*yb);
//		printf("\n");
	}
	yflu=ytry-tt*log(ran1(&idum));
	if (yflu < *yhi) {
		y[ihi]=ytry;
		*yhi=yflu;
		for (j=1;j<=ndim;j++) {
			psum[j] += ptry[j]-p[ihi][j];
			p[ihi][j]=ptry[j];
		}
	}
	free_dvector(ptry,1,ndim);
	return yflu;
}
#undef NRANSI
/* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */
