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

LAST MODIFIED: 15 June, 2004

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
          --ntrace1:   GENERIC data file 1
          --ntrace2:   GENERIC data file 2 up to --ntrace8


OPTIONAL INPUTS:
          --factor     factor for values of tau: 1, factor, factor^2, 
                       factor^3 (default 1.25)
	  --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).
          --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) 
          --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)
          --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).
          --topt4:     trace option 4. If specified, then assumes 
                       tracefile is traceopt=4 (default:none)
          --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.

          Users cannot input two tracefile numbers of different type
          with the same (n)trace number, i.e.

          --trace1=... --ntrace1 is not allowed, only --trace1=... --ntrace2=...


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


*=== Create trace reading routine (readtrace.gms)
$onecho > %gams.scrdir%readtrace.gms
$eolcom ,#
display '++++ TRACEOPT=3';


Set modelname, modeltype, solvername, nlpdef, mipdef;
Set col /julian, direction,eqnum,varnum,dvarnum,nz,nlnz, optfile,modelstat,solvestat,obj,objest,res,iter,dom,nodes/


Alias(u1,u2,u3,u4,u5,*);
Set domain(*,*,*,*,*);


Table tracedata(*,*,*,*,*,col) 
$ondelim
modelname, modeltype, solvername, nlpdef, mipdef, julian, direction,eqnum,varnum,dvarnum,nz,nlnz, optfile,modelstat,solvestat,obj,objest,res,iter,dom,nodes
$offlisting
$include %tracefile%
$onlisting
$offdelim
;


domain(u1,u2,u3,u4,u5)$sum( col, tracedata(u1,u2,u3,u4,u5,col)) = yes;

loop(domain(u1,u2,u3,u4,u5),
  modelname(u1)  = yes;
  modeltype(u2)  = yes;
  solvername(u3) = yes;
  nlpdef(u4)     = yes;
  mipdef(u5)     = yes;
);
$offecho 



*=== Convert non-GAMS data files to GAMS trace files
$onecho > %gams.scrdir%tconvawk 
{print $1 "," $2 "," $3 ",none,none,-1", $4, "-1,-1,-1,-1,-1,-1," $5 "," $6 "," $7 ",-1," $8 "," $9 ",-1,-1,#"} 
$offecho
$if set ntrace1 $call 'awk -F, -f "%gams.scrdir%tconvawk" %ntrace1% > %gams.scrdir%C_%ntrace1%'
$if set ntrace1 $set trace1 %ntrace1%
$if set ntrace2 $call 'awk -F, -f "%gams.scrdir%tconvawk" %ntrace2% > %gams.scrdir%C_%ntrace2%'
$if set ntrace2 $set trace2 %ntrace2%
$if set ntrace3 $call 'awk -F, -f "%gams.scrdir%tconvawk" %ntrace3% > %gams.scrdir%C_%ntrace3%'
$if set ntrace3 $set trace1 %ntrace3%
$if set ntrace4 $call 'awk -F, -f "%gams.scrdir%tconvawk" %ntrace4% > %gams.scrdir%C_%ntrace4%'
$if set ntrace4 $set trace1 %ntrace4%
$if set ntrace5 $call 'awk -F, -f "%gams.scrdir%tconvawk" %ntrace5% > %gams.scrdir%C_%ntrace5%'
$if set ntrace5 $set trace1 %ntrace5%
$if set ntrace6 $call 'awk -F, -f "%gams.scrdir%tconvawk" %ntrace6% > %gams.scrdir%C_%ntrace6%'
$if set ntrace6 $set trace1 %ntrace6%
$if set ntrace7 $call 'awk -F, -f "%gams.scrdir%tconvawk" %ntrace7% > %gams.scrdir%C_%ntrace7%'
$if set ntrace7 $set trace1 %ntrace7%
$if set ntrace8 $call 'awk -F, -f "%gams.scrdir%tconvawk" %ntrace8% > %gams.scrdir%C_%ntrace8%'
$if set ntrace8 $set trace1 %ntrace8%

$echo " " > %gams.scrdir%empty
$if not set trace1 $set trace1 %gams.scrdir%empty
$if not set trace2 $set trace2 %gams.scrdir%empty
$if not set trace3 $set trace3 %gams.scrdir%empty 
$if not set trace4 $set trace4 %gams.scrdir%empty 
$if not set trace5 $set trace5 %gams.scrdir%empty
$if not set trace6 $set trace6 %gams.scrdir%empty
$if not set trace7 $set trace7 %gams.scrdir%empty
$if not set trace8 $set trace8 %gams.scrdir%empty


$if not set objmin     $set objmin 1e-1
$if not set delta      $set delta Inf
$if not set optonly    $set optonly 0
$if not set total      $set total 1

$if not set modelfile  $goto nomodelfile


*=== 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
/;

