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

LAST MODIFIED: 2008/11/07

FILENAME: square.gms

PURPOSE:  Creates a solve square to compare solve/fail
          results. User must input two tracefiles (one for each
          solver).

          Works for LP, RMIP, NLP, RMINLP, MIP, or 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.


REQUIRED INPUTS:
          --trace1:    GAMS trace file of solver 1
          --trace2:    GAMS trace file of solver 2

OPTIONAL INPUTS:
          --bnd:       relative objective function difference. Prints
                       "optimal" solution in boldface for HTML output
                       for one solver if the relative objective function
                       difference is greater than --bnd. (default 1e-5)
          --gaptol     up to which gap a model is counted as solved for a solver
                       (default: bnd)
          --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).
          --outdir:    output directory (default: no directory)
          --outfile:   root names of square output file (default is out
                       --> out_sqr.txt)
          --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.
          --resmin:    minimum resource time threshold. If both solvers
                       report a resource time below --resmin, then no
                       ratio is reported. In this case, both solver times
                       are considered the same (default 0.05).
          --useobjest  whether objective estimates should be used to
                       calculate gaps and solution quality (default: 0).


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.

          Discrete and continous models cannot be mixed in an analysis.

$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 option
$if not set outfile   $set  outfile out
$if not set bnd       $set bnd     1e-5
$if not set objmin    $set objmin  1e-1
$if not set resmin    $set resmin  5e-2
$if not set tsame     $set tsame   10
$if not set tfaster   $set tfaster 50
$if not set reslim    $set reslim -1
$if not set useobjest $set useobjest 0
$if not set gaptol    $set gaptol  %bnd%
$if not set uselog    $set uselog 0

$if not set modellib        $set modellib none
$if %modellib% == linlib    $set liburl "http://www.gamsworld.org/performance/plib"
$if %modellib% == minlplib  $set liburl "http://www.gamsworld.org/minlp/minlplib"
$if %modellib% == globallib $set liburl "http://www.gamsworld.org/global/globallib"


*== Create master trace file and read in results using readtrace.gms
$echo > %gams.scrdir%tracedata
$call 'awk -f readtrace.awk %trace1% >> %gams.scrdir%tracedata'
$call 'awk -f readtrace.awk %trace2% >> %gams.scrdir%tracedata'
$set tracefile %gams.scrdir%tracedata
$include readtrace.gms
$call 'rm -f %gams.scrdir%tracedata'

$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


*== If output directory specified then create directory and other hosuekeeping
$if not set outdir $goto nooutdir
$set outfile %outdir%/%outfile%
execute 'mkdir %outdir%'
execute 'cp style.css %outdir%'
$label nooutdir


Set bool  / no, yes /;
Parameter list_text /no 1, yes 2/;
Parameter list_html/no 1, yes 2/;


Parameter bound        "relative objc func difference threshold"
                        / %bnd% /
          resourcemin  "min resource time threshold for ratio"
                        / %resmin%/
;


*=== Possible solver outcomes
Set all_results / opt    "optimal",
                  lopt   "locally optimal"
                  feas   "feasible",
		  infeas "infeasible",
		  unbnd  "unbounded",
		  f      "fail"
                  nd     "no data"
/;


Set result_cont(all_results)  "possible results for cont. models"
    result_disc (all_results) "possible results for discrete models"
;

result_cont(all_results) = yes;
result_disc(all_results) = yes;
*result_disc("lopt")      = no;


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

Set result(*), resultA(*), resultB(*);

if(is_discrete,
   result(result_disc)  = yes;
   resultA(result_disc) = yes;
   resultB(result_disc) = yes;
else
   result(result_cont)  = yes;
   resultA(result_cont) = yes;
   resultB(result_cont) = yes;
);

* compute a reslim if not provided by user
Parameter reslim / %reslim% /;
if (reslim=-1,
  reslim=smax((modelname(m),modeltype,solvername,nlpdef,mipdef)$(not mapval(tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'res'))), tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'res'))
);
display reslim;

* check if we have objective estimates, so we can compute gaps
Parameter have_objest;
if (not %useobjest%,
  have_objest=0
else
  have_objest = sum((modelname(m),modeltype,solvername,nlpdef,mipdef)$(not mapval(tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'objest'))), 1);
);


Alias(solvers,solvername)

Parameter square(*,*,*)   GAMS solver results for square;
square(modelname(m),solvername,result) = 0;

Parameter modstat /0/;
Parameter solstat /0/;
Parameter mstat(*,*,*), sstat(*,*,*);
Parameter i /0/;


*=== Compute totals for square
loop( (modelname(m),modeltype,solvername,nlpdef,mipdef),
   i = i+1;

   modstat =
       tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'modelstat');
   solstat =
       tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'solvestat');

*=== Assign DISCRETE model square only if data exists
   if(is_discrete,
      if (modstat and solstat,
         if( (modstat=1 or modstat=15) and solstat=1,
            square(modelname,solvername,"opt") = yes;
            mstat(modelname,solvername,"opt") = modstat;
            sstat(modelname,solvername,"opt") = solstat;
         elseif(modstat=2),
            square(modelname,solvername,"lopt") = yes;
            mstat(modelname,solvername,"lopt") = modstat;
            sstat(modelname,solvername,"lopt") = solstat;
         elseif( (modstat=8 or modstat=16) and
                 (solstat=1 or solstat=2 or solstat=3 or solstat=4 or solstat=5 or solstat=8)),
            square(modelname,solvername,"feas") = yes;
            mstat(modelname,solvername,"feas") = modstat;
            sstat(modelname,solvername,"feas") = solstat;
         elseif( (modstat=4 or modstat=5 or modstat=10 or modstat=19)and solstat=1),
            square(modelname,solvername,"infeas") = yes;
            mstat(modelname,solvername,"infeas") = modstat;
            sstat(modelname,solvername,"infeas") = solstat;
         elseif( (modstat=3 or modstat=18) and solstat=1  ),
            square(modelname,solvername,"unbnd") = yes;
            mstat(modelname,solvername,"unbnd") = modstat;
            sstat(modelname,solvername,"unbnd") = solstat;
         else
            square(modelname,solvername,"f") = yes;
            mstat(modelname,solvername,"f") = modstat;
            sstat(modelname,solvername,"f") = solstat;
         );
      );
*=== Assign CONTINUOUS model square only if data exists
   else
      if (modstat and solstat,
         if( (modstat=1 or modstat=15) and solstat=1,
            square(modelname,solvername,"opt") = yes;
            mstat(modelname,solvername,"opt") = modstat;
            sstat(modelname,solvername,"opt") = solstat;
         elseif( (modstat=2)),
            square(modelname,solvername,"lopt") = yes;
            mstat(modelname,solvername,"lopt") = modstat;
            sstat(modelname,solvername,"lopt") = solstat;
         elseif( (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)),
            square(modelname,solvername,"feas") = yes;
            mstat(modelname,solvername,"feas") = modstat;
            sstat(modelname,solvername,"feas") = solstat;
         elseif( (modstat=4 or modstat=5 or modstat=19) and solstat=1),
            square(modelname,solvername,"infeas") = yes;
            mstat(modelname,solvername,"infeas") = modstat;
            sstat(modelname,solvername,"infeas") = solstat;
         elseif( (modstat=3 or modstat=18) and solstat=1),
            square(modelname,solvername,"unbnd") = yes;
            mstat(modelname,solvername,"unbnd") = modstat;
            sstat(modelname,solvername,"unbnd") = solstat;
         else
            square(modelname,solvername,"f") = yes;
            mstat(modelname,solvername,"f") = modstat;
            sstat(modelname,solvername,"f") = solstat;
         );
      );
   );
);



