/*
**             MINEPS
**
**             Version 1.0
**
**             Estimate smallest epsilon so that epsilon-core is not empty
**             in two-dimensional space
**
**             (c) 2003 by Thomas Braeuninger, University of Konstanz,
**             Germany; email: thomas.braeuninger@uni-konstanz.de
**
**
**             citation:
**             Thomas Braeuninger (2007) 'Stability in Spatial Voting Games WITH
**             Restricted Preference Maximizing", Journal of Theoretical
**             Politics 19(2): 173-91.
**
**
** Purpose:    MINEPS computes the yolk and/or approximates the minimal epsilon that 
**             guarantees that the epsilon-core is not empty in a two-dimensional space
**             Ideal points of committee members are provided by the user or
**             selected randomly (uniform or normal distribution)
**             
**             
**
** Input:      (optional)
**             dataset of actors' ideal points
**                m x 2 matrix:    ideal points of m actors in two dimensions
**
**
** Output:     mineps01.log
**             case number of experiment (simul), committee size (n)
**
**             mineps02.log
**             case number of experiment (simul), number of iterations for computing minimal epsilon (it), committee size (n),
**             majority rule threshold (thres), yolk center x (ycx), yolk center y (ycy)
**             yolk radius (yr), epsilon-core center x (ccx), epslion core center y (ccy), minimal epsilon (eps),
**             ratio between minimal epsilon and yolk radius (ratio), approximate area of epsilon-core (area)
**
**             mineps03.log
**             case number of experiment (simul), mean of ideal points in x (mean_x), mean of ideal points in y (mean_y),            
**             standard deviation of ideal points in x (std_x), standard deviation of ideal points in y (std_y),            
**             minimum of ideal points x (min_x), minimum of ideal points in y (min_y),            
**             maximum of ideal points x (max_x), maximum of ideal points in y (max_y)            
**
**             mineps04.log
**             case number of experiment (simul), ideal point in x (x), ideal point in y (y) 
**
**             mineps05.log
**             error report
*/







/*
** setbegin - define globals and open output files
*/

   new;

   @ variables @
   simul = 4000;              @ number of simulations @
   committee_min = 3;         @ minimal size of committee @
   committee_max = 101;       @ maximal size of committee @
   space_size = 10;           @ two-dimensional space is [0,space_size] x [0,space_size] (for uniform distribution, only) @
   logprec = 4;               @ precision of epsilon relative to yolk radius (in log) @
   maxiteration = 150;        @ max. number of iterations @
   decrease = 1;              @ factor to decrease epsilon by incircle radius @
   rose = 16;                 @ number of directions of grid star @

   @ constants @
   precfac = 10^logprec;
   maxepsilon = 1/(2+sqrt(3));
   rdeflator = 1;
   sim = 0;
   T=0;
   it=0;
   @ open output files @
   show_tmpvar = 0;
   outwidth 256;
   format /rd 10,7;
   name1 = "mineps01.log";
   name2 = "mineps02.log";
   name3 = "mineps03.log";
   name4 = "mineps04.log";
   name5 = "mineps05.log";
   output file = ^name1 on;
   "     simul          n";
   output file = ^name2 on;
   "     simul         it          n      thres        ycx        ycy";;
   "         yr        ccx        ccy        eps      ratio           area";;
   "       time        eyc";
   output file = ^name3 on;
   "     simul     mean_x     mean_y      std_x      std_y      min_x      min_y";;
   "      max_x      max_y";
   output file = ^name4 on;
   "     simul          x          y";
   output file = ^name5 on;
   "error report";
   output off;

goto mainproc;





/*
**
** PROC datadraw
**
** Creates uniform or normal distribution of ideal points of an odd number
** of actors (committee_min < m < committee_max) in a two-dimensional
** space [0,space_size] x [0,space_size]
**
** FORMAT
**         x     = datadraw(un) 
** 
** INPUT
**         un    = scalar, type of distribution (0=uniform, 1=normal)
** OUTPUT
**         x     = N x 1 matrix, ideal points in R2
** GLOBALS
**         committee_min = scalar, minimum number of actors
**         committee_max = scalar, maximum number of actors
**         space_size    = scalar, size of space (must be positive)
**
*/

