$ontext ===============================================================

LAST MODIFIED: 2008/11/07

FILENAME: pproffile.gms

PURPOSE:  Creates a performance profile for a number solvers (up to
          eight) to compare their relative performance. Only optimal
          solutions are considered succesful.

	  User must input one tracefiles for each solver. Currently
          works for LP, RMIP, NLP, RMINLP, MIP, and MINLP  type models.

	  Each trace file can have model/solve information for only 1
          solver. For example trace1.trc can have info for all LPs in
          PerformanceLib for solver A, and trace2.trc can have
          information for all LPs using solver B.

PROFILE DESCRIPTION:
          Suppose that we define a performance ratio:

          rho(p,s) =    resource time for problem p for solver s
                      -------------------------------------------------
                      min resource time for problem p over all solvers

          which can be thought of the relative solver time for a problem
          with respect to the best possible time.

          Then the performance profile p_s(tau) for a solver s is the
          probability that a performance ratio rho is within a factor of
          tau of the best possible ratio.

INTERPRETATION OF PERFORMANCE PROFILE:
          If we are only interested in maximizing the number of "wins"
          (i.e. maximizing the number of models solved the fastest),
          then we need only consider p_s(1). If we are interested only
          in determining if a solver has a high probability of success
          we are interested in

                   lim        p_s(tau)
                   tau_>inf

REFERENCE:
          Dolan, E. and More, J., "Benchmarking optimiation software
          with performance profiles", Mathematical Programming, Vol 91
          (2), pp.201-213, 2002.

REQUIRED INPUTS:
	  --trace1:    GAMS trace file of solver 1
	  --trace2:    GAMS trace file of solver 2 etc up to --trace8


OPTIONAL INPUTS:
          --delta:     considered optimal if objective value is within
                       --delta of best objective value, where --delta
                       is the relative objective function error
                       (default: Inf, i.e. ignore)
          --factor     factor for values of tau: 1, factor, factor^2,
                       factor^3 (default 1.25)
          --gaptol:    considered optimal if gap is smaller-equal --gaptol
                       (default: Inf, i.e. ignore)
          --modelfile: filename containing subset of models to consider.
                       Useful if only a subset of models from the
                       tracefile are to be used in the solver square
                       comparison.
          --numpts:    number of data points for each solver for
                       performance profile (default: 30)
          --objmin:    minimum obj. threshhold. If both solvers report
		       an obj. less than the threshhold, the objdiff
                       is computed as absolute difference, otherwise
                       use relative difference (default 0.1).
          --optonly:   Plots only if model is truly optimal (not only
                       feasible). Set to1 if only optimal models.
                       Otherwise considers both feasible and truly optimal
                       models. (default: 0)
          --outdir:    output directory (default: no directory)
          --outfile:   output filename (default: profile.txt)
	  --resmin:    minimum resource time threshold. If a solver
                       reports a resource time below --resmin, then
                       it is increased to resmin before being used
                       in ratios. (default 0.05).
          --total:     indicator if total models to be solved by any
                       solver should be displayed. Use --total=0
                       if to be ignored.
                       (default: 1)


REMARKS:  Each tracefile can contain information from only one
          (1) solver for a given set of models. For example,
          tracefiles containing solver A and solver B will cause
          error message upon compilation.

$offtext ==============================================================

*=== Error checks
$if not set trace1     $goto errors_notrace1
$if not set trace2     $goto errors_notrace2
$if not exist %trace1% $goto errors_trace1notexist
$if not exist %trace2% $goto errors_trace2notexist


*=== Set default options
$if not set objmin     $set objmin 1e-1
$if not set delta      $set delta Inf
$if not set gaptol     $set gaptol Inf
$if not set optonly    $set optonly 0
$if not set total      $set total 1
$if not set numpts     $set numpts  30
$if not set outfile    $set outfile "profile.txt"
$if not set factor     $set factor  1.2675
$if not set resmin     $set resmin  5e-2


$if not set colselect $setglobal colselect res
$setglobal tracecols    direction,modelstat,solvestat,obj,objest,%colselect%
$setglobal alltracecols modelname,modeltype,solvername,nlpdef,mipdef,%tracecols%