*=== Assign square for model/solver with no trace data
loop( (modelname(m),solvername),
   if(not sum[(modeltype,nlpdef,mipdef), tracedata(modelname,modeltype,solvername,nlpdef,mipdef,'modelstat')],
      square(modelname,solvername,"nd") = yes;
   );
);




*=== Compute result totals
parameter resulttotals(*,*) GAMS solver totals for square;
          resulttotals(resultA, resultB) = 0;


Set solverA, solverB;
Parameter i;
i=1;
loop(solvername,
   if(i eq 1,
      solverA(solvername) = solvername(solvername);
      i = 2;
   elseif(i eq 2),
      solverB(solvername) = solvername(solvername);
   );
);



Set solversused(*);
solversused(solvername) = solverA(solvername) + solverB(solvername);


*=== Result totals for square
loop(resultA,
   loop(resultB,
      resulttotals(resultA, resultB) =  sum((m,solverA,solverB),
          square(m,solverA,resultA) and square(m, solverB,resultB) );
   );
);



*==== Extract resource times
Parameter resA(*)   "resource time of solver A",
          resB(*)   "resource time of solver B"
          objA(*)   "objective value using solver A"
          objB(*)   "objective value using solverB"
          objEstA(*) "objective estimate using solver A"
          objEstB(*) "objective estimate using solver B"
          gapA(*)   "gap using solver A"
          gapB(*)   "gap using solver B"
          dir(*)    "direction of optimization: 0=min 1=max"
          modstatA(*) "model status of solver A"
          solstatA(*) "solver status of solver A"
          modstatB(*) "model status of solver B"
          solstatB(*) "solver status of solver B";

loop( (modelname(m),modeltype,solverA,nlpdef,mipdef),

   modstat = tracedata(modelname,modeltype,solverA,nlpdef,mipdef,'modelstat');
   solstat = tracedata(modelname,modeltype,solverA,nlpdef,mipdef,'solvestat');

   if (modstat and solstat,
      resA(modelname) =
          tracedata(modelname,modeltype,solverA,nlpdef,mipdef,'res');
      objA(modelname) =
          tracedata(modelname,modeltype,solverA,nlpdef,mipdef,'obj');
      objEstA(modelname) =
          tracedata(modelname,modeltype,solverA,nlpdef,mipdef,'objest');
      dir(modelname)  =
          tracedata(modelname,modeltype,solverA,nlpdef,mipdef,'direction');
      modstatA(modelname) =
          tracedata(modelname,modeltype,solverA,nlpdef,mipdef,'modelstat');
      solstatA(modelname) =
          tracedata(modelname,modeltype,solverA,nlpdef,mipdef,'solvestat');
    );
);


loop( (modelname(m),modeltype,solverB,nlpdef,mipdef),

   modstat = tracedata(modelname,modeltype,solverB,nlpdef,mipdef,'modelstat');
   solstat = tracedata(modelname,modeltype,solverB,nlpdef,mipdef,'solvestat');

   if (modstat and solstat,
      resB(modelname) =
         tracedata(modelname,modeltype,solverB,nlpdef,mipdef,'res');
      objB(modelname) =
         tracedata(modelname,modeltype,solverB,nlpdef,mipdef,'obj');
      objEstB(modelname) =
         tracedata(modelname,modeltype,solverB,nlpdef,mipdef,'objest');
      modstatB(modelname) =
         tracedata(modelname,modeltype,solverB,nlpdef,mipdef,'modelstat');
      solstatB(modelname) =
         tracedata(modelname,modeltype,solverB,nlpdef,mipdef,'solvestat');
    );
);




parameter solvernum /0/;
parameter doheader   "dummy variable to check if header already written"
          ratio      "ratio of resource times"
          gapA       "gap for solver A"
          gapB       "gap for solver B"
          inconsistency   "whether there is some inconsistency in data (negative gap, ...)"
          objdiff(*) "relative difference in obj values"
          objestdiff(*)  "relative difference in objective estimates"
;

*===========================================================================
*=== Output results to _sqr.htm file
*===========================================================================

file fhtm /"%outfile%_sqr.htm"/;
put fhtm;


*=== Output square

put / '<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">'
    / '<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">'
    / '<head>'
    / '   <link rel="stylesheet" type="text/css" href="style.css" media="screen,projection" />'
    / '   <title>Solver Square Comparison</title>'
    / '</head>'
   // '<body>'
   // '<div id="container">'
   // '<div id="main">'

put / '<h2>Solver Square Comparison: Considers all models. </h2>' /;

put / '<a name="top"></a>';

put / '<p class="block">Results created on %system.date% %system.time% '
       / 'using the PAVER Server/Performance tools available at'
       / '<a href="http://www.gamsworld.eu/performance/paver">'
       / 'http://www.gamsworld.eu/performance/paver</a></p>' /

put // '<p>Compares all solver return outcomes (for example optimal, '
      'locally optimal, infeasible, unbounded,'
    / 'fail) of one solver with all return outcomes of '
      'another solver. Interrupt'
    / 'denotes resource or iteration limit has been reached.</p>' //;

solvernum = 0;
loop(solvers,
   solvernum = solvernum+1;
   put$(solvernum eq 1) / '<p>Solver ' solvers.tl:0 ' is represented on '
                          'the left (rows) and ';
   put$(solvernum eq 2) / 'solver ' solvers.tl:0 ' on top (columns).</p>' //
);

put / '<p>See the <a href="#SquareDef">solver return definitions</a> for return codes.'
    / 'Models having trace data only in one trace file are listed in'
    / 'the "no data" column of the other.</p>';

put / '<table>'
    / '   <tr>'
    / '      <th>Input File Number</th>'
    / '      <th>Filename</th>'
    / '   </tr>'
    / '   <tr>'
    / '      <td align=left><b>Tracefile 1:</b></td>'
    / '      <td align=right><a href="%trace1%">%trace1%</a></td>'
    / '   </tr>'
    / '   <tr>'
    / '      <td align=left><b>Tracefile 2:</b></td>'
    / '      <td align=right><a href="%trace2%">%trace2%</a></td>'
    / '   </tr>'
    / '   <tr>'
    / '      <td align=left><b>Solvers used :</b></td>';

