/* glpscl.c (problem scaling routines) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*  Copyright (C) 2000-2013 Free Software Foundation, Inc.
*  Written by Andrew Makhorin <mao@gnu.org>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "env.h"
#include "misc.h"
#include "prob.h"

/***********************************************************************
*  min_row_aij - determine minimal |a[i,j]| in i-th row
*
*  This routine returns minimal magnitude of (non-zero) constraint
*  coefficients in i-th row of the constraint matrix.
*
*  If the parameter scaled is zero, the original constraint matrix A is
*  assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
*
*  If i-th row of the matrix is empty, the routine returns 1. */

static double min_row_aij(glp_prob *lp, int i, int scaled)
{     GLPAIJ *aij;
      double min_aij, temp;
      xassert(1 <= i && i <= lp->m);
      min_aij = 1.0;
      for (aij = lp->row[i]->ptr; aij != NULL; aij = aij->r_next)
      {  temp = fabs(aij->val);
         if (scaled) temp *= (aij->row->rii * aij->col->sjj);
         if (aij->r_prev == NULL || min_aij > temp)
            min_aij = temp;
      }
      return min_aij;
}

/***********************************************************************
*  max_row_aij - determine maximal |a[i,j]| in i-th row
*
*  This routine returns maximal magnitude of (non-zero) constraint
*  coefficients in i-th row of the constraint matrix.
*
*  If the parameter scaled is zero, the original constraint matrix A is
*  assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
*
*  If i-th row of the matrix is empty, the routine returns 1. */

static double max_row_aij(glp_prob *lp, int i, int scaled)
{     GLPAIJ *aij;
      double max_aij, temp;
      xassert(1 <= i && i <= lp->m);
      max_aij = 1.0;
      for (aij = lp->row[i]->ptr; aij != NULL; aij = aij->r_next)
      {  temp = fabs(aij->val);
         if (scaled) temp *= (aij->row->rii * aij->col->sjj);
         if (aij->r_prev == NULL || max_aij < temp)
            max_aij = temp;
      }
      return max_aij;
}

/***********************************************************************
*  min_col_aij - determine minimal |a[i,j]| in j-th column
*
*  This routine returns minimal magnitude of (non-zero) constraint
*  coefficients in j-th column of the constraint matrix.
*
*  If the parameter scaled is zero, the original constraint matrix A is
*  assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
*
*  If j-th column of the matrix is empty, the routine returns 1. */

static double min_col_aij(glp_prob *lp, int j, int scaled)
{     GLPAIJ *aij;
      double min_aij, temp;
      xassert(1 <= j && j <= lp->n);
      min_aij = 1.0;
      for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next)
      {  temp = fabs(aij->val);
         if (scaled) temp *= (aij->row->rii * aij->col->sjj);
         if (aij->c_prev == NULL || min_aij > temp)
            min_aij = temp;
      }
      return min_aij;
}

/***********************************************************************
*  max_col_aij - determine maximal |a[i,j]| in j-th column
*
*  This routine returns maximal magnitude of (non-zero) constraint
*  coefficients in j-th column of the constraint matrix.
*
*  If the parameter scaled is zero, the original constraint matrix A is
*  assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
*
*  If j-th column of the matrix is empty, the routine returns 1. */

static double max_col_aij(glp_prob *lp, int j, int scaled)
{     GLPAIJ *aij;
      double max_aij, temp;
      xassert(1 <= j && j <= lp->n);
      max_aij = 1.0;
      for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next)
      {  temp = fabs(aij->val);
         if (scaled) temp *= (aij->row->rii * aij->col->sjj);
         if (aij->c_prev == NULL || max_aij < temp)
            max_aij = temp;
      }
      return max_aij;
}

/***********************************************************************
*  min_mat_aij - determine minimal |a[i,j]| in constraint matrix
*
*  This routine returns minimal magnitude of (non-zero) constraint
*  coefficients in the constraint matrix.
*
*  If the parameter scaled is zero, the original constraint matrix A is
*  assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
*
*  If the matrix is empty, the routine returns 1. */

