cuiksystem.c

Go to the documentation of this file.

#include "cuiksystem.h"

#include "error.h"
#include "boolean.h"
#include "defines.h"
#include "box_list.h"
#include "box_heap.h"
#include "simplex.h"
#include "random.h"
#include "filename.h"
#include "geom.h"

#include <stdlib.h>
#include <gsl/gsl_linalg.h>

#include <math.h>
#include <string.h>
#include <sys/time.h>
#include <signal.h>
#if (_USE_MPI)
  #include <mpi.h>
  #include <sys/resource.h>
  #include <time.h>
  #include <unistd.h>
#endif



/* Internal functions */

void DummifyAndAddEquation(Tparameters *p,Tequation *eq,TCuikSystem *cs);

void DummifyCuikSystem(Tparameters *p,TCuikSystem *cs);


double EvaluateEqMin(void *b,void *cs);

unsigned int ReduceBoxEquationWise(Tparameters *p,Tbox *b,TCuikSystem *cs);

unsigned int ReduceBox(Tparameters *p,unsigned int varMask,Tbox *b,TCuikSystem *cs);

void PostProcessBox(Tparameters *p,unsigned int c,
                    FILE *f_out,
                    Tlist *sol,
                    Theap *boxes,
                    Tbox *b,
                    TCuikSystem *cs);

void CSRemoveUnusedVars(Tparameters *p,TCuikSystem *cs);


boolean CSRemoveVarsWithCtRange(Tparameters *p,
                                boolean *replaced,TLinearConstraint *lc,
                                Tbox *borig,TCuikSystem *cs);

boolean CSRemoveLCVars(Tparameters *p,unsigned int level,boolean *changed,
                       boolean *replaced,TLinearConstraint *lc,
                       Tbox *borig,TCuikSystem *cs);
boolean SimplifyCuikSystem(Tparameters *p,TCuikSystem *cs);

void ComputeSimpCuikSystemJacobian(TCuikSystem *cs);

void DeleteSimpCuikSystemJacobian(TCuikSystem *cs);

boolean UpdateCuikSystem(Tparameters *p,TCuikSystem *cs);

void UnUpdateCuikSystem(TCuikSystem *cs);

unsigned int ComputeSplitDimInt(Tparameters *p,Tbox *b,TCuikSystem *cs);

void SaveCSState(char *fname,Tlist *lb,TCuikSystem *cs);

void LoadCSState(char *fname,Tlist *lb,TCuikSystem *cs);


/************************************************************************/
/************************************************************************/
/************************************************************************/

/*
  The equation is dummified and the resulting equations and new variables
  are added to the cuik system.

  Dummification is performed in the last step of the simplification and, thus,
  this function operates in the sets of simplified equations and variables

  Dummification consists in changing all monomials of type x^2 or x*y by new
  variables 'z' and adding the corresponding equations (z=x^2 or z=x*y) to
  the system.

  The level of dummification is controlled by the 'dummify' parameter.

  This internal function is hidden to the user. In this way dummification is
  applied only after the system simplification.
  See AddEquation2CS and SimplifyCuikSystem.
*/
void DummifyAndAddEquation(Tparameters *p,Tequation *eq,TCuikSystem *cs)
{
  if (eq!=NULL)
    {
      unsigned int nf;
      unsigned int i,nv,mv;
      Tvariable_set *vars;
      unsigned int dummify;
      double epsilon;
      unsigned int eqType;

      dummify=(unsigned int)(GetParameter(CT_DUMMIFY,p));
      epsilon=GetParameter(CT_EPSILON,p);

      vars=GetEquationVariables(eq);
      eqType=GetEquationType(eq);

      nv=VariableSetSize(vars);
      mv=NVariables(&(cs->variables));
      for(i=0;i<nv;i++)
        {
          if (GetVariableN(i,vars)>=mv)
            Error("Unknown variable in AddEquation2CS");
        }

      nf=GetNumMonomials(eq);
      if (nf>0)
        {
          boolean keep; /* True if the equation is kept in its original form */
          
          /* Dummy equations are not dummified again */
          /*dummify equations with just one monomial introduce a
            too simple variable x=ct that is not simplified since
            the dummifications is performed after the simplification. */
          if ((eqType==DUMMY_EQ)||(GetNumMonomials(eq)==1))
            keep=TRUE; 
          else
            {
              switch (dummify)
                {
                case 0:
                  keep=TRUE;
                  break;
                case 1:
                  keep=((LinearEquation(eq))||
                        (ParabolaEquation(eq))|| 
                        (SaddleEquation(eq))|| 
                        (CircleEquation(eq))|| 
                        (SphereEquation(eq)));
                  break;
                case 2:
                  keep=((LinearEquation(eq))||
                        (ParabolaEquation(eq))|| 
                        (SaddleEquation(eq)));
                  break;
                case 3:
                  keep=((LinearEquation(eq))||
                        (ParabolaEquation(eq))|| 
                        (SaddleEquation(eq))|| 
                        (CircleEquation(eq))|| 
                        (SphereEquation(eq))||
                        (GetEquationCmp(eq)!=EQU)); /*Inequalities are kept in original form*/
            
                  break;
                case 4:
                  keep=((LinearEquation(eq))||
                        (ParabolaEquation(eq))|| 
                        (SaddleEquation(eq))||
                        (GetEquationCmp(eq)!=EQU)); /*Inequalities are kept in original form*/
                  break;
                default:
                  Error("Undefined value for DUMMIFY parameter");
                }
            }

          if (keep)
            AddEquation(eq,&(cs->equations));
          else
            {
              unsigned int j,p,vid1,vid2,vid3;
              Tmonomial *f,fnew;
              Tequation eqnew; /* Bilinear and quadratic monomials produce new auxiliary equations */
              Tequation eqlineal; /* The original equation is transformed into a lineal equation */
              char *vname;
              Tinterval range;
              Tvariable v;
              char *n1,*n2;

              InitEquation(&eqlineal);
              
              for(j=0;j<nf;j++) /* For all monomials in the equation to be dummified */
                {

                  f=GetMonomial(j,eq);
                  vars=GetMonomialVariables(f);
                  nv=VariableSetSize(vars);

                  switch(nv)
                    {
                    case 1: /*A single variable in the monomial*/
                      p=GetVariablePowerN(0,vars);

                      switch(p)
                        {
                        case 1: /*One variable with degree 1*/
                          /* A linear monomial: just add it to the linear equation*/
                      
                          AddMonomial(f,&eqlineal); /* A copy of monomial 'f' is stored in eqlineal*/
                          break;

                        case 2:/*One varialbe With degree 2*/                 
                          /* Add a linear monomial to eqlinear with a dummy variable
                             and add the corresponding parabola equation*/

                          vid1=GetVariableN(0,vars);

                          /* Generate a new variable dummy that will stand for  x^2 */
                          n1=GetVariableName(GetVariable(vid1,&(cs->variables)));
                          NEW(vname,strlen(n1)+20,char);
                          sprintf(vname,"dummy_%s_2",n1);
                                                    
                          vid2=GetVariableID(vname,&(cs->variables));

                          if (vid2==NO_UINT)
                            {
                              /* if the dummy variable is a new one, add the corresponding
                                 dummy equation to the system */

                              NewVariable(DUMMY_VAR,vname,&v);
                              IntervalPow(GetVariableInterval(GetVariable(vid1,&(cs->variables))),2,&range);
                              SetVariableInterval(&range,&v);
                              
                              vid2=AddVariable(&v,&(cs->variables));
                              DeleteVariable(&v);
                              
                              /* Generate a new equation x^2-dummy=0 */
                              InitEquation(&eqnew);
                              
                              InitMonomial(&fnew);
                              AddVariable2Monomial(vid1,2,&fnew);
                              AddCt2Monomial( 1.0,&fnew);
                              AddMonomial(&fnew,&eqnew);
                              DeleteMonomial(&fnew);
                              
                              InitMonomial(&fnew);
                              AddVariable2Monomial(vid2,1,&fnew);
                              AddCt2Monomial(-1.0,&fnew);
                              AddMonomial(&fnew,&eqnew);
                              DeleteMonomial(&fnew);
                              
                              SetEquationCmp(EQU,&eqnew);
                              SetEquationValue(0.0,&eqnew);
                              /*The derive equations is of the same type as the origiona one*/
                              SetEquationType(DUMMY_EQ,&eqnew);
                              
                              AddEquation(&eqnew,&(cs->equations));
                              DeleteEquation(&eqnew);
                            }
                          free(vname);

                          /*Now add the lineal monomial using the new dummy variable to eqlineal*/
                          InitMonomial(&fnew);
                          AddVariable2Monomial(vid2,1,&fnew);
                          AddCt2Monomial(GetMonomialCt(f),&fnew);
                          AddMonomial(&fnew,&eqlineal);
                          DeleteMonomial(&fnew);
                          break;

                        default:
                          Error("The input system to DummifyCuikSystem is not in canonical form");
                          break;
                        }
                      break;

                    case 2: /*Two variables in the monomial*/
                      if ((GetVariablePowerN(0,vars)==1)&&
                          (GetVariablePowerN(1,vars)==1)) /*The two of them with degree 1*/
                        {
                          /*Get the two variables involved in the product we have to rewrite*/
                          vid1=GetVariableN(0,vars);
                          vid2=GetVariableN(1,vars);
                          
                          /*Generate a new dummy variable*/
                          n1=GetVariableName(GetVariable(vid1,&(cs->variables)));
                          n2=GetVariableName(GetVariable(vid2,&(cs->variables)));
                          NEW(vname,strlen(n1)+strlen(n2)+20,char);
                          sprintf(vname,"dummy_%s_%s",n1,n2);
                                        
                          vid3=GetVariableID(vname,&(cs->variables));
                          if (vid3==NO_UINT)
                            {

                              NewVariable(DUMMY_VAR,vname,&v);
                              IntervalProduct(GetVariableInterval(GetVariable(vid1,&(cs->variables))),
                                              GetVariableInterval(GetVariable(vid2,&(cs->variables))),
                                              &range);
                              SetVariableInterval(&range,&v);
                              
                              vid3=AddVariable(&v,&(cs->variables));
                              DeleteVariable(&v);

                              /* Generate the dummy equation X*Y-dummy=0*/
                              InitEquation(&eqnew);
                              
                              InitMonomial(&fnew);
                              AddVariable2Monomial(vid1,1,&fnew);
                              AddVariable2Monomial(vid2,1,&fnew);
                              AddCt2Monomial( 1.0,&fnew);
                              AddMonomial(&fnew,&eqnew);
                              DeleteMonomial(&fnew);
                              
                              InitMonomial(&fnew);
                              AddVariable2Monomial(vid3,1,&fnew);
                              AddCt2Monomial(-1.0,&fnew);
                              AddMonomial(&fnew,&eqnew);
                              DeleteMonomial(&fnew);
                              
                              SetEquationCmp(EQU,&eqnew);
                              SetEquationValue(0.0,&eqnew);
                              /*The derive equations is of the same type as the origiona one*/
                              SetEquationType(DUMMY_EQ,&eqnew);

                              AddEquation(&eqnew,&(cs->equations));
                              DeleteEquation(&eqnew);
                            }
                          free(vname);

                          /*Now add the lineal monomial using the new dummy variable to eqlienal*/
                          InitMonomial(&fnew);
                          AddVariable2Monomial(vid3,1,&fnew);
                          AddCt2Monomial(GetMonomialCt(f),&fnew);
                          AddMonomial(&fnew,&eqlineal);
                          DeleteMonomial(&fnew);
                        }
                      else
                        /* One of the 2 variables in the monomial has degree >1 */
                        Error("The input system to DummifyCuikSystem is not in canonical form");
                      break;

                    default:
                      /* A monomial with more than 2 variables */
                      Error("The input system to DummifyCuikSystem is not in canonical form");
                    }

                } /* Proceed with next monomial in the under-dummification equation*/

              /* The linearized original equation is added to the system */

              SetEquationCmp(GetEquationCmp(eq),&eqlineal);
              SetEquationValue(GetEquationValue(eq),&eqlineal);
              /*The resulting lineal equation is of the same type as the original equation*/
              SetEquationType(eqType,&eqlineal);

              AddEquation(&eqlineal,&(cs->equations));
              DeleteEquation(&eqlineal);
            }
        }
    }
}

void DummifyCuikSystem(Tparameters *p,TCuikSystem *cs)
{
  unsigned int i,neq;
  Tequations eqs;
  unsigned int dummify;
  double epsilon;
  
  dummify=(unsigned int)(GetParameter(CT_DUMMIFY,p));
  epsilon=GetParameter(CT_EPSILON,p);

  CopyEquations(&eqs,&(cs->equations));
  DeleteEquations(&(cs->equations));
  InitEquations(&(cs->equations));

  neq=NEquations(&eqs);
  for(i=0;i<neq;i++)
    DummifyAndAddEquation(p,GetEquation(i,&eqs),cs);

  DeleteEquations(&eqs);
}