solvernum = 0;
loop(solvers,
   solvernum = solvernum+1;
   put$(solvernum eq 1)
    / '      <td align=right>' solvers.tl:0 '</td>'
    / '   </tr>';
   put$(solvernum eq 2)
    / '   <tr>'
    / '      <td></td>'
    / '      <td align=right>' solvers.tl:0 '</td>'
    / '   </tr>' /;
);

put / '   <tr>'
    / '      <td align=left><b>Modeltype(s)</b></td>'
    / '      <td align=right>';

loop(modeltype,
   put '<tt> &nbsp ' modeltype.tl:0 '</tt>'
);
put / '     </td>'
    / '   </tr>'
    / '</table>';

put / '<hr>';

scalar tmp/1/;
solvernum = 0;
put / '<table>'
    / '<tr>';
loop(solvers,
   solvernum = solvernum+1;
   if(solvernum = 1,
      put /  '   <td></td>';
      loop(result,
         put /  '   <td align=center bgcolor="#cccccc"> ';
         tmp=1;
         loop(solvername,
            if(tmp=2,
               put  solvername.tl:0 ': <br>';
            else
               tmp=tmp+1;
            );
         );
         put /  '<b>' all_results.te(result):0
                '  </b></td>';
      );
      put /  '   <td align=right><b>total ' solvers.tl:0 '</b></td>';
   );
);
put / '</tr>';

loop(resultA,
   put / '<tr>'
       / '   <td>';
   solvernum = 1;
   loop(solvername,
      if(solvernum=1,
         put solvername.tl:0 ': <br>';
      );
      solvernum=solvernum+1;
   );
   put '<b>' all_results.te(resultA):0 '</b></td>';
   loop(resultB,
     put   '   <td align=right bgcolor="#eeeeee">'
     if(resulttotals(resultA,resultB),
        put  '<b><a href="#' resultA.tl:0 resultB.tl:0   '"><tt>'
             resulttotals(resultA, resultB):0:0 '</tt></b></a>';
     else
        put  '-';
     );
     put   '</td>' /;
   );

   if ( sum(resultB, resulttotals(resultA, resultB ) ),
      put    '   <td align=right><tt><b><a href="#solverA'
             resultA.tl:0  '">'
             sum(resultB, resulttotals(resultA, resultB ) ):0:0
             '</a></b></tt></td>';
   else
      put   '   <td align=right><tt>-</tt></td>';
   )

   put / '</tr>';
);

solvernum = 0;
loop(solvers,
   solvernum = solvernum+1;
   put$(solvernum eq 2) '<tr>' / '   <td bgcolor="#cccccc"><b>total '
       solvers.tl:0 '</b></td>' ;
);

loop(resultB,
   if ( sum(resultA, resulttotals(resultA, resultB ) ),
      put / '   <td align=right bgcolor="#cccccc"><tt><b><a href="#solverB'
         resultB.tl:0  '">'
         sum(resultA, resulttotals(resultA, resultB ) ):0:0
         '</a></b></tt></td>';
   else
      put / '   <td align=right bgcolor="#cccccc"><tt>-</tt></td>';
   )
);

put / '   <td align=right bgcolor="#eeeeee"><a href="#all"><tt><b>' sum( (resultA,resultB),
          resulttotals(resultA,resultB)  ):0:0 '</b></tt></a></td>' ;
put / '</tr>'
    / '</table>' //;


put / '<hr>'

* Generate model/solver return code definitions
put / '<h3><a name="SquareDef"></a>Solver return definitions:</h3>'
    / '<table>';
if(is_discrete,
    put / '<tr>'
        / '   <th align=left><b>Outcome</b></th>'
        / '   <th align=right><b><a href="status_codes.htm">Model Status</a></b></th>'
        / '   <th align=right><b><a href="status_codes.htm">Solver Status</a></b></th>'
        / '</tr>'
        / '<tr>'
        / '   <td align=left>optimal</td><td align=right> 1 or 15 </td>'
        / '   <td align=right> 1 </td>'
        / '</tr>'
        / '<tr>'
        / '   <td align=left>locally optimal</td>'
        / '   <td align=right> 2</td>'
        / '   <td align=right> any </td>'
        / '</tr>'
        / '<tr>'
        / '   <td align=left>feasible</td>'
        / '   <td align=right> 8 or 16</td>'
        / '   <td align=right> 1 or 2 or 3 or 4 or 5 </td>'
        / '</tr>'
        / '<tr>'
        / '   <td align=left>infeasible</td>'
        / '   <td align=right> 4 or 5 or 10 or 19 </td>'
        / '   <td align=right> 1 </td>'
        / '</tr>'
        / '<tr>'
        / '   <td align=left>unbounded</td><td align=right> 3 or 18</td><td align=right> 1 </td>'
        / '</tr>'
        / '<tr>'
        / '   <td align=left>fail</td>'
        / '   <td align=right> all other</td>'
        / '   <td align=right> all other </td>'
        / '</tr>'
else
    put / '<tr>'
        / '   <td align=left><b>Outcome</b></td>'
        / '   <td align=right><b><a href="status_codes.htm">Model Status</a></b></td>'
        / '   <td align=right><b><a href="status_codes.htm">Solver Status</a></b></td>'
        / '</tr>'
        / '<tr>'
        / '   <td align=left>optimal</td><td align=right> 1 or 15</td>'
        / '   <td align=right> 1 </td>'
        / '</tr>'
        / '<tr>'
        / '   <td align=left>locally optimal</td>'
        / '   <td align=right> 2</td>'
        / '   <td align=right> any </td>'
        / '</tr>'
        / '<tr>'
        / '   <td align=left>feasible</td>'
        / '   <td align=right>7 or 16 or 17 </td>'
        / '   <td align=right> 1 or 2 or 3 or 4 or 5</td>'
        / '</tr>'
        / '<tr>'
        / '   <td align=left>infeasible</td>'
        / '   <td align=right>4 or 5  or 19</td>'
        / '   <td align=right> 1 </td>'
        / '</tr>'
        / '<tr>'
        / '   <td align=left>unbounded</td>'
        / '   <td align=right>3 or 18 </td>'
        / '   <td align=right> 1 </td>'
        / '</tr>'
        / '<tr>'
        / '   <td align=left>fail</td>'
        / '   <td align=right> all other</td>'
        / '   <td align=right> all other </td>'
        / '</tr>'
);
put / '</table>' //


*=== Output models

put / '<hr>';

put / '<h3>Solver Resource Times</h3>'
    / '<ul>'
    / '<li>Models for each solver pair outcome. Listed are the '
    / 'solver resource times <tt>TIME(.)</tt> in seconds,'
    / 'as well as the ratio <tt>RATIO(./.)</tt> of'
    / 'resource times for the two solvers if both solved optimally.'

    / '<li>Also listed are the objective values <tt>OBJ(.)</tt> using both solvers.'
    / 'The <b>better solution</b> found is listed in boldface. A solution'
    / 'is considered better, if the relative objective function '
    / 'difference is greater than ';

