/* Program to accompany the article
   "How to find a feasible point for a constrained optimization in SAS"
   by Rick Wicklin, published 19JUN2017 on The DO Loop blog,
   http://blogs.sas.com/content/iml/2017/06/19/find-feasible-point-optimization-sas.html
*/

proc iml;
con = {  0   0   .   .,    /* param min */
        10  10   .   .,    /* param max */
         3  -2  -1  10,    /* 3*x1 + -2*x2 LE 10 */
         5  10  -1  56,    /* 5*x1 + 10*x2 LE 56 */
         4   2   1   7 };  /* 4*x1 +  2*x2 GE  7 */

guess = {0 0};                /* arbitrary p-dimensional point */
call nlpfea(z, guess, con);  /* x0 is feasible point */
print guess[c={"x" "y"}], z[c={"Tx" "Ty"}];

start func(x);
   x1 = x[,1];   x2 = x[,2];
   return ( -(x1-3)##4 + -(x2-2)##2 + 0.1*sin(x1#x2));
finish;

opt = {1,    /* find maximum of function     */
       3};   /* print a little bit of output */

x0 = {0 0};
call nlpnra(rc, x_Opt, "func", x0, opt) blc=con;
print rc;

/* write feasible points to SAS data set */

NumPts = 36;
twopi = 2*constant('pi');
x=.; y=.; Tx=.; Ty=.;
create feasible var {"x" "y" "Tx" "Ty"};
do i = 1 to NumPts;
   x = 2.5 + 5*cos((i-1)/NumPts * twopi);  /* guess on circle */
   y = 2.5 + 5*sin((i-1)/NumPts * twopi);
   x = 3 + 4.2*cos((i-1)/NumPts * twopi);  /* guess on circle */
   y = 3 + 4.2*sin((i-1)/NumPts * twopi);
   call nlpfea(feasPt, x||y, con);         /* transform into feasible */
   Tx = feasPt[1]; Ty = feasPt[2];         
   append;
end;
close;
QUIT;


/* visualize the feasible region and feasible points */
data polygon;
ID = 1;
input PolyX PolyY;
datalines;
1.75   0
3.3333 0
5.3    2.95
0      5.6
0      3.5
;

data All;
set polygon feasible;
run;

ods graphics / width=500px height=480px;
title "Transform Guess into the Feasible Region";
proc sgplot data=All aspect=1;
   polygon id=id x=PolyX y=PolyY /lineattrs=(thickness=3) fill outline 
                                  name="F" legendlabel="Feasible Region";
   scatter x=x y=y / markerattrs=(symbol=CircleFilled size=6) name="G" legendlabel="Guess";
   scatter x=Tx y=Ty / markerattrs=GraphData2(symbol=CircleFilled size=10) name="P"  legendlabel="Transformed" ;
   vector x=Tx y=Ty / xorigin=x yorigin=y lineattrs=(color=black);
   xaxis display=(nolabel);
   yaxis display=(nolabel);
   keylegend "F" "G" "P";
run;