unsigned int ReduceBoxEquationWise(Tparameters *p,Tbox *b,TCuikSystem *cs)
{
  boolean empty;
  boolean small;
  boolean significantReduction;
  Tbox b_orig;
  unsigned int i,j;
  double rho,smallSigma,epsilon;
  double mi_r;
  double reduction;

  rho=GetParameter(CT_RHO,p);
  epsilon=GetParameter(CT_EPSILON,p);
  smallSigma=GetParameter(CT_SMALL_SIGMA,p);

  /*default values*/
  empty=FALSE;
  significantReduction=FALSE;
  small=FALSE;

  #if (_DEBUG>2)
    printf("    Croping box: ");
    PrintBox(stdout,b);
    printf("      Size in: %g\n",GetBoxSize(cs->systemVar,b));
  #endif

  do {
    CopyBox(&b_orig,b);

    for(j=0;((!empty)&&(j<cs->nequations));j++)
      {
        #if (_DEBUG>5)
          printf("          Croping with equation %u: \n",j);
          printf("            Input box: "); 
          PrintBox(stdout,b);
          printf("            s: %g  v:%g\n",GetBoxSize(NULL,b),GetBoxVolume(NULL,b));
        #endif
        if (CropEquation(j,epsilon,rho,b,&(cs->equations))==EMPTY_BOX)
          empty=TRUE;
        #if (_DEBUG>5)
          if (empty)
            printf("            Empty output box\n");
          else
            {
              printf("            Output box: "); 
              PrintBox(stdout,b);
              printf("            s: %g  v:%g\n",GetBoxSize(NULL,b),GetBoxVolume(NULL,b));
            }
        #endif
      }

    if (!empty)
      {
        mi_r=1.0;
        for(i=0;i<cs->nvariables;i++)
          {
            /*variables of all types are affected by crop*/
            reduction=IntervalSize(GetBoxInterval(i,b))/IntervalSize(GetBoxInterval(i,&b_orig));
            if (reduction<mi_r) mi_r=reduction; 
          }

        significantReduction=(mi_r<rho);
        small=(GetBoxSize(cs->systemVar,b)<=smallSigma);

        #if (_DEBUG>2)
          printf("      Crop reduction: %f\n",mi_r);
          if (small)
            printf("        Tiny boxes are not reduced any more\n");
          else
            {
              if (significantReduction) 
                printf("        Large reduction. Cropping again\n");
            }
        #endif
      }
    DeleteBox(&b_orig);
    
  } while ((!empty)&&(!small)&&(significantReduction));
  
  #if (_DEBUG>2)
    if (empty)
      printf("    Empty in Crop\n");
    else
      {
        printf("    Croped box: ");
        PrintBox(stdout,b);
        printf("      Size out: %g\n",GetBoxSize(cs->systemVar,b));
      }
  #endif
    
  if (empty)
    return(EMPTY_BOX);
  else
    return(REDUCED_BOX);
}

/*
  This is the basic function of the solver. This function includes the
  ShrinkBox procedure described in our papers plus a loop interating
  ShrinkBox as far as the box is significantly reduced (and still not empty).

  Thus, ReduceBox takes an input box 'b' and reduces it as much as possible by removing
  parts that do not include any solution point.
  It returns information about what happened in the box shrinking process

      ERROR_IN_PROCESS: When one of the simplex failed due to numerical inestabilities
                        The better the simplex implementation the less likely this
                        output.

      EMPTY_BOX: The box includes no solution. In this case the output box is
                 undefined (it includes empty intervals).

      REDUCED_BOX: The box was reduced and it has to be processed further to determine
                   if it includes any solution.

      REDUCED_BOX_WITH_SOLUTION: The box was been reduced and one of the reduction
                                 steps was contractive. This indicates the box
                                 includes a solution point.
 */
unsigned int ReduceBox(Tparameters *p,unsigned int varMask,Tbox *b,TCuikSystem *cs)
{
  /*The box can only be reduced if we have constraints*/
  if (cs->empty)
    return(REDUCED_BOX);
  else
    {
      boolean empty;
      boolean significantReduction;
      boolean small;
      boolean hasSolution;
      boolean err;
      boolean boxOK;
      unsigned int e;
      double mi_r,ma_s;
      double reduction,size;
      unsigned int i,j;
      Tbox b_orig;
      TSimplex lp;
      double epsilon,smallSigma,sigma,rho,lr2tm_q,lr2tm_s;
      unsigned int safe_simplex;
      boolean repeatLoop;

      if (!cs->updated)
        Error("Non-updated CS in ReduceBox");

      epsilon=GetParameter(CT_EPSILON,p);
      sigma=GetParameter(CT_SIGMA,p);
      smallSigma=GetParameter(CT_SMALL_SIGMA,p);
      rho=GetParameter(CT_RHO,p);
      lr2tm_q=GetParameter(CT_LR2TM_Q,p);
      lr2tm_s=GetParameter(CT_LR2TM_S,p);
      safe_simplex=(unsigned int)(GetParameter(CT_SAFE_SIMPLEX,p));

      /*default values*/
      empty=FALSE;
      hasSolution=FALSE;
      err=FALSE;                
      boxOK=TRUE;
      small=FALSE;
      significantReduction=FALSE;
      repeatLoop=FALSE;

      do { 
        #if (_DEBUG>2)
          printf("\n  Go for box reduction via crop+simplex\n");        
          printf("  InputBox: "); PrintBox(stdout,b);
          printf("    s: %g  v:%g\n",GetBoxSize(NULL,b),GetBoxVolume(NULL,b));
        #endif

        /*****************************************************************/
        /* One step of box reduction */

        /*#if ((!_USE_MPI)&&(_DEBUG==1))*/
        #if ((!_USE_MPI)&&(_DEBUG>1))
          fprintf(stderr,".");fflush(stderr);
        #endif

        NewBoxReduction(&(cs->st));

        /* First we do an equation wise reduction of the ranges -> Crop equation 
           repeated for each equation as far as we get a significant box reduction */

        empty=(ReduceBoxEquationWise(p,b,cs)==EMPTY_BOX);

        /* If the equation-wise reduction gives a full box, proceed with the global
           test based on linear programming */

        if (empty)
          repeatLoop=FALSE; /*There is no need to iterate*/
        else
          {
            SimplexCreate(epsilon,cs->nvariables,&lp);

            /* Add the linear constraints bounding the equations. This include
               croping the input ranges if necessary */

            for(j=0;((!empty)&&(j<cs->nequations));j++)
              {
                if (!AddEquation2Simplex(j,
                                         lr2tm_q,lr2tm_s,
                                         epsilon,rho,
                                         b,&(cs->variables),
                                         &lp,&(cs->equations)))
                  empty=TRUE;
              }

            boxOK=TRUE;
            err=FALSE;
            small=FALSE;
            significantReduction=FALSE;

            /* Use the simplex to reduce the ranges for the variables */
            if ((!empty)&&(SimplexNRows(&lp)>0)) /* In extreme cases the simplex tableau can be empty 
                                                    and then nothings can be reduced*/
              { 
                /* Store the original box to check the actual reduction */
                CopyBox(&b_orig,b);

                #if (_SIMPLEX_ENGINE!=_CLP) /* in CLP Reset Simplex has no effect */
                  ResetSimplex(&lp);
                #endif

                for(i=0;((boxOK)&&(i<cs->nvariables));i++)
                  {
                    /* The crop takes care of reducing the ranges of the dummy 
                       variables in function of the non dummy ones (system,
                       secondary and cartesian).
                    */

                    if ((cs->varType[i]&varMask)&&(!SimplexColEmpty(i,&lp)))
                      {
                        /*Reduce the ranges for variable 'i'*/
                        e=ReduceRange(epsilon,safe_simplex,i,b,&lp);
                        boxOK=(e==REDUCED_BOX); 
                      }             
                  }

                /* Analize the output out of the range reduction */
                if (boxOK)
                  {
                    /* Compute the minimum reduction ratio for all ranges. The smaller
                       the reductin ratio, the maximal the box shrink for that dimension*/
                    mi_r=1.0;
                    ma_s=0.0;
                    for(i=0;i<cs->nvariables;i++)
                      {
                        if (cs->notDummyVar[i])
                          {
                            size=IntervalSize(GetBoxInterval(i,b));
                            if (size>ma_s) ma_s=size;

                            reduction=size/IntervalSize(GetBoxInterval(i,&b_orig));
                            if (reduction<mi_r) mi_r=reduction; 
                          }
                      }

                    /* If the output box is fully included in the input one this means
                       the box includes a solution (for sure) */
                    if ((ma_s<lr2tm_q)&&(ma_s<lr2tm_s)&&
                        (BoxInclusion(cs->systemVar,b,&b_orig)))
                      {
                        hasSolution=TRUE;
                        #if(!_USE_MPI)
                          fprintf(stderr,"(>!<)");
                        #endif
                      }

                    small=(ma_s<=smallSigma);
                    significantReduction=(mi_r<rho); 
                  }
                else
                  {
                    empty=(e==EMPTY_BOX);
                    err=(e==ERROR_IN_PROCESS);
                  }

                DeleteBox(&b_orig);
              }

            repeatLoop=((!empty)&&(!err)&&(!small)&&(significantReduction));
            
            SimplexDelete(&lp);

            #if (_DEBUG>2)
            if (empty)
              printf("    Empty\n");
            else
              {
                if (err)
                  printf("    Simplex Error\n");
                else
                  {
                    printf("  OutputBox: "); PrintBox(stdout,b);
                    printf("  Reduction ->  %f\n",mi_r);
                  
                    printf("  Size after reduction: %g\n ",GetBoxSize(cs->systemVar,b));
                    printf("  Volume after reduction: %g\n ",GetBoxVolume(cs->systemVar,b));
                  
                    if (significantReduction)
                      printf("    Large enough reduction-> Try again all equations\n");
                    else
                      printf("    Not enough reduction.\n");
                  }
              }
            #endif 

          } /* end if not empty after */

      } while(repeatLoop);
      
      if (err)
        e=ERROR_IN_PROCESS;
      else
        {
          if (empty)
            e=EMPTY_BOX;
          else
            {
              /* If at least on step was a contraction this means
                 there is a solution inside the box*/
              if (hasSolution)
                e=REDUCED_BOX_WITH_SOLUTION; 
              else
                e=REDUCED_BOX;
            }
        }

      return(e);
    }
}

/*
  This procedure is called just after using ReduceBox to analyzed what
  happened with the box and what to do next.ç

  'c' is the return code out of Reduce Box

  If the box is small (max side size smaller than sigma), we consider it 
  a solution and we store it in the list 'sol' (if defined), we save 
  the solution into 'f_out' (if defined).
  
  Otherwise, we split the box in two sub-boxes and we store them in list 'boxes'
  The dimension along which the box is splitte is determined  according to 
  parameter 'error_split'. The two sub-boxes are added to 'boxes' in the first
  of in the last position according to parameter 'depth_first'.
*/
void PostProcessBox(Tparameters *p,
                    unsigned int c,
                    FILE *f_out,
                    Tlist *sol,
                    Theap *boxes,Tbox *b,
                    TCuikSystem *cs)
{
  boolean empty,small;
  Tbox b_orig;
  double sigma;

  empty=(c==EMPTY_BOX);
        
  if (empty)
    {
      NewEmptyBox(&(cs->st));
      #if (_DEBUG>1)
        printf("  The box is empty -> deleted\n");
      #endif
      #if (_DEBUG>0)
        fprintf(stderr,"e\n");
      #endif
    }
  else
    {
      /* c==REDUCED_BOX or c==REDUCED_BOX_WITH_SOLUTION or c==ERROR_IN_PROCESS*/
      #if (_DEBUG>1)
        printf("  Resulting box: ");
        PrintBox(stdout,b);
      #endif
      #if (_DEBUG>0)
        fprintf(stderr,"<%g %g>",
                GetBoxVolume(cs->systemVar,b),
                GetBoxSize(cs->systemVar,b));
      #endif  
        
      /* We initialize a box in the original system and set the default ranges
         (ranges in case the variable in the original system is not in the
         simplified one... this is not likely to occur here but in the system
         to sybsystem relation but in this way the interface is the same)*/
      BoxFromVariables(&b_orig,&(cs->orig_variables)); 
      UpdateOriginalFromSimple(b,&b_orig,cs->orig2sd);
      #if (_DEBUG>1)
        printf("  Original box: ");
        PrintBox(stdout,&b_orig);
      #endif

      sigma=GetParameter(CT_SIGMA,p);

      /*The original system can have variables not included in the simplified one. For instance, 
        variables not in the system equations due to the particular constants of the problem under
        analysis are eliminated, even though they appear in dummy equations (this is common
        when ROT_REP=1, see link.h). Those removed variables are never reduced and thus,
        the original box never reduces below sigma. This is why we use the simplified box
        (i.e., the one we actually reduce) to assess whether or not we have a solution. */
      small=(GetBoxSize(cs->systemVar,b)<=sigma);

      if (small)
        {

          #if (_DEBUG>1)
            printf("  The box is a solution (size %g vs [%g]) ->add to the list of solutions\n",
                   GetBoxSize(cs->orig_systemVar,&b_orig),sigma);
          #endif

          NewSolutionBox((c==REDUCED_BOX_WITH_SOLUTION),
                         GetBoxVolume(cs->orig_systemVar,&b_orig),
                         GetBoxDiagonal(cs->orig_systemVar,&b_orig),
                         &(cs->st));
          
          if(c==ERROR_IN_PROCESS)
            NewRBError(&(cs->st));
 
          #if (_DEBUG>0)
            if(c==REDUCED_BOX_WITH_SOLUTION)
              fprintf(stderr,"S\n");
            else
              fprintf(stderr,"s\n");
          #endif        
              
          
          if (f_out!=NULL)
            {
              if(c==ERROR_IN_PROCESS)
                fprintf(f_out," E");
              
              fprintf(f_out," %3u (err:%g",
                      GetNSolutionBoxes(&(cs->st)),
                      ErrorInSolution(&b_orig,cs));
              
              if (cs->searchMode==MINIMIZATION_SEARCH)
                fprintf(f_out," min:%g",EvaluateEqMin((void *)b,(void *)cs));

              fprintf(f_out," tm:%g):",GetElapsedTime(&(cs->st)));

              if (c==REDUCED_BOX_WITH_SOLUTION)
                fprintf(f_out," <S> ");

              /*The output does not include the dummy variables (they are 
                not relevant any more)*/
              PrintBoxSubset(f_out,cs->orig_notDummyVar,cs->orig_varNames,&b_orig); 
              fflush(f_out);          
            }
          if (sol!=NULL)
            {
              Tbox *sb;

              NEW(sb,1,Tbox); 

              CopyBox(sb,&b_orig);

              AddLastElement((void *)sb,sol);
            }

          /*If we want to ad box to a solution list, we have to copy 'b' 
            somewhere else ('b' is delete after this function !)*/
        }
      else
        {
          /* Non empty large box or error in ReduceBox */

          Tbox b1,b2;
          unsigned int d;

          #if (_DEBUG>1)
            printf("  Box size : %g vs. %g\n",GetBoxSize(cs->orig_systemVar,&b_orig),sigma);
            printf("  The box has to be splitted (and deleted)\n");
          #endif

          NewSplittedBox(&(cs->st));

          d=ComputeSplitDimInt(p,b,cs);

          SplitBox(d,CUT_POINT,&b1,&b2,b);

          #if (_DEBUG>1)
            printf("  Splitting along dimension %u (internal)\n",d);
            if (cs->searchMode==MINIMIZATION_SEARCH)
              printf("      Box penalgy : %g\n",EvaluateEqMin((void *)&b1,(void *)cs));
            printf("        SubBox:"); PrintBox(stdout,&b1);
            if (cs->searchMode==MINIMIZATION_SEARCH)
              printf("      Box penalty : %g\n",EvaluateEqMin((void *)&b2,(void *)cs));
            printf("        SubBox:"); PrintBox(stdout,&b2);
          #endif

          AddBox2HeapOfBoxes(&b1,boxes);
          AddBox2HeapOfBoxes(&b2,boxes);
          DeleteBox(&b1);
          DeleteBox(&b2);
         
          if (c==ERROR_IN_PROCESS)
            {
              #if (_DEBUG>0)
                fprintf(stderr,"E[%u]\n",d);
              #endif
              #if (_DEBUG>1)
                printf("  Splitted due to simplex error\n");
              #endif
              NewRBError(&(cs->st));
            }
          #if (_DEBUG>0)
          else
            fprintf(stderr,"b[%u]\n",d);
          #endif
        }

      DeleteBox(&b_orig);  
    }
}

