/*      program epscore5.c    Craig Tovey  August 2009-June 2010 */
/*       this program finds the minimum value of epsilon		 */
/*      for which the epsilon core exists, and the corresponding */
/*      point, for 3 or 5 or 7 voters in the plane.                   */
/*       This is a modification of version epscore4.c            */
/*			which is a small modification of					*/
/*      epscore1.c that outputs some data to a file, checks      */
/*      if the epsilon function is convex, and allows you to     */
/*      add an offset to the starting point and decrease stepsize*/
/*      to help search for multiple local optima.  epscore4 only */
/*		handled 3 or 5 voters					                 */
/*      Last modified June 6 2010                                */

#include <stdio.h>
#include <math.h>

FILE *f1;

#define N  7    /* max # of ideal points  */
#define M  40   /* max size grid to check for convexity */
#define Stepsilon 0.000000005 /* accuracy level desired = least step size */
#define Stepmin 0.1 /* min value of initial step size */
#define convepsilon 0.001 /* slack permitted in convexity check */
#define pi 3.14159


int acquire_data();
void describe_epscore(int n);
void convexity_check(int n); 
void compute_epscore();
void report_core();
long double epsvalue(long double x,long double y);
int improve(long double stepsize);
int incore(long double x, long double y,long double eps);
int threecircles(int i1, int i2,int i3,	long double r1,long double r2,long double r3);
void report_finagle_radius();


long double xov[N+1];  /* x coordinates of voters */
long double yov[N+1];  /* y coordinates of voters */
long double dovv[N+1][N+1]; /* distance from one voter to another voter */
long double eoc;      /* epsilon value of epsilon core */
long double xoc,yoc;  /* x,y coordinates of core */
int    nov;       /* number of voters */

main()
{ 
  int err;   long double v;

  f1 = fopen("epcore.7.txt","w");
  err = acquire_data();
  err = (err || (f1 == NULL));
  if (err==0)
  {
/*	  v = epsvalue(0.3,0.3);
	  printf("epsilon value at .3,.3 is %lf\n",v);*/
	  v = epsvalue(0.83,-3.51);
	  printf("epsilon value at 0,0 is %lf\n",v);
	compute_epscore();     /* call the functions you want to call */
	report_core();
	report_finagle_radius();
/*	describe_epscore(12);*/
	convexity_check(40);
  }
  else
  { 
    printf("Error %d opening file or acquiring data\n",err);
  };
  printf("bye \n");
}

/*  This function sets the value of nov to 3 or 5 and sets the voter
    ideal point coordinates xov and yov.  From these it computes
	the inter-voter distances dovv.                                 */
int acquire_data()     
{
  int error;
  int i,j;  
  error=0;
 nov=3;  /* set this as you want */
 if (nov==3)
 {
  i=1;  /* set these as you want */
  xov[i]=-8.16; yov[i]=0.34;
  i++;
  xov[i]=10.31; yov[i]=0.48;
  i++;
  xov[i]=0.8; yov[i]=2.06;

  
 }                                  
 else if (nov==5)
 {
		 xov[4]=4.0; xov[3]=6.0; xov[5]=11.0;            xov[2]=1.0;
         yov[4]=8.0; yov[3]=3.0; yov[5]=14.0;  yov[2]=5.0;
		 xov[1]=9.0; yov[1]=2.5; 
             /* set these as you want; above is my nonconvex case */
/*	    xov[1]=0.0; yov[1]=0.0;
		xov[2]=1.0; yov[2]=1.0;
		xov[3]=-1.0; yov[3]=1.0;
		xov[4]=0.6; yov[4]=2.0;
		xov[5]=-0.6; yov[5]=2.0;*/

		 if (1)
		 {
			xov[1]=1.0; yov[1]=0.0; /* regular pentagon */
			xov[2]=cos(2.0*pi/5.0); yov[2]=sin(2.0*pi/5.0);
			xov[3]=cos(4.0*pi/5.0); yov[3]=sin(4.0*pi/5.0);
			xov[4]=xov[3]; yov[4] = -yov[3];
			xov[5]=xov[2]; yov[5] = -yov[2];
		 };
 }
 else 
 {				/* regular heptagon */
	xov[1]=1.0; yov[1]=0.0; 
	xov[2]=cos(2.0*pi/7.0); yov[2]=sin(2.0*pi/7.0);
	xov[3]=cos(4.0*pi/7.0); yov[3]=sin(4.0*pi/7.0);
	xov[4]=cos(6.0*pi/7.0); yov[4]=sin(6.0*pi/7.0);
	xov[5]=xov[4]; yov[5]=-yov[4];
	xov[6]=xov[3]; yov[6]=-yov[3];
	xov[7]=xov[2]; yov[7]=-yov[2];

 };


  for (i=1; i < nov; i++)
  {
	for (j=i+1; j<=nov; j++)
	{
		dovv[i][j] = sqrt( ((xov[i]-xov[j])*(xov[i]-xov[j])) + ((yov[i]-yov[j])*(yov[i]-yov[j])) );
		dovv[j][i] = dovv[i][j];
	}
  }
		
  return error;
}