static double min_mat_aij(glp_prob *lp, int scaled)
{     int i;
      double min_aij, temp;
      min_aij = 1.0;
      for (i = 1; i <= lp->m; i++)
      {  temp = min_row_aij(lp, i, scaled);
         if (i == 1 || min_aij > temp)
            min_aij = temp;
      }
      return min_aij;
}

/***********************************************************************
*  max_mat_aij - determine maximal |a[i,j]| in constraint matrix
*
*  This routine returns maximal magnitude of (non-zero) constraint
*  coefficients in the constraint matrix.
*
*  If the parameter scaled is zero, the original constraint matrix A is
*  assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
*
*  If the matrix is empty, the routine returns 1. */

static double max_mat_aij(glp_prob *lp, int scaled)
{     int i;
      double max_aij, temp;
      max_aij = 1.0;
      for (i = 1; i <= lp->m; i++)
      {  temp = max_row_aij(lp, i, scaled);
         if (i == 1 || max_aij < temp)
            max_aij = temp;
      }
      return max_aij;
}

/***********************************************************************
*  eq_scaling - perform equilibration scaling
*
*  This routine performs equilibration scaling of rows and columns of
*  the constraint matrix.
*
*  If the parameter flag is zero, the routine scales rows at first and
*  then columns. Otherwise, the routine scales columns and then rows.
*
*  Rows are scaled as follows:
*
*                         n
*     a'[i,j] = a[i,j] / max |a[i,j]|,  i = 1,...,m.
*                        j=1
*
*  This makes the infinity (maximum) norm of each row of the matrix
*  equal to 1.
*
*  Columns are scaled as follows:
*
*                         m
*     a'[i,j] = a[i,j] / max |a[i,j]|,  j = 1,...,n.
*                        i=1
*
*  This makes the infinity (maximum) norm of each column of the matrix
*  equal to 1. */

static void eq_scaling(glp_prob *lp, int flag)
{     int i, j, pass;
      double temp;
      xassert(flag == 0 || flag == 1);
      for (pass = 0; pass <= 1; pass++)
      {  if (pass == flag)
         {  /* scale rows */
            for (i = 1; i <= lp->m; i++)
            {  temp = max_row_aij(lp, i, 1);
               glp_set_rii(lp, i, glp_get_rii(lp, i) / temp);
            }
         }
         else
         {  /* scale columns */
            for (j = 1; j <= lp->n; j++)
            {  temp = max_col_aij(lp, j, 1);
               glp_set_sjj(lp, j, glp_get_sjj(lp, j) / temp);
            }
         }
      }
      return;
}

/***********************************************************************
*  gm_scaling - perform geometric mean scaling
*
*  This routine performs geometric mean scaling of rows and columns of
*  the constraint matrix.
*
*  If the parameter flag is zero, the routine scales rows at first and
*  then columns. Otherwise, the routine scales columns and then rows.
*
*  Rows are scaled as follows:
*
*     a'[i,j] = a[i,j] / sqrt(alfa[i] * beta[i]),  i = 1,...,m,
*
*  where:
*                n                        n
*     alfa[i] = min |a[i,j]|,  beta[i] = max |a[i,j]|.
*               j=1                      j=1
*
*  This allows decreasing the ratio beta[i] / alfa[i] for each row of
*  the matrix.
*
*  Columns are scaled as follows:
*
*     a'[i,j] = a[i,j] / sqrt(alfa[j] * beta[j]),  j = 1,...,n,
*
*  where:
*                m                        m
*     alfa[j] = min |a[i,j]|,  beta[j] = max |a[i,j]|.
*               i=1                      i=1
*
*  This allows decreasing the ratio beta[j] / alfa[j] for each column
*  of the matrix. */