void CSRemoveUnusedVars(Tparameters *p,TCuikSystem *cs)
{
  unsigned int i,vt;
  double epsilon;

  epsilon=GetParameter(CT_EPSILON,p);

  /*We proceed back to front because removing a variable change the
    indexes of variables above it.*/
  i=NVariables(&(cs->variables));
  while(i>0)
    {
      i--;
      vt=GetVariableTypeN(i,&(cs->variables));
      //if (((vt==DUMMY_VAR)&&(!UsedVarInNonDummyEquations(i,&(cs->equations))))||
      //  ((vt!=DUMMY_VAR)&&(!UsedVarInEquations(i,&(cs->equations)))))
      if (((vt==DUMMY_VAR)&&(!UsedVarInNonDummyEquations(i,&(cs->equations))))||
          ((vt!=SYSTEM_VAR)&&(!UsedVarInEquations(i,&(cs->equations)))))
        {
          /* ... Therefore, we can still have dummy equations with the non-used variables*/
          RemoveEquationsWithVar(epsilon,i,&(cs->equations)); 
          RemoveVariable(i,&(cs->variables));
        }
    }
}

boolean CSRemoveVarsWithCtRange(Tparameters *p,
                                boolean *replaced,TLinearConstraint *lc,
                                Tbox *borig,TCuikSystem *cs)
{
  boolean found,hold;
  unsigned int nv,origID1;
  unsigned int id1;
  char *name1;
  double epsilon;
  double ct2;
  Tinterval *range;
  
  epsilon=GetParameter(CT_EPSILON,p);

  hold=TRUE;

  do {
    found=FALSE;

    /* A variable with a ct range can be removed */
    nv=NVariables(&(cs->variables));
    id1=0;
    while((!found)&&(id1<nv))
      {
        name1=GetVariableName(GetVariable(id1,&(cs->variables)));
        origID1=GetVariableID(name1,&(cs->orig_variables));
        if (origID1==NO_UINT)
          Error("Undefined ID1 in CSRemoveVarsWithCtRange");

        range=GetBoxInterval(origID1,borig);
        if (IntervalSize(range)==0)
          {
            ct2=LowerLimit(range); /* ==UpperLimit(range) */
            found=TRUE;
          }
        else
          id1++;
      }

    /*if we found something to simplify, apply it to the rest of equations
      (this creates a new set of equations that replaces the previous one)*/
    if (found)
      {
        replaced[origID1]=TRUE;
        AddCt2LinearConstraint(ct2,&(lc[origID1]));

        #if (_DEBUG>5) 
          printf("Ct Range Replacement: %s [%u-%u] ->%.12g\n",name1,origID1,id1,ct2);
        #endif
            
        if (hold)
          {
            TLinearConstraint lc;

            InitLinearConstraint(&lc);
            AddCt2LinearConstraint(ct2,&lc);
            hold=ReplaceVariableInEquations(epsilon,id1,&lc,&(cs->equations));
            DeleteLinearConstraint(&lc);
            RemoveVariable(id1,&(cs->variables));
            
            #if (_DEBUG>6)
            if(hold)
              {
                char **varNames;
                
                
                printf("%%=======================================================\n");
                printf("%% One less variable (ct range)\n"); 
                printf("%%****************************************\n");
                PrintVariables(stdout,&(cs->variables));
                
                nv=NVariables(&(cs->variables));
                NEW(varNames,nv,char *);
                GetVariableNames(varNames,&(cs->variables));
                PrintEquations(stdout,varNames,&(cs->equations));
                free(varNames);
              }
            #endif
          }
      }
  } while((hold)&&(found));

  return(hold);
}

boolean CSRemoveLCVars(Tparameters *p,unsigned int nTerms,boolean *changed,
                       boolean *replaced,TLinearConstraint *lc,
                       Tbox *borig,TCuikSystem *cs)
{
  boolean found,hold;
  unsigned int m,i,j,nv,neq,origID,origID1;
  Tequation *eq;
  unsigned int id;
  char *name,*name1;
  double epsilon;
  Tvariable *v;
  TLinearConstraint lcr,lcc;
  boolean *systemVars;
  Tinterval error;
  unsigned int simplificationLevel;
  #if (_DEBUG>5) 
    char **varNamesO;

    nv=NVariables(&(cs->orig_variables));
    NEW(varNamesO,nv,char *);
    GetVariableNames(varNamesO,&(cs->orig_variables));
  #endif

  nv=NVariables(&(cs->variables));
  NEW(systemVars,nv,boolean);
  for(i=0;i<nv;i++)
    systemVars[i]=IsSystemVariable(i,&(cs->variables));

  epsilon=GetParameter(CT_EPSILON,p);
  simplificationLevel=(unsigned int)(GetParameter(CT_SIMPLIFICATION_LEVEL,p));

  *changed=FALSE;
  hold=TRUE;

  found=FALSE;
  i=0;
  neq=NEquations(&(cs->equations));
  while((!found)&&(i<neq))
    {
      eq=GetEquation(i,&(cs->equations));

      found=IsSimplificable(simplificationLevel,nTerms,systemVars,&id,&lcr,eq); /* 't' is the type of simplification*/
      if (!found)
        i++;
    }

  /*if we found something to simplify, apply it to the rest of equations
    (this creates a new set of equations that replaces the previous one)*/
  if (found)
    {
      *changed=TRUE;

      /* Translate the variables identifiers from local to global */
      v=GetVariable(id,&(cs->variables));
      name=GetVariableName(v);
      origID=GetVariableID(name,&(cs->orig_variables));
      if (origID==NO_UINT)
        Error("Undefined ID1 in CSRemoveLCVars");

      /* Translate the 'lcr' linear constraint that refers to the cuiksystem
         in the current state of simplification to the original system */
      m=GetNumTermsInLinearConstraint(&lcr);
      for(i=0;i<m;i++)
        {
          v=GetVariable(GetLinearConstraintVariable(i,&lcr),&(cs->variables));
          name1=GetVariableName(v);
          origID1=GetVariableID(name1,&(cs->orig_variables));
          if (origID1==NO_UINT)
            Error("Undefined ID in CSRemoveLCVars");
          AddTerm2LinearConstraint(origID1,
                                   GetLinearConstraintCoefficient(i,&lcr),
                                   &(lc[origID]));
        }
      GetLinearConstraintError(&error,&lcr);
      if (IntervalSize(&error)>ZERO)
        Error("Linear constraint used for variable replacement must have ct right-hand");
      AddCt2LinearConstraint(-IntervalCenter(&error),&(lc[origID]));

      replaced[origID]=TRUE;

      #if (_DEBUG>5) 
        printf("Var Replacement: %s[%u-%u]->",name,origID,id);
        PrintLinearConstraint(stdout,FALSE,varNamesO,&(lc[origID]));
      #endif

      CopyLinearConstraint(&lcc,&(lc[origID]));
      AddTerm2LinearConstraint(origID,-1.0,&lcc);
      hold=CropLinearConstraint(ZERO,borig,&lcc);
      if (!hold)
        printf("\nNon intersecting ranges\n");
      DeleteLinearConstraint(&lcc);

      /*If origID was used for the replacement of another variable
        we have to propagate the replacement of origID1 to this other
        variable*/
      nv=NVariables(&(cs->orig_variables));
      for(j=0;((hold)&&(j<nv));j++)
        {
          if ((replaced[j])&&(LinearConstraintIncludes(origID,&(lc[j]))))
            {
              TLinearConstraint ltmp;
              double ct;

              ct=RemoveTermFromLinearConstraint(origID,&(lc[j]));
              CopyLinearConstraint(&ltmp,&(lc[origID]));
              ScaleLinearConstraint(ct,&ltmp);
              AddLinearConstraints(&ltmp,&(lc[j]));
              
              CopyLinearConstraint(&lcc,&(lc[j]));
              AddTerm2LinearConstraint(j,-1.0,&lcc);
              hold=CropLinearConstraint(ZERO,borig,&lcc);
              if (!hold) 
                printf("\nNon intersecting ranges\n");
              DeleteLinearConstraint(&lcc);
            }
        }

      if (hold)
        {
          hold=ReplaceVariableInEquations(epsilon,id,&lcr,&(cs->equations));
          RemoveVariable(id,&(cs->variables));

          #if (_DEBUG>6)
            if(hold)
              {
                char **varNames;
                
                printf("%%=======================================================\n");
                printf("%% One less variable (replaced)\n"); 
                printf("%%****************************************\n");
                PrintVariables(stdout,&(cs->variables));
                
                nv=NVariables(&(cs->variables));
                NEW(varNames,nv,char *);
                GetVariableNames(varNames,&(cs->variables));
                PrintEquations(stdout,varNames,&(cs->equations));
                free(varNames);
              }
          #endif
        }

      DeleteLinearConstraint(&lcr);
    }
  #if (_DEBUG>5)
    free(varNamesO);
  #endif
  free(systemVars);
  return(hold);
}