$echo > %gams.scrdir%tracedata
$if not set ntrace1 $echo '$include %trace1%' >> %gams.scrdir%tracedata
$if not set ntrace2 $echo '$include %trace2%' >> %gams.scrdir%tracedata
$if not set ntrace3 $echo '$include %trace3%' >> %gams.scrdir%tracedata
$if not set ntrace4 $echo '$include %trace4%' >> %gams.scrdir%tracedata
$if not set ntrace5 $echo '$include %trace5%' >> %gams.scrdir%tracedata
$if not set ntrace6 $echo '$include %trace6%' >> %gams.scrdir%tracedata
$if not set ntrace7 $echo '$include %trace7%' >> %gams.scrdir%tracedata
$if not set ntrace8 $echo '$include %trace8%' >> %gams.scrdir%tracedata
$if set ntrace1 $echo '$include %gams.scrdir%C_%ntrace1%' >> %gams.scrdir%tracedata
$if set ntrace2 $echo '$include %gams.scrdir%C_%ntrace2%' >> %gams.scrdir%tracedata
$if set ntrace3 $echo '$include %gams.scrdir%C_%ntrace3%' >> %gams.scrdir%tracedata
$if set ntrace4 $echo '$include %gams.scrdir%C_%ntrace4%' >> %gams.scrdir%tracedata
$if set ntrace5 $echo '$include %gams.scrdir%C_%ntrace5%' >> %gams.scrdir%tracedata
$if set ntrace6 $echo '$include %gams.scrdir%C_%ntrace6%' >> %gams.scrdir%tracedata
$if set ntrace7 $echo '$include %gams.scrdir%C_%ntrace7%' >> %gams.scrdir%tracedata
$if set ntrace8 $echo '$include %gams.scrdir%C_%ntrace8%' >> %gams.scrdir%tracedata
$set tracefile %gams.scrdir%tracedata
$include %gams.scrdir%readtrace
$goto donemodelfile


*=== All models in trace files are to be analyzed
$label nomodelfile

$echo > %gams.scrdir%tracedata
$if not set ntrace1 $echo '$include %trace1%' >> %gams.scrdir%tracedata
$if not set ntrace2 $echo '$include %trace2%' >> %gams.scrdir%tracedata
$if not set ntrace3 $echo '$include %trace3%' >> %gams.scrdir%tracedata
$if not set ntrace4 $echo '$include %trace4%' >> %gams.scrdir%tracedata
$if not set ntrace5 $echo '$include %trace5%' >> %gams.scrdir%tracedata
$if not set ntrace6 $echo '$include %trace6%' >> %gams.scrdir%tracedata
$if not set ntrace7 $echo '$include %trace7%' >> %gams.scrdir%tracedata
$if not set ntrace8 $echo '$include %trace8%' >> %gams.scrdir%tracedata
$if set ntrace1 $echo '$include %gams.scrdir%C_%ntrace1%' >> %gams.scrdir%tracedata
$if set ntrace2 $echo '$include %gams.scrdir%C_%ntrace2%' >> %gams.scrdir%tracedata
$if set ntrace3 $echo '$include %gams.scrdir%C_%ntrace3%' >> %gams.scrdir%tracedata
$if set ntrace4 $echo '$include %gams.scrdir%C_%ntrace4%' >> %gams.scrdir%tracedata
$if set ntrace5 $echo '$include %gams.scrdir%C_%ntrace5%' >> %gams.scrdir%tracedata
$if set ntrace6 $echo '$include %gams.scrdir%C_%ntrace6%' >> %gams.scrdir%tracedata
$if set ntrace7 $echo '$include %gams.scrdir%C_%ntrace7%' >> %gams.scrdir%tracedata
$if set ntrace8 $echo '$include %gams.scrdir%C_%ntrace8%' >> %gams.scrdir%tracedata
$set tracefile %gams.scrdir%tracedata
$include %gams.scrdir%readtrace


Set m(*);
m(modelname) = yes;

$label donemodelfile


*=== 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,
  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 /;     
     );
  );
);



*=== Set default options
$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


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"  
          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)$sum( (modeltype,solvername,nlpdef,mipdef),
                 tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'direction') )     
                      = yes; 

Set discrete_mods / MIP, MINLP/;
loop(modeltype,
   is_discrete = sum(discrete_mods, sameas(modeltype, discrete_mods));
);

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,solvername),
   valid_obj(modelname,solvername) = Inf; 
   cansolve(modelname,solvername) = NO;
)
scalar objval;

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,'res');
   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) ), 
         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) ),
         cansolve(modelname,solvername) = yes;
      );
   );
   if(cansolve(modelname,solvername),
      valid_obj(modelname,solvername) = objval;
   );
);



*=== Now get actual best objective value
loop( modelname, 
   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;
restime(modelname,solvername) = Inf;

loop( (modelname(m),modeltype,solvername,nlpdef,mipdef)$
                              tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'modelstat'),
   i = i+1;
   objdiff_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,'res');
   maximize(modelname) = 
      tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'direction');
   objval = tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'obj');
  


*=== 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; 
   );

 
   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,
	    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) and objdiff_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,
	    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) and objdiff_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) = smin(solvername, restime(modelname,solvername));
resmin(modelname)$(resmin(modelname)=0) = %resmin%;
resrel(modelname,solvername)$(resmin(modelname) lt Inf) = 
       restime(modelname,solvername) / resmin(modelname);
resrel(modelname,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 = card(modelname);

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$(
            resrel(modelname,solvername) le cut(point)), 1)/nmod;
pprof(solvername,'Solved') = 100*sum(modelname$(
            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,
   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 ---"
);