static void gm_scaling(glp_prob *lp, int flag)
{     int i, j, pass;
      double temp;
      xassert(flag == 0 || flag == 1);
      for (pass = 0; pass <= 1; pass++)
      {  if (pass == flag)
         {  /* scale rows */
            for (i = 1; i <= lp->m; i++)
            {  temp = min_row_aij(lp, i, 1) * max_row_aij(lp, i, 1);
               glp_set_rii(lp, i, glp_get_rii(lp, i) / sqrt(temp));
            }
         }
         else
         {  /* scale columns */
            for (j = 1; j <= lp->n; j++)
            {  temp = min_col_aij(lp, j, 1) * max_col_aij(lp, j, 1);
               glp_set_sjj(lp, j, glp_get_sjj(lp, j) / sqrt(temp));
            }
         }
      }
      return;
}

/***********************************************************************
*  max_row_ratio - determine worst scaling "quality" for rows
*
*  This routine returns the worst scaling "quality" for rows of the
*  currently scaled constraint matrix:
*
*              m
*     ratio = max ratio[i],
*             i=1
*  where:
*                 n              n
*     ratio[i] = max |a[i,j]| / min |a[i,j]|,  1 <= i <= m,
*                j=1            j=1
*
*  is the scaling "quality" of i-th row. */

static double max_row_ratio(glp_prob *lp)
{     int i;
      double ratio, temp;
      ratio = 1.0;
      for (i = 1; i <= lp->m; i++)
      {  temp = max_row_aij(lp, i, 1) / min_row_aij(lp, i, 1);
         if (i == 1 || ratio < temp) ratio = temp;
      }
      return ratio;
}

/***********************************************************************
*  max_col_ratio - determine worst scaling "quality" for columns
*
*  This routine returns the worst scaling "quality" for columns of the
*  currently scaled constraint matrix:
*
*              n
*     ratio = max ratio[j],
*             j=1
*  where:
*                 m              m
*     ratio[j] = max |a[i,j]| / min |a[i,j]|,  1 <= j <= n,
*                i=1            i=1
*
*  is the scaling "quality" of j-th column. */

static double max_col_ratio(glp_prob *lp)
{     int j;
      double ratio, temp;
      ratio = 1.0;
      for (j = 1; j <= lp->n; j++)
      {  temp = max_col_aij(lp, j, 1) / min_col_aij(lp, j, 1);
         if (j == 1 || ratio < temp) ratio = temp;
      }
      return ratio;
}

/***********************************************************************
*  gm_iterate - perform iterative geometric mean scaling
*
*  This routine performs iterative geometric mean scaling of rows and
*  columns of the constraint matrix.
*
*  The parameter it_max specifies the maximal number of iterations.
*  Recommended value of it_max is 15.
*
*  The parameter tau specifies a minimal improvement of the scaling
*  "quality" on each iteration, 0 < tau < 1. It means than the scaling
*  process continues while the following condition is satisfied:
*
*     ratio[k] <= tau * ratio[k-1],
*
*  where ratio = max |a[i,j]| / min |a[i,j]| is the scaling "quality"
*  to be minimized, k is the iteration number. Recommended value of tau
*  is 0.90. */

static void gm_iterate(glp_prob *lp, int it_max, double tau)
{     int k, flag;
      double ratio = 0.0, r_old;
      /* if the scaling "quality" for rows is better than for columns,
         the rows are scaled first; otherwise, the columns are scaled
         first */
      flag = (max_row_ratio(lp) > max_col_ratio(lp));
      for (k = 1; k <= it_max; k++)
      {  /* save the scaling "quality" from previous iteration */
         r_old = ratio;
         /* determine the current scaling "quality" */
         ratio = max_mat_aij(lp, 1) / min_mat_aij(lp, 1);
#if 0
         xprintf("k = %d; ratio = %g\n", k, ratio);
#endif
         /* if improvement is not enough, terminate scaling */
         if (k > 1 && ratio > tau * r_old) break;
         /* otherwise, perform another iteration */
         gm_scaling(lp, flag);
      }
      return;
}