/*
  Simplifies a CuikSystem removing variables have constant value and those
  that depend linearly on other variables. Moreover, a primitive (but numerically
  save) form of Gaussian elimination is performed to remove as more variables
  and equations as possible.
  More elaborate forms of system simplifications are not implemented because they
  result in numerical errors (that can change the system solutions).
  The output mappings give information about how the original and the simplified
  systems are related.

  This is part of the UpdateCuikSystem.
*/
boolean SimplifyCuikSystem(Tparameters *p,TCuikSystem *cs)
{

  unsigned int i,nv,nvn,newVarID;
  boolean *replaced;
  TLinearConstraint *lc,lcS;
  boolean changed,anyChange,hold;
  char *name1;
  double epsilon;
  Tbox borig,bsimp;
  unsigned int nTerms,m,k;
  Tinterval error;
  unsigned int simplificationLevel;

  epsilon=GetParameter(CT_EPSILON,p);
  simplificationLevel=(unsigned int)(GetParameter(CT_SIMPLIFICATION_LEVEL,p));

  #if(_DEBUG>5)
    unsigned int iteration;
  #endif

  /* Get a copy of the original cuiksystem */
  CopyVariables(&(cs->variables),&(cs->orig_variables));
  CopyEquations(&(cs->equations),&(cs->orig_equations));

  nv=NVariables(&(cs->orig_variables)); /* number of variables in the original system */

  NEW(replaced,nv,boolean);
  NEW(lc,nv,TLinearConstraint);
  
  BoxFromVariables(&borig,&(cs->orig_variables)); /*original ranges to be possibly reduced during
                                                    simplification */
  for(i=0;i<nv;i++)
    {
      replaced[i]=FALSE;
      InitLinearConstraint(&(lc[i]));
    }

  hold=TRUE; /* in the simplification process we can discover that the input system
                is inconsitent (i.e., that it does not hold)*/
 
  if (simplificationLevel>0) /* For detailed debug systems are not simplified */
    { 
      #if(_DEBUG>5)
        iteration=0;
      #endif

      CSRemoveUnusedVars(p,cs);
      hold=CSRemoveVarsWithCtRange(p,replaced,lc,&borig,cs);
    
      anyChange=TRUE; /* Set to TRUE to trigger the simplification process */
      while((hold)&&(anyChange))
        {
          anyChange=FALSE; /* Will be set to TRUE if at least one variable is removed
                              in the loop below */
          nTerms=0; /* First we try to remove ct variables, then scaled ones,i.e., we
                       progressivelly increase the number of terms used in the
                       replacements. */
          while((hold)&&(nTerms<=MAX_TERMS_SIMP)) /*hard-coded maximum complexity of the simplifications*/
            {
              hold=CSRemoveLCVars(p,nTerms,&changed,replaced,lc,&borig,cs);
              anyChange=((anyChange)||(changed));
              if (changed)
                nTerms=0;
              else
                nTerms++;
            }

          /* At this point no more ct or scaled varibles can be removed. Try to
             simplify the problem via Gaussian elimination*/
          if ((hold)&&(anyChange))
            {
              anyChange=FALSE;
              do {
                changed=GaussianElimination(&(cs->equations));
                anyChange=((anyChange)||(changed));
              } while (changed);

              #if (_DEBUG>6)
              {
                char **varNames;
                unsigned int nvg;

                printf("%%=======================================================\n");
                printf("%% After gaussian elimination\n"); 
                printf("%%****************************************\n");
                PrintVariables(stdout,&(cs->variables));
              
                nvg=NVariables(&(cs->variables));
                NEW(varNames,nvg,char *);
                GetVariableNames(varNames,&(cs->variables));
                PrintEquations(stdout,varNames,&(cs->equations));
                free(varNames);
              }
              #endif
            }
        } 
    }

  #if(_DEBUG>4)
    if (hold)
      {
        char **varNames;
        unsigned int nvg;

        printf("%%=======================================================\n");
        printf("%%FINAL SIMPLIFIED SYSTEM\n"); 
        printf("%%****************************************\n");
        PrintVariables(stdout,&(cs->variables));
        
        nvg=NVariables(&(cs->variables));
        NEW(varNames,nvg,char *);
        GetVariableNames(varNames,&(cs->variables));
        PrintEquations(stdout,varNames,&(cs->equations));
        free(varNames);
      }
    else
      printf("Inconsistent input system\n");
  #endif

  if (hold)
    {
      /* Get a copy of the simplified but not still dummified cuiksystem */
      CopyVariables(&(cs->simp_variables),&(cs->variables));
      CopyEquations(&(cs->simp_equations),&(cs->equations));

      /*The final step is to dummify the cuiksystem*/
      DummifyCuikSystem(p,cs);

      /*define the mappings from the gathered information*/

      NEW(cs->orig2sd,1,Tmapping);
      NEW(cs->orig2s,1,Tmapping);

      InitMapping(&(cs->orig_variables),&(cs->variables),cs->orig2sd);
      InitMapping(&(cs->orig_variables),&(cs->simp_variables),cs->orig2s);


      /* now we complement the default initialization with the replacements computed above */
      for(i=0;i<nv;i++)
        {
          if (replaced[i])
            {
              m=GetNumTermsInLinearConstraint(&(lc[i]));
              InitLinearConstraint(&lcS);
              for(k=0;k<m;k++)
                {
                  name1=GetVariableName(GetVariable(GetLinearConstraintVariable(k,&(lc[i])),&(cs->orig_variables)));
                  newVarID=GetVariableID(name1,&(cs->variables));
                  
                  AddTerm2LinearConstraint(newVarID,GetLinearConstraintCoefficient(k,&(lc[i])),&lcS);
                }
              GetLinearConstraintError(&error,&(lc[i]));
              AddCt2LinearConstraint(-IntervalCenter(&error),&lcS);
              
              SetOriginalVarRelation(i,&lcS,cs->orig2sd);
              SetOriginalVarRelation(i,&lcS,cs->orig2s);

              DeleteLinearConstraint(&lcS);
            }
        }

      SimpleFromOriginal(&borig,&bsimp,cs->orig2sd);
      nvn=NVariables(&(cs->variables)); 
      for(i=0;i<nvn;i++)
        SetVariableInterval(GetBoxInterval(i,&bsimp),GetVariable(i,&(cs->variables)));
      DeleteBox(&bsimp);
    }

  DeleteBox(&borig);

  /* free the information used for the mapping */
  for(i=0;i<nv;i++)
    DeleteLinearConstraint(&lc[i]);
  free(replaced);
  free(lc);

  return(hold);
}

/*
  Removes the information stored in the cuiksystem during the update
*/
void UnUpdateCuikSystem(TCuikSystem *cs)
{
  
  if (cs->updated)
    {
      /*Removes cached information that is used during solving*/
     
      /*Firt the box sorting information*/
      if (cs->searchMode==MINIMIZATION_SEARCH)
        DeleteEquation(&(cs->eqMin));

      /*Now the information about the simplified+dummified system*/
      if (cs->orig2sd!=NULL)
        {
          DeleteMapping(cs->orig2sd);
          free(cs->orig2sd);
          cs->orig2sd=NULL;
        }

      cs->nequations=0;
      cs->nvariables=0;  
      if (cs->systemVar!=NULL) 
        free(cs->systemVar);
      cs->systemVar=NULL;
      if (cs->notDummyVar!=NULL)
        free(cs->notDummyVar);
      cs->notDummyVar=NULL;
      if (cs->varType!=NULL)
        free(cs->varType);
      cs->varType=NULL;

      DeleteEquations(&(cs->equations));
      DeleteVariables(&(cs->variables));

      /*Now remove the information about the simplified system*/
      if (cs->orig2s!=NULL)
        {
          DeleteMapping(cs->orig2s);
          free(cs->orig2s);
          cs->orig2s=NULL;
        }
      DeleteSimpCuikSystemJacobian(cs);

      cs->simp_nequations=0;
      cs->simp_nvariables=0;;
      cs->simp_nee=0;

      DeleteEquations(&(cs->simp_equations));
      DeleteVariables(&(cs->simp_variables));

      /*Finally remove the information about the original system*/ 
      cs->orig_nequations=0;
      cs->orig_nvariables=0;  
      if (cs->orig_systemVar!=NULL) 
        free(cs->orig_systemVar);
      cs->orig_systemVar=NULL;
      if (cs->orig_notDummyVar!=NULL)
        free(cs->orig_notDummyVar);
      cs->orig_notDummyVar=NULL;

      free(cs->orig_varNames);
      cs->orig_varNames=NULL;

      /*Mark the system as not updated*/
      cs->updated=FALSE;
    }
}

/*
  Simplifies the CuikSystem and caches information that is useful during 
  solving (cached->faster access).
*/
boolean UpdateCuikSystem(Tparameters *p,TCuikSystem *cs)
{
  boolean hold=TRUE;

  if (!cs->updated)
    { 
      unsigned int j;

      /* we update even if we have no equations to be able to deal
         with systems without equations (non constrained problems)
         In this case cuik is a standard approximated cell decomposition 
         approach
       */
      hold=SimplifyCuikSystem(p,cs);
      
      if (hold)
        {
          /* Update the information in the original system */
          cs->orig_nequations=NEquations(&(cs->orig_equations));
          cs->orig_nvariables=NVariables(&(cs->orig_variables));

          if (cs->orig_nvariables==0)
            {
              cs->orig_systemVar=NULL;
              cs->orig_notDummyVar=NULL;
              cs->orig_varNames=NULL;
            }
          else
            {
              NEW(cs->orig_systemVar,cs->orig_nvariables,boolean);
              NEW(cs->orig_notDummyVar,cs->orig_nvariables,boolean);
          
              for(j=0;j<cs->orig_nvariables;j++)
                {
                  cs->orig_systemVar[j]=IsSystemVariable(j,&(cs->orig_variables));
                  cs->orig_notDummyVar[j]=!IsDummyVariable(j,&(cs->orig_variables));
                }

              NEW(cs->orig_varNames,cs->orig_nvariables,char *);
              GetVariableNames(cs->orig_varNames,&(cs->orig_variables));
            }

          /* Update the simplifed cuiksystem */
          cs->simp_nequations=NEquations(&(cs->simp_equations));
          cs->simp_nvariables=NVariables(&(cs->simp_variables));
          cs->simp_nee=NEqualityEquations(&(cs->simp_equations));

          ComputeSimpCuikSystemJacobian(cs);

          /* Update the information in the simplified+dummified cuiksystem*/
          cs->nequations=NEquations(&(cs->equations));
          cs->nvariables=NVariables(&(cs->variables));

          cs->empty=((cs->nvariables==0)||(cs->nequations==0));

          if ((cs->nequations>0)&&(cs->nvariables==0))
            Error("System with equations but without variables !!");

          if (cs->nvariables==0)
            {
              cs->systemVar=NULL;
              cs->notDummyVar=NULL;
              cs->varType=NULL;
            }
          else
            {
              NEW(cs->systemVar,cs->nvariables,boolean);
              NEW(cs->notDummyVar,cs->nvariables,boolean);
              NEW(cs->varType,cs->nvariables,unsigned int);
          
              for(j=0;j<cs->nvariables;j++)
                {
                  cs->systemVar[j]=IsSystemVariable(j,&(cs->variables));
                  cs->notDummyVar[j]=!IsDummyVariable(j,&(cs->variables));
                  cs->varType[j]=GetVariableTypeN(j,&(cs->variables));
                }
            }

          /* Mark the system as updated */
          cs->updated=TRUE;
        }
      
      /* Update the information about the sorting criteria for boxes */
      if (cs->searchMode==MINIMIZATION_SEARCH)
        {
          if (GetNumMonomials(&(cs->orig_eqMin))==0)
            {
              ResetEquation(&(cs->orig_eqMin));
              cs->searchMode=DEPTH_FIRST_SEARCH;
            }
          else
            RewriteEquation(GetParameter(CT_EPSILON,p),
                            cs->orig2sd,&(cs->eqMin),&(cs->orig_eqMin));
        }
    }

  return(hold);
}

/*
  Selects one dimension along which to split box 'b'. The dimension can
  be selected according to the size or according to the error that each
  variable induce in each equation.
  The output is a index in the simplified system
 */
unsigned int ComputeSplitDimInt(Tparameters *p,Tbox *b,TCuikSystem *cs)
{
  
  UpdateCuikSystem(p,cs);

  if (cs->nvariables==0)
    {
      Error("Splitting an empty cuiksystem");
      return(NO_UINT);
    }
  else
    {
      unsigned int i;
      Tinterval *is;
      double *splitDim;
      unsigned int *d;
      unsigned int n;
      double m,s;
      unsigned int splitType;

      splitType=(unsigned int)(GetParameter(CT_SPLIT_TYPE,p));

      NEW(splitDim,cs->nvariables,double);
      NEW(d,cs->nvariables,unsigned int);

      is=GetBoxIntervals(b);

      if ((cs->nequations>0)&&(splitType==1))
        {
          for(i=0;i<cs->nvariables;i++)
            splitDim[i]=0.0;
      
          /*Add the error contribution*/
          for(i=0;i<cs->nequations;i++)
            UpdateSplitWeight(i,splitDim,b,&(cs->equations));

          /*Do not split along unbounded variables*/
          for(i=0;i<cs->nvariables;i++)
            {
              if(IntervalSize(&(is[i]))>=INF)
                splitDim[i]=0.0;
            }
        }
      else
        {
          if (splitType!=2) 
            {
              /*size-based split or error-based split without any equation*/
              for(i=0;i<cs->nvariables;i++)
                {
                  s=IntervalSize(&(is[i]));
                  if (s<INF)
                    splitDim[i]=s;
                  else
                    splitDim[i]=0.0;
                }
            }
          else
            {
              /*splitType==2 -> random selection of the split dimension */ 
              double epsilon;

              epsilon=GetParameter(CT_EPSILON,p);

              for(i=0;i<cs->nvariables;i++)
                {
                  s=IntervalSize(&(is[i]));
                  if ((s>10*epsilon)&&(s<INF))
                    splitDim[i]=1.0;
                  else
                    splitDim[i]=0.0;
                }
              
            }
        }

      m=-1.0;
      n=0;
      for(i=0;i<cs->nvariables;i++)
        {
          if (cs->systemVar[i])
            {
              if (splitDim[i]>m)
                {
                  m=splitDim[i];
                  d[0]=i;
                  n=1;
                }
              else
                {
                  if (splitDim[i]==m)
                    {
                      d[n]=i;
                      n++;
                    }
                }
            }
        
        }
 
      if (n>0) 
        i=d[randomMax(n-1)];
      else
        {
          Warning("There is no way where to bisect the box");
          i=randomMax(cs->nvariables-1);
        }

      free(d);
      free(splitDim);

      return(i);
    }
}

void SaveCSState(char *fname,Tlist *lb,TCuikSystem *cs)
{  
  FILE *f;
  char *tmpName;
  double tm;

  /*We initially save to a temporary file and then we rename it
    This increases the atomiticy of this operation (if the state
    file is corrupted the recovery operation would fail) */
  NEW(tmpName,strlen(fname)+10,char);
  sprintf(tmpName,"%s.tmp",fname);
  
  /*open in binary mode*/
  f=fopen(tmpName,"wb");
  if (!f)
    Error("Could not open file for saving CSState");

  tm=GetElapsedTime(&(cs->st));
  fwrite(&tm,sizeof(double),1,f);
  SaveStatistics(f,&(cs->st));
  SaveListOfBoxes(f,lb);

  fclose(f);

  remove(fname);
  rename(tmpName,fname);
  free(tmpName);
}