fhtm.nr = 2;
put / bound:0 '.';

put / 'If both solutions are less than %objmin%, we use the absolute'
    / 'difference.';
put / '<li>Solver resource time ratios for a particular model are '
      'listed only if one solver has resource greater than '
       resourcemin:0  '.';

put / '</ul>'
    / '<p><hr><p>';

fhtm.nr = 1;


*=== Output square results for solver outcome pairs
loop( (resultA,resultB),
   doheader = yes;

   loop( (m(modelname),solverA,solverB),
       if( square(modelname,solverA,resultA) and
                 square(modelname,solverB,resultB),
         if(doheader,
            put / '<a name="' resultA.tl:0 resultB.tl:0 '"></a>';
            put / '<h3>'solverA.tl:0 ': '
                   all_results.te(resultA):>25 ' -- ';
            put    solverB.tl:0 ': '
                   all_results.te(resultB):>25
                  '&nbsp&nbsp&nbsp&nbsp<a href="#top"></h3><p>Back to top</a></p>'/;

            put / '<table width="100%" border=1>'
            put / '<tr>'

            put /  '   <th align=left><b>Modelname</b></th>'
                /  '   <th align=right><b>Time (' solverA.tl:0 ')</b></th>'
                /  '   <th align=right><b>Time (' solverB.tl:0 ')</b></th>'
                /  '   <th align=center><b>Ratio ('solverA.tl:0 '/' solverB.tl:0 ') </b></th>'

            if( not (sameas(resultA,"opt") or sameas(resultA,"lopt") or sameas(resultA,"feas") ),
               put /  '   <th align=right><b>Status ('solverA.tl:0 ') </b></th>';
            elseif(sameas(resultA,"nd")),
               put /  '   <th align=right><tt>--</tt></th>';
            else
               put /  '   <th align=right><b>Obj ('solverA.tl:0 ') </b></th>';
            );

            if( not ( sameas(resultB,"opt") or sameas(resultB,"lopt") or sameas(resultB,"feas") ),
               put /  '   <th align=right><b>Status ('solverB.tl:0 ') </b></th>';
            elseif(sameas(resultB,"nd")),
               put /  '   <th align=right><tt>--</tt></th>';
            else
               put /  '   <th align=right><b>Obj ('solverB.tl:0 ') </b></th>';
            );
            put    /  '</tr>';

            doheader = no;
         );

         if( resA(modelname) > 1e-5 and resB(modelname) < 1e-5,
            ratio = inf;
         elseif( resA(modelname) < 1e-5 and resB(modelname) < 1e-5),
            ratio =1;
         else
            ratio = resA(modelname) / resB(modelname);
         );

         put / '<tr>'
             / '   <td align=left><tt>' modelname.tl:0
               '</tt></td>';
         if (not sameas(resultA,"nd"),
            put / '   <td align=right><tt>' resA(modelname):15:4
                  '</tt></td>';
         else
            put / '   <td align=right><tt>--</tt></td>';
         )
         if (not sameas(resultB,"nd"),
            put / '   <td align=right><tt>' resB(modelname):15:4
                  '</tt></td>';
         else
            put / '   <td align=right><tt>--</tt></td>';
         )

*        === we compute relative objdiff if both obj > objmin,
*        === otherwise use absolute difference
         if( max( abs(objA(modelname)), abs(objB(modelname))) > %objmin%,
            objdiff(modelname) = ( objA(modelname)-objB(modelname) )
             / max( abs(objA(modelname)), abs(objB(modelname)), %objmin%);
         else
            objdiff(modelname) = ( objA(modelname)-objB(modelname) );
         );

*        === Output ratios of resource times
         if( (sameas(resultA,"opt") or
              sameas(resultA,"lopt") or
              sameas(resultA,"feas")
             ) and
             (sameas(resultB,"opt") or
              sameas(resultB,"lopt") or
              sameas(resultB,"feas") ),

            if( resA(modelname) < resourcemin and
                resB(modelname) < resourcemin,
               put / '   <td align=right><tt>---</tt></td>';
            elseif (ratio = inf),
               put / '   <td align=right><tt>' ratio:16 '.000</tt></td>';
            else
               put / '   <td align=right><tt>' ratio:0:3 '</tt></td>';
            );

*        === Output objective values (maximizing)
            if( dir(modelname),
               if (objdiff(modelname) > bound,
                  put / '   <td align=right><b><tt>' objA(modelname):40:8
                        '</tt></b></td>'
                      / '   <td align=right><tt>'  objB(modelname):20:8
                        '</tt></td>';
               elseif (-objdiff(modelname) > bound),
                  put / '   <td align=right><tt>' objA(modelname):40:8
                        '</tt></td>'
                      / '   <td align=right><b><tt>' objB(modelname):>20:8
                        '</tt></b></td>';
               else
                  put / '   <td align=right><tt>' objA(modelname):40:8
                        '</tt></td>'
                      / '   <td align=right><tt>' objB(modelname):20:8
                        '</tt></td>' ;
               );
*           === Output objective values (minimizing)
            else
               if (-objdiff(modelname) > bound,
                  put / '   <td align=right><b><tt>' objA(modelname):40:8
                        '</tt></b></td>'
                      / '   <td align=right><tt>'  objB(modelname):20:8
                        '</tt></td>'/;
               elseif (objdiff(modelname) > bound),
                  put / '   <td align=right><tt>' objA(modelname):40:8
                        '</tt></td>'
                      / '   <td align=right><b><tt>'  objB(modelname):20:8
                        '</tt></b></td>';
               else
                  put / '   <td align=right><tt>' objA(modelname):40:8
                        '</tt></td>'
                      / '   <td align=right><tt>' objB(modelname):20:8
                        '</tt></td>';
               );
            );

         elseif( (sameas(resultA,"opt") or
                  sameas(resultA,"lopt") or
                  sameas(resultA,"feas") ) ),
            put / '   <td align=right>---</td>'
                / '   <td align=right><tt>' objA(modelname):0:8 '</tt></td>';
            if( not sameas(resultB,"nd"),
               put / '   <td align=right><tt>'
                           'mstat(' modstatB(modelname):2:0 ')'
                           'sstat(' solstatB(modelname):2:0  ')</tt>'
                      '</td>';
            else
               put / '   <td align=right><tt> ---</tt></td>';
            );
         elseif( (sameas(resultB,"opt") or
                  sameas(resultB,"lopt") or
                  sameas(resultB,"feas") ) ),
            put / '   <td align=right>---</td>';
            if( not sameas(resultA,"nd"),
               put / '   <td align=right><tt>'
                           'mstat(' modstatA(modelname):2:0 ')'
                           'sstat(' solstatA(modelname):2:0  ')</tt>'
                      '</td>';
            else
               put / '   <td align=right><tt> ---</tt></td>';
            );
            put / '   <td align=right><tt>' objB(modelname):0:8 '</tt></td>';
         else
            put / '   <td align=right>---</td>';
                if(not sameas(resultA,"nd"),
                   put / '   <td align=right><tt>'
                              'mstat(' modstatA(modelname):2:0 ')'
                              'sstat(' solstatA(modelname):2:0  ')</tt>'
                         '</td>';
                else
                   put / '   <td align=right><tt> --- </tt></td>';
                )
                if(not sameas(resultB,"nd"),
                   put / '   <td align=right><tt>'
                              'mstat(' modstatB(modelname):2:0 ')'
                              'sstat(' solstatB(modelname):2:0  ')</tt>'
                         '</td>';
                else
                   put / '   <td align=right><tt> --- </tt></td>';
                )
         );
         put / '</tr>';
      );
   );

   if(not doheader,
      put / '</table>'
      put / '<p><a href="#top">Back to top</a></p>'
      put / '<hr>' ///;
   );

);