/***********************************************************************
*  NAME
*
*  scale_prob - scale problem data
*
*  SYNOPSIS
*
*  #include "glpscl.h"
*  void scale_prob(glp_prob *lp, int flags);
*
*  DESCRIPTION
*
*  The routine scale_prob performs automatic scaling of problem data
*  for the specified problem object. */

static void scale_prob(glp_prob *lp, int flags)
{     static const char *fmt =
         "%s: min|aij| = %10.3e  max|aij| = %10.3e  ratio = %10.3e\n";
      double min_aij, max_aij, ratio;
      xprintf("Scaling...\n");
      /* cancel the current scaling effect */
      glp_unscale_prob(lp);
      /* report original scaling "quality" */
      min_aij = min_mat_aij(lp, 1);
      max_aij = max_mat_aij(lp, 1);
      ratio = max_aij / min_aij;
      xprintf(fmt, " A", min_aij, max_aij, ratio);
      /* check if the problem is well scaled */
      if (min_aij >= 0.10 && max_aij <= 10.0)
      {  xprintf("Problem data seem to be well scaled\n");
         /* skip scaling, if required */
         if (flags & GLP_SF_SKIP) goto done;
      }
      /* perform iterative geometric mean scaling, if required */
      if (flags & GLP_SF_GM)
      {  gm_iterate(lp, 15, 0.90);
         min_aij = min_mat_aij(lp, 1);
         max_aij = max_mat_aij(lp, 1);
         ratio = max_aij / min_aij;
         xprintf(fmt, "GM", min_aij, max_aij, ratio);
      }
      /* perform equilibration scaling, if required */
      if (flags & GLP_SF_EQ)
      {  eq_scaling(lp, max_row_ratio(lp) > max_col_ratio(lp));
         min_aij = min_mat_aij(lp, 1);
         max_aij = max_mat_aij(lp, 1);
         ratio = max_aij / min_aij;
         xprintf(fmt, "EQ", min_aij, max_aij, ratio);
      }
      /* round scale factors to nearest power of two, if required */
      if (flags & GLP_SF_2N)
      {  int i, j;
         for (i = 1; i <= lp->m; i++)
            glp_set_rii(lp, i, round2n(glp_get_rii(lp, i)));
         for (j = 1; j <= lp->n; j++)
            glp_set_sjj(lp, j, round2n(glp_get_sjj(lp, j)));
         min_aij = min_mat_aij(lp, 1);
         max_aij = max_mat_aij(lp, 1);
         ratio = max_aij / min_aij;
         xprintf(fmt, "2N", min_aij, max_aij, ratio);
      }
done: return;
}

/***********************************************************************
*  NAME
*
*  glp_scale_prob - scale problem data
*
*  SYNOPSIS
*
*  void glp_scale_prob(glp_prob *lp, int flags);
*
*  DESCRIPTION
*
*  The routine glp_scale_prob performs automatic scaling of problem
*  data for the specified problem object.
*
*  The parameter flags specifies scaling options used by the routine.
*  Options can be combined with the bitwise OR operator and may be the
*  following:
*
*  GLP_SF_GM      perform geometric mean scaling;
*  GLP_SF_EQ      perform equilibration scaling;
*  GLP_SF_2N      round scale factors to nearest power of two;
*  GLP_SF_SKIP    skip scaling, if the problem is well scaled.
*
*  The parameter flags may be specified as GLP_SF_AUTO, in which case
*  the routine chooses scaling options automatically. */

void glp_scale_prob(glp_prob *lp, int flags)
{     if (flags & ~(GLP_SF_GM | GLP_SF_EQ | GLP_SF_2N | GLP_SF_SKIP |
                    GLP_SF_AUTO))
         xerror("glp_scale_prob: flags = 0x%02X; invalid scaling option"
            "s\n", flags);
      if (flags & GLP_SF_AUTO)
         flags = (GLP_SF_GM | GLP_SF_EQ | GLP_SF_SKIP);
      scale_prob(lp, flags);
      return;
}

/* eof */