void LoadCSState(char *fname,Tlist *lb,TCuikSystem *cs)
{
  FILE *f;
  double tm;
  Titerator it;
  
  f=fopen(fname,"rb");
  if (!f)
    Error("Could not open file for loading CSState");

  fread(&tm,sizeof(double),1,f);
  LoadStatistics(f,&(cs->st));
  LoadListOfBoxes(f,lb);

  fclose(f);
  
  InitIterator(&it,lb);
  First(&it);
  if (cs->nvariables!=GetBoxNIntervals((Tbox *)GetCurrent(&it)))
    Error("The saved .state and the problem dimensionality do not match");

  /*Pretend that the process was started tm seconds before now*/
  SetInitialTime(GetTime(&(cs->st))-tm,&(cs->st));
}


/************************************************************************/
/************************************************************************/
/************************************************************************/

/*
  Warms up the cuiksystem structure
*/
void InitCuikSystem(TCuikSystem *cs)
{
  InitVariables(&(cs->orig_variables));
  InitEquations(&(cs->orig_equations));

  /* cs->variables and cs->equations are initialized
     via Copy when simplifying the system*/ 
  cs->empty=TRUE;

  cs->updated=FALSE;

  cs->nequations=0;
  cs->nvariables=0;

  cs->systemVar=NULL;
  cs->notDummyVar=NULL;
  cs->varType=NULL;

  cs->orig2sd=NULL;

  cs->orig_nequations=0;
  cs->orig_nvariables=0;

  cs->orig_systemVar=NULL;
  cs->orig_notDummyVar=NULL;

  cs->searchMode=DEPTH_FIRST_SEARCH;
  InitEquation(&(cs->orig_eqMin));
}

/*
  Copies one cuiksystem structure into another (that is assumed as
  initialized but empty).
*/
void CopyCuikSystem(TCuikSystem *cs_dst,TCuikSystem *cs_src)
{
  cs_dst->empty=cs_src->empty; 

  CopyEquations(&(cs_dst->orig_equations),&(cs_src->orig_equations));
  CopyVariables(&(cs_dst->orig_variables),&(cs_src->orig_variables));

  CopyStatistics(&(cs_dst->st),&(cs_src->st));

  cs_dst->searchMode=cs_src->searchMode;
  CopyEquation(&(cs_dst->orig_eqMin),&(cs_src->orig_eqMin));

  cs_dst->updated=cs_src->updated; 
  if (cs_dst->updated)
    {
      CopyEquations(&(cs_dst->equations),&(cs_src->equations));
      CopyVariables(&(cs_dst->variables),&(cs_src->variables));

      NEW(cs_dst->orig2sd,1,Tmapping);
      CopyMapping(cs_dst->orig2sd,cs_src->orig2sd);

      cs_dst->nequations=cs_src->nequations;
      cs_dst->nvariables=cs_src->nvariables;
      if (cs_dst->nvariables==0)
        {
          cs_dst->systemVar=NULL;
          cs_dst->notDummyVar=NULL;
          cs_dst->varType=NULL;
        }
      else
        {
          NEW(cs_dst->systemVar,cs_dst->nvariables,boolean);
          memcpy(cs_dst->systemVar,cs_src->systemVar,sizeof(boolean)*cs_dst->nvariables);
          
          NEW(cs_dst->notDummyVar,cs_dst->nvariables,boolean);
          memcpy(cs_dst->notDummyVar,cs_src->notDummyVar,sizeof(boolean)*cs_dst->nvariables);
          
          NEW(cs_dst->varType,cs_dst->nvariables,unsigned int);
          memcpy(cs_dst->varType,cs_src->varType,sizeof(unsigned int)*cs_dst->nvariables);
        }

      cs_dst->orig_nequations=cs_src->orig_nequations;
      cs_dst->orig_nvariables=cs_src->orig_nvariables;
      if (cs_dst->orig_nvariables==0)
        {
          cs_dst->orig_systemVar=NULL;
          cs_dst->orig_notDummyVar=NULL;
        }
      else
        {
          NEW(cs_dst->orig_systemVar,cs_dst->orig_nvariables,boolean);
          memcpy(cs_dst->orig_systemVar,cs_src->orig_systemVar,sizeof(boolean)*cs_dst->orig_nvariables);

          NEW(cs_dst->orig_notDummyVar,cs_dst->orig_nvariables,boolean);
          memcpy(cs_dst->orig_notDummyVar,cs_src->orig_notDummyVar,sizeof(boolean)*cs_dst->orig_nvariables);
        }
        
      if (cs_dst->searchMode==MINIMIZATION_SEARCH)
        CopyEquation(&(cs_dst->eqMin),&(cs_src->eqMin));
    }
  else
    {
      cs_dst->orig2sd=NULL;

      cs_dst->nequations=0;
      cs_dst->nvariables=0;
      cs_dst->systemVar=NULL;  
      cs_dst->notDummyVar=NULL;
      cs_dst->varType=NULL;

      cs_dst->orig_nequations=0;
      cs_dst->orig_nvariables=0;
      cs_dst->orig_systemVar=NULL;  
      cs_dst->orig_notDummyVar=NULL;
    }
}

void CuikSystemMerge(Tparameters *p,TCuikSystem *cs1,TCuikSystem *cs2,TCuikSystem *cs)
{
  unsigned int neq2,nvar1,nvar2;
  unsigned int i;

  InitCuikSystem(cs);

  cs->empty=((cs1->empty)&&(cs2->empty)); 

  cs->updated=FALSE; 

  cs->orig2sd=NULL;
  
  cs->nequations=0;
  cs->nvariables=0;
  cs->systemVar=NULL;  
  cs->notDummyVar=NULL;
  cs->varType=NULL;
  
  cs->orig_nequations=0;
  cs->orig_nvariables=0;
  cs->orig_systemVar=NULL;  
  cs->orig_notDummyVar=NULL;

  /* One of the sets of variables is supposed to be a sub-set of the
     other.
     If nvar1 is larger than the larger set was in cs1 and nothing has
     to be done*/
  CopyVariables(&(cs->orig_variables),&(cs1->orig_variables));
  nvar1=NVariables(&(cs1->orig_variables));
  nvar2=NVariables(&(cs2->orig_variables));
  for(i=nvar1;i<nvar2;i++)
    AddVariable2CS(GetVariable(i,&(cs2->orig_variables)),cs);

  CopyEquations(&(cs->orig_equations),&(cs1->orig_equations));
  neq2=NEquations(&(cs2->orig_equations));
  for(i=0;i<neq2;i++)
    AddEquation2CS(p,GetEquation(i,&(cs2->orig_equations)),cs);

  cs->searchMode=cs1->searchMode;
  CopyEquation(&(cs->orig_eqMin),&(cs1->orig_eqMin));
  AccumulateEquations(&(cs2->orig_eqMin),1,&(cs->orig_eqMin));
}

double EvaluateEqMin(void *b,void *cs)
{
  double v=0;

  if (!((TCuikSystem *)cs)->updated)
    Error("The CuikSystem must be updated before using EvaluateEqMin");

  if (((TCuikSystem *)cs)->searchMode==MINIMIZATION_SEARCH)
    {      
      #if (0)
        double ct;
        unsigned int nv;
        double *p;
        unsigned int i;
        
        nv=NVariables(&(((TCuikSystem *)cs)->variables));
        
        ct=GetEquationValue(&(((TCuikSystem *)cs)->eqMin));
        
        NEW(p,nv,double);
        
        for(i=0;i<nv;i++)
          p[i]=IntervalCenter(GetBoxInterval(i,(Tbox *)b));

        v=EvaluateEquation(p,&(((TCuikSystem *)cs)->eqMin))-ct; 

        free(p);
      #else
        Tinterval ie;
        double ct;

        ct=GetEquationValue(&(((TCuikSystem *)cs)->eqMin));
        
        EvaluateEquationInt(GetBoxIntervals((Tbox *)b),&ie,&(((TCuikSystem *)cs)->eqMin));
        IntervalOffset(&ie,-ct,&ie);

        /*v=LowerLimit(&ie);*/
        /*v=UpperLimit(&ie);*/
        v=IntervalCenter(&ie);
      #endif
    }
  return(v);
}

boolean CmpBoxesEquation(void *b1,void *b2,void *cs)
{
  return(EvaluateEqMin(b1,cs)<EvaluateEqMin(b2,cs));
}

void SetCSSearchMode(unsigned int sm,Tequation *eqMin,TCuikSystem *cs)
{
  switch(sm)
    {
    case DEPTH_FIRST_SEARCH:
      cs->searchMode=DEPTH_FIRST_SEARCH;
      break;
    case BREADTH_FIRST_SEARCH:
      cs->searchMode=BREADTH_FIRST_SEARCH;
      break;
    case MINIMIZATION_SEARCH:
      cs->searchMode=MINIMIZATION_SEARCH;
      DeleteEquation(&(cs->orig_eqMin));
      CopyEquation(&(cs->orig_eqMin),eqMin);
      break;
    default:
      Error("Unkonwn search mode in SetCSSearchMode");
    }
}

void AddTerm2SearchCriterion(double w,unsigned int v,double val,TCuikSystem *cs)
{
  if (w!=0.0)
    {
      Tequation eq;
      Tmonomial m;
      
      InitEquation(&eq);
      SetEquationCmp(EQU,&eq);
      InitMonomial(&m);
      
      AddCt2Monomial(w,&m);
      AddVariable2Monomial(v,2,&m);
      AddMonomial(&m,&eq);
      ResetMonomial(&m);
      
      AddCt2Monomial(-w*2*val,&m);
      AddVariable2Monomial(v,1,&m);
      AddMonomial(&m,&eq);
      ResetMonomial(&m);
      
      AddCt2Monomial(w*val*val,&m);
      AddMonomial(&m,&eq);
      ResetMonomial(&m);
      
      if (cs->searchMode==MINIMIZATION_SEARCH)
        AccumulateEquations(&eq,1,&(cs->orig_eqMin));
      else
        SetCSSearchMode(MINIMIZATION_SEARCH,&eq,cs);
      
      DeleteMonomial(&m);
      DeleteEquation(&eq);
    }
}

/*
  Adds a equation to the CuikSystem. 
  Equations as added by the user are not dummified. This helps
  when simplifying the system (dummify variables can prevent some
  simplifications or make them more "obscure")
  The dummificitaion is applied after the simplification of the
  system (see SimplifyCuikSystem)
*/
void AddEquation2CS(Tparameters *p,Tequation *eq,TCuikSystem *cs)
{
  if (cs->updated)
    UnUpdateCuikSystem(cs);

  if (eq!=NULL)
    AddEquation(eq,&(cs->orig_equations));
}

/*
  Checks if the system has a variable with the given name.
  If so, we return the variable identifier.
  If not, we create a new variable
 */
unsigned int AddVariable2CS(Tvariable *v,TCuikSystem *cs)
{
  unsigned int id;

  if (cs->updated)
    UnUpdateCuikSystem(cs);

  id=GetVariableID(GetVariableName(v),&(cs->orig_variables));

  if (id==NO_UINT)
    id=AddVariable(v,&(cs->orig_variables));

  return(id);
}

/*
  Returns a copy of the variables stored in the CuikSystem
*/
void GetCSVariables(Tvariables *vs,TCuikSystem *cs)
{
  CopyVariables(vs,&(cs->orig_variables));
}

/*
  Returns the number of variables in the CuikSystem
*/
unsigned int GetCSNumVariables(TCuikSystem *cs)
{
  return(NVariables(&(cs->orig_variables)));
}

/*
  Returns the number of system variables in the CuikSystem
*/
unsigned int GetCSNumNonDummyVariables(TCuikSystem *cs)
{
  return(NVariables(&(cs->orig_variables))-GetNumDummyVariables(&(cs->orig_variables)));
}

/* 
   Gets a copy of variable with ID 'n'
*/
void GetCSVariable(unsigned int n,Tvariable *v,TCuikSystem *cs)
{
  CopyVariable(v,GetVariable(n,&(cs->orig_variables)));
}

/*
  Sets a new range for variable with ID 'n'
*/
void SetCSVariableRange(unsigned int n,Tinterval *r,TCuikSystem *cs)
{ 
  if (cs->updated)
    UnUpdateCuikSystem(cs); /* Variables with ct range are simplified in
                               the UpdateCuikSystem */

  SetVariableInterval(r,GetVariable(n,&(cs->orig_variables)));
}

/*
  Gets the ID of the variable with the given name
*/
unsigned int GetCSVariableID(char *name,TCuikSystem *cs)
{
  return(GetVariableID(name,&(cs->orig_variables)));
}

/*
  Returns an array of booleans with sv[i]=TRUE if variable
  'i' is a system variable.
  Also returns the number of variables in the CuikSystem
*/
unsigned int GetCSSystemVars(boolean **sv,TCuikSystem *cs)
{
  unsigned int i,n;

  n=NVariables(&(cs->orig_variables));
  NEW(*sv,n,boolean);

  for(i=0;i<n;i++)
    (*sv)[i]=((IsSystemVariable(i,&(cs->orig_variables)))||
              (IsSecondaryVariable(i,&(cs->orig_variables))));
  return(n);
}

/*
  Returns a copy of the equations in the cuiksystem
*/
void GetCSEquations(Tequations *eqs,TCuikSystem *cs)
{
  CopyEquations(eqs,&(cs->orig_equations));
}

/*
  Returns a copy of equation 'n' in the CuikSystem
*/
void GetCSEquation(unsigned int n,Tequation *eq,TCuikSystem *cs)
{
  CopyEquation(eq,GetEquation(n,&(cs->orig_equations)));
}

/*
  Returns the number of equations in the CuikSytem
*/
unsigned int GetCSNumEquations(TCuikSystem *cs)
{
  return(NEquations(&(cs->orig_equations)));
}