*== Create master trace file and read in results using readtrace.gms
$echo > %gams.scrdir%tracedata
$if set trace1 $call 'awk -f readtrace.awk -v gamscolumns="%alltracecols%" %trace1% >> %gams.scrdir%tracedata'
$if set trace2 $call 'awk -f readtrace.awk -v gamscolumns="%alltracecols%" %trace2% >> %gams.scrdir%tracedata'
$if set trace3 $call 'awk -f readtrace.awk -v gamscolumns="%alltracecols%" %trace3% >> %gams.scrdir%tracedata'
$if set trace4 $call 'awk -f readtrace.awk -v gamscolumns="%alltracecols%" %trace4% >> %gams.scrdir%tracedata'
$if set trace5 $call 'awk -f readtrace.awk -v gamscolumns="%alltracecols%" %trace5% >> %gams.scrdir%tracedata'
$if set trace6 $call 'awk -f readtrace.awk -v gamscolumns="%alltracecols%" %trace6% >> %gams.scrdir%tracedata'
$if set trace7 $call 'awk -f readtrace.awk -v gamscolumns="%alltracecols%" %trace7% >> %gams.scrdir%tracedata'
$if set trace8 $call 'awk -f readtrace.awk -v gamscolumns="%alltracecols%" %trace8% >> %gams.scrdir%tracedata'
$set tracefile %gams.scrdir%tracedata
$include readtrace.gms
$call 'rm -f %tracefile%'

$ifthen set modelfile
*=== If only the subset of models as listed in "modelfile" is to be analyzed
Set m(*) the subset of models to consider /
$onempty
$include %modelfile%
$offempty
/;
$else
*=== If all models in trace files are to be analyzed
Set m(*);
m(modelname) = yes;
$endif


*=== Check to see if all models are represented in all trace files.
*=== If not, send warning message
execute 'rm -f status.log'
file fput /status.log /;
put fput;

Parameter warning_models /0/;
loop(modelname(m),
  loop(solvername,
     if(sum( (modeltype,nlpdef,mipdef),
       tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'modelstat'))=0,
         warning_models = 1;
         put "WARNING:  model " modelname.tl:0 " doesn't exist for solver "
              solvername.tl:0 /;
     );
  );
);


$if not set outdir $goto nooutdir
$set outfile %outdir%/%outfile%
execute 'mkdir %outdir%'
execute 'cp style.css %outdir%'
$label nooutdir


Parameter resourcemin  "min resource time threshold for ratio" / %resmin%/;

Alias(solvers,solvername)

Parameter restime(*,*)	  "resource usage"
          bestobj(*)      "best objective value"
          objdiff(*,*)    "relative obj diff wrt to best obj"
          gap(*,*)        "gap (obj estimate vs obj)"
          modstat         "model status"            /0/
          solstat         "solver status"           /0/
          resusd          "resource time used"      /0/
          resmax          "max resource time used"  /0/
          i               "index"                   /0/;


*=== Check if discrete models exist
Scalar    is_discrete/0/;
Parameter maximize(*) "indicator if maximization problem" ;
maximize(modelname(m))$sum( (modeltype,solvername,nlpdef,mipdef),
                 tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'direction') )
                      = yes;

Set discrete_mods / MIP, MINLP, MIQCP/;
is_discrete = sum((modeltype,discrete_mods)$sameas(modeltype, discrete_mods), 1);

display is_discrete;


*=== Get best obejctive value over all solvers
Set       cansolve(*,*)     "indicator if solver can solve particular model at all";
Parameter valid_obj(*,*)    "objective value for solver, if optimal or feasible"
loop((modelname(m),solvername),
   valid_obj(modelname,solvername) = Inf;
   cansolve(modelname,solvername) = NO;
)
scalar objval;
scalar bound;

loop( (modelname(m),modeltype,solvername,nlpdef,mipdef)$
                        tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'modelstat'),
   i = i+1;

   modstat =
      tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'modelstat');
   solstat =
      tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'solvestat');
   resusd  =
      tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'%colselect%');
   maximize(modelname) =
      tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'direction');
   objval = tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'obj');


*  === Check if the model can be solved at all (and by which solver). if so, assign
*  === objvalue for that solver for that model
   if(is_discrete,
      if( (modstat=1 or modstat=15) and solstat=1,
         cansolve(modelname,solvername) = yes;
      elseif( (modstat=2 or modstat=8 or modstat=16)
               and (solstat=1 or solstat=2 or solstat=3 or solstat=4 or solstat=5 or solstat=8) ),
         cansolve(modelname,solvername) = yes;
      );
   else
      if( (modstat=1 or modstat=15) and solstat=1,
         cansolve(modelname,solvername) = yes;
      elseif( (modstat=2 or modstat=7 or modstat=16 or modstat=17)
               and (solstat=1 or solstat=2 or solstat=3 or solstat=4 or solstat=5 or solstat=8) ),
         cansolve(modelname,solvername) = yes;
      );
   );
   if(cansolve(modelname,solvername),
      valid_obj(modelname,solvername) = objval;
   );
);