*=== Output result totals for solver A outcome
loop( resultA,
   doheader = yes;

   loop( (m(modelname),solverA),
      if ( square(modelname,solverA,resultA),
         if(doheader,
            put / '<a name="solverA' resultA.tl:0  '"></a>' ;
            put / '<h3>'solverA.tl:0 ': ' all_results.te(resultA):0 ':</h3>'
                  '<p><a href="#top">Back to top</a></p>'/;

            put / '<table>';
            put / '<tr>'

            put /  '   <th align=left><b>Modelname</b></th>'
                /  '   <th align=right><b>Time (' solverA.tl:0 ')</b></th>'
            if( not (sameas(resultA,"opt") or
                     sameas(resultA,"lopt") or
                     sameas(resultA,"feas")),
               put /  '<th align=right><b>Status ('solverA.tl:0 ') </b></th>';
            else
               put /  '<th align=right><b>Obj ('solverA.tl:0 ') </b></th>';
            );
            put / '</tr>';

            doheader = no;
         );

         put / '<tr>'
             / '   <td align=left><tt>' modelname.tl:0 '</tt></td>';
         if(sameas(resultA,"nd") or
            sameas(resultA,"infeas") or
            sameas(resultA,"f"),
               put /  '   <td align=right><tt> -- </tt></td>';
         else
            put / '   <td align=right><tt>' resA(modelname):15:4
                  '</tt></td>';
         );


         if( sameas(resultA,"opt") or
             sameas(resultA,"lopt") or
             sameas(resultA,"feas"),
            put / '   <td align=right><tt>' objA(modelname):0:8
                      '</tt></td>';
         elseif( sameas(resultA,"nd")),
            put / '   <td align=right><tt> -- </tt></td>';
         else
            put / '   <td align=right><tt>mstat(' modstatA(modelname):2:0
                  ') sstat(' solstatA(modelname):2:0  ')</tt></td>';
         );
         put / '</tr>';
      );
   );

   if(not doheader,
      put / '</tr></table>'
      put / '<p align=right><a href="#top">Back to top</a></p>'
      put / '<hr>' ///;
   );

);



*=== Output result totals for solver B outcome
display square, resB, objB;

loop( resultB,
   doheader = yes;

   loop( (m(modelname),solverB),
      if ( square(modelname,solverB,resultB),
         if(doheader,
            put / '<a name="solverB' resultB.tl:0  '"></a>' ;
            put / '<h3>'solverB.tl:0 ': '
                  all_results.te(resultB):0 ':</h3>'
                  '<p><a href="#top">Back to top</a></p>'/;

            put / '<table>';
            put / '<tr>'

            put /  '   <th align=left><b>Modelname</b></th>'
                /  '   <th align=right><b>Time (' solverB.tl:0 ')</b></th>'
            if( not (sameas(resultB,"opt") or
                     sameas(resultB,"lopt") or
                     sameas(resultB,"feas") ),
               put /  '   <th align=right><b>Status ('solverB.tl:0 ') </b></th>';
            else
               put /  '   <th align=right><b>Obj ('solverB.tl:0 ') </b></th>';
            );
            put    /  '</tr>';

            doheader = no;
         );

         put / '<tr>'
             / '   <td align=left><tt>' modelname.tl:0 '</tt></td>';
         if(sameas(resultB,"nd") or
            sameas(resultB,"infeas") or
            sameas(resultB,"f"),
               put /  '   <td align=right><tt> -- </tt></td>';
         else
            put / '   <td align=right><tt>' resB(modelname):15:4
                  '</tt></td>';
         );


         if( sameas(resultB,"opt") or
             sameas(resultB,"lopt") or
             sameas(resultB,"feas"),
            put / '   <td align=right><tt>' objB(modelname):0:8 '</tt></td>';
         elseif( sameas(resultB,"nd")),
            put / '   <td align=right><tt> -- </tt></td>';
         else
            put / '   <td align=right><tt>'
                        'mstat(' modstatB(modelname):2:0 ')'
                        'sstat(' solstatB(modelname):2:0  ')</tt>'
                  '</td>';
         );
         put / '</tr>';
      );
   );

   if(not doheader,
      put / '</tr></table>'
      put / '<p align=right><a href="#top">Back to top</a></p>';
   );

);


* === Output results of all solves with additional info (shifted geometric mean
* === and % unsolved)
Parameter resultcount / 0 /;

loop( (m(modelname),solverA,solverB),
   loop( (resultA,resultB)$(square(modelname,solverA,resultA) and
                            square(modelname,solverB,resultB)
       ),

       resultcount = resultcount+1;
       if( not (sameas(resultA,"opt") or
                sameas(resultA,"lopt") or
                sameas(resultA,"feas") ),
          objA(modelname) = Inf;
          objestA(modelname) = Inf;
       );
       if( not (sameas(resultB,"opt") or
                sameas(resultB,"lopt") or
                sameas(resultB,"feas") ),
          objB(modelname) = Inf;
          objestB(modelname) = Inf;
       );

       if (objestA(modelname) = NA, objestA(modelname) = Inf);
       if (objestB(modelname) = NA, objestB(modelname) = Inf);

*   === we compute relative objdiff if both obj > objmin,
*   === otherwise use absolute difference
       if( max( abs(objA(modelname)), abs(objB(modelname))) > %objmin%,
          objdiff(modelname) = ( objA(modelname)-objB(modelname) ) /
              max( abs(objA(modelname)), abs(objB(modelname)), %objmin%);
       else
          objdiff(modelname) = ( objA(modelname)-objB(modelname) );
       );

*   === we compute relative objestdiff if both objest > objmin,
*   === otherwise use absolute difference
       if (objestA(modelname) < Inf and objestB(modelname) < Inf,
          if( max( abs(objestA(modelname)), abs(objestB(modelname))) > %objmin%,
             objestdiff(modelname) =
               ( objestA(modelname)-objestB(modelname) ) /
               max( abs(objestA(modelname)), abs(objestB(modelname)), %objmin%);
          else
             objestdiff(modelname) = ( objestA(modelname)-objestB(modelname) );
          );
       else
          objestdiff(modelname) = 0;
       );

*   === we compute relative gap if obj and objest > objmin,
*   === otherwise use absolute difference
       if( abs(objestA(modelname)) >= Inf,
          gapA(modelname) = Inf;
          objestA(modelname) = Inf;
       else
          if( abs(objestA(modelname)) > %objmin%,
             gapA(modelname) = (objA(modelname) - objestA(modelname)) /
                 abs(objestA(modelname));
          else
             gapA(modelname) = objA(modelname) - objestA(modelname);
          );
*     === negate if maximizing
          if( dir(modelname), gapA(modelname) = -gapA(modelname) );
       );

       if( abs(objestB(modelname)) >= Inf,
          gapB(modelname) = Inf;
          objestB(modelname) = Inf;
       else
          if( abs(objestB(modelname)) > %objmin%,
             gapB(modelname) = (objB(modelname) - objestB(modelname)) /
                     abs(objestB(modelname));
          else
             gapB(modelname) = objB(modelname) - objestB(modelname);
          );
*      === negate if maximizing
          if( dir(modelname), gapB(modelname) = -gapB(modelname) );
       );
   );
);