void GetCuikSystemJacobian(Tequation ***J,TCuikSystem *cs)
{
  unsigned int i,j,neqs,nvars;

  neqs=NEquations(&(cs->orig_equations));
  nvars=NVariables(&(cs->orig_variables));

  NEW((*J),neqs,Tequation *);
  for(i=0;i<neqs;i++)
    {
      NEW((*J)[i],nvars,Tequation);
      for(j=0;j<nvars;j++)
        DeriveEquation(j,&((*J)[i][j]),GetEquation(i,&(cs->orig_equations)));
    }
}

void DeleteCuikSystemJacobian(Tequation **J,TCuikSystem *cs)
{
  unsigned int i,j,neqs,nvars;
  
  neqs=NEquations(&(cs->orig_equations));
  nvars=NVariables(&(cs->orig_variables));

  for(i=0;i<neqs;i++)
    {
      for(j=0;j<nvars;j++)
        DeleteEquation(&(J[i][j]));
      free(J[i]);
    }
  free(J);
}


void ComputeSimpCuikSystemJacobian(TCuikSystem *cs)
{
  unsigned int i,j,k;
  Tequation *eq;

  NEW(cs->J,cs->simp_nee,Tequation *);
  k=0;
  for(i=0;i<cs->simp_nee;i++)
    {
      eq=GetEquation(i,&(cs->simp_equations));
      if (GetEquationCmp(eq)==EQU)
        {
          NEW(cs->J[k],cs->simp_nvariables,Tequation);
          for(j=0;j<cs->simp_nvariables;j++)
            DeriveEquation(j,&(cs->J[k][j]),eq);
          k++;
        }
    }
}

void DeleteSimpCuikSystemJacobian(TCuikSystem *cs)
{
  unsigned int i,j;
  
  if (!cs->updated)
    Error("DeleteSimpCuikSystemJacobian can only be used on updated cuiksystems");

  for(i=0;i<cs->simp_nee;i++)
    {
      for(j=0;j<cs->simp_nvariables;j++)
        DeleteEquation(&(cs->J[i][j]));
      free(cs->J[i]);
    }
  free(cs->J);
}

/*
  Reduces a box as much as possible (according to parameters in 'p').
  If the set of parameters is undefined, we use a default set.
  This procedure is used in CuikPlan.
  
  This is basically another flavor of ReduceBox where parameters
  are given in one structure instead of as different inputs

  The input box and output boxes are refered to the original
  cuiksystem (the non-simplified one)
  
  The output can be EMPTY_BOX
                    REDUCED_BOX
                    REDUCED_BOX_WITH_SOLUTION
                    ERROR_IN_PROCESS
  
*/
unsigned int MaxReduction(Tparameters *p,unsigned int varMask,Tbox *b,TCuikSystem *cs)
{
  unsigned int e;
  
  if (UpdateCuikSystem(p,cs))
    {
      Tbox bsimp;
      
      SimpleFromOriginal(b,&bsimp,cs->orig2sd);
      e=ReduceBox(p,(varMask==0?~DUMMY_VAR:varMask),&bsimp,cs);
      if (e!=EMPTY_BOX)
        {
          /* ERROR_IN_PROCESS, REDUCED_BOX, REDUCED_BOX_WITH_SOLUTION*/
          UpdateOriginalFromSimple(&bsimp,b,cs->orig2sd);
        }
      DeleteBox(&bsimp);
    }
  else
    e=EMPTY_BOX;

  return(e);
}


boolean SampleCuikSystem(Tparameters *p,char *fname,Tlist *sb,
                         unsigned int nsamples,unsigned int ntries,
                         unsigned int ndof,TCuikSystem *cs)
{
  Tbox init_box;
  boolean haveSample;

  GenerateInitialBox(&init_box,cs);

  haveSample=SampleCuikSystemInBox(p,fname,sb,nsamples,ntries,ndof,&init_box,cs);
  
  DeleteBox(&init_box);
  
  return(haveSample);
}

boolean SampleCuikSystemInBox(Tparameters *p,char *fname,Tlist *sb,
                              unsigned int nsamples,unsigned int ntries,
                              unsigned int ndof,
                              Tbox *init_box,TCuikSystem *cs)
{
  Tbox orig_ranges,b,*bAux,*box;
  Tlist solutions;
  unsigned int k,d,ns,nv,ne,i,nt;
  double v;
      
  Tfilename fsols_box;
  Tfilename fsols_point;

  FILE *f_out_box;  
  FILE *f_out_point;  

  Titerator it;
  boolean *systemVars;
  Tinterval *r;
  unsigned int spt;
  double ms;

  boolean ok;
  boolean end;

  /* free var = variable not in equations -> must be fixed always */
  unsigned int *freeVars;
  unsigned int nFreeVars;

  /*Get a copy of the ranges in 'cs' for all the variables*/
  GenerateInitialBox(&orig_ranges,cs);
  
  /*Set the ranges for the 'cs' variables to those in the given box (init_box)*/
  if (cs->updated) 
    UnUpdateCuikSystem(cs);

  nv=NVariables(&(cs->orig_variables));
  for(k=0;k<nv;k++)
    SetVariableInterval(GetBoxInterval(k,init_box),GetVariable(k,&(cs->orig_variables)));

  NEW(freeVars,nv,unsigned int);
  nFreeVars=0;
  for(k=0;k<nv;k++)
    {
      if ((GetVariableTypeN(k,&(cs->orig_variables))==SYSTEM_VAR)&&
          (!UsedVarInEquations(k,&(cs->orig_equations))))
        {
          freeVars[nFreeVars]=k;
          nFreeVars++;
        }
    }

  /*Change the SPLIT_TYPE parameter to random split*/
  spt=(unsigned int)GetParameter(CT_SPLIT_TYPE,p);
  ChangeParameter(CT_SPLIT_TYPE,2,p); /* select split dimentions at random */
  /*Change the SMALL_SIGMA paramter to ensure it is small enough.
    SMALL_SIGMA controls the precission of the isolated solutions and, thus, of the
    sampled points. */
  ms=GetParameter(CT_SMALL_SIGMA,p);
  ChangeParameter(CT_SMALL_SIGMA,100*GetParameter(CT_EPSILON,p),p);

  if (fname==NULL)
    {
      f_out_box=NULL;
      f_out_point=NULL;
    }
  else
    {
      /* Check if we can create the output files */
      CreateFileName(NULL,fname,"_samples",SOL_EXT,&fsols_box);
      f_out_box=fopen(GetFileFullName(&fsols_box),"w"); 
      if (!f_out_box)
        Error("Could not open box sample output file");
      
      CreateFileName(NULL,fname,NULL,SAMPLE_EXT,&fsols_point);
      f_out_point=fopen(GetFileFullName(&fsols_point),"w"); 
      if (!f_out_point)
        Error("Could not open box sample output file");
    }

  if (sb!=NULL)
    InitListOfBoxes(sb);

  nv=GetCSNumVariables(cs);
  ne=GetCSNumEquations(cs);

  if (ndof==NO_UINT)
    {
      /* try to get the number of dof from the system (vars-equations) */
      if (ne>=nv)
        Error("Over/Well-constrained system. Please indicate the number of dof in the call");
      ndof=nv-ne;
    }

  if (MaxReduction(p,~DUMMY_VAR,init_box,cs)==EMPTY_BOX)
    end=TRUE;
  else
    end=FALSE;

  if (nFreeVars>ndof)
    Error("More free vars than degrees of freedom");

  ns=0; /*samples already determined*/
  nt=0;
  while((!end)&&(ns<nsamples))
    {
      do {
        ok=TRUE;

        /*fix 'ndof' varibles*/
        CopyBox(&b,init_box);
        for(k=0;k<nv;k++)
          SetCSVariableRange(k,GetBoxInterval(k,&b),cs);
        
        for(k=0;k<nFreeVars;k++)
          {     
            r=GetBoxInterval(freeVars[k],&b);
            v=randomInInterval(r);
            
            #if (_DEBUG>1)
              printf("[CS]->fixing %u to %f\n",freeVars[k],v);
            #endif
            
            NewInterval(v,v,r);
            
            SetCSVariableRange(freeVars[k],r,cs);
          }

        for(k=nFreeVars;((ok)&&(k<ndof));k++)
          {
            /*fix 'ndof' varibles*/
            d=ComputeSplitDim(p,&b,cs);

            if (d==NO_UINT)
              ok=FALSE;
            else
              {
                r=GetBoxInterval(d,&b);
                v=randomInInterval(r);
                
                #if (_DEBUG>1)
                  printf("[CS]->fixing %u to %f\n",d,v);
                #endif

                NewInterval(v,v,r);
                
                SetCSVariableRange(d,r,cs);
              }
          }
        if (ok)
          {           
            #if (_DEBUG>1)
              printf("[CS]");PrintBox(stdout,&b);
            #endif
          }
        else
          {
            #if (_DEBUG>1)
              printf("Inconsistent system while sampling\n");
            #endif
          }
      
        DeleteBox(&b);
     
        nt++;
        if ((ntries!=NO_UINT)&&(nt==ntries))
          end=TRUE;
      } while((!ok)&&(!end));

      if (!end)
        {
          
          InitListOfBoxes(&solutions);

          /*See if there are solutions for the fixed ranges*/
          SolveCuikSystem(p,nsamples-ns,FALSE,NULL,f_out_box,&solutions,cs);
          
          /*We have to print the samples*/
          GetCSSystemVars(&systemVars,cs);
          InitIterator(&it,&solutions);
          First(&it);
          while((!EndOfList(&it))&&(ns<nsamples))
            {
              bAux=(Tbox *)GetCurrent(&it);

              if (f_out_point!=NULL)
                {
                  /* Print the center of the box into the point file */
                  /* We only print values for the system vars */
                  for(i=0;i<nv;i++)
                    {
                      if (systemVars[i])
                        fprintf(f_out_point," %.12g",IntervalCenter(GetBoxInterval(i,bAux)));
                    }
                  fprintf(f_out_point,"\n");
                }
              if (sb!=NULL)
                {
                  NEW(box,1,Tbox); 
            
                  CopyBox(box,bAux);
            
                  AddLastElement((void *)box,sb);
                }
        
              ns++;
              Advance(&it);
            }
          free(systemVars);

          fflush(f_out_point);
          fflush(f_out_box);

          /*Remove all boxes from the list and the current box*/
          DeleteListOfBoxes(&solutions);
        }
    } 

  if (f_out_box!=NULL)
    {
      DeleteFileName(&fsols_box);
      fclose(f_out_box);
    }
  if (f_out_point!=NULL)
    {
      fclose(f_out_point);
      DeleteFileName(&fsols_point);
    }

  /*Restore the user given split_type and min_sigma*/
  ChangeParameter(CT_SPLIT_TYPE,spt,p);
  ChangeParameter(CT_SMALL_SIGMA,ms,p);

  /*Return the original ranges to the 'cs' variables*/
  if (cs->updated) 
    UnUpdateCuikSystem(cs);

  nv=NVariables(&(cs->orig_variables));
  for(k=0;k<nv;k++)
    SetVariableInterval(GetBoxInterval(k,&orig_ranges),GetVariable(k,&(cs->orig_variables)));

  free(freeVars);

  DeleteBox(&orig_ranges);
  
  return(!end);
}

/*
  Takes a cuik systema and returns all solutions. 
  Solutions are stored in file 'f_out' (if defined) and in list 'sol'
  (if defined), or both.

  This procedure is used for single-cpu versions of Cuik.
 */