*=== Now get actual best objective value
loop( modelname(m),
   if( maximize(modelname),
      bestobj(modelname) = smax(solvername$(valid_obj(modelname,solvername)<>Inf), valid_obj(modelname,solvername));
   else
      bestobj(modelname) = smin(solvername$(valid_obj(modelname,solvername)<>Inf), valid_obj(modelname,solvername));

   );
);




*=== Get resource times for each (model,solver) pair
*=== Initially set time to infinity, in case not all trace files
*=== have data for all models.
scalar objdiff_ok;
scalar gap_ok;
restime(modelname(m),solvername) = Inf;

loop( (modelname(m),modeltype,solvername,nlpdef,mipdef)$
                              tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'modelstat'),
   i = i+1;
   objdiff_ok = NO;
   gap_ok = NO;

   modstat =
      tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'modelstat');
   solstat =
      tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'solvestat');
   resusd  =
      tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'%colselect%');
   maximize(modelname) =
      tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'direction');
   objval = tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'obj');
   bound = tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'objest');
   bound$(bound=NA) = Inf$(maximize(modelname)) - Inf$(not maximize(modelname));

*=== Use relative objective value error for |best obj|>objmin, but absolute
*=== error otherwise
   objdiff(modelname,solvername)$(abs[bestobj(modelname)] > %objmin%) =
            abs(objval-bestobj(modelname)) / max(abs(bestobj(modelname)),%objmin%);
   objdiff(modelname,solvername)$(abs[bestobj(modelname)] < %objmin%) =
            abs(objval-bestobj(modelname));
   if(%delta% < Inf,
      objdiff_ok$(objdiff(modelname,solvername) <= %delta%) = YES;
   else
      objdiff_ok$(objdiff(modelname,solvername) < Inf) = YES;
   );

*=== Use relative gap for |best bound|>objmin, but absolute gap otherwise
   if (cansolve(modelname, solvername),
     gap(modelname,solvername)$(abs[bound] >= Inf) = Inf;
     gap(modelname,solvername)$(abs[bound] < Inf and abs[bound] >  %objmin%) = (objval-bound) / max(abs(bound),%objmin%);
     gap(modelname,solvername)$(abs[bound] < Inf and abs[bound] <  %objmin%) =  objval-bound;
     gap(modelname,solvername)$(abs[bound] < Inf and maximize(modelname)) = -gap(modelname,solvername);
     if (gap(modelname,solvername) < -%gaptol%,
       put "WARNING: solver " solvername.tl " report a negative gap (" gap(modelname,solvername) ") for model " modelname.tl /;
       warning_models = warning_models + 2;
       gap(modelname,solvername) = Inf;
     );
    if(%gaptol% < Inf,
        gap_ok$(gap(modelname,solvername) <= %gaptol%) = YES;
     else
*       === if we have no gap tolerance, than we also accept results from solvers that do not report and objestimate
        gap_ok = YES;
*        gap_ok$(gap(modelname,solvername) < Inf) = YES;
     );
*=== a model is not counted as solved if the gap is too large
     cansolve(modelname,solvername)$(not gap_ok) = NO;
   else
     gap(modelname,solvername) = Inf;
     gap_ok = NO
   );


   if(is_discrete,
*=== Assign resource data only if data exists
      if (modstat and solstat,
         if( (modstat=1 or modstat=15) and solstat=1 and objdiff_ok and gap_ok,
	    restime(modelname,solvername) = max( resusd, resourcemin);
         elseif(not %optonly% and (modstat=2 or modstat=8 or modstat=16)
                 and (solstat=1 or solstat=2 or solstat=3 or solstat=4 or solstat=5 or solstat=8) and objdiff_ok and gap_ok),
            restime(modelname,solvername) = max( resusd, resourcemin);
*=== If model is not optimal, mark it with Infinite restime.
         else
	    restime(modelname,solvername) = Inf;
         );
      );
   else
*=== Assign resource data only if data exists
      if (modstat and solstat,
         if( (modstat=1 or modstat=15) and solstat=1 and objdiff_ok and gap_ok,
	    restime(modelname,solvername) = max( resusd, resourcemin);
         elseif(not %optonly% and (modstat=2 or modstat=7 or modstat=16 or modstat=17)
                  and (solstat=1 or solstat=2 or solstat=3 or solstat=4 or solstat=5 or solstat=8) and objdiff_ok and gap_ok),
            restime(modelname,solvername) = max( resusd, resourcemin);
*=== If model is not optimal, mark it with Infinite restime.
         else
	    restime(modelname,solvername) = Inf;
         );
      );
   );
);