Scalar
   resAavg      average of resource times over all models (as shifted geom. mean)
   resBavg      average of resource times over all models (as shifted geom. mean)
   solvedA      number of models solved by solver A
   solvedB      number of models solved by solver B
   gapAbetter   whether gap of solver A is better than gap of solver B
   gapBbetter   whether gap of solver B is better than gap of solver A
;
Scalar timeshift / 10 /;

if (have_objest,
*=== compute shifted geometric mean of resA or reslim; take reslim if smaller or not solved or solved above gap tolerance
 resAavg =
   prod( (m(modelname),solverA,solverB),
    (  sum( (resultA,resultB)$
            (square(modelname,solverA,resultA) and
             square(modelname,solverB,resultB) and
             not mapval(resA(modelname)) and gapA(modelname) <= %gaptol%
            ), min(resA(modelname),reslim))
     + sum( (resultA,resultB)$
            (square(modelname,solverA,resultA) and
             square(modelname,solverB,resultB) and
             (   mapval(resA(modelname))  or gapA(modelname) >  %gaptol%)
            ), reslim)
     + timeshift
    )**(1/resultcount)
   ) - timeshift;

*=== compute shifted geometric mean of resA or reslim; take reslim if smaller or not solved or solved above gap tolerance
 resBavg =
   prod( (m(modelname),solverA,solverB),
    (  sum( (resultA,resultB)$
            (square(modelname,solverA,resultA) and
             square(modelname,solverB,resultB) and
             not mapval(resB(modelname)) and gapB(modelname) <= %gaptol%
            ), min(resB(modelname),reslim))
     + sum( (resultA,resultB)$
            (square(modelname,solverA,resultA) and
             square(modelname,solverB,resultB) and
             (   mapval(resB(modelname))  or gapB(modelname) >  %gaptol%)
            ), reslim)
     + timeshift
    )**(1/resultcount)
   ) - timeshift;

*=== compute number of models solved by solver A (= gap below tolerance)
 solvedA =
  sum( (m(modelname),solverA,solverB),
  sum( (resultA,resultB)$
       (square(modelname,solverA,resultA) and
        square(modelname,solverB,resultB) and
        not mapval(resB(modelname)) and gapB(modelname) <= %gaptol%
       ), 1));

*=== compute number of models solved by solver B (= gap below tolerance)
 solvedB =
  sum( (m(modelname),solverA,solverB),
  sum( (resultA,resultB)$
       (square(modelname,solverA,resultA) and
        square(modelname,solverB,resultB) and
        not mapval(resB(modelname)) and gapB(modelname) <= %gaptol%
       ), 1));
);

loop((solverA,solverB),
            put / '<a name="all"></a>'
                / '<h3>' solverA.tl:0 ': all -- ' solverB.tl:0 ': all</h3>'
                / '<a href="#top">Back to top</a><br/><br/>';
);

put / '<p>A "<tt>(max)</tt>" behind a model names indicates that this model is a maximization model.</p>';

if (have_objest,
   put / '<p>Rows with <font class="redproblem">red background</font> in the model name '
       / 'mark <b>inconsistencies</b>, caused by negative gaps or objective values of one solver '
       / 'that are better than the objestive estimate of the other solver.</p>'
       / '<p>A model is counted as <i>solved</i> if the gap reported by the solver is at '
       / 'most %gaptol%. The gap is in this case displayed with a <font class="greenopt">green background</font>.</p>'
       / '<p>A shifted geometric mean of the solver time is computed, where unsolved '
       /  'instances and instances above the timelimit are counted with the timelimit ('
          reslim:0 ' seconds). The shift is 10s.</p>';
);

loop((solverA,solverB),
   put / '<table>'
       / '<tr>'
       /  '   <td align=left><b>Modelname</b></td>'
       /  '   <td align=right><b>Time (' solverA.tl:0 ')</b></td>'
       /  '   <td align=right><b>Time (' solverB.tl:0 ')</b></td>'
       /  '   <td align=center><b>Ratio ('solverA.tl:0 '/ ' solverB.tl:0 ') </b></td>'
       /  '   <td align=right><b>Obj ('solverA.tl:0 ') </b></td>'
       /  '   <td align=right><b>Obj ('solverB.tl:0 ') </b></td>';
 if (have_objest,
   put /  '   <td align=right><b>ObjEst ('solverA.tl:0 ') </b></td>'
       /  '   <td align=right><b>ObjEst ('solverB.tl:0 ') </b></td>'
       /  '   <td align=right><b>Gap ('solverA.tl:0') </b></td>'
       /  '   <td align=right><b>Gap ('solverB.tl:0') </b></td>'
    );
   put /  '</tr>'
);

Scalar logexists     / 0 /;
Scalar logbz2exists  / 0 /;
Scalar loggzexists   / 0 /;
Scalar logzipexists  / 0 /;

file fx / '%gams.scrdir%file_exists.cmd' /;
$if     %system.filesys% == UNIX $setglobal prefix 'source '
$if not %system.filesys% == UNIX $setglobal prefix ' '