void describe_epscore(int n)   /* computes eps values on an n by n grid around minimum */
{
  int i,j;     /* indices */
  long double s;  /* size of grid */
  long double x,y; /* current grid location */
  long double e;

  fprintf(f1,"x and y coordinates of ideal points are:\n");
  for (i=1; i <= nov; i++)
  {
	  fprintf(f1,"%lf , %lf\n",xov[i],yov[i]);
  }
  compute_epscore();
  s=0.0;
  for (i=1; i<nov; i++)        /* compute average internode distance */
  {
	  for (j=i+1; j<= nov; j++)
	  {
		  s += dovv[i][j];
	  };
  };
  s = s / ((nov-1.0) * nov);


  y = yoc - s/2.0;  x = xoc - s/2.0;   /* write out grid */
  fprintf(f1,"Lower left corner of grid is at %lf ,  %lf\n",x,y);
  fprintf(f1,"step size is %lf \n",s/n);
  for (i=1; i <= n; i++)
  {
	  x = xoc - s/2.0;
	  for (j=1; j <= n; j++)
	  {
		 e=epsvalue(x,y);
		 fprintf(f1,"  %lf  ",e);
		 x += s / n ;
	  };
	  fprintf(f1,"\n");
	  y += s / n;
  };
  

}



void convexity_check(int n)   /* checks for convexity of eps values on an n by n grid around minimum */
{
  int i,j;     /* indices */
  long double s;  /* size of grid */
  long double x,y; /* current grid location */
  long double e;
  long double epsarray[M+1][M+1];

  fprintf(f1,"x and y coordinates of ideal points are:\n");
  for (i=1; i <= nov; i++)
  {
	  fprintf(f1,"%lf , %lf\n",xov[i],yov[i]);
  }
  compute_epscore();
  s=0.0;
  for (i=1; i<nov; i++)        /* compute average internode distance */
  {
	  for (j=i+1; j<= nov; j++)
	  {
		  s += dovv[i][j];
	  };
  };
  s = s / ((nov-1.0) * nov);


  y = yoc - s/2.0;  x = xoc - s/2.0;   /* write out grid and store in array */
  fprintf(f1,"Lower left corner of grid is at %lf ,  %lf\n",x,y);
  fprintf(f1,"step size is %lf \n",s/n);
  for (i=1; i <= n; i++)
  {
	  x = xoc - s/2.0;
	  for (j=1; j <= n; j++)
	  {
		 e=epsvalue(x,y);
		 fprintf(f1,"  %lf  ",e);
		 /*if (e >= 100.0) {fprintf(f1,"at %lf,%lf\n",x,y); return;};*/
		 epsarray[i][j] = e;
		 x += s / n ;
	  };
	  fprintf(f1,"\n");
	  y += s / n;
  };
  

  /* now look for nonconvexities */
  for (i=2; i<n; i++)
  {
	  for (j=2;j<n;j++)
	  {
		  if (epsarray[i][j] - convepsilon > (epsarray[i-1][j] + epsarray[i+1][j])/2.0) 
			  fprintf(f1,"nonconvexity at %d,%d [+-1,0]   e=%lf\n",i,j,epsarray[i][j]);
			if (epsarray[i][j]  - convepsilon > (epsarray[i][j-1] + epsarray[i][j+1])/2.0)
			  fprintf(f1,"nonconvexity at %d,%d [0,+-1]\n",i,j);
			if (epsarray[i][j]  - convepsilon > (epsarray[i-1][j-1] + epsarray[i+1][j+1])/2.0)
			  fprintf(f1,"nonconvexity at %d,%d  +-[1,1]\n",i,j);
			if (epsarray[i][j]  - convepsilon > (epsarray[i-1][j+1] + epsarray[i+1][j-1])/2.0)
			  fprintf(f1,"nonconvexity at %d,%d  +-[1,-1]\n",i,j);
	  };
  };
}