void SolveCuikSystem(Tparameters *p,unsigned int nsols,
                     boolean restart,char *fstate,
                     FILE *f_out,Tlist *sol,TCuikSystem *cs)
{
  Theap boxes;
  #if (_DEBUG>0)
    double box_volume;
  #endif
  #if (_DEBUG>1)
    unsigned int nbox=0;
  #endif 

  Tbox b;
  unsigned int c; /* output of ReduceBox */
  unsigned int current_level; /* level of the current box in the search three */
  Tbox init_box;
  boolean done;
  unsigned int nb; /*boxes for recovery mode*/
  unsigned int statePeriod;

  /************************************************************************************/
  /************************************************************************************/
  /************************************************************************************/

  if (UpdateCuikSystem(p,cs))
    {
      
      switch(cs->searchMode)
        {
        case DEPTH_FIRST_SEARCH:
          InitHeapOfBoxes(CmpBoxDepthFirst,NULL,&boxes);
          break;
        case BREADTH_FIRST_SEARCH:
          InitHeapOfBoxes(CmpBoxBreadthFirst,NULL,&boxes);
          break;
        case MINIMIZATION_SEARCH:
          InitHeapOfBoxes(CmpBoxesEquation,(void *)cs,&boxes);
          break;
        }

      if ((restart)&&(fstate!=NULL))
        {
          Tlist pboxes;

          LoadCSState(fstate,&pboxes,cs);
          AddList2Heap(&pboxes,&boxes);
          DeleteListOfBoxes(&pboxes);
        }
      else
        {
          BoxFromVariables(&init_box,&(cs->variables));      
          AddBox2HeapOfBoxes(&init_box,&boxes);
          DeleteBox(&init_box);

          InitStatistics(1,HeapOfBoxesVolume(cs->systemVar,&boxes),&(cs->st));
        }

      current_level=0;
      nb=0;
      statePeriod=GetParameter(CT_STATE_PERIOD,p);

      done=FALSE;

      /*the cuik solve process itself*/
      do {
        if(!HeapEmpty(&boxes))
          {
            /*Get a new box*/

            if ((statePeriod>0)&&(fstate!=NULL))
              {
                nb++;
                if (nb==statePeriod)
                  {
                    Tlist pboxes;

                    Heap2List(&pboxes,&boxes);
                    SaveCSState(fstate,&pboxes,cs);
                    DeleteListOfBoxes(&pboxes);
                    nb=0;
                  }
              }

            NewBoxProcessed(&(cs->st));

            /*The box is extracted to avoid another process 
              to get the same box*/
            ExtractMinElement(&b,&boxes); 

            current_level=GetBoxLevel(&b);

            NewMaxLevel(current_level,&(cs->st));

            #if (_DEBUG>0)
              fprintf(stderr,"<%g %g>[%u]",
                      GetBoxVolume(cs->systemVar,&b),
                      GetBoxSize(cs->systemVar,&b),
                      current_level);
              box_volume=GetBoxVolume(cs->systemVar,&b);
            #endif

            #if (_DEBUG>1)
              printf("\n\n%u: Analizing box: ",nbox);
              nbox++;
              PrintBox(stdout,&b);
              printf("      Box level   : %u\n",GetBoxLevel(&b));
              printf("      Box volume  : %g\n",GetBoxVolume(cs->systemVar,&b));
              printf("      Box size    : %g\n",GetBoxSize(cs->systemVar,&b));
              if (cs->searchMode==MINIMIZATION_SEARCH)
                printf("      Box penalty : %g\n",EvaluateEqMin((void *)&b,(void *)cs));
            #endif

            /* Reduce the box by all but dummy vars */;
            c=ReduceBox(p,~DUMMY_VAR,&b,cs);
        
            #if (_DEBUG>0)
              fprintf(stderr,"  -> ");
            #endif  

            PostProcessBox(p,c,f_out,sol,&boxes,&b,cs);

            if (nsols>0)
              done=(GetNSolutionBoxes(&(cs->st))==nsols);

            DeleteBox(&b);
          }
      } while((!done)&&(!HeapEmpty(&boxes))); /*Until everything is completed*/

      DeleteHeap(&boxes); /*for sanity we delete the empty list*/

      PrintStatistics((f_out==NULL?stdout:f_out),&(cs->st));
    }
}


#if (_USE_MPI)
/*
  In multi-cpu versions we use a MPI version of CuikSolve.
  One of the CPUs takes the role of scheduler in charge of managing the
  lists of boxes pending to be processed and classifing boxes out of
  reduction process as solutions of boxes to be split.
  The rest of CPUs take the charge of reducing boxes (i.e., of executing
  ReduceBox to the boxes received from the scheduler) and of returning
  the resulting boxes (with the associated result code) to the
  scheduler.

  The scheduler has to take into account that some child nodes can
  'die' while reducing boxes. So it has to keep track of when a box was
  sent to a child and prediodically check whether a timeout was been reached.
  If so, the corresponding processor is considered as 'dead' the
  box sent is recovered and added to the list of boxes to be processed
  (i.e., send to other child nodes still alive).

  The scheduler should considere the extreme case where all child process die
  and should take care of "awakening" dead processors.
  This is not implemented.... till now.
*/
void MPI_SolveCuikSystem(Tparameters *p,
                         boolean restart,char *fstate,
                         FILE *f_out,TCuikSystem *cs)
{
  Theap boxes;  /*queued boxes*/
  Tbox b;
  unsigned int c;
  unsigned int nr;
  unsigned int current_level;

  /*Variables to set up the parallelism*/
  unsigned int in_process; /* number of boxes send to child processors and still
                              in process */
  unsigned int available; /* number of child processors available */
  signed int np; /* number of processors used (np-1) is the number of childs (process
                    0 is the scheduler) */
  MPI_Status status; /* Output status of a MPI command */
  MPI_Request *request; /* Communication request. When the scheduller sends out a box
                           is starts a communication request */
  unsigned int buffer_size; /* Size of a buffer of doubles where we store the 
                               information of a box */
  double **buffers_out; /* boxes send to the child processors */
  double **buffers_in;  /* boxes returned by the child processors */
  signed int i,completed; /* MPI_Testany returns completed=1 if any of the pending 
                             communications requests has been replied (then, 'i' 
                             is the number of the compleated request)*/
  boolean *lost_packet; /* true if a box was never returned */
  time_t *time_stamp; /* time when we send out a box */
  time_t deadline; /* deadline (in seconds) used to compute the maximum expected
                      return time for a box*/
  Tbox init_box;  /* initial search space */
  unsigned int nb;
  unsigned int statePeriod;

  /*lowest possible priority but not as low as that of the children processors*/
  /*setpriority(PRIO_PROCESS,0,PRIO_MAX-1); */

  MPI_Comm_size(MPI_COMM_WORLD,&np); /*total number of processes (take into 
                                       account that this routine is process number 0)*/
  
  if (UpdateCuikSystem(p,cs))
    {
      switch(cs->searchMode)
        {
        case DEPTH_FIRST_SEARCH:
          InitHeapOfBoxes(CmpBoxDepthFirst,NULL,&boxes);
          break;
        case BREADTH_FIRST_SEARCH:
          InitHeapOfBoxes(CmpBoxBreadthFirst,NULL,&boxes);
          break;
        case MINIMIZATION_SEARCH:
          InitHeapOfBoxes(CmpBoxesEquation,(void *)cs,&boxes);
          break;
        }

      BoxFromVariables(&init_box,&(cs->variables));
      
      /* Reserve space for the boxes sent/received and associated information */
      /*All boxes are the same size as the initial one*/        
      buffer_size=GetBoxBufferSize(&init_box);

      if ((restart)&&(fstate!=NULL))
        {
          Tlist pboxes;

          LoadCSState(fstate,&pboxes,cs);
          AddList2Heap(&pboxes,&boxes);
          DeleteListOfBoxes(&pboxes);
        }
      else
        {
          AddBox2HeapOfBoxes(&init_box,&boxes);
          InitStatistics(np,HeapOfBoxesVolume(cs->systemVar,&boxes),&(cs->st));
        }
      DeleteBox(&init_box);

      current_level=0;

      NEW(buffers_out,np,double*);
      NEW(buffers_in,np,double*);
    
      NEW(request,np,MPI_Request);
      NEW(lost_packet,np,boolean);
      NEW(time_stamp,np,time_t);

      for(i=1;i<np;i++)
        {
          request[i]=MPI_REQUEST_NULL;
          NEW(buffers_out[i],buffer_size,double);
          NEW(buffers_in[i],buffer_size,double);
          lost_packet[i]=FALSE;
        }

      request[0]=MPI_REQUEST_NULL; /* Processor 0 is the scheduler, do not use it as a child */
      in_process=0;
      available=(np-1);
        
      nb=0;
      statePeriod=GetParameter(CT_STATE_PERIOD,p);

      /*the cuik solve process itself*/
      do {

        if ((!HeapEmpty(&boxes))&& /* there is something to be processes */
            (available>0)) /* and someone ready to process it */
          {
            /* Search for a free child (at least one must exists) */
            i=1;
            while ((i<np)&&(request[i]!=MPI_REQUEST_NULL)) i++;
                  
            if (i==np)
              Error("wrong number of available processors");
            else
              {
                if ((statePeriod>0)&&(fstate!=NULL))
                  {
                    nb++;
                    if (nb==statePeriod)
                      {
                        Tbox btmp;
                        unsigned int ctmp,ntmp;
                        Tlist pboxes;
                        unsigned int kk;

                        Heap2List(&pboxes,&boxes);
                        for(kk=1;kk<np;kk++)
                          {
                            if (request[kk]!=MPI_REQUEST_NULL)
                              { 
                                InitBox(cs->nvariables,NULL,&btmp);
                                Buffer2Box(&ctmp,&ntmp,buffers_out[kk],&btmp);
                                AddFirstElement((void *)&btmp,&pboxes);
                              }
                          }
                        SaveCSState(fstate,&pboxes,cs);
                        DeleteListOfBoxes(&pboxes);
                        nb=0;
                      }
                  }

                /*Get a new box*/
                ExtractMinElement(&b,&boxes); 
            
                /* Transform the box into a buffer of doubles */
                Box2Buffer(REDUCED_BOX,0,buffers_out[i],&b);
            
                /* Send it to the iddle child processors */
                /* and start waiting for the reply (i.e., the reduced box) */
                if (MPI_Send(buffers_out[i],buffer_size,MPI_DOUBLE,i,20,MPI_COMM_WORLD)==MPI_SUCCESS)
                  {
                    /* The box is sent and now we start waiting for it */
                    if (MPI_Irecv(buffers_in[i],buffer_size,MPI_DOUBLE,i,20,MPI_COMM_WORLD,&(request[i]))==MPI_SUCCESS)
                      {

                        current_level=GetBoxLevel(&b);
                        
                        NewBoxProcessed(&(cs->st));
                        NewMaxLevel(current_level,&(cs->st));
                        
                        /* One more box is in process */
                        in_process++;
                        
                        /* One less processor is available */
                        available--;
                        
                        #if (_DEBUG>0)
                          /* Inform about the box launchig */
                          fprintf(stderr,"b<v:%g, s:%g, l:%u>-> %u\n",
                                  GetBoxVolume(cs->systemVar,&b),
                                  GetBoxSize(cs->systemVar,&b),
                                  current_level,i);
                          /*
                          fprintf(stderr,"   A:%d   P:%d  Q:%u\n", 
                                  available,in_process,HeapSize(&boxes));
                          */
                        #endif  
                        
                        /* Store the time at which we sent out the box */
                        time_stamp[i]=time(NULL);
                        /* and, mark the sent box as not lost */
                        lost_packet[i]=FALSE;
                        
                        #if (_DEBUG>1)
                          printf("Sending box to processors : %u (iddle -> %u box_in_list-> %u)\n",
                                 i,np-in_process-1,HeapSize(&boxes));
                          printf("  Box: ");
                          PrintBox(stdout,&b);
                        #endif
                        
                        DeleteBox(&b);
                      }
                    else
                      {
                        /*The box will never be back: add it to the list and free the processor */
                        #if (_DEBUG>0)
                          fprintf(stderr,"RF-> %u\n",i);
                        #endif
                        
                        AddBox2HeapOfBoxes(&b,&boxes);
                        request[i]=MPI_REQUEST_NULL;
                      }
                  }
                else
                  {
                    /* We didn't manage to send out the box, just return it to the boxes
                       to be processes*/ 
                    #if (_DEBUG>0)
                      fprintf(stderr,"SF-> %u\n",i);
                    #endif

                    AddBox2HeapOfBoxes(&b,&boxes);
                  }
              }
          }

        /* see if we already got any reply from the slaves :-) */
        if ((MPI_Testany(np,request,&i,&completed,&status)==MPI_SUCCESS)&&
            (completed))
          {
            /* if 'completed' 'i' is the request completed */
            #if (_DEBUG>1)
              printf("Receiving box from processors : %u\n",i);
            #endif

            available++; /* either if lost or not, the processor becomes available from now  on */
            request[i]=MPI_REQUEST_NULL; 

            if (!lost_packet[i])
              {
                /* If so, reconstruct the output box and analyze it */
                InitBox(cs->nvariables,NULL,&b);

                Buffer2Box(&c,&nr,buffers_in[i],&b);

                AddNBoxReductions(nr,&(cs->st));

                in_process--; /* lost ones are already discounted when classified as lost 
                                 (to be able to finish even if lost ones exits)*/

                #if (_DEBUG>0)
                  fprintf(stderr,"  %u->b<v:%g s:%g l:%u>(t=%lu): ",
                          i,
                          GetBoxVolume(cs->systemVar,&b),
                          GetBoxSize(cs->systemVar,&b),
                          GetBoxLevel(&b),
                          time(NULL)-time_stamp[i]);
                  /*
                    fprintf(stderr,"   A:%d   P:%d  Q:%u\n", 
                    available,in_process,HeapSize(&boxes)); 
                  */
                #endif

                PostProcessBox(p,c,f_out,NULL,&boxes,&b,cs);

                DeleteBox(&b);
              }
            #if (_DEBUG>0)
            else
              {
                /* We receive a box that was thought to be lost. The only think we
                   do is to mark the child as available again by setting request[i]=MPI_REQUEST_NULL;
                   below */
                fprintf(stderr,"  lb[%u](t=%lu > %u)\n",i,
                        time(NULL)-time_stamp[i],MPI_TREAT_BOX_TIMEOUT(cs));
              }
            #endif
          }


        /* See if one of the send boxes is lost.
           All boxes sent before the deadline are considered lost.
           Lost boxes are recovered and added to the list of pending boxes.
           The processor to which the box was assigned is discarted for future
           jobs.
        */
        deadline=time(NULL)-(time_t)MPI_TREAT_BOX_TIMEOUT(cs);
        for(i=1;i<np;i++)
          {
            if ((!lost_packet[i])&& /*not already lost*/
                (request[i]!=MPI_REQUEST_NULL)&& /*and we are waiting for a reply*/
                (time_stamp[i]<deadline)) /*for too long actually*/
              {
                /* we considere the sub-task as lost, we add the box for further process by 
                   another computer and we mark the current message as lost */

                lost_packet[i]=TRUE; /* mark the packet as lost */

                InitBox(cs->nvariables,NULL,&b); /*Recover it and put it again in the list of 
                                                   boxes to be processed*/
                Buffer2Box(&c,&nr,buffers_out[i],&b); /* Error code (c) and additional information
                                                         (nr) are empty (this is a box whose process 
                                                         was not concluded)*/
            
                NewLostBox(&(cs->st));
                #if (_DEBUG>0)
                  fprintf(stderr,"  l[%u] (elapsed %lu)\n",i,time(NULL)-time_stamp[i]);
                #endif

                /* a time out tipically means a error in the simplex */
                PostProcessBox(p,ERROR_IN_PROCESS,f_out,NULL,&boxes,&b,cs);
            
                DeleteBox(&b);

                in_process--; /* one less box is in process (Observe that the processors is
                                 not marked as available) */
              }
          }

        fflush(stderr);

      } while((!HeapEmpty(&boxes))||(in_process>0)); /*Until everything is completed (nothing in the
                                                       is to be processed and all boxes returned
                                                       from children)*/

      /* Send a magic packed to the child so that they kill themselves */
      
      for(i=1;i<np;i++)
        {
          if (request[i]!=MPI_REQUEST_NULL)
            {
              /*we have a pending communication with this processor -> cancel it*/
              MPI_Cancel(&(request[i]));
            }
          else
            { 
              /*nothing pending with this processors, just order it to die*/
              buffers_out[i][0]=-1.0; 
              MPI_Send(buffers_out[i],buffer_size,MPI_DOUBLE,i,20,MPI_COMM_WORLD); 
            }
        }

      free(request);
      free(lost_packet);
      free(time_stamp);

      DeleteHeap(&boxes);

      if (f_out!=NULL)
        {
          fprintf(f_out,"\n\nCuik executed in %u (%u) child processors\n",available,np-1);
          PrintStatistics(f_out,&(cs->st));
        }

      DeleteStatistics(&(cs->st));

      for(i=1;i<np;i++)
        {
          free(buffers_in[i]);
          free(buffers_out[i]);
        }
        
      free(buffers_in);
      free(buffers_out);
    }
}