loop( (m(modelname),solverA,solverB),
  loop( (resultA,resultB)$(square(modelname,solverA,resultA) and square(modelname,solverB,resultB)),

*        === check that gaps are positive and that not one solver reports a solution better than the others bound
         inconsistency = no;
         if (have_objest,
           inconsistency$(gapA(modelname) < -bound or gapB(modelname) < -bound) = yes;
           if( (sameas(resultA,"opt") or sameas(resultA,"lopt") or sameas(resultA,"feas") ) and
               (sameas(resultB,"opt") or sameas(resultB,"lopt") or sameas(resultB,"feas") ),

             if (dir(modelname),
               inconsistency$(objA(modelname) > objestB(modelname) + bound * max(abs(objestB(modelname)), %objmin%)) = yes;
               inconsistency$(objB(modelname) > objestA(modelname) + bound * max(abs(objestA(modelname)), %objmin%)) = yes;
             else
               inconsistency$(objestB(modelname) < Inf and objA(modelname) < objestB(modelname) - bound * max(abs(objestB(modelname)), %objmin%)) = yes;
               inconsistency$(objestA(modelname) < Inf and objB(modelname) < objestA(modelname) - bound * max(abs(objestA(modelname)), %objmin%)) = yes;
             );
           );
         );

         if( resA(modelname) > 1e-5 and resB(modelname) < 1e-5,
            ratio = inf;
         elseif( resA(modelname) < 1e-5 and resB(modelname) < 1e-5),
            ratio = 1;
         else
            ratio = resA(modelname) / resB(modelname);
         );
        put / '<tr>';
         if (inconsistency,
           put '   <td align=left class="redproblem">';
         else
           put '   <td align=left>';
         );
         put '<tt>';
$if set liburl
         put '<a href="%liburl%/' modelname.tl:0 '.htm">' modelname.tl:0 '</a>';
$if not set liburl
         put modelname.tl:0;

         if (dir(modelname), put ' (max)'; );
         put '</tt></td>';
$ifthen not %uselog% == 0
* -- check whether logfile for solverA exists
           putclose fx 'test -f logs/%trace1%/' modelname.tl:0 '.log';
           execute '%prefix% %gams.scrdir%file_exists.cmd';
           logexists = not errorlevel;
         if (not logexists,
           putclose fx 'test -f logs/%trace1%/' modelname.tl:0 '.log.bz2';
           execute '%prefix% %gams.scrdir%file_exists.cmd';
           logbz2exists = not errorlevel;
         if (not logbz2exists,
           putclose fx 'test -f logs/%trace1%/' modelname.tl:0 '.log.gz';
           execute '%prefix% %gams.scrdir%file_exists.cmd';
           loggzexists = not errorlevel;
         if (not logzipexists,
           putclose fx 'test -f logs/%trace1%/' modelname.tl:0 '.log.zip';
           execute '%prefix% %gams.scrdir%file_exists.cmd';
           logzipexists = not errorlevel;
         )));
         put fhtm;
$endif

         put / '   <td align=right><tt>'
         if (not sameas(resultA,"nd"),
            if (logexists, put / '<a href="logs/%trace1%/' modelname.tl:0 '.log">';
            else if (loggzexists, put '<a href="logs/%trace1%/' modelname.tl:0 '.log.gz">';
            else if (logbz2exists, put '<a href="logs/%trace1%/' modelname.tl:0 '.log.bz2">';
            else if (logzipexists, put '<a href="logs/%trace1%/' modelname.tl:0 '.log.zip">'; ))))
            put resA(modelname):15:4 ;
            if (logexists or loggzexists or logbz2exists or logzipexists, put '</a>');
         else
            put 'no data';
         )
         put '</tt></td>';
$ifthen not %uselog% == 0
*        -- check whether logfile for solverB exists
           putclose fx 'test -f logs/%trace2%/' modelname.tl:0 '.log';
           execute '%prefix% %gams.scrdir%file_exists.cmd';
           logexists = not errorlevel;
         if (not logexists,
           putclose fx 'test -f logs/%trace2%/' modelname.tl:0 '.log.bz2';
           execute '%prefix% %gams.scrdir%file_exists.cmd';
           logbz2exists = not errorlevel;
         if (not logbz2exists,
           putclose fx 'test -f logs/%trace2%/' modelname.tl:0 '.log.gz';
           execute '%prefix% %gams.scrdir%file_exists.cmd';
           loggzexists = not errorlevel;
         if (not logzipexists,
           putclose fx 'test -f logs/%trace2%/' modelname.tl:0 '.log.zip';
           execute '%prefix% %gams.scrdir%file_exists.cmd';
           logzipexists = not errorlevel;
         )));
         put fhtm;
$endif

         put / '   <td align=right><tt>'
         if (not sameas(resultB,"nd"),
            if (logexists, put / '<a href="logs/%trace2%/' modelname.tl:0 '.log">';
            else if (loggzexists, put '<a href="logs/%trace2%/' modelname.tl:0 '.log.gz">';
            else if (logbz2exists, put '<a href="logs/%trace2%/' modelname.tl:0 '.log.bz2">';
            else if (logzipexists, put '<a href="logs/%trace2%/' modelname.tl:0 '.log.zip">'; ))))
            put resB(modelname):15:4 ;
            if (logexists or loggzexists or logbz2exists or logzipexists, put '</a>');
         else
            put 'no data';
         )
         put '</tt></td>';

         if( (sameas(resultA,"opt") or sameas(resultA,"lopt") or sameas(resultA,"feas") ) and
             (sameas(resultB,"opt") or sameas(resultB,"lopt") or sameas(resultB,"feas") ),

*           === Output ratios of resource times
            if( resA(modelname) < resourcemin and resB(modelname) < resourcemin,
               put / '   <td align=right><tt>---</tt></td>';
            elseif (ratio = inf),
               put / '   <td align=right><tt>' ratio:16 '.000</tt></td>';
            else
               put / '   <td align=right><tt>' ratio:0:3 '</tt></td>';
            );
*           === Output objective values (maximizing) and objective estimates (if have_objest)
            if( dir(modelname),
               if (objdiff(modelname) > bound,
                  put / '   <td align=right><b><tt>' objA(modelname):40:8 '</tt></b></td>'
                      / '   <td align=right><tt>'  objB(modelname):20:8 '</tt></td>';
               elseif (-objdiff(modelname) > bound),
                  put / '   <td align=right><tt>' objA(modelname):40:8 '</tt></td>'
                      / '   <td align=right><b><tt>' objB(modelname):>20:8 '</tt></b></td>';
               else
                  put / '   <td align=right><tt>' objA(modelname):40:8 '</tt></td>'
                      / '   <td align=right><tt>' objB(modelname):20:8 '</tt></td>' ;
               );
               if (have_objest,
                 if (objestdiff(modelname) > bound,
                    put / '   <td align=right><tt>' objestA(modelname):40:8 '</tt></td>'
                        / '   <td align=right><b><tt>'  objestB(modelname):20:8 '</tt></b></td>';
                 elseif (-objestdiff(modelname) > bound),
                    put / '   <td align=right><b><tt>' objestA(modelname):40:8 '</tt></b></td>'
                        / '   <td align=right><tt>' objestB(modelname):>20:8 '</tt></td>';
                else
                    put / '   <td align=right><tt>' objestA(modelname):40:8 '</tt></td>'
                        / '   <td align=right><tt>' objestB(modelname):20:8 '</tt></td>' ;
                 );
               )
*           === Output objective values (minimizing) and objective estimates (if have_objest)
            else
               if (-objdiff(modelname) > bound,
                  put / '   <td align=right><b><tt>' objA(modelname):40:8 '</tt></b></td>'
                      / '   <td align=right><tt>'  objB(modelname):20:8 '</tt></td>'/;
               elseif (objdiff(modelname) > bound),
                  put / '   <td align=right><tt>' objA(modelname):40:8 '</tt></td>'
                      / '   <td align=right><b><tt>'  objB(modelname):20:8 '</tt></b></td>';
               else
                  put / '   <td align=right><tt>' objA(modelname):40:8 '</tt></td>'
                      / '   <td align=right><tt>' objB(modelname):20:8 '</tt></td>';
               );
               if (have_objest,
                 if (-objestdiff(modelname) > bound,
                   put / '   <td align=right><tt>' objestA(modelname):40:8 '</tt></td>'
                       / '   <td align=right><b><tt>'  objestB(modelname):20:8 '</tt></b></td>';
                 elseif (objestdiff(modelname) > bound),
                   put / '   <td align=right><b><tt>' objestA(modelname):40:8 '</tt></b></td>'
                       / '   <td align=right><tt>' objestB(modelname):>20:8 '</tt></td>';
                 else
                   put / '   <td align=right><tt>' objestA(modelname):40:8 '</tt></td>'
                       / '   <td align=right><tt>' objestB(modelname):20:8 '</tt></td>' ;
                 );
               )
            );
*           === Output gaps
            if (have_objest,
              gapAbetter = gapB(modelname) - gapA(modelname) > bound and gapB(modelname) < Inf and gapA(modelname) < Inf;
              gapBbetter = gapA(modelname) - gapB(modelname) > bound and gapA(modelname) < Inf and gapB(modelname) < Inf;

              put / '   <td align=right ';
              if (gapA(modelname) <= %gaptol%, put 'class="greenopt"');
              put '>';
              if (gapAbetter, put '<b>');
              if (gapA(modelname) < Inf, put '<tt>' gapA(modelname):40:4 '</tt>', else put '<tt>NA<tt>' );
              if (gapAbetter, put '</b>');
              put '</td>';

              put / '   <td align=right ';
              if (gapB(modelname) <= %gaptol%, put 'class="greenopt"');
              put '>';
              if (gapBbetter, put '<b>');
              if (gapB(modelname) < Inf, put '<tt>' gapB(modelname):40:4 '</tt>', else put '<tt>NA<tt>' );
              if (gapBbetter, put '</b>');
              put '</td>';
            )
         elseif( (sameas(resultA,"opt") or sameas(resultA,"lopt") or sameas(resultA,"feas") ) ),
            put / '   <td align=right>---</td>'
                / '   <td align=right><b><tt>' objA(modelname):0:8 '</tt></b></td>';
            if( not sameas(resultB,"nd"),
               put / '   <td align=right><tt> mstat(' modstatB(modelname):2:0 ') sstat(' solstatB(modelname):0:0  ')</tt></td>';
            else
               put / '   <td align=right><tt> no data</tt></td>';
            );
            if (have_objest,
              put / '   <td align=right><tt>' objEstA(modelname):40:8 '</tt></td>'
                  / '   <td align=right>---</td>'
                  / '   <td align=right ';
              if (gapA(modelname) <= %gaptol%, put 'class="greenopt"');
              put   '><tt>' gapA(modelname):40:4 '</tt></td>'
                  / '   <td align=right><tt>---</tt></td>';
            )
         elseif( (sameas(resultB,"opt") or sameas(resultB,"lopt") or sameas(resultB,"feas")  )),
            put / '   <td align=right>---</td>';
            if( not sameas(resultA,"nd"),
               put / '   <td align=right><tt> mstat(' modstatA(modelname):2:0 ') sstat(' solstatA(modelname):0:0  ')</tt></td>';
            else
               put / '   <td align=right><tt> no data</tt></td>';
            );
            put / '   <td align=right><b><tt>' objB(modelname):0:8 '</tt></b></td>';
            if (have_objest,
              put / '   <td align=right>---</td>'
                  / '   <td align=right><tt>' objEstB(modelname):40:8 '</tt></td>'
                  / '   <td align=right><tt>---</tt></td>'
                  / '   <td align=right '
              if (gapB(modelname) <= %gaptol%, put 'class="greenopt"');
              put '><tt>' gapB(modelname):40:4 '</tt></td>'
            );
         else
            put / '   <td align=right>---</td>';
                if(not sameas(resultA,"nd"),
                   put / '   <td align=right><tt> mstat(' modstatA(modelname):2:0 ') sstat(' solstatA(modelname):0:0  ')</tt></td>';
                else
                   put / '   <td align=right><tt> no data </tt></td>';
                )
                if(not sameas(resultB,"nd"),
                   put / '   <td align=right><tt> mstat(' modstatB(modelname):2:0 ') sstat(' solstatB(modelname):0:0  ')</tt></td>';
                else
                   put / '   <td align=right><tt> no data </tt></td>';
                )
            if (have_objest,
              put / '   <td align=right>---</td>'
                  / '   <td align=right>---</td>'
                  / '   <td align=right>---</td>'
                  / '   <td align=right>---</td>';
            )
         );
         put / '</tr>';
      );
);