void compute_epscore()   /* finds a local minimum of the epsilon function */
{
  long double stepsize,rn,e;
  int i,k;
  xoc=0.0; yoc=0.0;
  for (i=1; i<=nov;i++)
  {
	  xoc += xov[i]; yoc += yov[i];
  };
  rn = (long double) nov;
  xoc =  xoc/rn; yoc =  yoc/rn; /* initialize to center of voter mass */
  /* now add offset to help look for other local optima */
  xoc += 5.0; yoc += 7.0;   /* choose whatever offset or starting point you want here */
  stepsize = fabs(xoc - xov[1]) + fabs(yoc - yov[1]); /* heuristic initial step size */
  if (stepsize <  Stepmin) {stepsize=Stepmin;};
  /* optionally change stepsize to small value to help find new local optima */
   stepsize = 0.1; 
  while (stepsize > Stepsilon)
  {
    k=1;
	e = epsvalue(xoc,yoc);
	printf("Starting local search from %lf,%lf where the epsilon value is %lf\n",xoc,yoc,e);
	while (k) {k=improve(stepsize);};
	printf("moved to %lf,%lf with step size %lf \n",xoc,yoc,stepsize);
	stepsize = stepsize / 4.0;
  };
}



int improve(long double stepsize) /* start at xoc,yoc improve to best neighbor w.r.t. stepsize */
                                  /* does not assume that eoc is the epsilon value at xoc,yoc */
{
	long double xoc1,yoc1,eoc1,xoc2,yoc2,eoc2;
	int improved,jx,jy;

	long double bestxoc,bestyoc,besteoc;
	xoc1=xoc; yoc1=yoc; /* local copies of xoc,yoc */
	eoc1 = epsvalue(xoc1,yoc1);
	improved = 0;
	besteoc = eoc1; bestxoc = xoc1; bestyoc =yoc1; /* set initial point as incumbent */
	for (jx=-2; jx<=2; jx+=1)
	{
		for (jy=-2; jy<=2; jy+=1)
		{
			xoc2 = xoc1+(jx * stepsize); yoc2=yoc1+(jy * stepsize);
			eoc2 = epsvalue(xoc2,yoc2);
			if (eoc2 < besteoc) {besteoc=eoc2; bestxoc = xoc2; bestyoc = yoc2; improved=1;};
		}  /* may later need to use separate x and y stepsizes */
	}
	xoc=bestxoc; yoc=bestyoc;eoc=besteoc;
	return(improved);
}


long double epsvalue(long double x,long double y)  /* returns least epsilon s.t. x,y is in epsilon-core */
{

	long double e;
	long double delta;
	long double lower;
	long double upper;
	long double middle;
	int l;
	delta=0.01*Stepsilon; /* delta is desired level of accuracy for returned value */
	lower=0.0;
	upper= 500.0; /* should compute not just guess upper bound */
	while (upper - lower > delta)
	{
        middle = (upper+lower)/2.0;
		l = incore(x,y,middle);
		if (l) {upper= middle;}
		else {lower= middle;}
	};
	e = upper;
	return(e);
}