Parameter resmin(*)    "Best resource usage for all solvers"
          resrel(*,*)  "Resource usage relative to the best = res / resmin" ;

resmin(modelname(m)) = smin(solvername, restime(modelname,solvername));
resmin(modelname(m))$(resmin(modelname)=0) = %resmin%;
resrel(modelname(m),solvername)$(resmin(modelname) lt Inf) =
       restime(modelname,solvername) / resmin(modelname);
resrel(modelname(m),solvername)$(resmin(modelname) eq Inf) = Inf;


display resrel;



*=== Define performance profile sets and parameters
Set point / p1*p%numpts% /;
Parameter cutfact    "factor for tau: tau_i+1 = tau_i*cutfact"  / %factor% /,
          cutcur     "initial tau for performance probability function" / 1 /,
          cut(point) "variable tau for performance probability function",
          nmod       "number of models";
nmod = sum(modelname(m),1);

Parameter pprof(*,*)  "performance profile (in percent)";


* === Compute variable tau for prob. function
loop(point,
    cut(point) = cutcur;
    cutcur = cutcur * cutfact
);


*=== Compute pprof: performance profile probability function
pprof('cutpoint',point) = cut(point);
pprof(solvername,point) = 100*sum(modelname(m)$(
            resrel(modelname,solvername) le cut(point)), 1)/nmod;
pprof(solvername,'Solved') = 100*sum(modelname(m)$(
            resrel(modelname,solvername) lt Inf), 1)/nmod;

*=== Compute supersolver result: percentage of models that can be solved
*=== by any solver
scalar is_solved
       numsolved /0/;

numsolved = 0;
loop(modelname(m),
   is_solved = 0;
   if(sum( solvername$cansolve(modelname,solvername),1),
      is_solved = 1;
   );
   numsolved = numsolved+is_solved;
);
display "++++", numsolved;

pprof('supersolver',point) = 100*numsolved / nmod;




*=== Write table to profile file
$if     set minratio $goto dominratio

file fout / "%outfile%" /;
put fout;

put fout "   Cutpoint":>12;
loop(solvername,
   put solvername.tl:>15;
);
if(%total%,
   put ' CAN_SOLVE':>15;
);
put /;

loop(point,
   put cut(point):>12:3;
   loop(solvername,
      put pprof(solvername,point):>15:3;
   );
   if(%total%,
      put pprof('supersolver',point):>15:3;
   );
   put /;
);

putclose fout;




*=== Create data in form for plotting
Set item / cut, profile /;
Parameter plotdata(*, point, item);
plotdata(solvername,point,"cut") = cut(point);
plotdata(solvername,point,"profile") = pprof(solvername,point);
if(%total%,
   plotdata('CAN_SOLVE',point,"cut") = cut(point);
   plotdata('CAN_SOLVE',point,"profile") = pprof('supersolver',point);
);

display plotdata;
$goto done


*=== Write table to profile file for pprof values > --minratio
$label dominratio
$if not set minratio $set minratio 0


file fout / %outfile%.txt /;
put fout;
fout.pw = 1000;


put 'Min ratio = %minratio%'/;

put 'Reformulation/Solver':27 'Ratio_1':15 'Cutpoint':10 'Ratio_inf':15 /;

loop((point,solvername)$(ord(point)=%numpts% and
          pprof(solvername,point)>%minratio% ),
   put solvername.tl:23  pprof(solvername,'p1'):12:5  cut(point):15
       pprof(solvername,point):10:5 /;
)


putclose fout;


$label done



*=== If trace data does not exist for all models for all solvers, then
*=== give warning.
if(warning_models,
   execute "echo ---"
   execute "echo --- WARNING: trace data does not exist for all models for"
   execute "echo ---          all solvers. See the file status.log for details."
   execute "echo ---"
);

$goto noerrors


*=== Error messages
$label errors_notrace1
$log ---   ABORTING --trace1 option required: Trace file not specified..."
$abort "Aborting because --trace1 option not specified"

$label errors_notrace2
$log ---   ABORTING --trace2 option required: Trace file not specified..."
$abort   "Aborting because --trace2 option not specified"

$label errors_trace1notexist
$log "---     ABORTING --trace1 file not found..."
$abort "Aborting because trace file 1 not found"

$label errors_trace2notexist
$log ---     ABORTING --trace2 file not found...
$abort "Aborting because trace file 2 not found"


$label noerrors