*=== Output mean values and percentage of unsolved instances
if (have_objest,
put / '</tr>'
    / '   <td><b>Shifted Geometric Mean<br/>(Solved)</b></td>'
    / '   <td align=right>' resAavg '</td>'
    / '   <td align=right>' resBavg '</td>'
    / '   <td align=right>&nbsp;</td>'
    / '   <td align=right>&nbsp;</td>'
    / '   <td align=right>&nbsp;</td>'
    / '   <td align=right>&nbsp;</td>'
    / '   <td align=right>&nbsp;</td>'
    / '   <td align=right>&nbsp;</td>'
    / '   <td align=right>&nbsp;</td>';
put / '</tr><tr>'
    / '   <td><b>% Unsolved</b></td>'
    / '   <td align=right>' (100*(resultcount-solvedA)/resultcount) '%</td>'
    / '   <td align=right>' (100*(resultcount-solvedB)/resultcount) '%</td>'
    / '<td>&nbsp;</td><td>&nbsp;</td><td>&nbsp;</td>'
    / '<td>&nbsp;</td><td>&nbsp;</td><td>&nbsp;</td><td>&nbsp;</td>';
);

put / '</tr>'
    / '</table>'
    / '<p align=right><a href="#top">Back to top</a></p>'
    / '<hr>';


put / '</div> <!-- end div=main -->'
    / '<div id="footer">'
    / '<p>Generated using the PAVER Server/Performance Tools.</a></p>'
   // '</div> '
   // '</div> <!-- end div=container -->'
   // '</body>'
   // '</html>'

putclose fhtm;

$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


*=== Warning message if more than two solvers defined in trace file
file fput / "%gams.scrdir%solvers.txt" /;
put fput "--- ";

if (card(solvers) > 2,
  loop(solvers,
     put fput solvers.tl:0 "  ";
  );
  putclose fput;
  execute "echo ---"
  execute "echo --- WARNING: more than two solvers defined in trace files."
  execute "echo ---          Results may not be what you intend."
  execute "echo ---          Solvers defined are:"
  execute "head %gams.scrdir%solvers.txt"
  execute "echo ---"
);

putclose fput;