simplex_lpsolve.c

Go to the documentation of this file.

#include "simplex.h"
#include "error.h"
#include "geom.h"
#include "defines.h"

#include <math.h>
#include <string.h>

void SimplexCreate(double epsilon,unsigned int ncols,TSimplex *s)
{
  double ct;

  /* LP: Create an empty LP problem */
  s->lp=make_lp(0,ncols);
  if (s->lp==NULL)
    Error("Error creating the simplex (simplex_lpsolve)");

  /* LP: Adjust the verbosity of the LP problem */
  #if (_DEBUG<5)
    set_verbose(s->lp,NEUTRAL);
  #endif

  s->inf=get_infinite(s->lp);

  /* LP: Adjust LP internal parameters */
  set_presolve(s->lp,PRESOLVE_NONE,0);
  set_scaling(s->lp,SCALE_NONE);
  
  ct=get_epsb(s->lp);
  if (ct>(epsilon*1e-2))
    set_epsb(s->lp,ct*1e-2);

  ct=get_epsd(s->lp);
  if (ct>(epsilon*1e-2))
    set_epsd(s->lp,ct*1e-2);
  
  ct=get_epsel(s->lp);
  if (ct>(epsilon*1e-2))
    set_epsel(s->lp,ct*1e-2);

  ct=get_epspivot(s->lp);
  if (ct>(epsilon*1e-4))
    set_epspivot(s->lp,ct*1e-4);
  
  /* LP: Minimize */
  set_minim(s->lp);

  /* LP: Timeout in seconds*/
  set_timeout(s->lp,ncols*SIMPLEX_TIMEOUT); 
}

/* LP: Reset the LP internal information */
void ResetSimplex(TSimplex *s)
{
  default_basis(s->lp);
}

unsigned int SimplexNColumns(TSimplex *s)
{
  return(get_Ncolumns(s->lp));
}

unsigned int SimplexNRows(TSimplex *s)
{
  return(get_Nrows(s->lp));
}

/* LP: set bounds for a given column */
void SimplexSetColBounds(unsigned int ncol,Tinterval *i,TSimplex *s)
{
  double l,u;
  unsigned int ncolInt;

  ncolInt=ncol+1;

  l=LowerLimit(i);
  u=UpperLimit(i);

  if (l==-INF) l=-s->inf;
  if (l== INF) l= s->inf;
  if (u==-INF) u=-s->inf;
  if (u== INF) u= s->inf;
  
  set_lowbo(s->lp,ncolInt,l);
  set_upbo(s->lp,ncolInt,u);
}

/* LP: get bounds for a given column */
void SimplexGetColBounds(unsigned int ncol,Tinterval *i,TSimplex *s)
{
  double l,u;
  unsigned int ncolInt;

  ncolInt=ncol+1;

  l=get_lowbo(s->lp,ncolInt);
  u=get_upbo(s->lp,ncolInt);

  if (l==-s->inf) l=-INF;
  if (l== s->inf) l= INF;
  if (u==-s->inf) u=-INF;
  if (u== s->inf) u= INF;  

  NewInterval(l,u,i);
}

void SimplexGetColConstraint(unsigned int ncol,TLinearConstraint *lc,TSimplex *s)
{
  unsigned int m,i,k;
  double *val;
  signed int *ind;
  Tinterval error;

  m=SimplexNRows(s);

  NEW(val,m+1,double);
  NEW(ind,m+1,signed int);

  k=get_columnex(s->lp,ncol+1,val,ind);

  InitLinearConstraint(lc);
  for(i=0;i<k;i++)
    AddTerm2LinearConstraint(ind[i]-1,val[i],lc);
 
  SimplexGetColBounds(ncol,&error,s);
  SetLinearConstraintError(&error,lc);

  free(val);
  free(ind);
}

boolean SimplexColEmpty(unsigned int ncol,TSimplex *s)
{  
  unsigned int m;
  double *val;
  signed int *ind;
  signed int k;

  m=SimplexNRows(s);

  NEW(val,m+1,double);
  NEW(ind,m+1,signed int);

  k=get_columnex(s->lp,ncol+1,val,ind);

  free(val);
  free(ind);

  return(k==0);
}

double SimplexGetColPrimal(unsigned int ncol,TSimplex *s)
{
  return(get_var_primalresult(s->lp,1+SimplexNRows(s)+ncol));
}

double SimplexGetColDual(unsigned int ncol,TSimplex *s)
{
  return(get_var_dualresult(s->lp,1+SimplexNRows(s)+ncol));
}

/* LP: set bounds for a given row */
void SimplexSetRowBounds(unsigned int nrow,Tinterval *i,TSimplex *s)
{
  double l,u;
  int t;  
  unsigned int nrowInt;

  nrowInt=nrow+1;

  t=get_constr_type(s->lp,nrowInt);

  l=LowerLimit(i);
  u=UpperLimit(i);

  if (l==-INF)
    {
      if (u==INF)
        {
          set_constr_type(s->lp,nrowInt,GE);
          set_rh(s->lp,nrowInt,-s->inf);
          set_rh_range(s->lp,nrowInt,s->inf);
        }
      else
        {
          set_constr_type(s->lp,nrowInt,LE);
          set_rh(s->lp,nrowInt,u);
        }
    }
  else
    {
      if (u==INF)
        {
          set_constr_type(s->lp,nrowInt,GE);
          set_rh(s->lp,nrowInt,l);
        }
      else
        {
          if (u==l)
            {
              set_constr_type(s->lp,nrowInt,EQ);
              set_rh(s->lp,nrowInt,l);
            }
          else  
            {
              set_constr_type(s->lp,nrowInt,GE);
              set_rh(s->lp,nrowInt,l);
              set_rh_range(s->lp,nrowInt,IntervalSize(i));
            }
        }
    }
}