int incore(long double x, long double y,long double eps)
{
	/* returns true if x,y is in the epsilon core with epsilon == eps */
	/*  i.e. x,y is not eps-dislodged by any point in the plane       */
    /*  only handles 3 or 5 or 7 ideal points                              */
int l;
int i,j,k;
long double di,dj,dk;
long double dov[N+1]; /* distance from voter v to the point x,y */
int hyper[N+1][N+1][N+1];  /* hypergraph of which 3-tuples of circles intersect */

for (i=1; i<=nov; i++)
{								/* calculate dov[] */
	dov[i] = sqrt( ((xov[i]-x)*(xov[i]-x)) + ((yov[i]-y)*(yov[i]-y)) );
};
if (nov == 3)          /* handle 3 voter case separately */
{

    l = (dovv[1][2] >= dov[1] + dov[2] - eps -eps);
	l = l && (dovv[1][3] >= dov[1] + dov[3] - eps - eps);
	l = l && (dovv[2][3] >= dov[2] + dov[3] - eps - eps);
	return(l);
}
else
if (nov == 5)         /* handle 5 voter case separately */
{
	for (i=1; i<=3; i++)     /* does any triple i,j,k of voters prefer
							    some y to x?                     */
	{
		di = dov[i] - eps;
		for (j=i+1; j<=4; j++)
		{
			dj = dov[j] - eps;
			for (k=j+1; k<=5; k++)
			{
				dk = dov[k] - eps;    /* y exists if the three circles have
									     an intersection */
				if (threecircles(i,j,k,di,dj,dk) ) return(0);
			};
		};
	};
	return(2);  /* no such y exists for any triple i,j,k */
}
else 
if (nov == 7)  /* handle 7 voter case separately */
{
	/* build hypergraph of 3-edges */
	/* hyper[i][j][k] == 1 if vi vj vk have nonempty intersection */
	for (i=1; i<=nov-2; i++)
	{
		di = dov[i] - eps;
		for (j=i+1; j<=nov-1; j++)
		{
			dj = dov[j] - eps;
			for (k=j+1; k<= nov; k++)
			{
				dk = dov[k] - eps; 
				l = threecircles(i,j,k,di,dj,dk);
				hyper[i][j][k]=l;
				hyper[i][k][j]=l;
				hyper[j][i][k]=l;
				hyper[j][k][i]=l;
				hyper[k][i][j]=l;
				hyper[k][j][i]=l;
			};
		};
	};

	/* does hypergraph have a clique of size 4? If so not in epsilon-core */
	for (i=1; i<=nov-3; i++)
	{	for (j=i+1; j<=nov-2; j++)
		{	for (k=j+1; k<=nov-1; k++)
			{	for (l=k+1; l<=nov; l++)
				{
					if (hyper[i][j][k] && hyper[i][j][l] && hyper[i][k][l] && hyper[j][k][l])
					return(0);
				};
			};
		};
	};
	return(1);
}
else
/* general case is stubbed for now                       */
				/* algorithm must build the hypergraph of 3-edges        */
				/* and find if the max clique has more than nov/2 voters */
{
	l = 1;
};
return(l);
}



int threecircles(int i1, int i2,int i3,	long double r1,long double r2,long double r3)