proc datadraw(un);

   local m,x;

   m = 0;
   do until m/2 /= round(m/2) AND m>=committee_min AND m<=committee_max;
      m = round(rndu(1,1)*(committee_max-committee_min)+committee_max);
   endo;
   if un == 0;
      x = round(rndu(m,2)*space_size*1000000)/1000000;
   else;
      x = round(rndn(m,2)*space_size*1000000)/1000000;
   endif;
   retp(x);

endp;




/*
**
** PROC incircle_app
**
** Approximates center and radius of the incircle of a convex polygon 
** by searching along the rays of a star at the center of gravity
**
** FORMAT
**         {c,r} = incircle_app(expts) 
** 
** INPUT
**         expts = K x 2 matrix, extremal points of convex polygon in R2
** OUTPUT
**         c     = 2 x 1 column vector, center of incircle
**         r     = scalar, radius of incircle
**
*/

proc(2) = incircle_app(excore);

   local n,z,r,z_new,i,ii,test,rtest;

   n = rows(excore);
   z = sumc(excore)./n;
   r = mindisthull(z',excore);
   z_new = z;
   i = 0;
   do until i == n;
      i = i + 1;
      ii = 0;
      do until ii == 10;
         ii = ii + 1;
         test = z' + ii*(excore[i,.]-z')/10;
         rtest = mindisthull(test,excore);
         if rtest > r;
            r = rtest;
            z_new = test';
         endif;
      endo;
   endo;
   retp(z_new,r);

endp;




/*
**
** PROC mindisthull
**
** Calculates minimal distance between elements of a convex polygon in R2
** to the boundary of the polygon
**
** FORMAT
**         d     = mindisthull(point,expts) 
** 
** INPUT
**         point = N x 2 matrix, N points in R2
**         expts = K x 2 matrix, extremal points of convex hull in R2
** OUTPUT
**         d     = N x 2 matrix, minimal distance of points to boundary
**
*/