/* LP: get bounds for a given row */
void SimplexGetRowBounds(unsigned int nrow,Tinterval *i,TSimplex *s)
{
  signed int t;
  double r,g;
  unsigned int nrowInt;

  nrowInt=nrow+1;

  t=get_constr_type(s->lp,nrowInt);
  r=get_rh(s->lp,nrowInt);
  g=get_rh_range(s->lp,nrowInt);

  switch(t)
    {
    case LE:
      if (g==s->inf)
        NewInterval(-INF,r,i);
      else
        NewInterval(r-g,r,i);
      break;
    case GE:
      if (g==s->inf)
        {
          if (r==-s->inf)
            NewInterval(-INF,INF,i);
          else
            NewInterval(r,INF,i);
        }
      else
        NewInterval(r,r+g,i);
      break;
    case EQ:
      if (g==s->inf)
        NewInterval(r,r,i);
      else
        NewInterval(r,r+g,i);
      break;
    }
}

void SimplexGetRowConstraint(unsigned int nrow,TLinearConstraint *lc,TSimplex *s)
{
  unsigned int i,n,m;
  unsigned int nrowInt;
  double *val;
  Tinterval error;

  nrowInt=nrow+1;

  InitLinearConstraint(lc);

  n=SimplexNColumns(s);
  m=SimplexNRows(s);

  NEW(val,m+1,double);

  for(i=1;i<=n;i++)
    {
      get_column(s->lp,i,val);
      AddTerm2LinearConstraint(i-1,val[nrowInt],lc);
    }
  free(val);

  SimplexGetRowBounds(nrow,&error,s);
  SetLinearConstraintError(&error,lc);
}


double SimplexGetRowPrimal(unsigned int nrow,TSimplex *s)
{
  return(get_var_primalresult(s->lp,1+nrow));
}

double SimplexGetRowDual(unsigned int nrow,TSimplex *s)
{
  return(get_var_dualresult(s->lp,1+nrow));
}

void SimplexAddNewConstraintRaw(TLinearConstraint *lc,TSimplex *s)
{
  Tinterval bounds;
  unsigned int nrow;
  unsigned int i;
  signed int *ind;

  /* LP: Add a row to the tableau */
  nrow=SimplexNRows(s); /* number of the row to be added next */

  NEW(ind,n,signed int);
  for(i=0;i<n;i++)
    ind[i]=(signed int)GetLinearConstraintVariable(i,lc)+1;
        
  /* LP: set the row information  */
  add_constraintex(s->lp,n,
                   GetLinearConstraintCoefficients(lc),
                   ind,
                   EQ,0);

  free(ind);

  GetLinearConstraintError(&bounds,lc);
  SimplexSetRowBounds(nrow,&bounds,s);

#if (_DEBUG>4)
  printf("          Setting Row %u with range ",nrow);
  PrintInterval(stdout,&bounds);
  printf("\n            ");
  PrintLinearConstraint(stdout,TRUE,NULL,lc);
#endif
}

void SimplexSetOptimizationFunction(TLinearConstraint *obj,TSimplex *s)
{
  unsigned int i,n;
  signed int *ind;

  #if (_DEBUG>4)
    printf("Setting cost function: ");
  #endif

  n=GetNumTermsInLinearConstraint(obj);

  #if (_DEBUG>4) 
    for(i=0;i<n;i++)
      {
        printf("+%f*v[%u]",
               GetLinearConstraintCoefficient(i,obj),
               GetLinearConstraintVariable(i,obj));
      }
  #endif

  NEW(ind,n,signed int);
  for(i=0;i<n;i++)
    ind[i]=(signed int)GetLinearConstraintVariable(i,obj)+1;

  set_obj_fnex(s->lp,n,
               GetLinearConstraintCoefficients(obj),
               ind);
  free(ind);

  #if (_DEBUG>4)
    printf("\n");
  #endif
}

void SimplexGetOptimizationFunction(TLinearConstraint *obj,TSimplex *s)
{
  unsigned int i,n,m;
  double *val;

  InitLinearConstraint(obj);

  n=SimplexNColumns(s);
  m=SimplexNRows(s);

  NEW(val,m+1,double);

  for(i=1;i<=n;i++)
    {
      get_column(s->lp,i,val);
      AddTerm2LinearConstraint(0,i-1,val[0],obj);
    }
  free(val);
}

double SimplexGetOptimalValueRaw(TSimplex *s)
{
  return(get_objective(s->lp));
}

unsigned int SimplexSolve(TSimplex *s)
{
  unsigned int r;

  /*print_lp(s->lp); exit(0);*/

  r=solve(s->lp);

  if (r==OPTIMAL)
    return(REDUCED_BOX);
  else
    {
      if (r==INFEASIBLE)
        return(EMPTY_BOX);
      else
        {
          if (r==UNBOUNDED)
            return(UNBOUNDED_BOX);
          else
            return(ERROR_IN_PROCESS);
        }
    }
}

void SimplexDelete(TSimplex *s)
{
  delete_lp(s->lp);
}