/*
  This is the process executed by each child process
*/
void MPI_TreatBox(Tparameters *p,TCuikSystem *cs)
{
  if (UpdateCuikSystem(p,cs))
    {
      boolean end;
      unsigned int n;
      double *buffer;
      unsigned int c,nr;
      Tbox b;
      MPI_Status status;

      InitBox(cs->nvariables,NULL,&b);
      n=GetBoxBufferSize(&b);
      NEW(buffer,n,double);

      /* A child can be in execution in computers not only devoted to cuik.
         We reduce the priority of the cuik process to allow other processes
         (those own by the owners of the machines) to run smoothly. 
         This way we can enlarge the cluster of procesors to almost all
         computers in our Lab.
         In the future we have to develop programs to check which computers
         are on and to add them into the machine list (i.e., the members of
         the cluster)
      */
      /*setpriority(PRIO_PROCESS,0,PRIO_MAX); */

      end=FALSE; /* While the kill buffer is not received */
         
      while(!end)
        {
          /* Get new box from the master */
          if (MPI_Recv(buffer,n,MPI_DOUBLE,0,20,MPI_COMM_WORLD,&status)==MPI_SUCCESS)
            {
              /* The first double is used to send control messages between the 
                 scheduler and the slave. If negative, the slave dies*/
              if (buffer[0]<0) 
                end=TRUE;
              else
                {
                  /* If not the kill buffer, recover the box out of the buffer */
                  Buffer2Box(&c,&nr,buffer,&b);

                  ResetNBoxReductions(&(cs->st)); /* Let's count the box reductions for this 
                                                     box shrink only*/

                  #if (_DEBUG>1)
                    printf("Box from main processor:");
                    PrintBox(stdout,&b);
                    fflush(stdout);
                  #endif

                  /* And reduce it */
                  c=ReduceBox(p,~DUMMY_VAR,&b,cs);

                  #if (_DEBUG>1)
                    printf("Box out of ReduceBox (%u):",c);
                    PrintBox(stdout,&b);
                    fflush(stdout);
                  #endif

                  /* Then pack the result into another buffer */
                  Box2Buffer(c,GetNBoxReductions(&(cs->st)),buffer,&b);

                  /* and send it back to the master */
              
                  #if (_DEBUG>1)
                    if (MPI_Send(buffer,n,MPI_DOUBLE,0,20,MPI_COMM_WORLD)!=MPI_SUCCESS)
                      printf("Sending to master failed\n");
                  #else
                    MPI_Send(buffer,n,MPI_DOUBLE,0,20,MPI_COMM_WORLD);
                  #endif
                
                }
            }
          #if (_DEBUG>1)
          else
            printf("Receive from master failed\n");
          #endif
        }

      DeleteBox(&b);
      free(buffer);
    }
}
#endif



/*
  Returns a box defined by the ranges for the variables given by the user in the
  input file. This is the initial search space.
 */
void GenerateInitialBox(Tbox *box,TCuikSystem *cs)
{
  BoxFromVariables(box,&(cs->orig_variables));
}

void ReGenerateOriginalBox(Tparameters *p,Tbox *boxS,Tbox *boxO,TCuikSystem *cs)
{
  if (UpdateCuikSystem(p,cs))
    {
      BoxFromVariables(boxO,&(cs->orig_variables)); 
      UpdateOriginalFromSimple(boxS,boxO,cs->orig2sd); 
    }
  else
    Error("Inconsisten input cuiksystem");
}

/*
  Selects one dimension along which to split box 'b'. The dimension can
  be selected according to the size or according to the error that each
  variable induce in each equation.
  This is similar to ComputeSplitDimInt but the output is an index in 
  the original system and not in the simplified one
 */
unsigned int ComputeSplitDim(Tparameters *p,Tbox *b,TCuikSystem *cs)
{
  Tbox bsimp;
  unsigned int d;

  if (UpdateCuikSystem(p,cs))
    {
      SimpleFromOriginal(b,&bsimp,cs->orig2sd);
      d=GetVarIDInOriginal(ComputeSplitDimInt(p,&bsimp,cs),cs->orig2sd);
      DeleteBox(&bsimp);
    }
  else
    d=NO_UINT;

  return(d);
}

/*
  Returns TRUE if point stored in vector 'v' is included in box 'b'
  Dummy variables are not taken into account in this procedure.
*/
boolean PointInSystemBox(Tvector *v,Tbox *b,TCuikSystem *cs)
{
  unsigned int i,k,nv;
  double *val;
  boolean in;

  nv=NVariables(&(cs->orig_variables));

  if (nv!=GetBoxNIntervals(b))
    Error("Wrong box size in PointInSystemBox");
  
  k=0;
  in=TRUE;
  for(i=0;((in)&&(i<nv));i++)
    {
      if (!IsDummyVariable(i,&(cs->orig_variables)))
        {
          val=(double *)GetVectorElement(k,v);
          if (val==NULL) 
            Error("Vector with incorrect number of elements in PointInSystemBox");
          in=IsInside(*val,GetBoxInterval(i,b));
          k++;
        }
    }
  return(in);
}

/*
  Takes the center of the box and computes the error of this
  point when replaced into the system of equations.
  We use this function to check the validity of the outputs
  of Cuik.
 */
double ErrorInSolution(Tbox *b,TCuikSystem *cs)
{
  double *p; /*central point of the box*/
  unsigned int i;
  double e,max_error;
  Tequation *eq;
  unsigned int nv,ne;

  nv=NVariables(&(cs->orig_variables));

  if (nv!=GetBoxNIntervals(b))
    Error("Wrong box size in ErrorInSolution");

  NEW(p,nv,double);

  for(i=0;i<cs->orig_nvariables;i++)
    p[i]=IntervalCenter(GetBoxInterval(i,b));

  max_error=0.0;
  
  ne=NEquations(&(cs->orig_equations));
  for(i=0;i<ne;i++)
    {
      eq=GetEquation(i,&(cs->orig_equations));
 
      if ((GetEquationType(eq)==SYSTEM_EQ)&&
          (GetEquationCmp(eq)==EQU))
        {
          e=fabs(EvaluateEquation(p,eq)-GetEquationValue(eq)); 

          if (e>max_error) max_error=e;
        }
    }

  free(p);
  
  return(max_error);
}


double ErrorInInequalities(Tbox *b,TCuikSystem *cs)
{
  unsigned int i,ne,nv,cmp;
  Tequation *eq;
  double e,maxError,r,v;
  double *p;

  nv=NVariables(&(cs->orig_variables));

  if (nv!=GetBoxNIntervals(b))
    Error("Wrong box size in ErrorInSolution");

  NEW(p,nv,double);

  for(i=0;i<cs->orig_nvariables;i++)
    p[i]=IntervalCenter(GetBoxInterval(i,b));

  ne=NEquations(&(cs->orig_equations));
  maxError=0;
  i=0;
  while (i<ne)
    {
      eq=GetEquation(i,&(cs->orig_equations));

      cmp=GetEquationCmp(eq);
      switch (cmp)
        {
        case LEQ:
          r=EvaluateEquation(p,eq);
          v=GetEquationValue(eq);
          e=(r<=v?0:r-v);
          break;
        case GEQ:
          r=EvaluateEquation(p,eq);
          v=GetEquationValue(eq);
          e=(r>=v?0:v-r);
          break;
        default:
          e=0;
          break;
        }
      if (e>maxError)
        maxError=e;
      i++;
    }

  free(p);

  return(maxError);
}

boolean CoordInequalitiesHold(Tbox *b,TCuikSystem *cs)
{
  unsigned int i,ne,nv,cmp;
  Tinterval *is;
  boolean hold;
  Tequation *eq;
  Tinterval result;

  nv=NVariables(&(cs->orig_variables));

  if (nv!=GetBoxNIntervals(b))
    Error("Wrong box size in ErrorInSolution");

  is=GetBoxIntervals(b);
  
  ne=NEquations(&(cs->orig_equations));
  hold=TRUE;
  i=0;
  while((i<ne)&&(hold))
    {
      eq=GetEquation(i,&(cs->orig_equations));
      if (GetEquationType(eq)==COORD_EQ)
        {
          cmp=GetEquationCmp(eq);
          switch (cmp)
            {
            case LEQ:
              EvaluateEquationInt(is,&result,eq);
              hold=(UpperLimit(&result)<=GetEquationValue(eq));
              break;
            case GEQ:
              EvaluateEquationInt(is,&result,eq);
              hold=(LowerLimit(&result)>=GetEquationValue(eq));
              break;
            }
        }
      i++;
    }
  return(hold);
}


/*
  Prints the information stored in the cuiksystem
  It prints it in a form than can be parsed again.
*/
void PrintCuikSystem(Tparameters *p,FILE *f,TCuikSystem *cs)
{
  unsigned int nv;
  char **varNames;

  PrintVariables(f,&(cs->orig_variables));
  
  nv=NVariables(&(cs->orig_variables));
  NEW(varNames,nv,char *);
  GetVariableNames(varNames,&(cs->orig_variables));
  PrintEquations(f,varNames,&(cs->orig_equations));
  if (cs->searchMode==MINIMIZATION_SEARCH)
    {
      fprintf(f,"\n[SEARCH]\n\n MIN ");
      PrintMonomials(f,varNames,&(cs->orig_eqMin));
    }
  free(varNames);  
  fprintf(f,"\n\n");
}

/*
  Prints the information stored in the cuiksystem together with the simplified
  CuikSystem.
  It prints it in a form than can be parsed again.
*/
void PrintCuikSystemWithSimplification(Tparameters *p,FILE *f,TCuikSystem *cs)
{
  unsigned int nv;
  char **varNames;

  #if (_DEBUG>1)
    fprintf(f,"%%****************************************\n");
    fprintf(f,"%% Original system \n");
    PrintVariables(f,&(cs->orig_variables));
    
    nv=NVariables(&(cs->orig_variables));
    NEW(varNames,nv,char *);
    GetVariableNames(varNames,&(cs->orig_variables));
    PrintEquations(f,varNames,&(cs->orig_equations));
    if (cs->searchMode==MINIMIZATION_SEARCH)
      {
        fprintf(f,"\n[SEARCH]\n\n MIN ");
        PrintMonomials(f,varNames,&(cs->orig_eqMin));
      }
    fprintf(f,"\n\n");
    free(varNames);
  #endif

  if (UpdateCuikSystem(p,cs))
    {
      fprintf(f,"%%****************************************\n");
      fprintf(f,"%% Simplified system \n");
      fprintf(f,"%% SIMPLIFICATION_LEVEL: %u\n",
              (unsigned int)(GetParameter(CT_SIMPLIFICATION_LEVEL,p)));
      PrintMapping(f,cs->orig2sd);
      PrintVariables(f,&(cs->variables));

      nv=NVariables(&(cs->variables));
      NEW(varNames,nv,char *);
      GetVariableNames(varNames,&(cs->variables));
      PrintEquations(f,varNames,&(cs->equations));
      if (cs->searchMode==MINIMIZATION_SEARCH)
        {
          fprintf(f,"\n[SEARCH]\n\n MIN ");
          PrintMonomials(f,varNames,&(cs->eqMin));
        }
      free(varNames);
    }
  else
    fprintf(f,"INCONSISTENT INPUT CUIK SYSTEM\n");
}

void SaveCuikSystemSimplification(Tparameters *p,FILE *f,TCuikSystem *cs)
{
  if (UpdateCuikSystem(p,cs))
    SaveMapping(f,cs->orig2sd);
  else
    fprintf(f,"INCONSISTENT INPUT CUIK SYSTEM\n");
}

/*
  Destructor of the cuiksystem structure.
*/
void DeleteCuikSystem(TCuikSystem *cs)
{ 
  DeleteVariables(&(cs->orig_variables));
  DeleteEquations(&(cs->orig_equations));  
  DeleteStatistics(&(cs->st)); 

  UnUpdateCuikSystem(cs);

  cs->empty=TRUE;
}