{
	/* returns true if the interiors of the three circles centered at voters i1,i2,i3 with 
	   respective radii r1,r2,r3  have nonempty intersection                            */				 
	long double x1,x2,x3,x4;
	long double y1,y2,y3,y4;
	long double d,cosalpha,sinalpha;
	/* first first check for nonpositive radii */
	if (r1<0.0) return(0); if (r2<0.0) return(0); if (r3<0.0) return(0);
	/* first check for pairwise intersections */
	if (dovv[i1][i2] >= r1 + r2) return(0);
	if (dovv[i1][i3] >= r1 + r3) return(0);
	if (dovv[i2][i3] >= r2 + r3) return(0);
				/* next, check if any circle encloses another circle             */
	if (dovv[i1][i2] + r1 <= r2) return(1);
	if (dovv[i1][i2] + r2 <= r1) return(1);
	if (dovv[i1][i3] + r1 <= r3) return(1);
	if (dovv[i1][i3] + r3 <= r1) return(1);
	if (dovv[i3][i2] + r3 <= r2) return(1);
	if (dovv[i3][i2] + r2 <= r3) return(1);

				/* next, return yes if any center is within the other two circles */
	if ((dovv[i1][i2] < r2) && (dovv[i1][i3] < r3) ) return(1);
	if ((dovv[i2][i1] < r1) && (dovv[i2][i3] < r3) ) return(2);
	if ((dovv[i3][i1] < r1) && (dovv[i3][i2] < r2) ) return(3);

	/* compute endpoints of intersection of circles 1 and 2 */
    d = dovv[i1][i2];
	cosalpha = (r2*r2 + d*d - r1*r1) / (2.0*r2*d);  /* alpha is angle at c2, facing r1 */
	sinalpha = sqrt(1 - (cosalpha*cosalpha));
    x1 = (xov[i1] - xov[i2]) * ( r2/d);
	y1 = (yov[i1] - yov[i2]) * ( r2/d);/* x1,y1 is the radius vector of circle at voter i2 */
	                                      /* that points towards voter i1                  */
	x2 = xov[i2]; y2 = yov[i2];              /* start at voter i2 */
	x3 = x2;      y3 = y2;
	x2 += cosalpha * x1 + sinalpha*y1;     /* add vector rotated by alpha to get endpoints */
	y2 += cosalpha * y1 - sinalpha*x1;
	x3 += cosalpha * x1 - sinalpha*y1;
	y3 += cosalpha * y1 + sinalpha*x1;
 
						/* return yes if either endpoint is within 3rd circle */
	if ( (x2-xov[i3])*(x2-xov[i3]) + (y2-yov[i3])*(y2-yov[i3]) < r3*r3 ) return(4);
	if ( (x3-xov[i3])*(x3-xov[i3]) + (y3-yov[i3])*(y3-yov[i3]) < r3*r3 ) return(5);

						/* now check if a point on the boundary between the endpoints */
						/* is close enough											  */
 
	x4 = xov[i2] +  ( xov[i3]-xov[i2] )*r2/dovv[i2][i3];
	y4 = yov[i2] +  ( yov[i3]-yov[i2] )*r2/dovv[i2][i3];  /* x4,y4 by definition is   */
														 /* in circles 2 and 3 since it */
	                                                    /* is known they intersect      */							
	if ( (x4-xov[i1])*(x4-xov[i1]) + (y4-yov[i1])*(y4-yov[i1]) < r1*r1 ) return(6);
														/* now repeat for circle 1 */
	x4 = xov[i1] +  ( xov[i3]-xov[i1] )*r1/dovv[i1][i3];
	y4 = yov[i1] +  ( yov[i3]-yov[i1] )*r1/dovv[i1][i3];  /* x4,y4 by definition is   */
														 /* in circles 1 and 3 since it */
	                                                    /* is known they intersect      */							
	if ( (x4-xov[i2])*(x4-xov[i2]) + (y4-yov[i2])*(y4-yov[i2]) < r2*r2 ) return(7);




	return(0);
}




void report_core()
{
 
	 printf("Epsilon core point at coordinates %lf, %lf with epsilon = %lf\n",xoc,yoc,eoc);
     fprintf(f1,"Epsilon core local minimum at coordinates %lf, %lf with epsilon = %lf\n",xoc,yoc,eoc);
}




void report_finagle_radius()  /* computes finagle radius, yolk radius, and strong point */
							  /* of the triangle defined by xov[i],yov[i] i=1,2,3.      */
							  /* When nov==3 provides examples of strong point and yolk */
							  /* not being same as epsilon-core                         */
{
	long double s,a,b,c,alph,bet,gam,r,x,y,area;
    a = dovv[2][3];
	b = dovv[1][3];
	c = dovv[1][2];
	s = (a+b+c)/2.0; /* semiperimeter */
	alph = 1.0/(s-a);
	bet  = 1.0/(s-b);
	gam  = 1.0/(s-c);
	r    = 1.0/(alph+bet+gam+2*sqrt(alph*bet + alph*gam + bet*gam));
	printf("finagle radius of first three points = %lf\n",r);
	fprintf(f1,"finagle radius of first three points = %lf\n",r);
	/* now compute coordinates of center of yolk of first three points */
	x = (a * xov[1] + b * xov[2] + c * xov[3])/(a+b+c);
	y = (a * yov[1] + b * yov[2] + c * yov[3])/(a+b+c);
	fprintf(f1,"yolk center of first three points has coordinates %lf,%lf\n",x,y);
	r = sqrt((s-a)*(s-b)*(s-c)/s);
	fprintf(f1,"yolk radius of first three points = %lf\n",r);
	area = r * s;
	alph = 2.0 * area / (b*c);
	alph = asin(alph);
	bet  = 2.0 * area / (a*c);
	bet  = asin(bet);
	gam  = 2.0 * area / (a*b);
	gam  = asin(gam);
	x = (alph * xov[1] + bet * xov[2] + gam * xov[3])/pi;
	y = (alph * yov[1] + bet * yov[2] + gam * yov[3])/pi;
	fprintf(f1,"strong point of first three points has coordinates %lf,%lf\n",x,y);
}