proc mindisthull(point,excore);

   local n,k,abstand,excore1,ii,i,v,dv,r,mindist,lambda,t,y,minabs;

   n=rows(point);
   k=rows(excore);
   @ matrix of distances to edges and vertices @
   abstand=zeros(n,2*k)+999;

   ii=0;
   do until ii==n;
      ii=ii+1;
      excore1=excore|excore[1,.];
      i=0;
      do until i==k;
         i=i+1;
         @ distance to edges @
         t=(excore1[i+1,.]-excore1[i,.])';
         if (t'*t) NE 0;
         lambda= (t' * (point[ii,.]-excore1[i,.])') / (t' * t);
            if ((lambda GE 0) AND (lambda le 1));
               y=(excore1[i,.]+lambda*(excore1[i+1,.]-excore1[i,.]))';
               abstand[ii,i]=sqrt((point[ii,.]'-y)'*(point[ii,.]'-y));
            endif;
         endif;
         @ distance to vertices @
         abstand[ii,i+k]=
           sqrt((point[ii,.]'-excore[i,.]')'*(point[ii,.]'-excore[i,.]'));
      endo;
   endo;
   minabs=minc(abstand');
   retp(minabs);

endp;




/*
**
** PROC winempty_e
**
** Tests whether or not the epsilon-winset is empty by checking
** for all intersections of any two indifference curves
**
** FORMAT
**         dummy = winempty_e(sq,x,T,eta) 
** 
** INPUT
**         sq    = 1 x 2 vector, status quo
**         x     = N x 2 matrix, ideal points
**         T     = scalar, voting threshold
**         eta   = scalar, eta
** OUTPUT
**         dummy = scalar, denoting wether winset is empty or not 
**		   (0=empty, 1=not)
**
*/

proc winempty_e(sq,x,T,eta);

   local m,dxsq,i,j,vij,dij,nij,z,h,a,dist,votes;

   m = rows(x);
   dxsq = sqrt(sumc((sq-x)'.*(sq-x)'))-eta;
   i = 0;
   do until i ==m;
      i = i+1;
      if dxsq[i] < 0;
         dxsq[i] = 0;
      endif;
   endo;
   i=0;
   do until i == m;
      i = i + 1;
      j = 0;
      do until j == m;
         j = j + 1;
         vij = x[j,.] - x[i,.];
         dij = sqrt(vij*vij');
         if ((dij == 0) OR (dxsq[i] + dxsq[j] - dij <= 0));
            goto stepnext;
         endif;
         @ normal vector from xi to xj @
         nij = vij ./ dij;
         @ distances of intersection to xi @
         z = ( (dxsq[i]^2 - dxsq[j]^2) / dij + dij ) / 2;
         h = sqrt(dxsq[i]^2 - z^2);
         if iscplx(h) == 1;
            if abs(imag(h)) > .0000000001;
               output file = ^name5 on;
               screen off;
               "error, complex h ";;
               h;
               screen on;
               output off;
            endif;
            h = real(h);
         endif;
         @ a is right-turn intersection of preferred-to-sets of i and j @
         a = x[i,.] + nij .* z + ((nij[2]~nij[1]) .* (1~-1)) .* h;
         dist = dxsq - sqrt( sumc((x-a)'.*(x-a)') );
         votes = sumc(dist .>= -.00000000001);
         if votes >= T;
             retp(1);
         endif;
         stepnext:
      endo;
   endo;
   retp(0);

endp;





/*
** orient - gives orientation of point with respect to hyperplane H
*/

proc orient(H,point);

   local a,b,c,dis,sign,factor;

   factor=10^10;
   a=H[.,1];
   b=H[.,2];
   c=H[.,3];
   dis = (b-a)~(c-a)~(point-a);
   dis = det(dis);
   dis = round(dis*factor) / factor;
@show
"det(dis)   sign";
@
   sign = (dis>=0)-(dis<=0);
@format /rd 20,15;
dis~sign;
@
   retp (sign);

endp;



/*
** lire0 - counts unicameral weighted votes on hyperplane H,
** in H+ and in H-
*/

proc lire0(H,x,g);

   local i,z,zz,m;

   m = rows(x);
   zz=zeros(3,1);
@show
" ";
"this is lire0";
@
   i=0;
   do until i==m;
      i=i+1;
@show
"i";;
i;
@
      z=orient(H,x[i,.]');
@show
"z";;
z;
@
      zz[z+2,.]=zz[z+2,.]+g[i,.];
   endo;
@show
"this is end lire0";
@
   retp(zz);

endp;



/*
** bm2_un - limiting median hyperplane for two dimensions
**
** (finds median also if three or more points are on the line,
** procedure is for unweighted votes, only)
*/

proc bm2_un(number,th);

   local z,zli,zre;

   z=0;
   @ votes in H- and on H @
   zli=number[1,.]+number[2,.];
   @ votes in H+ and on H @
   zre=number[3,.]+number[2,.];
   if ( (zli>=th) AND (zre>=th) );
      z=9;
   endif;
   retp(z);

endp;






/*
** hesse - calculates normal vector of hyperplane
*/

proc(2) = hesse(H);

   local a,b,n,p,hn;

   a = H[.,1] - H[.,3];
   b = H[.,2] - H[.,3];
   n = {0,0,0};
   n[1] = a[2]*b[3] - a[3]*b[2];
   n[2] = a[3]*b[1] - a[1]*b[3];
   n[3] = a[1]*b[2] - a[2]*b[1];
   hn = (-1) .* n ./ (sqrt(n'*n));
   p = (-1) .* (hn' * H[.,3]);
   retp (hn,p);

endp;



/*
**
** PROC triangle
**
** Calculate vertices of triangle given three lines in 2-point form
**
** Formula used is taken from mathworld.wolfram.com, keyword 
** "line-line intersection". 
** Procedure checked 11.3.03: works
**
** FORMAT
**         {A,B,C,err} = triangle(a,b,c,d,e,f) 
** 
** INPUT
**         a     = 1 x 2 row vector (or 2 x 1 column vector), 
**                 first point of first line
**         b     = 1 x 2 row vector (or 2 x 1 column vector), 
**                 second point of first line
**         c     = 1 x 2 row vector (or 2 x 1 column vector), 
**                 first point of second line
**         d       etc.
**         e       etc.
**         f       etc.
** OUTPUT
**         A,B,C = 1 x 2 row vectors (or 2 x 1 column vectors), 
**                 first, second and third vertex
**         err   = scalar, error code
**                          0 = no error, three distinct vertices
**                         -2 = all lines intersect in one point,
**                              {A,A,(-9999,-9999),-2} is returned,
**                         -1 = no intersection, at least two lines 
**                              are parallel, {(-9999,-9999),(-9999,
**                              -9999),(-9999,-9999),-2} is returned
*/

proc(4) = triangle(a,b,c,d,e,f);

   local o,p1,p2,vert,i,ii,denom,D1,D2,AA,BB,CC,errorcode;

   @ makes sure that points are horizontally oriented @
   if rows(a) == 2;
      o = 1;
      a = a';
      b = b';
      c = c';
      d = d';
      e = e';
      f = f';
   else;
      o = 0;
   endif;
   p1 = a|c|e;
   p2 = b|d|f;
   vert = zeros(3,2)-9999;
   errorcode = 0;
   i = 0;
   do until i == 2;
      i = i+1;
      ii = i;
      do until ii == 3;
         ii = ii+1;
         denom = (p1[i,.]-p2[i,.]) | (p1[ii,.]-p2[ii,.]);
         denom = det(denom);
         if denom /= 0;
            D1 = det(p1[i,.]|p2[i,.]);
            D2 = det(p1[ii,.]|p2[ii,.]);
            vert[i+ii-2,1] =
              det((D1~(p1[i,1]-p2[i,1]))|(D2~(p1[ii,1]-p2[ii,1]))) /denom;
            vert[i+ii-2,2] =
              det((D1~(p1[i,2]-p2[i,2]))|(D2~(p1[ii,2]-p2[ii,2]))) /denom;
         else;
            @ at least two lines are parallel @
            errorcode = -1;
         endif;
      endo;
   endo;
   AA = vert[1,.];
   BB = vert[2,.];
   CC = vert[3,.];
   @ three lines intersect in one point @
   if AA == BB AND BB == CC AND errorcode == 0;
      errorcode = -2;
   endif;
   if ((errorcode == -1) AND (AA[1]==-9999 OR BB[1]==-9999 OR CC[1]==-9999));
      if AA[1]==-9999;
         AA = CC;
         CC = {-9999 -9999};
      endif;
      if BB[1]==-9999;
         BB = CC;
         CC = {-9999 -9999};
      endif;
   endif;
   if o == 1;
      retp(AA',BB',CC',errorcode);
   else;
      retp(AA,BB,CC,errorcode);
   endif;

endp;



/*
**
** PROC yolk
**
** Calculates center and radius of yolk
**
** FORMAT
**         {Yc,Yr} = yolk(x,thres) 
** 
** INPUT
**         x     = N x 2 matrix, ideal points
**         T     = scalar, voting threshold
** OUTPUT
**         Yc    = 1 x 2 row vector, yolk center
**         Yr    = scalar, yolk radius
**
*/

proc (2) = yolk(x,thres);

   local m,g,Hplane,i,ii,iii,a,b,c,v,w,H,number,g1,g2,mm,
         lh,test,vekt,offset,med,nmed,Area,radi,radimax,
         radimaxtriangle,errorcode,
         AA,BB,CC,DD,s,Ycenter,A1,B1,C1;

   @ median lines @
   m = rows(x);
   x = x~zeros(m,1);
   g = ones(m,1);
   Hplane = zeros(1,4);
   med = zeros(1,2);
   i=0;
   do until i==m;
      i=i+1;
      ii=i;
      do until ii==m;
         ii=ii+1;
         @ select two points a and c, b is arbitrary @
         a=x[i,.]';
         b=x[i,.]'; b[3]=1;
         c=x[ii,.]';
         v=a-c;
         w=b-c;
         H=a~b~c;
         @ count actors left, right and on the line @
         number=lire0(H,x,g);
         @ decide whether a line is a median line @
         mm=bm2_un(number,thres);
         lh = {0,0,0};
         lh[1] = v[2]*w[3] - v[3]*w[2];
         lh[2] = v[3]*w[1] - v[1]*w[3];
         lh[3] = v[1]*w[2] - v[2]*w[1];
         test=sumc(abs(lh));
    	 if (test>0) AND (mm/=0);
            if (mm==1) OR (mm==-1);
               {vekt,offset}=hesse(H);
               offset=mm.*offset;
               vekt=mm.*vekt;
               Hplane = Hplane | (vekt'~offset);
               med = med | (i~ii);
            endif;
            if (mm==9);
               {vekt,offset}=hesse(H);
               Hplane = Hplane | (vekt'~offset);
               med = med | (i~ii);
            endif;
         endif;
      endo;
   endo;
   med = med[2:rows(med),.];   @ median lines (number of two actors) @
   nmed = rows(med);           @ number of median lines @
   x = x[.,1 2];

   @ largest circle @
   radimax = 0;
   radimaxtriangle = {0 0 0};
   i=0;
   do until i==nmed;
      i=i+1;
      ii=i;
      do until ii==nmed;
         ii=ii+1;
         iii=ii;
         do until iii==nmed;
            iii=iii+1;
            @ AA,BB,CC are vertices of triangle of three median lines @
            {AA,BB,CC,errorcode} = triangle(x[med[i,1],.],x[med[i,2],.],
              x[med[ii,1],.],x[med[ii,2],.],x[med[iii,1],.],x[med[iii,2],.]);
            if errorcode /= 0;
               if AA==BB AND CC==-9999 AND radimax == 0;
                  radimaxtriangle = i~ii~iii;
               endif;
               break;
            endif;
            @ a,b,c are length of sides of triangle @
            a = sqrt(sumc( ((BB-CC).*(BB-CC))' ));
            b = sqrt(sumc( ((AA-CC).*(AA-CC))' ));
            c = sqrt(sumc( ((AA-BB).*(AA-BB))' ));
            s = (a+b+c)/2;
            @ area of triangle @
            Area = sqrt(s*(s-a)*(s-b)*(s-c));
            @ radius of incircle @
            radi = Area/s;
            if radi > radimax;
               @ save maximum radius and number of the three median lines @
               radimax = radi;
               radimaxtriangle = i~ii~iii;
            endif;
         endo;
      endo;
   endo;

   @ yolk center @
   @ consider triangle that determines yolk @
   @ i,i,iii are median line numbers @
   i   = radimaxtriangle[1];
   ii  = radimaxtriangle[2];
   iii = radimaxtriangle[3];
   {AA,BB,CC,errorcode} = triangle(x[med[i,1],.],x[med[i,2],.],
                                   x[med[ii,1],.],x[med[ii,2],.],
                                   x[med[iii,1],.],x[med[iii,2],.]);
   @ yolk has zero radius, all median lines intersect in one point @
   if AA==BB AND CC==-9999 AND radimax == 0;
      retp(AA,0);
   else;
   @ incenter is intersection of all angle-bisectors AA A1 etc @
   a = CC-BB;
   b = CC-AA;
   c = BB-AA;
   a = a ./ sqrt(a*a');
   b = b ./ sqrt(b*b');
   c = c ./ sqrt(c*c');
   A1 = (AA+b+AA+c)/2;
   B1 = (BB+a+BB-c)/2;
   C1 = (CC-a+CC-b)/2;
   {Ycenter,A1,B1,errorcode}=triangle(AA,A1,BB,B1,CC,C1);
   retp(Ycenter,radimax);
   endif;

endp;





/*
**
** PROC minepsilon
**
** Approximates minimum epsilon so that epsilon-core is not empty. Returns
** both epsilon and (unique) core element
**
** FORMAT
**         {Cc,eps,it,area,eyc} = minepsilon(x,T,Yc,Yr)
** 
** INPUT
**         x     = N x 2 matrix, ideal points
**         T     = scalar, voting threshold
**         Yc    = 1 x 2 row vector, yolk center
**         Yr    = scalar, yolk radius
** OUTPUT
**         Cc    = 1 x 2 row vector, unique core element of e-core
**         eps   = scalar, minimum epsilon
**         it    = scalar, number of iterations processed
**                         -1 = procedure did not converge
**                         -4 = error, yolk center is e-core with
**                              e==2*radius of yolk
**         area  = scalar, log of estimated area of e-core after last
**                 estimation
**         eyc   = scalar, minimum epsilon so that Yc is element of e-core
**
*/


proc (5) = minepsilon(x,T,Yc,Yr);

   local csry,m,prec,c,r,it,i,S16,incr1,test,cobor,Area,
         aa,bb,cc,s,rad,raddiv,tmp,startrad,eyc,
         errorcode,c_new,r_new,rose1,decrease1,problem,logArea;

   csry = csrlin;
   print "         N  iteration           core center  core area  increment";;
   print "    epsilon";
   print "---------- ---------- --------------------- ---------- ----------";;
   print " ----------";
   raddiv = 10;
   startrad = 2;

   step0:
   m = rows(x);
   prec = Yr/precfac;
   c = Yc;
   r = Yr*startrad;
   eyc = r;
   it = 0;
   Area = 0;
   rad = r/100;
   rose1 = rose;
   decrease1 = decrease;
   S16 = seqa(0,1,rose);
   S16 = sin((S16/rose)*2*pi)~cos((S16/rose)*2*pi);
   if winempty_e(c,x,T,r) == 1;
      output file = ^name5 on;
      screen off;
      "error, e-winset of yolk center is not empty! sim#/it#/epsilon ";;
      sim~it~r;
      screen on;
      output off;
      if startrad >= 10;
         it = -4;
         goto stepB;
      endif;
      startrad = startrad*2;
      it = -4;
      goto step0;
   endif;
   if winempty_e(Yc,x,T,r) == 0;
      eyc = r;
   endif;

   stepA:
   locate csry+2,1;
   m~it~c~Area~rad~r;
   it = it+1;
   @ estimate border of e-core all directions @
   cobor = zeros(rose1,2)+c;
   i=0;
   do until i == rose1;
      i = i+1;
      incr1 = r;
      do until incr1 < rad/raddiv;
         test = cobor[i,.] + incr1 .* S16[i,.];
         if winempty_e(test,x,T,r) ==0;
            cobor[i,.] = test;
         else;
            incr1 = incr1 /3;
         endif;
      endo;
   endo;
   @ approximate area of e-core @
   i = 0;
   Area = 0;
   cobor = cobor | cobor[1,.];
   do until i == rose1;
      i = i+1;
      aa = sqrt(sumc( ((cobor[i,.]-cobor[i+1,.])
           .*(cobor[i,.]-cobor[i+1,.]))' ));
      bb = sqrt(sumc( ((cobor[i,.]-c).*(cobor[i,.]-c))' ));
      cc = sqrt(sumc( ((c-cobor[i+1,.]).*(c-cobor[i+1,.]))' ));
      s = (aa+bb+cc)/2;
      @ area of triangle @
      Area = Area + sqrt(s*(s-aa)*(s-bb)*(s-cc));
   endo;
   cobor = cobor[1:rose1,.];
   @ end if Area is small @
   if Area == 0;
      raddiv = precfac;
      goto step0;
   endif;
   if Area < pi*prec^2 OR it >= maxiteration;
      goto stepB;
   endif;
   @ new estimate of core center @
   c_new = (sumc(cobor)./rows(cobor))';
   rad = mindisthull(c_new,cobor);
   r_new = r - rad*decrease1;
   do until winempty_e(c_new,x,T,r_new) == 0;
      r_new = (r + r_new)/2;
   endo;
   if r - r_new < rad/20;
      {c_new,rad} = incircle_app(cobor);
      c_new = c_new';
      r_new = r - rad*decrease1;
      do until winempty_e(c_new,x,T,r_new) == 0;
         r_new = (r + r_new)/2;
      endo;
   endif;
   if r - r_new < rad/1000;
      @ error report to file @
      output file = ^name5 on;
      screen off;
      "error, progress is less than r/1000, sim#/iter# ";;
      sim~it;
      screen on;
      output off;
      problem = problem + 1;
      it = it + .001;
      if problem >= 2;
         rose1 = rose1*2;
      endif;
      rdeflator = rdeflator*2;
      decrease1 = decrease1 * .8;
      if rose1 > 128;
         output file = ^name5 on;
         screen off;
         "error, e-winset of center is not empty! NO CONVERGENCE, sim#/iter#";;
         sim~it;
         screen on;
         output off;
         it = -1;
         goto stepB;
      endif;
      S16 = seqa(0,1,rose1);
      S16 = sin((S16/rose1)*2*pi)~cos((S16/rose1)*2*pi);
      goto stepA;
   else;
      problem = 0;
   endif;
   c = c_new;
   r = r_new;
   goto stepA;
 
   stepB:
   logArea = log(Area+10E-99);
   locate csry+2,1;
   m~it~c~Area~rad~r;
   retp(c,r,it,logArea,eyc);

endp;





/*
**
** main
**
*/

mainproc:

   screen on;
   cls;
   "MINEPS - Estimation of yolk and minimal epsilon of epsilon-core in two dimensions";
   "(c) 2003 Thomas Braeuninger, University of Konstanz, Germany";
   " ";   "Type of problem:";
   "(1) Compute yolk";
   "(2) Compute yolk and minimal epsilon";
   "(3) Simulate commitee majority voting and compute yolk";
   "(4) Simulate commitee majority voting and compute yolk and minimal epsilon";
   " ";   typ=con(1,1);
   screen off; 
   typ;
   screen on;
   cls;
   tstart = date;
   if (typ==1 OR typ==2);
      et = hsec;
      sim = 1;
      "Data file containing ideal points: ";
      name_xy=cons;
      name_xy;
      load x[] = ^name_xy;
      m = rows(x)/2;
      x=reshape(x,m,2);
      if round(m/2) /= round;
         "Error, committee size not odd!";
         end;
      endif;
      format /ld 1,0;
      "N =";;
      m;;
      format /rd 10,7;
      {Yc,Yr}=yolk(x,T);
      ", Yc = {";;
      Yc;;
      "} , Yr =";;
      Yr;
      if typ==1;
         {fp,fr,it,logA,eyc}=minepsilon(x,T,Yc,Yr);
         ratio = fr/Yr;
         "Yr/epsilon =";;
         ratio;
      endif;
      screen off;
      output file = ^name1 on;
      sim~m;
      output file = ^name2 on;
      if typ==1;
         sim~it~m~T~Yc~Yr;
      else;
         sim~it~m~T~Yc~Yr~fp~fr~ratio~logA~et~eyc;
      endif;
      output file = ^name3 on;
      sim~meanc(x)'~stdc(x)'~minc(x)'~maxc(x)';
      output file = ^name4 on;
      tmp = zeros(m,1)+sim;
      format /ma1;
      tmp~x;;
      format /mb1;
      output off;
      screen on;
      goto menueend;
   endif;
   if (typ==3 OR typ==4);
      "Number of experiments: ";
      simul=con(1,1);
      screen off; 
      simul;
      screen on;
      " ";
      "Minimal committee size: ";
      committee_min=con(1,1);
      screen off; 
      committee_min;
      screen on;
      " ";
      "Maximum committee size: ";
      committee_max=con(1,1);
      screen off; 
      committee_max;
      screen on;
      " ";
      "Uniform (0) or normal distribution (1) of ideal points: ";
      distrib=con(1,1);
      screen off; 
      distrib;
      screen on;
      if distrib==0;
        " ";
        "Size of policy space [0,?]x[0,?]: ";
        space_size=con(1,1);
        screen off; 
        space_size;
        screen on;
      endif;
      sim = 0;
      do until sim == simul;
        cls;
        sim = sim+1;
        et = hsec;
        x = datadraw(distrib);
        m = rows(x);
        T = m/2 + .5;
        format /ld 1,0;
        "sim =";;
        sim;;       
        ", N =";;
        m;;
        format /rd 10,7;
        {Yc,Yr}=yolk(x,T);
        ", Yc = {";;
        Yc;;
        "} , Yr =";;
        Yr;
        if typ==4;
          {fp,fr,it,logA,eyc}=minepsilon(x,T,Yc,Yr);
          ratio = fr/Yr;
          "Yr/epsilon =";;
          ratio;
        endif;
        et = hsec - et;
        screen off;
        output file = ^name1 on;
        sim~m;
        output file = ^name2 on;
        if typ==3;
          sim~it~m~T~Yc~Yr;
        else;
          sim~it~m~T~Yc~Yr~fp~fr~ratio~logA~et~eyc;
        endif;
        output file = ^name3 on;
        sim~meanc(x)'~stdc(x)'~minc(x)'~maxc(x)';
        output file = ^name4 on;
        tmp = zeros(m,1)+sim;
        format /ma1;
        tmp~x;;
        format /mb1;
        output off;
        screen on;
      endo;
      goto menueend;
   endif;
   menueend:
   tend = date;
   output file = ^name2 on;
   "elapsed time:  ";;
   etstr(ethsec(tstart,tend));;
   "\r";
   output off;

end;