cuiksystem.c
Go to the documentation of this file.
1 #include "cuiksystem.h"
2 
3 #include "error.h"
4 #include "boolean.h"
5 #include "defines.h"
6 #include "box_list.h"
7 #include "box_heap.h"
8 #include "simplex.h"
9 #include "random.h"
10 #include "filename.h"
11 #include "algebra.h"
12 
13 #include <stdlib.h>
14 
15 #include <math.h>
16 #include <string.h>
17 #include <sys/time.h>
18 #include <signal.h>
19 #if (_USE_MPI)
20  #include <mpi.h>
21  #include <sys/resource.h>
22  #include <time.h>
23  #include <unistd.h>
24 #endif
25 
26 
27 
37 /* Internal functions */
38 
56 
69 
70 
88 double EvaluateEqMin(void *b,void *cs);
89 
109 unsigned int ReduceBoxEquationWise(Tparameters *p,Tbox *b,TCuikSystem *cs);
110 
139 unsigned int ReduceBox(Tparameters *p,unsigned int varMask,Tbox *b,TCuikSystem *cs);
140 
162 void PostProcessBox(Tparameters *p,unsigned int c,
163  FILE *f_out,
164  Tlist *sol,
165  Theap *boxes,
166  Tbox *b,
167  TCuikSystem *cs);
168 
181 
182 
213  boolean *replaced,TLinearConstraint *lc,
214  Tbox *borig,TCuikSystem *cs);
215 
260 boolean CSRemoveLCVars(Tparameters *p,unsigned int simplificationLevel,
261  unsigned int level,boolean *changed,
262  boolean *replaced,TLinearConstraint *lc,
263  Tbox *borig,TCuikSystem *cs);
289 
307 void AddJacobianEquationsInt(Tparameters *p,boolean *selectedVars,TJacobian *J,
308  TCuikSystem *cs);
309 
327 
342 
362 unsigned int ComputeSplitDimInt(Tparameters *p,Tbox *b,TCuikSystem *cs);
363 
383 void SaveCSState(char *fname,Tlist *lb,TCuikSystem *cs);
384 
398 void LoadCSState(char *fname,Tlist *lb,TCuikSystem *cs);
399 
400 
401 /************************************************************************/
402 /************************************************************************/
403 /************************************************************************/
404 
405 /*
406  The equation is dummified and the resulting equations and new variables
407  are added to the cuik system.
408 
409  Dummification is performed in the last step of the simplification and, thus,
410  this function operates in the sets of simplified equations and variables
411 
412  Dummification consists in changing all monomials of type x^2 or x*y by new
413  variables 'z' and adding the corresponding equations (z=x^2 or z=x*y) to
414  the system.
415 
416  The level of dummification is controlled by the 'dummify' parameter.
417 
418  This internal function is hidden to the user. In this way dummification is
419  applied only after the system simplification.
420  See AddEquation2CS and SimplifyCuikSystem.
421 */
423 {
424  if (eq!=NULL)
425  {
426  unsigned int nf;
427  unsigned int i,nv,mv;
428  Tvariable_set *vars;
429  unsigned int dummify;
430  unsigned int eqType;
431 
432  dummify=(unsigned int)(GetParameter(CT_DUMMIFY,p));
433 
434  vars=GetEquationVariables(eq);
435  eqType=GetEquationType(eq);
436 
437  nv=VariableSetSize(vars);
438  mv=NVariables(&(cs->variables));
439  for(i=0;i<nv;i++)
440  {
441  if (GetVariableN(i,vars)>=mv)
442  Error("Unknown variable in AddEquation2CS");
443  }
444 
445  nf=GetNumMonomials(eq);
446  if (nf>0)
447  {
448  boolean keep=TRUE; /* True if the equation is kept in its original form */
449 
450  /* Dummy equations are not dummified again */
451  /*dummify equations with just one monomial introduce a
452  too simple variable x=ct that is not simplified since
453  the dummifications is performed after the simplification. */
454  if ((eqType==DUMMY_EQ)||(GetNumMonomials(eq)==1)||(GetEquationOrder(eq)>2))
455  keep=TRUE;
456  else
457  {
458  switch (dummify)
459  {
460  case 0:
461  keep=TRUE;
462  break;
463  case 1:
464  keep=((LinearEquation(eq))||
465  (ParabolaEquation(eq))||
466  (SaddleEquation(eq))||
467  (CircleEquation(eq))||
468  (SphereEquation(eq)));
469  break;
470  case 2:
471  keep=((LinearEquation(eq))||
472  (ParabolaEquation(eq))||
473  (SaddleEquation(eq)));
474  break;
475  case 3:
476  keep=((LinearEquation(eq))||
477  (ParabolaEquation(eq))||
478  (SaddleEquation(eq))||
479  (CircleEquation(eq))||
480  (SphereEquation(eq))||
481  (GetEquationCmp(eq)!=EQU)); /*Inequalities are kept in original form*/
482 
483  break;
484  case 4:
485  keep=((LinearEquation(eq))||
486  (ParabolaEquation(eq))||
487  (SaddleEquation(eq))||
488  (GetEquationCmp(eq)!=EQU)); /*Inequalities are kept in original form*/
489  break;
490  default:
491  Error("Undefined value for DUMMIFY parameter");
492  }
493  }
494 
495  if (keep)
496  AddEquation(eq,&(cs->equations));
497  else
498  {
499  unsigned int j,p,vid1,vid2,vid3;
500  Tmonomial *f,fnew;
501  Tequation eqnew; /* Bilinear and quadratic monomials produce new auxiliary equations */
502  Tequation eqlineal; /* The original equation is transformed into a lineal equation */
503  char *vname;
504  Tinterval range;
505  Tvariable v;
506  char *n1,*n2;
507 
508  InitEquation(&eqlineal);
509 
510  for(j=0;j<nf;j++) /* For all monomials in the equation to be dummified */
511  {
512 
513  f=GetMonomial(j,eq);
514  vars=GetMonomialVariables(f);
515  nv=VariableSetSize(vars);
516 
517  switch(nv)
518  {
519  case 1: /*A single variable in the monomial*/
520  if (GetVariableFunctionN(0,vars)!=NFUN)
521  AddMonomial(f,&eqlineal);
522  else
523  {
524  p=GetVariablePowerN(0,vars);
525 
526  switch(p)
527  {
528  case 1: /*One variable with degree 1*/
529  /* A linear monomial: just add it to the linear equation*/
530 
531  AddMonomial(f,&eqlineal); /* A copy of monomial 'f' is stored in eqlineal*/
532  break;
533 
534  case 2:/*One varialbe With degree 2*/
535  /* Add a linear monomial to eqlinear with a dummy variable
536  and add the corresponding parabola equation*/
537 
538  vid1=GetVariableN(0,vars);
539 
540  /* Generate a new variable dummy that will stand for x^2 */
541  n1=GetVariableName(GetVariable(vid1,&(cs->variables)));
542  NEW(vname,strlen(n1)+20,char);
543  sprintf(vname,"dummy_%s_2",n1);
544 
545  vid2=GetVariableID(vname,&(cs->variables));
546 
547  if (vid2==NO_UINT)
548  {
549  /* if the dummy variable is a new one, add the corresponding
550  dummy equation to the system */
551 
552  NewVariable(DUMMY_VAR,vname,&v);
553  IntervalPow(GetVariableInterval(GetVariable(vid1,&(cs->variables))),2,&range);
554  SetVariableInterval(&range,&v);
555 
556  vid2=AddVariable(&v,&(cs->variables));
557  DeleteVariable(&v);
558 
559  /* Generate a new equation x^2-dummy=0 */
560  InitEquation(&eqnew);
561 
562  InitMonomial(&fnew);
563  AddVariable2Monomial(NFUN,vid1,2,&fnew);
564  AddCt2Monomial( 1.0,&fnew);
565  AddMonomial(&fnew,&eqnew);
566  DeleteMonomial(&fnew);
567 
568  InitMonomial(&fnew);
569  AddVariable2Monomial(NFUN,vid2,1,&fnew);
570  AddCt2Monomial(-1.0,&fnew);
571  AddMonomial(&fnew,&eqnew);
572  DeleteMonomial(&fnew);
573 
574  SetEquationCmp(EQU,&eqnew);
575  SetEquationValue(0.0,&eqnew);
576  /*The derive equations is of the same type as the origiona one*/
577  SetEquationType(DUMMY_EQ,&eqnew);
578 
579  AddEquation(&eqnew,&(cs->equations));
580  DeleteEquation(&eqnew);
581  }
582  free(vname);
583 
584  /*Now add the lineal monomial using the new dummy variable to eqlineal*/
585  InitMonomial(&fnew);
586  AddVariable2Monomial(NFUN,vid2,1,&fnew);
587  AddCt2Monomial(GetMonomialCt(f),&fnew);
588  AddMonomial(&fnew,&eqlineal);
589  DeleteMonomial(&fnew);
590  break;
591 
592  default:
593  Error("The input system to DummifyCuikSystem is not in canonical form");
594  break;
595  }
596  }
597  break;
598 
599  case 2: /*Two variables in the monomial*/
600  if ((GetVariableFunctionN(0,vars)!=NFUN)||
601  (GetVariableFunctionN(1,vars)!=NFUN))
602  AddMonomial(f,&eqlineal);
603  else
604  {
605  if ((GetVariablePowerN(0,vars)==1)&&
606  (GetVariablePowerN(1,vars)==1)) /*The two of them with degree 1*/
607  {
608  /*Get the two variables involved in the product we have to rewrite*/
609  vid1=GetVariableN(0,vars);
610  vid2=GetVariableN(1,vars);
611 
612  /*Generate a new dummy variable*/
613  n1=GetVariableName(GetVariable(vid1,&(cs->variables)));
614  n2=GetVariableName(GetVariable(vid2,&(cs->variables)));
615  NEW(vname,strlen(n1)+strlen(n2)+20,char);
616  sprintf(vname,"dummy_%s_%s",n1,n2);
617 
618  vid3=GetVariableID(vname,&(cs->variables));
619  if (vid3==NO_UINT)
620  {
621 
622  NewVariable(DUMMY_VAR,vname,&v);
625  &range);
626  SetVariableInterval(&range,&v);
627 
628  vid3=AddVariable(&v,&(cs->variables));
629  DeleteVariable(&v);
630 
631  /* Generate the dummy equation X*Y-dummy=0*/
632  InitEquation(&eqnew);
633 
634  InitMonomial(&fnew);
635  AddVariable2Monomial(NFUN,vid1,1,&fnew);
636  AddVariable2Monomial(NFUN,vid2,1,&fnew);
637  AddCt2Monomial( 1.0,&fnew);
638  AddMonomial(&fnew,&eqnew);
639  DeleteMonomial(&fnew);
640 
641  InitMonomial(&fnew);
642  AddVariable2Monomial(NFUN,vid3,1,&fnew);
643  AddCt2Monomial(-1.0,&fnew);
644  AddMonomial(&fnew,&eqnew);
645  DeleteMonomial(&fnew);
646 
647  SetEquationCmp(EQU,&eqnew);
648  SetEquationValue(0.0,&eqnew);
649  /*The derive equations is of the same type as the origiona one*/
650  SetEquationType(DUMMY_EQ,&eqnew);
651 
652  AddEquation(&eqnew,&(cs->equations));
653  DeleteEquation(&eqnew);
654  }
655  free(vname);
656 
657  /*Now add the lineal monomial using the new dummy variable to eqlienal*/
658  InitMonomial(&fnew);
659  AddVariable2Monomial(NFUN,vid3,1,&fnew);
660  AddCt2Monomial(GetMonomialCt(f),&fnew);
661  AddMonomial(&fnew,&eqlineal);
662  DeleteMonomial(&fnew);
663  }
664  else
665  /* One of the 2 variables in the monomial has degree >1 */
666  /* Right now this should never happen since equations
667  with order>2 are directly added to the outpu system */
668  Error("The input system to DummifyCuikSystem is not in canonical form");
669  }
670  break;
671 
672  default:
673  /* A monomial with more than 2 variables */
674  /* Right now this should never happen since equations
675  with order>2 are directly added to the outpu system */
676  Error("The input system to DummifyCuikSystem is not in canonical form");
677  }
678 
679  } /* Proceed with next monomial in the under-dummification equation*/
680 
681  /* The linearized original equation is added to the system */
682 
683  SetEquationCmp(GetEquationCmp(eq),&eqlineal);
684  SetEquationValue(GetEquationValue(eq),&eqlineal);
685  /*The resulting lineal equation is of the same type as the original equation*/
686  SetEquationType(eqType,&eqlineal);
687 
688  AddEquation(&eqlineal,&(cs->equations));
689  DeleteEquation(&eqlineal);
690  }
691  }
692  }
693 }
694 
696 {
697  if (cs->scalar)
698  {
699  unsigned int i,neq;
700  Tequations eqs;
701 
702  CopyEquations(&eqs,&(cs->equations));
703  DeleteEquations(&(cs->equations));
704  InitEquations(&(cs->equations));
705 
706  neq=NEquations(&eqs);
707  for(i=0;i<neq;i++)
708  DummifyAndAddEquation(p,GetEquation(i,&eqs),cs);
709 
710  DeleteEquations(&eqs);
711  }
712 }
713 
715 {
716  boolean empty;
717  boolean small;
718  boolean significantReduction;
719  Tbox b_orig;
720  unsigned int i,j;
721  double rho,smallSigma,epsilon;
722  double mi_r;
723  double reduction;
724 
725  rho=GetParameter(CT_RHO,p);
726  epsilon=GetParameter(CT_EPSILON,p);
727  smallSigma=GetParameter(CT_SMALL_SIGMA,p);
728 
729  /*default values*/
730  empty=FALSE;
731  significantReduction=FALSE;
732  small=FALSE;
733 
734  #if (_DEBUG>2)
735  printf(" Croping box: ");
736  PrintBox(stdout,b);
737  printf(" Size in: %g\n",GetBoxSize(cs->systemVar,b));
738  #endif
739 
740  do {
741  CopyBox(&b_orig,b);
742 
743  for(j=0;((!empty)&&(j<cs->nequations));j++)
744  {
745  #if (_DEBUG>5)
746  printf(" Croping with equation %u: \n",j);
747  printf(" Input box: ");
748  PrintBox(stdout,b);
749  printf(" s: %g v:%g\n",GetBoxSize(NULL,b),GetBoxVolume(NULL,b));
750  #endif
751  if (CropEquation(j,ANY_TYPE_VAR,epsilon,rho,b,&(cs->variables),&(cs->equations))==EMPTY_BOX)
752  empty=TRUE;
753  #if (_DEBUG>5)
754  if (empty)
755  printf(" Empty output box\n");
756  else
757  {
758  printf(" Output box: ");
759  PrintBox(stdout,b);
760  printf(" s: %g v:%g\n",GetBoxSize(NULL,b),GetBoxVolume(NULL,b));
761  }
762  #endif
763  }
764 
765  if (!empty)
766  {
767  mi_r=1.0;
768  for(i=0;i<cs->nvariables;i++)
769  {
770  /*variables of all types are affected by crop*/
771  reduction=IntervalSize(GetBoxInterval(i,b))/IntervalSize(GetBoxInterval(i,&b_orig));
772  if (reduction<mi_r) mi_r=reduction;
773  }
774 
775  significantReduction=(mi_r<rho);
776  small=(GetBoxSize(cs->systemVar,b)<=smallSigma);
777 
778  #if (_DEBUG>2)
779  printf(" Crop reduction: %f\n",mi_r);
780  if (small)
781  printf(" Tiny boxes are not reduced any more\n");
782  else
783  {
784  if (significantReduction)
785  printf(" Large reduction. Cropping again\n");
786  }
787  #endif
788  }
789  DeleteBox(&b_orig);
790 
791  } while ((!empty)&&(!small)&&(significantReduction));
792 
793  #if (_DEBUG>2)
794  if (empty)
795  printf(" Empty in Crop\n");
796  else
797  {
798  printf(" Croped box: ");
799  PrintBox(stdout,b);
800  printf(" Size out: %g\n",GetBoxSize(cs->systemVar,b));
801  }
802  #endif
803 
804  if (empty)
805  return(EMPTY_BOX);
806  else
807  return(REDUCED_BOX);
808 }
809 
810 /*
811  This is the basic function of the solver. This function includes the
812  ShrinkBox procedure described in our papers plus a loop interating
813  ShrinkBox as far as the box is significantly reduced (and still not empty).
814 
815  Thus, ReduceBox takes an input box 'b' and reduces it as much as possible by removing
816  parts that do not include any solution point.
817  It returns information about what happened in the box shrinking process
818 
819  ERROR_IN_PROCESS: When one of the simplex failed due to numerical inestabilities
820  The better the simplex implementation the less likely this
821  output.
822 
823  EMPTY_BOX: The box includes no solution. In this case the output box is
824  undefined (it includes empty intervals).
825 
826  REDUCED_BOX: The box was reduced and it has to be processed further to determine
827  if it includes any solution.
828 
829  REDUCED_BOX_WITH_SOLUTION: The box was been reduced and one of the reduction
830  steps was contractive. This indicates the box
831  includes a solution point.
832  */
833 unsigned int ReduceBox(Tparameters *p,unsigned int varMask,Tbox *b,TCuikSystem *cs)
834 {
835  /*The box can only be reduced if we have constraints*/
836  if (cs->simp_nequations==0)
837  return(REDUCED_BOX);
838  else
839  {
840  boolean empty;
841  boolean significantReduction;
842  boolean small;
843  boolean hasSolution;
844  boolean err;
845  boolean boxOK;
846  unsigned int e;
847  double mi_r,ma_s;
848  double reduction,size;
849  unsigned int i,j;
850  Tbox b_orig;
851  TSimplex lp;
852  double epsilon,smallSigma,rho,lr2tm_q,lr2tm_s;
853  unsigned int safeSimplex;
854  boolean repeatLoop;
855 
856  if ((!cs->updated)||(!cs->consistent))
857  Error("Non-updated/consistent CS in ReduceBox");
858 
859  epsilon=GetParameter(CT_EPSILON,p);
860  smallSigma=GetParameter(CT_SMALL_SIGMA,p);
861  rho=GetParameter(CT_RHO,p);
862  lr2tm_q=GetParameter(CT_LR2TM_Q,p);
863  lr2tm_s=GetParameter(CT_LR2TM_S,p);
864  safeSimplex=(unsigned int)(GetParameter(CT_SAFE_SIMPLEX,p));
865 
866  /*default values*/
867  empty=FALSE;
868  hasSolution=FALSE;
869  err=FALSE;
870  boxOK=TRUE;
871  small=FALSE;
872  significantReduction=FALSE;
873  repeatLoop=FALSE;
874 
875  do {
876  #if (_DEBUG>2)
877  printf("\n Go for box reduction via crop+simplex\n");
878  printf(" InputBox: "); PrintBox(stdout,b);
879  printf(" s: %g v:%g\n",GetBoxSize(NULL,b),GetBoxVolume(NULL,b));
880  #endif
881 
882  /*****************************************************************/
883  /* One step of box reduction */
884 
885  /*#if ((!_USE_MPI)&&(_DEBUG==1))*/
886  #if ((!_USE_MPI)&&(_DEBUG>1))
887  fprintf(stderr,".");fflush(stderr);
888  #endif
889 
890  NewBoxReduction(&(cs->st));
891 
892  /* First we do an equation wise reduction of the ranges -> Crop equation
893  repeated for each equation as far as we get a significant box reduction */
894 
895  empty=(ReduceBoxEquationWise(p,b,cs)==EMPTY_BOX);
896 
897  /* If the equation-wise reduction gives a full box, proceed with the global
898  test based on linear programming */
899 
900  if (empty)
901  repeatLoop=FALSE; /*There is no need to iterate*/
902  else
903  {
904  SimplexCreate(epsilon,cs->nvariables,&lp);
905 
906  /* Add the linear constraints bounding the equations. This include
907  croping the input ranges if necessary */
908 
909  for(j=0;((!empty)&&(j<cs->nequations));j++)
910  {
911  if (!AddEquation2Simplex(j,
912  lr2tm_q,lr2tm_s,
913  epsilon,safeSimplex,rho,
914  b,&(cs->variables),
915  &lp,&(cs->equations)))
916  empty=TRUE;
917  }
918 
919  boxOK=TRUE;
920  err=FALSE;
921  small=FALSE;
922  significantReduction=FALSE;
923 
924  /* Use the simplex to reduce the ranges for the variables */
925  if ((!empty)&&(SimplexNRows(&lp)>0)) /* In extreme cases the simplex tableau can be empty
926  and then nothings can be reduced*/
927  {
928  /* Store the original box to check the actual reduction */
929  CopyBox(&b_orig,b);
930 
931  #if (_SIMPLEX_ENGINE!=_CLP) /* in CLP Reset Simplex has no effect */
932  ResetSimplex(&lp);
933  #endif
934 
935  for(i=0;((boxOK)&&(i<cs->nvariables));i++)
936  {
937  /* The crop takes care of reducing the ranges of the dummy
938  variables in function of the non dummy ones (system,
939  secondary and cartesian).
940  */
941 
942  if ((cs->varType[i]&varMask)&&(!SimplexColEmpty(i,&lp)))
943  {
944  /*Reduce the ranges for variable 'i'*/
945  e=ReduceRange(epsilon,safeSimplex,i,b,&lp);
946  boxOK=(e==REDUCED_BOX);
947  }
948  }
949 
950  /* Analize the output out of the range reduction */
951  if (boxOK)
952  {
953  /* Compute the minimum reduction ratio for all ranges. The smaller
954  the reductin ratio, the maximal the box shrink for that dimension*/
955  mi_r=1.0;
956  ma_s=0.0;
957  for(i=0;i<cs->nvariables;i++)
958  {
959  if (cs->notDummyVar[i])
960  {
961  size=IntervalSize(GetBoxInterval(i,b));
962  if (size>ma_s) ma_s=size;
963 
964  reduction=size/IntervalSize(GetBoxInterval(i,&b_orig));
965  if (reduction<mi_r) mi_r=reduction;
966  }
967  }
968 
969  /* If the output box is fully included in the input one this means
970  the box includes a solution (for sure) */
971  if ((ma_s<lr2tm_q)&&(ma_s<lr2tm_s)&&
972  (BoxInclusion(cs->systemVar,b,&b_orig)))
973  {
974  hasSolution=TRUE;
975  #if (!_USE_MPI)
976  fprintf(stderr,"(>!<)");
977  #endif
978  }
979 
980  small=(ma_s<=smallSigma);
981  significantReduction=(mi_r<rho);
982  }
983  else
984  {
985  empty=(e==EMPTY_BOX);
986  err=(e==ERROR_IN_PROCESS);
987  }
988 
989  DeleteBox(&b_orig);
990  }
991 
992  repeatLoop=((!empty)&&(!err)&&(!small)&&(significantReduction));
993 
994  SimplexDelete(&lp);
995 
996  #if (_DEBUG>2)
997  if (empty)
998  printf(" Empty\n");
999  else
1000  {
1001  if (err)
1002  printf(" Simplex Error\n");
1003  else
1004  {
1005  printf(" OutputBox: "); PrintBox(stdout,b);
1006  printf(" Reduction -> %f\n",mi_r);
1007 
1008  printf(" Size after reduction: %g\n ",GetBoxSize(cs->systemVar,b));
1009  printf(" Volume after reduction: %g\n ",GetBoxVolume(cs->systemVar,b));
1010 
1011  if (significantReduction)
1012  printf(" Large enough reduction-> Try again all equations\n");
1013  else
1014  printf(" Not enough reduction.\n");
1015  }
1016  }
1017  #endif
1018 
1019  } /* end if not empty after */
1020 
1021  } while(repeatLoop);
1022 
1023  if (err)
1024  e=ERROR_IN_PROCESS;
1025  else
1026  {
1027  if (empty)
1028  e=EMPTY_BOX;
1029  else
1030  {
1031  /* If at least on step was a contraction this means
1032  there is a solution inside the box*/
1033  if (hasSolution)
1035  else
1036  e=REDUCED_BOX;
1037  }
1038  }
1039 
1040  return(e);
1041  }
1042 }
1043 
1044 /*
1045  This procedure is called just after using ReduceBox to analyzed what
1046  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; unsigned int nNewton; double epsilon; 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 epsilon=GetParameter(CT_EPSILON,p); nNewton=(unsigned int)GetParameter(CT_MAX_NEWTON_ITERATIONS,p); if ((nNewton>0)&&(f_out!=NULL)&&(GetBoxSize(cs->systemVar,b)>epsilon)) { double *sol; Tbox b_sol; NEW(sol,cs->orig_nvariables,double); if ((CuikNewtonInBox(p,&b_orig,sol,&b_sol,cs)&(CONVERGED_IN_BOX|CONVERGED_IN_GLOBAL))>0) { fprintf(f_out,"Newton [s:%f v:%g l:%u t:%g]: ", GetBoxSize(cs->systemVar,b), GetBoxVolume(cs->systemVar,b), GetBoxLevel(b), GetElapsedTime(&(cs->st))); PrintBoxSubset(f_out,cs->orig_notDummyVar,cs->orig_varNames,&b_sol); } free(sol); DeleteBox(&b_sol); fflush(f_out); } 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 (_DEBUG>1) printf(" Solution number %u with erorr %g\n", GetNSolutionBoxes(&(cs->st)), ErrorInSolution(&b_orig,cs)); #endif 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 ((!VarIncluded(i,GetEquationVariables(&(cs->orig_eqMin))))&& (((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 simplificationLevel, 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; Tbox currentBox; #if (_DEBUG>5) char **varNamesO; nv=NVariables(&(cs->orig_variables)); NEW(varNamesO,nv,char *); GetVariableNames(varNamesO,&(cs->orig_variables)); #endif nv=NVariables(&(cs->variables)); InitBox(nv,NULL,&currentBox); NEW(systemVars,nv,boolean); for(i=0;i<nv;i++) { systemVars[i]=IsSystemVariable(i,&(cs->variables)); SetBoxInterval(i,GetVariableInterval(GetVariable(i,&(cs->variables))),&currentBox); } epsilon=GetParameter(CT_EPSILON,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,&currentBox,&id,&lcr,eq); /* 't' is the type of simplification*/ if (!found) i++; } DeleteBox(&currentBox); /*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 var in original system 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 variable in original system 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,ANY_TYPE_VAR,borig,&(cs->orig_variables),&lcc); if (!hold) printf("\nNon intersecting ranges (I)\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])); DeleteLinearConstraint(&ltmp); CopyLinearConstraint(&lcc,&(lc[j])); AddTerm2LinearConstraint(j,-1.0,&lcc); hold=CropLinearConstraint(ZERO,ANY_TYPE_VAR,borig,&(cs->orig_variables),&lcc); if (!hold) printf("\nNon intersecting ranges (II)\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; Tbox borig,bsimp; unsigned int nTerms,m,k; Tinterval error; unsigned int simplificationLevel; if (!cs->scalar) simplificationLevel=0; /* Matrix equations -> no simplification */ else { simplificationLevel=(unsigned int)(GetParameter(CT_SIMPLIFICATION_LEVEL,p)); if ((!PolynomialEquations(&(cs->orig_equations)))&&(simplificationLevel>1)) simplificationLevel=1; } #if (_DEBUG>6) { char **varNames; printf("%%=======================================================\n"); printf("%% Original system\n"); printf("%%****************************************\n"); PrintVariables(stdout,&(cs->orig_variables)); nv=NVariables(&(cs->orig_variables)); NEW(varNames,nv,char *); GetVariableNames(varNames,&(cs->orig_variables)); PrintEquations(stdout,varNames,&(cs->orig_equations)); free(varNames); } #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)*/ CSRemoveUnusedVars(p,cs); /* This is done even without simplification */ if (simplificationLevel>0) /* For detailed debug systems are not simplified */ { 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,simplificationLevel,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); } } #if(_DEBUG>4) if (hold) { char **varNames; unsigned int nvg; printf("%%=======================================================\n"); printf("%%FINAL fianl 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 SimpleFromOriginal(&borig,&bsimp,cs->orig2s); nvn=NVariables(&(cs->simp_variables)); for(i=0;i<nvn;i++) SetVariableInterval(GetBoxInterval(i,&bsimp),GetVariable(i,&(cs->simp_variables))); DeleteBox(&bsimp); 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; } if (cs->simp_nequations>0) DeleteJacobian(&(cs->J)); cs->simp_nequations=0; cs->simp_nvariables=0;; cs->simp_nee=0; DeleteEquations(&(cs->simp_equations)); DeleteVariables(&(cs->simp_variables)); if (cs->simp_tp!=NULL) free(cs->simp_tp); /*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; /*We assume the system to be consistent although this value is not used unless cs->update is TRUE*/ cs->consistent=TRUE; } } /* Simplifies the CuikSystem and caches information that is useful during solving (cached->faster access). */ boolean UpdateCuikSystem(Tparameters *p,TCuikSystem *cs) { if (!cs->updated) { unsigned int j; boolean found; /* 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 */ cs->consistent=SimplifyCuikSystem(p,cs); if (cs->consistent) { /* Update the information in the original system */ cs->orig_nequations=NEquations(&(cs->orig_equations)); cs->orig_nvariables=NVariables(&(cs->orig_variables)); if ((cs->orig_nequations>0)&&(cs->orig_nvariables==0)) Error("System with equations but without 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); /* Note that orig_systemVar include the secodary variables while systemVar (the one for the simplified/dummified system) keeps the difference between system and secondary varibles */ for(j=0;j<cs->orig_nvariables;j++) { cs->orig_systemVar[j]=((IsSystemVariable(j,&(cs->orig_variables)))|| (IsSecondaryVariable(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)); cs->simp_empty=((cs->simp_nvariables==0)&&(cs->simp_nequations==0)); if (cs->simp_empty) cs->simp_tp=NULL; else { GetVariablesTopology(&(cs->simp_tp),&(cs->simp_variables)); found=FALSE; j=0; while ((!found)&&(j<cs->simp_nvariables)) { found=(cs->simp_tp[j]!=TOPOLOGY_R); j++; } if (!found) { free(cs->simp_tp); cs->simp_tp=NULL; } } if (cs->simp_nequations>0) InitJacobian(&(cs->simp_variables),&(cs->simp_equations),&(cs->J)); /* 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->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)); } } /* 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 RewriteEquation2Simp(GetParameter(CT_EPSILON,p), cs->orig2sd,&(cs->eqMin),&(cs->orig_eqMin)); } /* Mark the system as updated )(althogh it might not be consistent) */ cs->updated=TRUE; } } return(cs->consistent); } /* 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) { if (!UpdateCuikSystem(p,cs)) Error("Inconsistent cuiksystem in ComputeSplitDimInt"); 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; double epsilon; epsilon=GetParameter(CT_EPSILON,p); 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 or too small variables*/ for(i=0;i<cs->nvariables;i++) { s=IntervalSize(&(is[i])); if ((s<10*epsilon)||(s>=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 */ 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) { #if (RANDOMNESS) i=d[randomMax(n-1)]; #else i=d[0]; #endif } 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) { InitConstants(&(cs->constants)); 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->simp_empty=TRUE; cs->updated=FALSE; cs->consistent=TRUE; /*this is not used unless updated=TRUE*/ cs->scalar=TRUE; cs->nequations=0; cs->nvariables=0; cs->systemVar=NULL; cs->notDummyVar=NULL; cs->varType=NULL; cs->orig2sd=NULL; cs->orig2s=NULL; cs->simp_tp=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)); } void VerifyCuikSystem(Tparameters *p,TCuikSystem *cs) { if (!UpdateCuikSystem(p,cs)) Error("Inconsistent input cuiksystem"); } /* Copies one cuiksystem structure into another (that is assumed as initialized but empty). */ void CopyCuikSystem(TCuikSystem *cs_dst,TCuikSystem *cs_src) { CopyConstants(&(cs_dst->constants),&(cs_src->constants)); cs_dst->updated=cs_src->updated; cs_dst->consistent=cs_src->consistent; cs_dst->empty=cs_src->empty; cs_dst->scalar=cs_src->scalar; CopyStatistics(&(cs_dst->st),&(cs_src->st)); cs_dst->searchMode=cs_src->searchMode; CopyEquation(&(cs_dst->orig_eqMin),&(cs_src->orig_eqMin)); CopyEquations(&(cs_dst->orig_equations),&(cs_src->orig_equations)); CopyVariables(&(cs_dst->orig_variables),&(cs_src->orig_variables)); if (cs_dst->updated) { /*The simplified dummified system*/ NEW(cs_dst->orig2sd,1,Tmapping); CopyMapping(cs_dst->orig2sd,cs_src->orig2sd); CopyEquations(&(cs_dst->equations),&(cs_src->equations)); CopyVariables(&(cs_dst->variables),&(cs_src->variables)); 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); } if (cs_dst->searchMode==MINIMIZATION_SEARCH) CopyEquation(&(cs_dst->eqMin),&(cs_src->eqMin)); /*The simplified system*/ NEW(cs_dst->orig2s,1,Tmapping); CopyMapping(cs_dst->orig2s,cs_src->orig2s); cs_dst->simp_empty=cs_src->simp_empty; CopyEquations(&(cs_dst->simp_equations),&(cs_src->simp_equations)); CopyVariables(&(cs_dst->simp_variables),&(cs_src->simp_variables)); cs_dst->simp_nequations=cs_src->simp_nequations; cs_dst->simp_nvariables=cs_src->simp_nvariables; cs_dst->simp_nee=cs_src->simp_nee; if (cs_src->simp_tp!=NULL) { NEW(cs_dst->simp_tp,cs_dst->simp_nvariables,unsigned int); memcpy(cs_dst->simp_tp,cs_src->simp_tp,sizeof(unsigned int)*cs_dst->simp_nvariables); } else cs_dst->simp_tp=NULL; CopyJacobian(&(cs_dst->J),&(cs_src->J)); /*Cached elements from the original system*/ 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; cs_dst->orig_varNames=NULL; } else { NEW(cs_dst->orig_systemVar,cs_dst->orig_nvariables,boolean); memcpy(cs_dst->orig_systemVar,cs_src->orig_systemVar,cs_dst->orig_nvariables*sizeof(boolean)); NEW(cs_dst->orig_notDummyVar,cs_dst->orig_nvariables,boolean); memcpy(cs_dst->orig_notDummyVar,cs_src->orig_notDummyVar,cs_dst->orig_nvariables*sizeof(boolean)); NEW(cs_dst->orig_varNames,cs_dst->orig_nvariables,char *); GetVariableNames(cs_dst->orig_varNames,&(cs_dst->orig_variables)); } } else { cs_dst->orig2sd=NULL; cs_dst->orig2s=NULL; cs_dst->nequations=0; cs_dst->nvariables=0; cs_dst->systemVar=NULL; cs_dst->notDummyVar=NULL; cs_dst->varType=NULL; cs_dst->simp_tp=NULL; cs_dst->orig_nequations=0; cs_dst->orig_nvariables=0; cs_dst->orig_systemVar=NULL; cs_dst->orig_notDummyVar=NULL; cs_dst->orig_varNames=NULL; } } void CuikSystemMerge(Tparameters *p,TCuikSystem *cs1,TCuikSystem *cs2,TCuikSystem *cs) { unsigned int nvar1,nvar2; unsigned int i; InitCuikSystem(cs); MergeConstants(&(cs1->constants),&(cs2->constants),&(cs->constants)); cs->empty=((cs1->empty)&&(cs2->empty)); cs->simp_empty=((cs1->simp_empty)&&(cs2->simp_empty)); cs->updated=FALSE; cs->consistent=TRUE; /*not used untill updated==TRUE*/ cs->scalar=((ScalarEquations(&(cs1->orig_equations)))&& (ScalarEquations(&(cs->orig_equations)))); 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)); MergeEquations(&(cs2->orig_equations),&(cs->orig_equations)); 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)->consistent) Error("The CuikSystem must be consistent before using EvaluateEqMin"); if (((TCuikSystem *)cs)->searchMode==MINIMIZATION_SEARCH) { #if (EQ_MIN_IN_CENTER) unsigned int nv; double *p; unsigned int i; nv=NVariables(&(((TCuikSystem *)cs)->variables)); NEW(p,nv,double); for(i=0;i<nv;i++) p[i]=IntervalCenter(GetBoxInterval(i,(Tbox *)b)); v=EvaluateWholeEquation(p,&(((TCuikSystem *)cs)->eqMin)); 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(NFUN,v,2,&m); AddMonomial(&m,&eq); ResetMonomial(&m); AddCt2Monomial(-w*2*val,&m); AddVariable2Monomial(NFUN,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)); } void AddMatrixEquation2CS(Tparameters *p,TMequation *eq,TCuikSystem *cs) { if (cs->updated) UnUpdateCuikSystem(cs); if (eq!=NULL) { cs->scalar=FALSE; if (!SimplifiedMEquation(eq)) Error("Adding a non-simplified matrix equation to the cuiksystem"); //SimplifyMEquation(eq); AddMatrixEquation(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)); } void GetCSVariableNames(char **varNames,TCuikSystem *cs) { /* We do not use the cached orig_varnames since the cs can be non-updated */ GetVariableNames(varNames,&(cs->orig_variables)); } /* Returns the number of variables in the CuikSystem */ unsigned int GetCSNumVariables(TCuikSystem *cs) { return(NVariables(&(cs->orig_variables))); } unsigned int GetCSNumSystemVariables(TCuikSystem *cs) { return(GetNumSystemVariables(&(cs->orig_variables))+ GetNumSecondaryVariables(&(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))); } char *GetCSVariableName(unsigned int id,TCuikSystem *cs) { return(VariableName(id,&(cs->orig_variables))); } boolean IsSystemVarInSimpCS(Tparameters *p,char *v,TCuikSystem *cs) { unsigned int n; if (!UpdateCuikSystem(p,cs)) Error("Inconsistent input cuiksystem"); n=GetVariableID(v,&(cs->simp_variables)); return((n!=NO_UINT)&& ((IsSystemVariable(n,&(cs->simp_variables)))|| (IsSecondaryVariable(n,&(cs->simp_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); } unsigned int GetCSVarTopology(unsigned int vID,TCuikSystem *cs) { return(GetVariableTopology(GetVariable(vID,&(cs->orig_variables)))); } /* 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) { if (!cs->scalar) Error("GetCSEquation only operates on scalar systems"); CopyEquation(eq,GetEquation(n,&(cs->orig_equations))); } boolean IsCSPolynomial(TCuikSystem *cs) { return(PolynomialEquations(&(cs->orig_equations))); } boolean IsCSScalar(TCuikSystem *cs) { return(cs->scalar); } /* Returns the number of equations in the CuikSytem */ unsigned int GetCSNumEquations(TCuikSystem *cs) { return(NEquations(&(cs->orig_equations))); } void GetCSJacobian(TJacobian *J,TCuikSystem *cs) { InitJacobian(&(cs->orig_variables),&(cs->orig_equations),J); } unsigned int GetSimpCSTopology(Tparameters *p, unsigned int **t,TCuikSystem *cs) { if (!UpdateCuikSystem(p,cs)) Error("Inconsistent input cuiksystem"); return(GetVariablesTopology(t,&(cs->simp_variables))); } void GetSimpCSJacobian(Tparameters *p,TJacobian *J,TCuikSystem *cs) { if (!cs->scalar) GetCSJacobian(J,cs); else { if (!UpdateCuikSystem(p,cs)) Error("Inconsistent input cuiksystem"); InitJacobian(&(cs->simp_variables),&(cs->simp_equations),J); } } void AddJacobianEquationsInt(Tparameters *p,boolean *selectedVars,TJacobian *J, TCuikSystem *cs) { Tinterval one_one,zero_one; unsigned int *nvl; Tvariable v; char vname[100]; Tequation eq,eq1; unsigned int i,j,nvars,neqs; if (!cs->scalar) Error("AddJacobianEquationsInt only operates on scalar systems"); GetJacobianSize(&neqs,&nvars,J); NewInterval(-1.0,1.0,&one_one); NewInterval( 0.0,1.0,&zero_one); /* From this point we start to modify cs and thus, all queries to the original cs must be done on csOrig*/ /*define the lambda variables*/ NEW(nvl,neqs,unsigned int); for(i=0;i<neqs;i++) { sprintf(vname,"_lambda_%u",i); NewVariable(SYSTEM_VAR,vname,&v); SetVariableInterval((i==0?&zero_one:&one_one),&v); nvl[i]=AddVariable2CS(&v,cs); DeleteVariable(&v); } /*Lambdas must be normalized*/ GenerateGeneralNormEquation(neqs,nvl,1.0,&eq); AddEquation2CS(p,&eq,cs); DeleteEquation(&eq); /*define the linear combination using the lambda*/ for(j=0;j<nvars;j++) { if ((selectedVars==NULL)||(selectedVars[j])) { InitEquation(&eq); SetEquationCmp(EQU,&eq); SetEquationType(SYSTEM_EQ,&eq); for(i=0;i<neqs;i++) { CopyEquation(&eq1,GetJacobianEquation(i,j,J)); VarScaleEquation(nvl[i],&eq1); AccumulateEquations(&eq1,1.0,&eq); DeleteEquation(&eq1); } AddEquation2CS(p,&eq,cs); DeleteEquation(&eq); } } } void AddJacobianEquations(Tparameters *p,boolean *selectedVars,TCuikSystem *cs) { unsigned int neqs,nvars,i,ns; TCuikSystem csOrig; TJacobian J; if (!UpdateCuikSystem(p,cs)) Error("Inconsistent input cuiksystem"); CopyCuikSystem(&csOrig,cs); GetCSJacobian(&J,&csOrig); GetJacobianSize(&neqs,&nvars,&J); if (selectedVars==NULL) ns=nvars; else { ns=0; for(i=0;i<nvars;i++) if (selectedVars[i]) ns++; } if ((neqs>0)&&(ns>0)) AddJacobianEquationsInt(p,selectedVars,&J,cs); DeleteJacobian(&J); UnUpdateCuikSystem(cs); DeleteCuikSystem(&csOrig); } void AddSimplifiedJacobianEquations(Tparameters *p,boolean *selectedVars,TCuikSystem *cs) { unsigned int neqs,nvars,i,ns; boolean *simpSelectedVars; unsigned int simpID,nt; TLinearConstraint lc; TCuikSystem csOrig; if (!cs->scalar) Error("AddSimplifiedJacobianEquations only operates on scalar systems"); if (!UpdateCuikSystem(p,cs)) Error("Inconsistent input cuiksystem"); CopyCuikSystem(&csOrig,cs); neqs=cs->simp_nee; nvars=cs->simp_nvariables; NEW(simpSelectedVars,nvars,boolean); if (selectedVars==NULL) ns=nvars; else { /*Translated the selected vars from original to simplified*/ for(i=0;i<nvars;i++) simpSelectedVars[i]=TRUE; ns=0; for(i=0;i<cs->orig_nvariables;i++) { if (!selectedVars[i]) { GetOriginalVarRelation(i,&lc,cs->orig2s); nt=GetNumTermsInLinearConstraint(&lc); /* if nt is zero the variable is set to constant and must be not taken into account*/ if (nt>0) { ns++; /*We use the first variable in the simplified system as equivalent to the original variable*/ simpID=GetLinearConstraintVariable(0,&lc); simpSelectedVars[simpID]=FALSE; #if (_DEBUG>0) printf(" Var %s -> %s\n",cs->orig_varNames[i], VariableName(simpID,&(cs->simp_variables))); #endif } #if (_DEBUG>0) else printf(" Var %s -> constant (not considered)\n",cs->orig_varNames[i]); #endif DeleteLinearConstraint(&lc); } } } if ((neqs>0)&&(ns>0)) { TJacobian J; Tequation *eq; unsigned int j; /* We have to translate the simplified Jacobian matrix to the original system */ InitJacobian(&(csOrig.simp_variables),&(csOrig.simp_equations),&J); for(i=0;i<neqs;i++) { for(j=0;j<nvars;j++) { eq=GetJacobianEquation(i,j,&J); RewriteEquation2Orig(csOrig.orig2s,eq,eq); } } #if (_DEBUG>1) { char **varNames; NEW(varNames,csOrig.orig_nvariables,char*); GetVariableNames(varNames,&(csOrig.orig_variables)); for(j=0;j<nvars;j++) { for(i=0;i<neqs;i++) { printf("J[%u,%u] ->",i,j); PrintEquation(stdout,varNames,GetJacobianEquation(i,j,&J)); } } free(varNames); } #endif AddJacobianEquationsInt(p,simpSelectedVars,&J,cs); DeleteJacobian(&J); } UnUpdateCuikSystem(cs); DeleteCuikSystem(&csOrig); } /* 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,double *reduction,Tbox *b,TCuikSystem *cs) { unsigned int e; if (!cs->scalar) Error("MaxReduction only operates on scalar systems"); if (UpdateCuikSystem(p,cs)) { Tbox bsimp; double sIn,sOut; SimpleFromOriginal(b,&bsimp,cs->orig2sd); sIn=GetBoxSize(cs->systemVar,&bsimp); e=ReduceBox(p,(varMask==0?~DUMMY_VAR:varMask),&bsimp,cs); if (e==EMPTY_BOX) *reduction=0; else { sOut=GetBoxSize(cs->systemVar,&bsimp); *reduction=sOut/sIn; /* 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; if (!cs->scalar) Error("SampleCuikSystem only operates on scalar systems"); 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; double reduction; boolean ok; boolean end; /* free var = variable not in equations -> must be fixed always */ unsigned int *freeVars; unsigned int nFreeVars; unsigned int nsols; if (!cs->scalar) Error("SampleCuikSystemInBox only operates on scalar systems"); /*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 parameter 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); /*keep the original nsols parameter*/ nsols=(unsigned int)GetParameter(CT_N_SOLUTIONS,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,LINKS_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,&reduction,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; /*number of tries so far*/ while(!end) { do { ok=TRUE; /*fix 'ndof' variables*/ 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' variables*/ /* This updates the system, including simplification */ 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 (_DEBUG>1) if (ok) {printf("[CS]");PrintBox(stdout,&b);} else 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*/ ChangeParameter(CT_N_SOLUTIONS,nsamples-ns,p); SolveCuikSystem(p,FALSE,NULL,NULL,f_out_box,&solutions,cs); /*We have to print the samples*/ GetCSSystemVars(&systemVars,cs); InitIterator(&it,&solutions); First(&it); while((!EndOfList(&it))&&(!end)) { 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++; if (ns<nsamples) end=TRUE; 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, min_sigma, and n_solutions*/ ChangeParameter(CT_SPLIT_TYPE,spt,p); ChangeParameter(CT_SMALL_SIGMA,ms,p); ChangeParameter(CT_N_SOLUTIONS,nsols,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); } boolean IncrementalSampleCuikSystem(Tparameters *p,char *fname,Tlist *sb, boolean *fixVars, unsigned int nsamples,unsigned int ntries, unsigned int ndof,TCuikSystem *cs) { Tbox init_box; boolean haveSample; if (!cs->scalar) Error("IncrementalSampleCuikSystem only operates on scalar systems"); GenerateInitialBox(&init_box,cs); haveSample=IncrementalSampleCuikSystemInBox(p,fname,sb,fixVars,nsamples,ntries,ndof,&init_box,cs); DeleteBox(&init_box); return(haveSample); } boolean IncrementalSampleCuikSystemInBox(Tparameters *p,char *fname,Tlist *sb, boolean *fixVars, unsigned int nsamples,unsigned int ntries, unsigned int ndof, Tbox *init_box,TCuikSystem *cs) { Tbox b,*bAux,*box; Tlist solutions; unsigned int k,d,ns,nv,ne,i,nt; double v; double smallSigma; Tfilename fsols_box; Tfilename fsols_point; FILE *f_out_box; FILE *f_out_point; Titerator it; boolean *systemVars; Tinterval *r,newr; Tbox *bl; double ms; double reduction; boolean ok; boolean end; boolean *fv; unsigned int *ind2fv; unsigned int nfv; unsigned int newtonIterations; double *sol; Tbox b_sol; Tbox orig_ranges; unsigned int level,*attempts,maxAttemptsPerLevel; unsigned int nsols; if (!cs->scalar) Error("IncrementalSampleCuikSystemInBox only operates on scalar systems"); /*Set the ranges for the 'cs' variables to those in the given box (init_box)*/ if (!UpdateCuikSystem(p,cs)) Error("Inconsistent cuiksystem in IncrementalSampleCuikSystemInBox"); BoxFromVariables(&orig_ranges,&(cs->orig_variables)); nv=NVariables(&(cs->orig_variables)); NEW(sol,nv,double); ms=GetParameter(CT_SMALL_SIGMA,p); smallSigma=100*GetParameter(CT_EPSILON,p); ChangeParameter(CT_SMALL_SIGMA,smallSigma,p); nsols=(unsigned int)GetParameter(CT_N_SOLUTIONS,p); newtonIterations=(unsigned int)GetParameter(CT_MAX_NEWTON_ITERATIONS,p); /* Array of variables that can be fixed */ NEW(fv,nv,boolean); NEW(ind2fv,nv,unsigned int); nfv=0; for(k=0;k<nv;k++) { if ((IntervalSize(GetBoxInterval(k,init_box))>smallSigma)&& (IsInSimple(k,cs->orig2sd))) { if (((fixVars!=NULL)&&(fixVars[k]))|| ((fixVars==NULL)&&(IsSystemVariable(k,&(cs->orig_variables))))) { fv[k]=TRUE; ind2fv[nfv]=k; nfv++; } else fv[k]=FALSE; } else fv[k]=FALSE; } 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,LINKS_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 (_DEBUG>1) printf("Initial reduction (s:%g)\n",GetBoxSumSide(NULL,init_box)); #endif if (MaxReduction(p,~DUMMY_VAR,&reduction,init_box,cs)==EMPTY_BOX) end=TRUE; else end=FALSE; #if (_DEBUG>1) if (!end) { printf("[CS] (s:%g)",GetBoxSumSide(NULL,init_box)); PrintBox(stdout,init_box); } else printf("Inconsistent system while sampling\n"); fflush(stdout); #endif NEW(attempts,ndof,unsigned int); NEW(bl,ndof,Tbox); maxAttemptsPerLevel=10; ns=0; /*samples already determined*/ nt=0; /*time we tried to sample so far*/ while((!end)&&(ns<nsamples)) { ok=TRUE; CopyBox(&b,init_box); /*fix ndof variables*/ level=0; for(k=0;k<ndof;k++) attempts[k]=0; while((level<ndof)&&(!end)) { /*Keep a copy of the box BEFORE fixing the new variable at the current 'level'. This is used when backtracking to recover the state before this assigment*/ CopyBox(&(bl[level]),&b); do { /*select one of the still not fixed variables*/ d=ind2fv[randomMax(nfv-1)]; r=GetBoxInterval(d,&b); } while (IntervalSize(r)<=smallSigma); #if (_DEBUG>1) printf("[CS](level %u/%u bs:%f)",level,ndof,GetBoxSumSide(NULL,&b)); #endif { Tinterval raux; double c,l; c=IntervalCenter(r); l=IntervalSize(r)/2; NewInterval(c-l,c+l,&raux); /*Now change the current box ,b*/ v=randomInInterval(&raux); } //v=randomInInterval(r); NewInterval(v,v,&newr); SetBoxInterval(d,&newr,&b); /* Change the variable ranges to those in b */ UnUpdateCuikSystem(cs); VariablesFromBox(&b,&(cs->orig_variables)); attempts[level]++; #if (_DEBUG>1) printf("(attempts %u/%u)->fixing %u to %f \n", attempts[level],maxAttemptsPerLevel,d,v); #endif /* And reduce the systems as much as possible */ ok=((UpdateCuikSystem(p,cs))&& (MaxReduction(p,~DUMMY_VAR,&reduction,&b,cs)!=EMPTY_BOX)); if (ok) { if (newtonIterations>0) { if ((CuikNewtonInBox(p,&b,sol,&b_sol,cs)&(CONVERGED_IN_GLOBAL|CONVERGED_IN_BOX))>0) { fprintf(f_out_box,"Newton (%g)",ErrorInSolution(&b_sol,cs)); PrintBoxSubset(f_out_box,cs->orig_notDummyVar,cs->orig_varNames,&b_sol); fflush(f_out_box); } DeleteBox(&b_sol); } level++; #if (_DEBUG>1) printf("[CS] (l:%u s:%g)",level,GetBoxSumSide(NULL,&b)); PrintBox(stdout,&b); #endif } else { #if (_DEBUG>1) printf("Inconsistent system while sampling. Backtracking\n"); #endif /*Undo the last assigment*/ DeleteBox(&b); CopyBox(&b,&(bl[level])); DeleteBox(&(bl[level])); if (attempts[level]>maxAttemptsPerLevel) { attempts[level]=0; if (level>0) { level--; /*move to the state just before the assigment at leve-1*/ DeleteBox(&b); CopyBox(&b,&(bl[level])); DeleteBox(&(bl[level])); } } } } if (ok) { /* Last degree of freedom is not fixed but we solve the cuiksystem in the box */ InitListOfBoxes(&solutions); /*See if there are solutions for the fixed ranges*/ ChangeParameter(CT_N_SOLUTIONS,nsamples-ns,p); SolveCuikSystem(p,FALSE,NULL,&b,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); } DeleteBox(&b); nt++; if ((ntries!=NO_UINT)&&(nt==ntries)) end=TRUE; } 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 small_sigma and number of solutions*/ ChangeParameter(CT_SMALL_SIGMA,ms,p); ChangeParameter(CT_N_SOLUTIONS,nsols,p); free(attempts); free(bl); free(fv); free(ind2fv); free(sol); UnUpdateCuikSystem(cs); VariablesFromBox(&orig_ranges,&(cs->orig_variables)); DeleteBox(&orig_ranges); return(!end); } unsigned int CuikNewtonSimp(Tparameters *p,double *x,TCuikSystem *cs) { boolean converged=FALSE; unsigned int out=DIVERGED; if (UpdateCuikSystem(p,cs)) { if (cs->simp_nequations==0) { out=CONVERGED_IN_BOX; } else { TNewton newton; double *Jx_d; double *b_d; double epsilon,nullSingularValue; double errorVal; int err; unsigned int nStep,nStepMax; #if (_DEBUG>2) { unsigned int i; printf("Starting Newton Iteration form: ["); for(i=0;i<cs->simp_nvariables;i++) printf("%f ",x[i]); printf("]\n"); } #endif InitNewton(cs->simp_nee,cs->simp_nvariables,&newton); Jx_d=GetNewtonMatrixBuffer(&newton); b_d=GetNewtonRHBuffer(&newton); epsilon=GetParameter(CT_EPSILON,p); nullSingularValue=epsilon/100; nStepMax=(unsigned int)GetParameter(CT_MAX_NEWTON_ITERATIONS,p); if (nStepMax==0) Error("Parameter MAX_NEWTON_ITERATIONS must be larger than 0 to use CuikNewton"); nStep=0; while((!converged)&&(nStep<nStepMax)) { EvaluateJacobianInVector(x,cs->simp_nee,cs->simp_nvariables,Jx_d,&(cs->J)); EvaluateEqualityEquations(FALSE,x,b_d,&(cs->simp_equations)); err=NewtonStep(nullSingularValue,x,&errorVal,&newton); if (err) nStep=nStepMax; else { #if (_DEBUG>2) { unsigned int i; printf(" Iteration: %u error: %g",nStep,errorVal); #if (_DEBUG>3) printf(" x: [\n"); for(i=0;i<cs->simp_nvariables;i++) printf("%f ",x[i]); printf("]"); #endif printf("\n"); } #endif /* stopCriterion test */ if(errorVal<epsilon) converged=TRUE; nStep++; } } if (cs->simp_tp!=NULL) ArrayPi2Pi(cs->simp_nvariables,cs->simp_tp,x); DeleteNewton(&newton); if ((converged)&&(ErrorInSimpCSEquations(p,x,cs)<epsilon)) out=CONVERGED_IN_GLOBAL; else out=DIVERGED; } } else Error("Inconsistent input cuiksystem"); return(out); } unsigned int CuikNewtonInBox(Tparameters *p,Tbox *bIn,double *sol,Tbox *b_sol,TCuikSystem *cs) { boolean converged=FALSE; unsigned int out=DIVERGED; if (UpdateCuikSystem(p,cs)) { if (cs->simp_nequations==0) { unsigned int i; Tbox ib; /* The system of equations is empty -> just pic a random value for the variables */ BoxFromVariables(&ib,&(cs->orig_variables)); for(i=0;i<cs->orig_nvariables;i++) sol[i]=randomInInterval(GetBoxInterval(i,&ib)); DeleteBox(&ib); InitBoxFromPoint(cs->orig_nvariables,sol,b_sol); out=CONVERGED_IN_BOX; } else { Tbox bsimp; Tinterval c,*r; unsigned int nv,ne,bs; unsigned int i,j,k,l; double *m,*h; double *x; double *Jx_d; double *b_d; /* Newton required variables */ TNewton newton; double epsilon; double nullSingularValue; double errorVal; unsigned int nStep,nStepMax; int err; Tbox box; /*randomSet(2341);*/ bs=GetBoxNIntervals(bIn); if (cs->orig_nvariables<bs) Error("Box size missmatch in CuikNewtonInBox"); else { if (cs->orig_nvariables>bs) { if (cs->orig_nvariables==(bs+GetNumDummyVariables(&(cs->orig_variables)))) { /* We are given a solution box without the dummy vars */ /* We need to compute the values for the dummy vars from the system ones*/ GenerateInitialBox(&box,cs); /*all variables to initial ranges*/ SetBoxSubset(cs->orig_systemVar,bIn,&box); /* only dummies to initial ranges */ RegenerateSolution(p,&box,cs); /* define dummies from system variables */ } else Error("Box size missmatch in CuikNewtonInBox"); } else CopyBox(&box,bIn); } nStepMax=(unsigned int)GetParameter(CT_MAX_NEWTON_ITERATIONS,p); if (nStepMax==0) Error("Parameter MAX_NEWTON_ITERATIONS must be larger than 0 to use CuikNewton"); SimpleFromOriginal(&box,&bsimp,cs->orig2s); /* We will add some variables and equations to deal with ranges */ #if (NEWTON_WITHIN_RANGES) nv=cs->simp_nvariables+cs->simp_nvariables; ne=cs->simp_nee+cs->simp_nvariables; #else nv=cs->simp_nvariables; ne=cs->simp_nee; #endif /* Space for the variable ranges */ NEW(m,cs->simp_nvariables,double); NEW(h,cs->simp_nvariables,double); for(i=0;i<cs->simp_nvariables;i++) { r=GetBoxInterval(i,&bsimp); m[i]=IntervalCenter(r); h[i]=IntervalSize(r)/2.0; } /* Matrices and arrays used in the Newton iteration */ NEW(x,nv,double); for(i=0;i<cs->simp_nvariables;i++) x[i]=randomInInterval(GetBoxInterval(i,&bsimp)); for(i=cs->simp_nvariables,j=0;i<nv;i++,j++) x[i]=asin((x[j]-m[j])/h[j]); #if (_DEBUG>2) { printf("Starting Newton Iteration form: x=["); for(i=0;i<nv;i++) printf("%f ",x[i]); printf("]\n"); } #endif InitNewton(ne,nv,&newton); Jx_d=GetNewtonMatrixBuffer(&newton); k=ne*nv; if (k!=(cs->simp_nee*cs->simp_nvariables)) { /* If we have limit constraints just init the full buffer to zero. */ for(i=0;i<k;i++) Jx_d[i]=0; } b_d=GetNewtonRHBuffer(&newton); epsilon=GetParameter(CT_EPSILON,p); nullSingularValue=epsilon/100; nStep=0; while((!converged)&&(nStep<nStepMax)) { EvaluateJacobianInVector(x,ne,nv,Jx_d,&(cs->J)); EvaluateEqualityEquations(FALSE,x,b_d,&(cs->simp_equations)); for(i=cs->simp_nvariables,j=0,l=cs->simp_nee;i<nv;i++,j++,l++) { /* 'j' is the index for the original variable 'v' */ /* 'i' is the index for the new variable 'q' */ NewtonSetMatrix(l,j,1,&newton); NewtonSetMatrix(l,i,-h[j]*cos(x[i]),&newton); NewtonSetRH(l,(x[j]-m[j]-h[j]*sin(x[i])),&newton); } errorVal=Norm(ne,b_d); //fprintf(stderr,"%u Newton Error: %g\n",nStep,errorVal); converged=(errorVal<epsilon); if (!converged) { err=NewtonStep(nullSingularValue,x,&errorVal,&newton); if (err) nStep=nStepMax; else { #if (_DEBUG>2) printf(" Iteration: %u error: %g",nStep,errorVal); #if (_DEBUG>3) printf(" x: [\n"); for(i=0;i<cs->simp_nvariables;i++) printf("%f ",x[i]); printf("]"); #endif printf("\n"); #endif nStep++; } } } DeleteNewton(&newton); if (cs->simp_tp!=NULL) ArrayPi2Pi(cs->simp_nvariables,cs->simp_tp,x); for(i=0;i<cs->simp_nvariables;i++) { NewInterval(x[i],x[i],&c); SetBoxInterval(i,&c,&bsimp); } CopyBox(b_sol,&box); UpdateOriginalFromSimple(&bsimp,b_sol,cs->orig2s); /* Some dummy variables in original system might not be present in the simplified system but their value can be trivially deduced from the system variables*/ RegenerateSolution(p,b_sol,cs); /*At this point b_sol must be a punctual box (fully defined) */ for(i=0;i<cs->orig_nvariables;i++) sol[i]=IntervalCenter(GetBoxInterval(i,b_sol)); if (!converged) //((!converged)||(ErrorInSolution(b_sol,cs)>epsilon)) out=DIVERGED; else { Tbox init_box; GenerateInitialBox(&init_box,cs); if ((ErrorInInequalities(b_sol,cs)>epsilon)|| (!PointInBox(cs->orig_notDummyVar,cs->orig_nvariables,sol,epsilon,&init_box))) out=CONVERGED_OUTSIDE_GLOBAL; else { if (PointInBox(cs->orig_systemVar,cs->orig_nvariables,sol,epsilon,&box)) out=CONVERGED_IN_BOX; else out=CONVERGED_IN_GLOBAL; } DeleteBox(&init_box); } DeleteBox(&bsimp); free(x); free(m); free(h); DeleteBox(&box); } } else Error("Inconsistent input cuiksystem"); return(out); } /* Find a solution of a cuiksystem set of equations using the Newton-Rhapson method. If the system is undetermined, the method is based on the Moor-Penrose generalised inverse. The method can diverge if the initial point is far from the solutions of the system. */ boolean CuikNewton(Tparameters *p,double *sol,Tbox *b_sol,TCuikSystem *cs) { Tbox init_box; boolean converged; GenerateInitialBox(&init_box,cs); converged=((CuikNewtonInBox(p,&init_box,sol,b_sol,cs)&(CONVERGED_IN_BOX|CONVERGED_IN_GLOBAL))>0); DeleteBox(&init_box); return(converged); } /* 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, boolean restart,char *fstate,Tbox *searchSpace, FILE *f_out,Tlist *sol,TCuikSystem *cs) { Theap boxes; #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; unsigned int nsols; if (!cs->scalar) Error("SolveCuikSystem only operates on scalar systems"); /************************************************************************************/ /************************************************************************************/ /************************************************************************************/ if (UpdateCuikSystem(p,cs)) { if (cs->simp_empty) { Tbox b; /* The system of equations is empty -> just pic a random value for the variables */ InitBox(cs->simp_nvariables,NULL,&b); PostProcessBox(p,REDUCED_BOX_WITH_SOLUTION,f_out,sol,&boxes,&b,cs); DeleteBox(&b); } else { 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 { if (searchSpace==NULL) { BoxFromVariables(&init_box,&(cs->variables)); AddBox2HeapOfBoxes(&init_box,&boxes); DeleteBox(&init_box); } else { SimpleFromOriginal(searchSpace,&init_box,cs->orig2sd); AddBox2HeapOfBoxes(&init_box,&boxes); } InitStatistics(1,HeapOfBoxesVolume(cs->systemVar,&boxes),&(cs->st)); } current_level=0; nb=0; statePeriod=(unsigned int)GetParameter(CT_STATE_PERIOD,p); nsols=(unsigned int)GetParameter(CT_N_SOLUTIONS,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); #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,Tbox *searchSpace, 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; unsigned int nsols; boolean done; if (!cs->scalar) Error("MPI_SolveCuikSystem only operates on scalar systems"); /*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 { if (searchSpace==NULL) AddBox2HeapOfBoxes(&init_box,&boxes); else { Tbox lb; SimpleFromOriginal(searchSpace,&lb,cs->orig2sd); AddBox2HeapOfBoxes(&lb,&boxes); DeleteBox(&lb); } 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); nsols=(unsigned int)GetParameter(CT_N_SOLUTIONS,p); done=FALSE; /*the cuik solve process itself*/ do { if ((!done)&& /*we don't have enough solutions*/ (!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); if (nsols>0) done=(GetNSolutionBoxes(&(cs->st))>=nsols); 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); /* Repeat until we don't have enough solutions and there are still some boxes from where we can get solutions. If we already have all the desired solutions but some boxes are in_process, we wait untill their process finishes (so that all processors are ready to receive the kill signal) */ } while(((!done)&&(!HeapEmpty(&boxes)))||(in_process>0)); /* 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 (!cs->scalar) Error("MPI_TreatBox only operates on scalar systems"); 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 GenerateSimpInitialBox(Tparameters *p,Tbox *box,TCuikSystem *cs) { if (!UpdateCuikSystem(p,cs)) Error("Inconsistent system in GenerateSimpInitialBox"); BoxFromVariables(box,&(cs->simp_variables)); } boolean RegenerateSolution(Tparameters *p,Tbox *b,TCuikSystem *cs) { boolean ok; if (UpdateCuikSystem(p,cs)) { unsigned int i; Tbox bInit; double epsilon,rho; Tinterval all; /* If we set the dummies to range [-INF,INF] they are not adjusted via crop. */ NewInterval(-INF/2,INF/2,&all); epsilon=GetParameter(CT_EPSILON,p); rho=GetParameter(CT_RHO,p); InitBox(cs->orig_nvariables,NULL,&bInit); /* System and secondary vars already have value in the input box but dummies and cartesian variables probably not. These are set to inf range to avoid any limit in their value when regenerated (boxes commming from Newton can have extrange values in dummy variables. */ for(i=0;i<cs->orig_nvariables;i++) { if ((IsDummyVariable(i,&(cs->orig_variables)))|| (IsCartesianVariable(i,&(cs->orig_variables)))) SetBoxInterval(i,&all,b); } DeleteBox(&bInit); ok=TRUE; /* First we define the value of the dummy variables from the system ones. They are related by dymmy equations and, given the system varibles, the dummy variables can be easily defined via crop. */ for(i=0;((ok)&&(i<cs->orig_nequations));i++) { if (IsDummyEquation(i,&(cs->orig_equations))) ok=(CropEquation(i,DUMMY_VAR,epsilon,rho,b,&(cs->orig_variables),&(cs->orig_equations))!=EMPTY_BOX); } /* Now we try to define the cartesian variables from the system and dummy ones. The separating plane cartesian variables can not be properly bounded this way but the ones for body vertexes can (and they are the relevant ones) */ for(i=0;((ok)&&(i<cs->orig_nequations));i++) { if (IsCoordEquation(i,&(cs->orig_equations))) ok=(CropEquation(i,CARTESIAN_VAR,epsilon,rho,b,&(cs->orig_variables),&(cs->orig_equations))!=EMPTY_BOX); } } else ok=FALSE; return(ok); } unsigned int RegenerateSolutionPoint(Tparameters *p,double *pt,double **rp,TCuikSystem *cs) { unsigned int i,k; Tbox b; if (!UpdateCuikSystem(p,cs)) Error("Inconsistent system in RegenerateSolutionPoint"); NEW(*rp,cs->orig_nvariables,double); k=0; for(i=0;i<cs->orig_nvariables;i++) { /*The input box only has values for the input variables*/ if (cs->orig_systemVar[i]) { (*rp)[i]=pt[k]; k++; } else (*rp)[i]=0.0; /*Non system variables are set to 0. This value is not used */ } InitBoxFromPoint(cs->orig_nvariables,*rp,&b); if (!RegenerateSolution(p,&b,cs)) Error("Invalid solution point in RegenerateSolutionPoint"); k=0; for(i=0;i<cs->orig_nvariables;i++) { /*The input box only has values for the input variables (system and secondary)*/ if (!cs->orig_systemVar[i]) (*rp)[i]=IntervalCenter(GetBoxInterval(i,&b)); } DeleteBox(&b); return(cs->orig_nvariables); } 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("Inconsistent input cuiksystem"); } unsigned int RegenerateOriginalPoint(Tparameters *p,double *s,double **o,TCuikSystem *cs) { if (UpdateCuikSystem(p,cs)) { PointFromVariables(o,&(cs->orig_variables)); UpdateOriginalPointFromSimple(s,*o,cs->orig2s); } else Error("Inconsistent input cuiksystem"); return(cs->orig_nvariables); } unsigned int GenerateSimplifiedPoint(Tparameters *p,double *o,double **s,TCuikSystem *cs) { if (UpdateCuikSystem(p,cs)) SimplePointFromOriginal(o,s,cs->orig2s); else Error("Inconsistent input cuiksystem"); return(cs->simp_nvariables); } /* 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 (!cs->scalar) Error("ComputeSplitDim only operates on scalar systems"); 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,0.0,GetBoxInterval(i,b)); k++; } } return(in); } void EvaluateCSEquations(double *p,double *r,TCuikSystem *cs) { EvaluateEqualityEquations(FALSE,p,r,&(cs->orig_equations)); } void EvaluateSimpCSEquations(Tparameters *pr,double *p,double *r,TCuikSystem *cs) { if (!UpdateCuikSystem(pr,cs)) Error("Inconsistent input cuiksystem"); EvaluateEqualityEquations(FALSE,p,r,&(cs->simp_equations)); } void EvaluateSubSetSimpCSEquations(Tparameters *pr,boolean *se,double *p,double *r,TCuikSystem *cs) { if (!UpdateCuikSystem(pr,cs)) Error("Inconsistent input cuiksystem"); EvaluateSubSetEqualityEquations(p,se,r,&(cs->simp_equations)); } void EvaluateCSJacobian(double *p,double ***J,TCuikSystem *cs) { TJacobian Ja; InitJacobian(&(cs->orig_variables),&(cs->orig_equations),&Ja); AllocateJacobianEvaluation(J,&Ja); EvaluateJacobian(p,*J,&Ja); DeleteJacobian(&Ja); } double ErrorInCSEquations(double *p,TCuikSystem *cs) { double *r,e; unsigned int neq; neq=NEqualityEquations(&(cs->orig_equations)); if (neq>0) { NEW(r,neq,double); EvaluateEqualityEquations(FALSE,p,r,&(cs->orig_equations)); e=Norm(neq,r); free(r); } else e=0.0; return(e); } double ErrorInSimpCSEquations(Tparameters *pr,double *p,TCuikSystem *cs) { double *r,e; unsigned int neq; if (!UpdateCuikSystem(pr,cs)) Error("Inconsistent input cuiksystem"); neq=NEqualityEquations(&(cs->simp_equations)); if (neq>0) { NEW(r,neq,double); EvaluateEqualityEquations(FALSE,p,r,&(cs->simp_equations)); e=Norm(neq,r); free(r); } else e=0.0; return(e); } double EvaluateCSCost(Tparameters *p,boolean simp,double *s,void *cs) { double v; if (simp) v=EvaluateWholeEquation(s,&(((TCuikSystem *)cs)->eqMin)); else v=EvaluateWholeEquation(s,&(((TCuikSystem *)cs)->orig_eqMin)); return(v); } /* 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 *r,maxError; //double e; unsigned int nv,neq; maxError=0; nv=NVariables(&(cs->orig_variables)); neq=NEqualityEquations(&(cs->orig_equations)); if ((nv>0)&&(neq>0)) { if (nv!=GetBoxNIntervals(b)) Error("Wrong box size in ErrorInSolution"); NEW(p,nv,double); for(i=0;i<nv;i++) p[i]=IntervalCenter(GetBoxInterval(i,b)); NEW(r,neq,double); EvaluateEqualityEquations(TRUE,p,r,&(cs->orig_equations)); maxError=Norm(neq,r); /* for(i=0;i<neq;i++) { e=fabs(r[i]); if (e>maxError) maxError=e; } */ free(r); free(p); } return(maxError); } double ErrorInInequalities(Tbox *b,TCuikSystem *cs) { unsigned int i,neq,nv; double maxError,*r,*p; nv=NVariables(&(cs->orig_variables)); if (nv!=GetBoxNIntervals(b)) Error("Wrong box size in ErrorInInequalities"); neq=NInequalityEquations(&(cs->orig_equations)); if (neq>0) { NEW(p,nv,double); for(i=0;i<cs->orig_nvariables;i++) p[i]=IntervalCenter(GetBoxInterval(i,b)); NEW(r,neq,double); EvaluateInequalityEquations(p,r,&(cs->orig_equations)); maxError=MaxVector(neq,r); free(r); free(p); } else maxError=0.0; return(maxError); } boolean InequalitiesHoldOnPoint(double *p,TCuikSystem *cs) { unsigned int neq; boolean hold; double *r; neq=NInequalityEquations(&(cs->orig_equations)); if (neq>0) { NEW(r,neq,double); EvaluateInequalityEquations(p,r,&(cs->orig_equations)); hold=(MaxVector(neq,r)==0.0); free(r); } else hold=TRUE; return(hold); } boolean SimpInequalitiesHoldOnPoint(Tparameters *pr,double *p,TCuikSystem *cs) { unsigned int neq; boolean hold; double *r; if (!UpdateCuikSystem(pr,cs)) Error("Inconsistent input cuiksystem"); neq=NInequalityEquations(&(cs->simp_equations)); if (neq>0) { NEW(r,neq,double); EvaluateInequalityEquations(p,r,&(cs->simp_equations)); hold=(MaxVector(neq,r)==0.0); free(r); } else hold=TRUE; return(hold); } double ErrorInSimpInequalitiesOnPoint(Tparameters *pr,double *p,TCuikSystem *cs) { unsigned int neq; double e,*r; if (!UpdateCuikSystem(pr,cs)) Error("Inconsistent input cuiksystem"); neq=NInequalityEquations(&(cs->simp_equations)); if (neq>0) { NEW(r,neq,double); EvaluateInequalityEquations(p,r,&(cs->simp_equations)); e=MaxVector(neq,r); free(r); } else e=0.0; return(e); } /* 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))); fprintf(f,"%% Variable reduction %u -> %u\n",cs->orig_nvariables, cs->simp_nvariables); fprintf(f,"%% Num syst+secu variables in original : %u \n", GetNumSystemVariables(&(cs->orig_variables))+ GetNumSecondaryVariables(&(cs->orig_variables))); fprintf(f,"%% Num syst+secu variables in simplified: %u \n", GetNumSystemVariables(&(cs->simp_variables))+ GetNumSecondaryVariables(&(cs->simp_variables))); 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) { DeleteConstants(&(cs->constants)); DeleteVariables(&(cs->orig_variables)); DeleteEquations(&(cs->orig_equations)); DeleteStatistics(&(cs->st)); DeleteEquation(&(cs->orig_eqMin)); UnUpdateCuikSystem(cs); cs->empty=TRUE; cs->simp_empty=TRUE; }
1047 
1048  'c' is the return code out of Reduce Box
1049 
1050  If the box is small (max side size smaller than sigma), we consider it
1051  a solution and we store it in the list 'sol' (if defined), we save
1052  the solution into 'f_out' (if defined).
1053 
1054  Otherwise, we split the box in two sub-boxes and we store them in list 'boxes'
1055  The dimension along which the box is splitte is determined according to
1056  parameter 'error_split'. The two sub-boxes are added to 'boxes' in the first
1057  of in the last position according to parameter 'depth_first'.
1058 */
1060  unsigned int c,
1061  FILE *f_out,
1062  Tlist *sol,
1063  Theap *boxes,Tbox *b,
1064  TCuikSystem *cs)
1065 {
1066  boolean empty,small;
1067  Tbox b_orig;
1068  double sigma;
1069  unsigned int nNewton;
1070  double epsilon;
1071 
1072  empty=(c==EMPTY_BOX);
1073 
1074  if (empty)
1075  {
1076  NewEmptyBox(&(cs->st));
1077  #if (_DEBUG>1)
1078  printf(" The box is empty -> deleted\n");
1079  #endif
1080  #if (_DEBUG>0)
1081  fprintf(stderr,"e\n");
1082  #endif
1083  }
1084  else
1085  {
1086  /* c==REDUCED_BOX or c==REDUCED_BOX_WITH_SOLUTION or c==ERROR_IN_PROCESS*/
1087  #if (_DEBUG>1)
1088  printf(" Resulting box: ");
1089  PrintBox(stdout,b);
1090  #endif
1091  #if (_DEBUG>0)
1092  fprintf(stderr,"<%g %g>",
1093  GetBoxVolume(cs->systemVar,b),
1094  GetBoxSize(cs->systemVar,b));
1095  #endif
1096 
1097  /* We initialize a box in the original system and set the default ranges
1098  (ranges in case the variable in the original system is not in the
1099  simplified one... this is not likely to occur here but in the system
1100  to sybsystem relation but in this way the interface is the same)*/
1101  BoxFromVariables(&b_orig,&(cs->orig_variables));
1102  UpdateOriginalFromSimple(b,&b_orig,cs->orig2sd);
1103  #if (_DEBUG>1)
1104  printf(" Original box: ");
1105  PrintBox(stdout,&b_orig);
1106  #endif
1107 
1108  epsilon=GetParameter(CT_EPSILON,p);
1109 
1110  nNewton=(unsigned int)GetParameter(CT_MAX_NEWTON_ITERATIONS,p);
1111  if ((nNewton>0)&&(f_out!=NULL)&&(GetBoxSize(cs->systemVar,b)>epsilon))
1112  {
1113  double *sol;
1114  Tbox b_sol;
1115 
1116  NEW(sol,cs->orig_nvariables,double);
1117 
1118  if ((CuikNewtonInBox(p,&b_orig,sol,&b_sol,cs)&(CONVERGED_IN_BOX|CONVERGED_IN_GLOBAL))>0)
1119  {
1120  fprintf(f_out,"Newton [s:%f v:%g l:%u t:%g]: ",
1121  GetBoxSize(cs->systemVar,b),
1122  GetBoxVolume(cs->systemVar,b),
1123  GetBoxLevel(b),
1124  GetElapsedTime(&(cs->st)));
1125  PrintBoxSubset(f_out,cs->orig_notDummyVar,cs->orig_varNames,&b_sol);
1126  }
1127 
1128  free(sol);
1129  DeleteBox(&b_sol);
1130  fflush(f_out);
1131  }
1132 
1133  sigma=GetParameter(CT_SIGMA,p);
1134 
1135  /*The original system can have variables not included in the simplified one. For instance,
1136  variables not in the system equations due to the particular constants of the problem under
1137  analysis are eliminated, even though they appear in dummy equations (this is common
1138  when ROT_REP=1, see link.h). Those removed variables are never reduced and thus,
1139  the original box never reduces below sigma. This is why we use the simplified box
1140  (i.e., the one we actually reduce) to assess whether or not we have a solution. */
1141  small=(GetBoxSize(cs->systemVar,b)<=sigma);
1142 
1143  if (small)
1144  {
1145  #if (_DEBUG>1)
1146  printf(" The box is a solution (size %g vs [%g]) ->add to the list of solutions\n",
1147  GetBoxSize(cs->orig_systemVar,&b_orig),sigma);
1148  #endif
1149 
1151  GetBoxVolume(cs->orig_systemVar,&b_orig),
1152  GetBoxDiagonal(cs->orig_systemVar,&b_orig),
1153  &(cs->st));
1154 
1155  if(c==ERROR_IN_PROCESS)
1156  NewRBError(&(cs->st));
1157 
1158  #if (_DEBUG>0)
1160  fprintf(stderr,"S\n");
1161  else
1162  fprintf(stderr,"s\n");
1163  #endif
1164 
1165 
1166  if (f_out!=NULL)
1167  {
1168  if(c==ERROR_IN_PROCESS)
1169  fprintf(f_out," E");
1170 
1171  fprintf(f_out," %3u (err:%g",
1172  GetNSolutionBoxes(&(cs->st)),
1173  ErrorInSolution(&b_orig,cs));
1174  #if (_DEBUG>1)
1175  printf(" Solution number %u with erorr %g\n",
1176  GetNSolutionBoxes(&(cs->st)),
1177  ErrorInSolution(&b_orig,cs));
1178  #endif
1179 
1181  fprintf(f_out," min:%g",EvaluateEqMin((void *)b,(void *)cs));
1182 
1183  fprintf(f_out," tm:%g):",GetElapsedTime(&(cs->st)));
1184 
1186  fprintf(f_out," <S> ");
1187 
1188  /*The output does not include the dummy variables (they are
1189  not relevant any more)*/
1190  PrintBoxSubset(f_out,cs->orig_notDummyVar,cs->orig_varNames,&b_orig);
1191  fflush(f_out);
1192  }
1193  if (sol!=NULL)
1194  {
1195  Tbox *sb;
1196 
1197  NEW(sb,1,Tbox);
1198 
1199  CopyBox(sb,&b_orig);
1200 
1201  AddLastElement((void *)sb,sol);
1202  }
1203 
1204  /*If we want to ad box to a solution list, we have to copy 'b'
1205  somewhere else ('b' is delete after this function !)*/
1206  }
1207  else
1208  {
1209  /* Non empty large box or error in ReduceBox */
1210 
1211  Tbox b1,b2;
1212  unsigned int d;
1213 
1214  #if (_DEBUG>1)
1215  printf(" Box size : %g vs. %g\n",GetBoxSize(cs->orig_systemVar,&b_orig),sigma);
1216  printf(" The box has to be splitted (and deleted)\n");
1217  #endif
1218 
1219  NewSplittedBox(&(cs->st));
1220 
1221  d=ComputeSplitDimInt(p,b,cs);
1222 
1223  SplitBox(d,CUT_POINT,&b1,&b2,b);
1224 
1225  #if (_DEBUG>1)
1226  printf(" Splitting along dimension %u (internal)\n",d);
1228  printf(" Box penalgy : %g\n",EvaluateEqMin((void *)&b1,(void *)cs));
1229  printf(" SubBox:"); PrintBox(stdout,&b1);
1231  printf(" Box penalty : %g\n",EvaluateEqMin((void *)&b2,(void *)cs));
1232  printf(" SubBox:"); PrintBox(stdout,&b2);
1233  #endif
1234 
1235  AddBox2HeapOfBoxes(&b1,boxes);
1236  AddBox2HeapOfBoxes(&b2,boxes);
1237  DeleteBox(&b1);
1238  DeleteBox(&b2);
1239 
1240  if (c==ERROR_IN_PROCESS)
1241  {
1242  #if (_DEBUG>0)
1243  fprintf(stderr,"E[%u]\n",d);
1244  #endif
1245  #if (_DEBUG>1)
1246  printf(" Splitted due to simplex error\n");
1247  #endif
1248  NewRBError(&(cs->st));
1249  }
1250  #if (_DEBUG>0)
1251  else
1252  fprintf(stderr,"b[%u]\n",d);
1253  #endif
1254  }
1255 
1256  DeleteBox(&b_orig);
1257  }
1258 }
1259 
1261 {
1262  unsigned int i,vt;
1263  double epsilon;
1264 
1265  epsilon=GetParameter(CT_EPSILON,p);
1266 
1267  /*We proceed back to front because removing a variable change the
1268  indexes of variables above it.*/
1269  i=NVariables(&(cs->variables));
1270  while(i>0)
1271  {
1272  i--;
1273  vt=GetVariableTypeN(i,&(cs->variables));
1274 
1275  if ((!VarIncluded(i,GetEquationVariables(&(cs->orig_eqMin))))&&
1276  (((vt==DUMMY_VAR)&&(!UsedVarInNonDummyEquations(i,&(cs->equations))))||
1277  ((vt!=SYSTEM_VAR)&&(!UsedVarInEquations(i,&(cs->equations))))))
1278  {
1279  /* ... Therefore, we can still have dummy equations with the non-used variables*/
1280  RemoveEquationsWithVar(epsilon,i,&(cs->equations));
1281  RemoveVariable(i,&(cs->variables));
1282  }
1283  }
1284 }
1285 
1287  boolean *replaced,TLinearConstraint *lc,
1288  Tbox *borig,TCuikSystem *cs)
1289 {
1290  boolean found,hold;
1291  unsigned int nv,origID1;
1292  unsigned int id1;
1293  char *name1;
1294  double epsilon;
1295  double ct2;
1296  Tinterval *range;
1297 
1298  epsilon=GetParameter(CT_EPSILON,p);
1299 
1300  hold=TRUE;
1301 
1302  do {
1303  found=FALSE;
1304 
1305  /* A variable with a ct range can be removed */
1306  nv=NVariables(&(cs->variables));
1307  id1=0;
1308  while((!found)&&(id1<nv))
1309  {
1310  name1=GetVariableName(GetVariable(id1,&(cs->variables)));
1311  origID1=GetVariableID(name1,&(cs->orig_variables));
1312  if (origID1==NO_UINT)
1313  Error("Undefined ID1 in CSRemoveVarsWithCtRange");
1314 
1315  range=GetBoxInterval(origID1,borig);
1316  if (IntervalSize(range)==0)
1317  {
1318  ct2=LowerLimit(range); /* ==UpperLimit(range) */
1319  found=TRUE;
1320  }
1321  else
1322  id1++;
1323  }
1324 
1325  /*if we found something to simplify, apply it to the rest of equations
1326  (this creates a new set of equations that replaces the previous one)*/
1327  if (found)
1328  {
1329  replaced[origID1]=TRUE;
1330  AddCt2LinearConstraint(ct2,&(lc[origID1]));
1331 
1332  #if (_DEBUG>5)
1333  printf("Ct Range Replacement: %s [%u-%u] ->%.12g\n",name1,origID1,id1,ct2);
1334  #endif
1335 
1336  if (hold)
1337  {
1338  TLinearConstraint lc;
1339 
1340  InitLinearConstraint(&lc);
1341  AddCt2LinearConstraint(ct2,&lc);
1342  hold=ReplaceVariableInEquations(epsilon,id1,&lc,&(cs->equations));
1344  RemoveVariable(id1,&(cs->variables));
1345 
1346  #if (_DEBUG>6)
1347  if(hold)
1348  {
1349  char **varNames;
1350 
1351 
1352  printf("%%=======================================================\n");
1353  printf("%% One less variable (ct range)\n");
1354  printf("%%****************************************\n");
1355  PrintVariables(stdout,&(cs->variables));
1356 
1357  nv=NVariables(&(cs->variables));
1358  NEW(varNames,nv,char *);
1359  GetVariableNames(varNames,&(cs->variables));
1360  PrintEquations(stdout,varNames,&(cs->equations));
1361  free(varNames);
1362  }
1363  #endif
1364  }
1365  }
1366  } while((hold)&&(found));
1367 
1368  return(hold);
1369 }
1370 
1371 boolean CSRemoveLCVars(Tparameters *p,unsigned int simplificationLevel,
1372  unsigned int nTerms,boolean *changed,
1373  boolean *replaced,TLinearConstraint *lc,
1374  Tbox *borig,TCuikSystem *cs)
1375 {
1376  boolean found,hold;
1377  unsigned int m,i,j,nv,neq,origID,origID1;
1378  Tequation *eq;
1379  unsigned int id;
1380  char *name,*name1;
1381  double epsilon;
1382  Tvariable *v;
1383  TLinearConstraint lcr,lcc;
1384  boolean *systemVars;
1385  Tinterval error;
1386  Tbox currentBox;
1387  #if (_DEBUG>5)
1388  char **varNamesO;
1389 
1390  nv=NVariables(&(cs->orig_variables));
1391  NEW(varNamesO,nv,char *);
1392  GetVariableNames(varNamesO,&(cs->orig_variables));
1393  #endif
1394 
1395  nv=NVariables(&(cs->variables));
1396  InitBox(nv,NULL,&currentBox);
1397  NEW(systemVars,nv,boolean);
1398  for(i=0;i<nv;i++)
1399  {
1400  systemVars[i]=IsSystemVariable(i,&(cs->variables));
1401  SetBoxInterval(i,GetVariableInterval(GetVariable(i,&(cs->variables))),&currentBox);
1402  }
1403 
1404  epsilon=GetParameter(CT_EPSILON,p);
1405 
1406  *changed=FALSE;
1407  hold=TRUE;
1408 
1409  found=FALSE;
1410  i=0;
1411  neq=NEquations(&(cs->equations));
1412  while((!found)&&(i<neq))
1413  {
1414  eq=GetEquation(i,&(cs->equations));
1415 
1416  found=IsSimplificable(simplificationLevel,nTerms,systemVars,&currentBox,&id,&lcr,eq); /* 't' is the type of simplification*/
1417  if (!found)
1418  i++;
1419  }
1420  DeleteBox(&currentBox);
1421 
1422  /*if we found something to simplify, apply it to the rest of equations
1423  (this creates a new set of equations that replaces the previous one)*/
1424  if (found)
1425  {
1426  *changed=TRUE;
1427 
1428  /* Translate the variables identifiers from local to global */
1429  v=GetVariable(id,&(cs->variables));
1430  name=GetVariableName(v);
1431  origID=GetVariableID(name,&(cs->orig_variables));
1432  if (origID==NO_UINT)
1433  Error("Undefined var in original system in CSRemoveLCVars");
1434 
1435  /* Translate the 'lcr' linear constraint that refers to the cuiksystem
1436  in the current state of simplification to the original system */
1438  for(i=0;i<m;i++)
1439  {
1441  name1=GetVariableName(v);
1442  origID1=GetVariableID(name1,&(cs->orig_variables));
1443  if (origID1==NO_UINT)
1444  Error("Undefined variable in original system in CSRemoveLCVars");
1445  AddTerm2LinearConstraint(origID1,
1447  &(lc[origID]));
1448  }
1449  GetLinearConstraintError(&error,&lcr);
1450  if (IntervalSize(&error)>ZERO)
1451  Error("Linear constraint used for variable replacement must have ct right-hand");
1452  AddCt2LinearConstraint(-IntervalCenter(&error),&(lc[origID]));
1453 
1454  replaced[origID]=TRUE;
1455 
1456  #if (_DEBUG>5)
1457  printf("Var Replacement: %s[%u-%u]->",name,origID,id);
1458  PrintLinearConstraint(stdout,FALSE,varNamesO,&(lc[origID]));
1459  #endif
1460 
1461  CopyLinearConstraint(&lcc,&(lc[origID]));
1462  AddTerm2LinearConstraint(origID,-1.0,&lcc);
1463  hold=CropLinearConstraint(ZERO,ANY_TYPE_VAR,borig,&(cs->orig_variables),&lcc);
1464  if (!hold)
1465  printf("\nNon intersecting ranges (I)\n");
1466  DeleteLinearConstraint(&lcc);
1467 
1468  /*If origID was used for the replacement of another variable
1469  we have to propagate the replacement of origID1 to this other
1470  variable*/
1471  nv=NVariables(&(cs->orig_variables));
1472  for(j=0;((hold)&&(j<nv));j++)
1473  {
1474  if ((replaced[j])&&(LinearConstraintIncludes(origID,&(lc[j]))))
1475  {
1476  TLinearConstraint ltmp;
1477  double ct;
1478 
1479  ct=RemoveTermFromLinearConstraint(origID,&(lc[j]));
1480  CopyLinearConstraint(&ltmp,&(lc[origID]));
1481  ScaleLinearConstraint(ct,&ltmp);
1482  AddLinearConstraints(&ltmp,&(lc[j]));
1483  DeleteLinearConstraint(&ltmp);
1484 
1485  CopyLinearConstraint(&lcc,&(lc[j]));
1486  AddTerm2LinearConstraint(j,-1.0,&lcc);
1487  hold=CropLinearConstraint(ZERO,ANY_TYPE_VAR,borig,&(cs->orig_variables),&lcc);
1488  if (!hold)
1489  printf("\nNon intersecting ranges (II)\n");
1490  DeleteLinearConstraint(&lcc);
1491  }
1492  }
1493 
1494  if (hold)
1495  {
1496  hold=ReplaceVariableInEquations(epsilon,id,&lcr,&(cs->equations));
1497  RemoveVariable(id,&(cs->variables));
1498 
1499  #if (_DEBUG>6)
1500  if(hold)
1501  {
1502  char **varNames;
1503 
1504  printf("%%=======================================================\n");
1505  printf("%% One less variable (replaced)\n");
1506  printf("%%****************************************\n");
1507  PrintVariables(stdout,&(cs->variables));
1508 
1509  nv=NVariables(&(cs->variables));
1510  NEW(varNames,nv,char *);
1511  GetVariableNames(varNames,&(cs->variables));
1512  PrintEquations(stdout,varNames,&(cs->equations));
1513  free(varNames);
1514  }
1515  #endif
1516  }
1517 
1518  DeleteLinearConstraint(&lcr);
1519  }
1520  #if (_DEBUG>5)
1521  free(varNamesO);
1522  #endif
1523  free(systemVars);
1524  return(hold);
1525 }
1526 
1527 /*
1528  Simplifies a CuikSystem removing variables have constant value and those
1529  that depend linearly on other variables. Moreover, a primitive (but numerically
1530  save) form of Gaussian elimination is performed to remove as more variables
1531  and equations as possible.
1532  More elaborate forms of system simplifications are not implemented because they
1533  result in numerical errors (that can change the system solutions).
1534  The output mappings give information about how the original and the simplified
1535  systems are related.
1536 
1537  This is part of the UpdateCuikSystem.
1538 */
1540 {
1541  unsigned int i,nv,nvn,newVarID;
1542  boolean *replaced;
1543  TLinearConstraint *lc,lcS;
1544  boolean changed,anyChange,hold;
1545  char *name1;
1546  Tbox borig,bsimp;
1547  unsigned int nTerms,m,k;
1548  Tinterval error;
1549  unsigned int simplificationLevel;
1550 
1551  if (!cs->scalar)
1552  simplificationLevel=0; /* Matrix equations -> no simplification */
1553  else
1554  {
1555  simplificationLevel=(unsigned int)(GetParameter(CT_SIMPLIFICATION_LEVEL,p));
1556  if ((!PolynomialEquations(&(cs->orig_equations)))&&(simplificationLevel>1))
1557  simplificationLevel=1;
1558  }
1559 
1560  #if (_DEBUG>6)
1561  {
1562  char **varNames;
1563 
1564  printf("%%=======================================================\n");
1565  printf("%% Original system\n");
1566  printf("%%****************************************\n");
1567  PrintVariables(stdout,&(cs->orig_variables));
1568 
1569  nv=NVariables(&(cs->orig_variables));
1570  NEW(varNames,nv,char *);
1571  GetVariableNames(varNames,&(cs->orig_variables));
1572  PrintEquations(stdout,varNames,&(cs->orig_equations));
1573  free(varNames);
1574  }
1575  #endif
1576 
1577 
1578  /* Get a copy of the original cuiksystem */
1579  CopyVariables(&(cs->variables),&(cs->orig_variables));
1580  CopyEquations(&(cs->equations),&(cs->orig_equations));
1581 
1582  nv=NVariables(&(cs->orig_variables)); /* number of variables in the original system */
1583 
1584  NEW(replaced,nv,boolean);
1585  NEW(lc,nv,TLinearConstraint);
1586 
1587  BoxFromVariables(&borig,&(cs->orig_variables)); /*original ranges to be possibly reduced during
1588  simplification */
1589  for(i=0;i<nv;i++)
1590  {
1591  replaced[i]=FALSE;
1592  InitLinearConstraint(&(lc[i]));
1593  }
1594 
1595  hold=TRUE; /* in the simplification process we can discover that the input system
1596  is inconsitent (i.e., that it does not hold)*/
1597 
1598  CSRemoveUnusedVars(p,cs); /* This is done even without simplification */
1599 
1600  if (simplificationLevel>0) /* For detailed debug systems are not simplified */
1601  {
1602  hold=CSRemoveVarsWithCtRange(p,replaced,lc,&borig,cs);
1603 
1604  anyChange=TRUE; /* Set to TRUE to trigger the simplification process */
1605  while((hold)&&(anyChange))
1606  {
1607  anyChange=FALSE; /* Will be set to TRUE if at least one variable is removed
1608  in the loop below */
1609  nTerms=0; /* First we try to remove ct variables, then scaled ones,i.e., we
1610  progressivelly increase the number of terms used in the
1611  replacements. */
1612  while((hold)&&(nTerms<=MAX_TERMS_SIMP)) /*hard-coded maximum complexity of the simplifications*/
1613  {
1614  hold=CSRemoveLCVars(p,simplificationLevel,nTerms,&changed,replaced,lc,&borig,cs);
1615  anyChange=((anyChange)||(changed));
1616  if (changed)
1617  nTerms=0;
1618  else
1619  nTerms++;
1620  }
1621 
1622  /* At this point no more ct or scaled varibles can be removed. Try to
1623  simplify the problem via Gaussian elimination*/
1624  if ((hold)&&(anyChange))
1625  {
1626  anyChange=FALSE;
1627  do {
1628  changed=GaussianElimination(&(cs->equations));
1629  anyChange=((anyChange)||(changed));
1630  } while (changed);
1631 
1632  #if (_DEBUG>6)
1633  {
1634  char **varNames;
1635  unsigned int nvg;
1636 
1637  printf("%%=======================================================\n");
1638  printf("%% After gaussian elimination\n");
1639  printf("%%****************************************\n");
1640  PrintVariables(stdout,&(cs->variables));
1641 
1642  nvg=NVariables(&(cs->variables));
1643  NEW(varNames,nvg,char *);
1644  GetVariableNames(varNames,&(cs->variables));
1645  PrintEquations(stdout,varNames,&(cs->equations));
1646  free(varNames);
1647  }
1648  #endif
1649  }
1650  }
1651  }
1652 
1653  #if(_DEBUG>4)
1654  if (hold)
1655  {
1656  char **varNames;
1657  unsigned int nvg;
1658 
1659  printf("%%=======================================================\n");
1660  printf("%%FINAL SIMPLIFIED SYSTEM\n");
1661  printf("%%****************************************\n");
1662  PrintVariables(stdout,&(cs->variables));
1663 
1664  nvg=NVariables(&(cs->variables));
1665  NEW(varNames,nvg,char *);
1666  GetVariableNames(varNames,&(cs->variables));
1667  PrintEquations(stdout,varNames,&(cs->equations));
1668  free(varNames);
1669  }
1670  else
1671  printf("Inconsistent input system\n");
1672  #endif
1673 
1674  if (hold)
1675  {
1676  /* Get a copy of the simplified but not still dummified cuiksystem */
1677  CopyVariables(&(cs->simp_variables),&(cs->variables));
1678  CopyEquations(&(cs->simp_equations),&(cs->equations));
1679 
1680  /*The final step is to dummify the cuiksystem*/
1681  DummifyCuikSystem(p,cs);
1682 
1683  /*define the mappings from the gathered information*/
1684 
1685  NEW(cs->orig2sd,1,Tmapping);
1686  NEW(cs->orig2s,1,Tmapping);
1687 
1688  InitMapping(&(cs->orig_variables),&(cs->variables),cs->orig2sd);
1689  InitMapping(&(cs->orig_variables),&(cs->simp_variables),cs->orig2s);
1690 
1691  /* now we complement the default initialization with the replacements computed above */
1692  for(i=0;i<nv;i++)
1693  {
1694  if (replaced[i])
1695  {
1696  m=GetNumTermsInLinearConstraint(&(lc[i]));
1697  InitLinearConstraint(&lcS);
1698  for(k=0;k<m;k++)
1699  {
1701  newVarID=GetVariableID(name1,&(cs->variables));
1702 
1703  AddTerm2LinearConstraint(newVarID,GetLinearConstraintCoefficient(k,&(lc[i])),&lcS);
1704  }
1705  GetLinearConstraintError(&error,&(lc[i]));
1706  AddCt2LinearConstraint(-IntervalCenter(&error),&lcS);
1707 
1708  SetOriginalVarRelation(i,&lcS,cs->orig2sd);
1709  SetOriginalVarRelation(i,&lcS,cs->orig2s);
1710 
1711  DeleteLinearConstraint(&lcS);
1712  }
1713  }
1714 
1715  #if(_DEBUG>4)
1716  if (hold)
1717  {
1718  char **varNames;
1719  unsigned int nvg;
1720 
1721  printf("%%=======================================================\n");
1722  printf("%%FINAL fianl SIMPLIFIED SYSTEM\n");
1723  printf("%%****************************************\n");
1724  PrintVariables(stdout,&(cs->variables));
1725 
1726  nvg=NVariables(&(cs->variables));
1727  NEW(varNames,nvg,char *);
1728  GetVariableNames(varNames,&(cs->variables));
1729  PrintEquations(stdout,varNames,&(cs->equations));
1730  free(varNames);
1731  }
1732  else
1733  printf("Inconsistent input system\n");
1734  #endif
1735 
1736  SimpleFromOriginal(&borig,&bsimp,cs->orig2s);
1737  nvn=NVariables(&(cs->simp_variables));
1738  for(i=0;i<nvn;i++)
1740  DeleteBox(&bsimp);
1741 
1742  SimpleFromOriginal(&borig,&bsimp,cs->orig2sd);
1743  nvn=NVariables(&(cs->variables));
1744  for(i=0;i<nvn;i++)
1746  DeleteBox(&bsimp);
1747  }
1748 
1749  DeleteBox(&borig);
1750 
1751  /* free the information used for the mapping */
1752  for(i=0;i<nv;i++)
1753  DeleteLinearConstraint(&lc[i]);
1754  free(replaced);
1755  free(lc);
1756 
1757  return(hold);
1758 }
1759 
1760 /*
1761  Removes the information stored in the cuiksystem during the update
1762 */
1764 {
1765 
1766  if (cs->updated)
1767  {
1768  /*Removes cached information that is used during solving*/
1769 
1770  /*Firt the box sorting information*/
1772  DeleteEquation(&(cs->eqMin));
1773 
1774  /*Now the information about the simplified+dummified system*/
1775  if (cs->orig2sd!=NULL)
1776  {
1777  DeleteMapping(cs->orig2sd);
1778  free(cs->orig2sd);
1779  cs->orig2sd=NULL;
1780  }
1781 
1782  cs->nequations=0;
1783  cs->nvariables=0;
1784  if (cs->systemVar!=NULL)
1785  free(cs->systemVar);
1786  cs->systemVar=NULL;
1787  if (cs->notDummyVar!=NULL)
1788  free(cs->notDummyVar);
1789  cs->notDummyVar=NULL;
1790  if (cs->varType!=NULL)
1791  free(cs->varType);
1792  cs->varType=NULL;
1793 
1794  DeleteEquations(&(cs->equations));
1795  DeleteVariables(&(cs->variables));
1796 
1797  /*Now remove the information about the simplified system*/
1798  if (cs->orig2s!=NULL)
1799  {
1800  DeleteMapping(cs->orig2s);
1801  free(cs->orig2s);
1802  cs->orig2s=NULL;
1803  }
1804 
1805  if (cs->simp_nequations>0)
1806  DeleteJacobian(&(cs->J));
1807 
1808  cs->simp_nequations=0;
1809  cs->simp_nvariables=0;;
1810  cs->simp_nee=0;
1811 
1814 
1815  if (cs->simp_tp!=NULL)
1816  free(cs->simp_tp);
1817 
1818  /*Finally remove the information about the original system*/
1819  cs->orig_nequations=0;
1820  cs->orig_nvariables=0;
1821  if (cs->orig_systemVar!=NULL)
1822  free(cs->orig_systemVar);
1823  cs->orig_systemVar=NULL;
1824  if (cs->orig_notDummyVar!=NULL)
1825  free(cs->orig_notDummyVar);
1826  cs->orig_notDummyVar=NULL;
1827 
1828  free(cs->orig_varNames);
1829  cs->orig_varNames=NULL;
1830 
1831  /*Mark the system as not updated*/
1832  cs->updated=FALSE;
1833 
1834  /*We assume the system to be consistent although this value is not
1835  used unless cs->update is TRUE*/
1836  cs->consistent=TRUE;
1837  }
1838 }
1839 
1840 /*
1841  Simplifies the CuikSystem and caches information that is useful during
1842  solving (cached->faster access).
1843 */
1845 {
1846  if (!cs->updated)
1847  {
1848  unsigned int j;
1849  boolean found;
1850 
1851  /* we update even if we have no equations to be able to deal
1852  with systems without equations (non constrained problems)
1853  In this case cuik is a standard approximated cell decomposition
1854  approach
1855  */
1856  cs->consistent=SimplifyCuikSystem(p,cs);
1857 
1858  if (cs->consistent)
1859  {
1860  /* Update the information in the original system */
1863 
1864  if ((cs->orig_nequations>0)&&(cs->orig_nvariables==0))
1865  Error("System with equations but without variables");
1866 
1867  if (cs->orig_nvariables==0)
1868  {
1869  cs->orig_systemVar=NULL;
1870  cs->orig_notDummyVar=NULL;
1871  cs->orig_varNames=NULL;
1872  }
1873  else
1874  {
1875  NEW(cs->orig_systemVar,cs->orig_nvariables,boolean);
1876  NEW(cs->orig_notDummyVar,cs->orig_nvariables,boolean);
1877 
1878  /* Note that orig_systemVar include the secodary variables while
1879  systemVar (the one for the simplified/dummified system) keeps
1880  the difference between system and secondary varibles */
1881  for(j=0;j<cs->orig_nvariables;j++)
1882  {
1883  cs->orig_systemVar[j]=((IsSystemVariable(j,&(cs->orig_variables)))||
1884  (IsSecondaryVariable(j,&(cs->orig_variables))));
1886  }
1887 
1888  NEW(cs->orig_varNames,cs->orig_nvariables,char *);
1890  }
1891 
1892  /* Update the simplifed cuiksystem */
1896  cs->simp_empty=((cs->simp_nvariables==0)&&(cs->simp_nequations==0));
1897 
1898  if (cs->simp_empty)
1899  cs->simp_tp=NULL;
1900  else
1901  {
1903  found=FALSE;
1904  j=0;
1905  while ((!found)&&(j<cs->simp_nvariables))
1906  {
1907  found=(cs->simp_tp[j]!=TOPOLOGY_R);
1908  j++;
1909  }
1910  if (!found)
1911  {
1912  free(cs->simp_tp);
1913  cs->simp_tp=NULL;
1914  }
1915  }
1916 
1917  if (cs->simp_nequations>0)
1918  InitJacobian(&(cs->simp_variables),&(cs->simp_equations),&(cs->J));
1919 
1920  /* Update the information in the simplified+dummified cuiksystem*/
1921  cs->nequations=NEquations(&(cs->equations));
1922  cs->nvariables=NVariables(&(cs->variables));
1923  cs->empty=((cs->nvariables==0)||(cs->nequations==0));
1924 
1925  if (cs->nvariables==0)
1926  {
1927  cs->systemVar=NULL;
1928  cs->notDummyVar=NULL;
1929  cs->varType=NULL;
1930  }
1931  else
1932  {
1933  NEW(cs->systemVar,cs->nvariables,boolean);
1934  NEW(cs->notDummyVar,cs->nvariables,boolean);
1935  NEW(cs->varType,cs->nvariables,unsigned int);
1936 
1937  for(j=0;j<cs->nvariables;j++)
1938  {
1939  cs->systemVar[j]=IsSystemVariable(j,&(cs->variables));
1940  cs->notDummyVar[j]=!IsDummyVariable(j,&(cs->variables));
1941  cs->varType[j]=GetVariableTypeN(j,&(cs->variables));
1942  }
1943  }
1944 
1945  /* Update the information about the sorting criteria for boxes */
1947  {
1948  if (GetNumMonomials(&(cs->orig_eqMin))==0)
1949  {
1950  ResetEquation(&(cs->orig_eqMin));
1952  }
1953  else
1955  cs->orig2sd,&(cs->eqMin),&(cs->orig_eqMin));
1956  }
1957 
1958  /* Mark the system as updated )(althogh it might not be consistent) */
1959  cs->updated=TRUE;
1960  }
1961  }
1962 
1963  return(cs->consistent);
1964 }
1965 
1966 /*
1967  Selects one dimension along which to split box 'b'. The dimension can
1968  be selected according to the size or according to the error that each
1969  variable induce in each equation.
1970  The output is a index in the simplified system
1971  */
1973 {
1974 
1975  if (!UpdateCuikSystem(p,cs))
1976  Error("Inconsistent cuiksystem in ComputeSplitDimInt");
1977 
1978  if (cs->nvariables==0)
1979  {
1980  Error("Splitting an empty cuiksystem");
1981  return(NO_UINT);
1982  }
1983  else
1984  {
1985  unsigned int i;
1986  Tinterval *is;
1987  double *splitDim;
1988  unsigned int *d;
1989  unsigned int n;
1990  double m,s;
1991  unsigned int splitType;
1992  double epsilon;
1993 
1994  epsilon=GetParameter(CT_EPSILON,p);
1995  splitType=(unsigned int)(GetParameter(CT_SPLIT_TYPE,p));
1996 
1997  NEW(splitDim,cs->nvariables,double);
1998  NEW(d,cs->nvariables,unsigned int);
1999 
2000  is=GetBoxIntervals(b);
2001 
2002  if ((cs->nequations>0)&&(splitType==1))
2003  {
2004  for(i=0;i<cs->nvariables;i++)
2005  splitDim[i]=0.0;
2006 
2007  /*Add the error contribution*/
2008  for(i=0;i<cs->nequations;i++)
2009  UpdateSplitWeight(i,splitDim,b,&(cs->equations));
2010 
2011  /*Do not split along unbounded or too small variables*/
2012  for(i=0;i<cs->nvariables;i++)
2013  {
2014  s=IntervalSize(&(is[i]));
2015  if ((s<10*epsilon)||(s>=INF))
2016  splitDim[i]=0.0;
2017  }
2018  }
2019  else
2020  {
2021  if (splitType!=2)
2022  {
2023  /*size-based split or error-based split without any equation*/
2024 
2025  for(i=0;i<cs->nvariables;i++)
2026  {
2027  s=IntervalSize(&(is[i]));
2028  if (s<INF)
2029  splitDim[i]=s;
2030  else
2031  splitDim[i]=0.0;
2032  }
2033  }
2034  else
2035  {
2036  /*splitType==2 -> random selection of the split dimension */
2037  for(i=0;i<cs->nvariables;i++)
2038  {
2039  s=IntervalSize(&(is[i]));
2040  if ((s>10*epsilon)&&(s<INF))
2041  splitDim[i]=1.0;
2042  else
2043  splitDim[i]=0.0;
2044  }
2045 
2046  }
2047  }
2048 
2049  m=-1.0;
2050  n=0;
2051  for(i=0;i<cs->nvariables;i++)
2052  {
2053  if (cs->systemVar[i])
2054  {
2055  if (splitDim[i]>m)
2056  {
2057  m=splitDim[i];
2058  d[0]=i;
2059  n=1;
2060  }
2061  else
2062  {
2063  if (splitDim[i]==m)
2064  {
2065  d[n]=i;
2066  n++;
2067  }
2068  }
2069  }
2070  }
2071 
2072  if (n>0)
2073  {
2074  #if (RANDOMNESS)
2075  i=d[randomMax(n-1)];
2076  #else
2077  i=d[0];
2078  #endif
2079  }
2080  else
2081  {
2082  Warning("There is no way where to bisect the box");
2083  i=randomMax(cs->nvariables-1);
2084  }
2085 
2086  free(d);
2087  free(splitDim);
2088 
2089  return(i);
2090  }
2091 }
2092 
2093 void SaveCSState(char *fname,Tlist *lb,TCuikSystem *cs)
2094 {
2095  FILE *f;
2096  char *tmpName;
2097  double tm;
2098 
2099  /*We initially save to a temporary file and then we rename it
2100  This increases the atomiticy of this operation (if the state
2101  file is corrupted the recovery operation would fail) */
2102  NEW(tmpName,strlen(fname)+10,char);
2103  sprintf(tmpName,"%s.tmp",fname);
2104 
2105  /*open in binary mode*/
2106  f=fopen(tmpName,"wb");
2107  if (!f)
2108  Error("Could not open file for saving CSState");
2109 
2110  tm=GetElapsedTime(&(cs->st));
2111  fwrite(&tm,sizeof(double),1,f);
2112  SaveStatistics(f,&(cs->st));
2113  SaveListOfBoxes(f,lb);
2114 
2115  fclose(f);
2116 
2117  remove(fname);
2118  rename(tmpName,fname);
2119  free(tmpName);
2120 }
2121 
2122 void LoadCSState(char *fname,Tlist *lb,TCuikSystem *cs)
2123 {
2124  FILE *f;
2125  double tm;
2126  Titerator it;
2127 
2128  f=fopen(fname,"rb");
2129  if (!f)
2130  Error("Could not open file for loading CSState");
2131 
2132  fread(&tm,sizeof(double),1,f);
2133  LoadStatistics(f,&(cs->st));
2134  LoadListOfBoxes(f,lb);
2135 
2136  fclose(f);
2137 
2138  InitIterator(&it,lb);
2139  First(&it);
2140  if (cs->nvariables!=GetBoxNIntervals((Tbox *)GetCurrent(&it)))
2141  Error("The saved .state and the problem dimensionality do not match");
2142 
2143  /*Pretend that the process was started tm seconds before now*/
2144  SetInitialTime(GetTime(&(cs->st))-tm,&(cs->st));
2145 }
2146 
2147 
2148 /************************************************************************/
2149 /************************************************************************/
2150 /************************************************************************/
2151 
2152 /*
2153  Warms up the cuiksystem structure
2154 */
2156 {
2157  InitConstants(&(cs->constants));
2158 
2159  InitVariables(&(cs->orig_variables));
2160  InitEquations(&(cs->orig_equations));
2161 
2162  /* cs->variables and cs->equations are initialized
2163  via Copy when simplifying the system*/
2164  cs->empty=TRUE;
2165  cs->simp_empty=TRUE;
2166 
2167  cs->updated=FALSE;
2168  cs->consistent=TRUE; /*this is not used unless updated=TRUE*/
2169 
2170  cs->scalar=TRUE;
2171 
2172  cs->nequations=0;
2173  cs->nvariables=0;
2174 
2175  cs->systemVar=NULL;
2176  cs->notDummyVar=NULL;
2177  cs->varType=NULL;
2178 
2179  cs->orig2sd=NULL;
2180 
2181  cs->orig2s=NULL;
2182  cs->simp_tp=NULL;
2183 
2184  cs->orig_nequations=0;
2185  cs->orig_nvariables=0;
2186 
2187  cs->orig_systemVar=NULL;
2188  cs->orig_notDummyVar=NULL;
2189 
2191  InitEquation(&(cs->orig_eqMin));
2192 }
2193 
2195 {
2196  if (!UpdateCuikSystem(p,cs))
2197  Error("Inconsistent input cuiksystem");
2198 }
2199 
2200 /*
2201  Copies one cuiksystem structure into another (that is assumed as
2202  initialized but empty).
2203 */
2205 {
2206  CopyConstants(&(cs_dst->constants),&(cs_src->constants));
2207 
2208  cs_dst->updated=cs_src->updated;
2209  cs_dst->consistent=cs_src->consistent;
2210  cs_dst->empty=cs_src->empty;
2211  cs_dst->scalar=cs_src->scalar;
2212 
2213  CopyStatistics(&(cs_dst->st),&(cs_src->st));
2214 
2215  cs_dst->searchMode=cs_src->searchMode;
2216  CopyEquation(&(cs_dst->orig_eqMin),&(cs_src->orig_eqMin));
2217 
2218  CopyEquations(&(cs_dst->orig_equations),&(cs_src->orig_equations));
2219  CopyVariables(&(cs_dst->orig_variables),&(cs_src->orig_variables));
2220 
2221  if (cs_dst->updated)
2222  {
2223  /*The simplified dummified system*/
2224  NEW(cs_dst->orig2sd,1,Tmapping);
2225  CopyMapping(cs_dst->orig2sd,cs_src->orig2sd);
2226 
2227  CopyEquations(&(cs_dst->equations),&(cs_src->equations));
2228  CopyVariables(&(cs_dst->variables),&(cs_src->variables));
2229 
2230  cs_dst->nequations=cs_src->nequations;
2231  cs_dst->nvariables=cs_src->nvariables;
2232 
2233  if (cs_dst->nvariables==0)
2234  {
2235  cs_dst->systemVar=NULL;
2236  cs_dst->notDummyVar=NULL;
2237  cs_dst->varType=NULL;
2238  }
2239  else
2240  {
2241  NEW(cs_dst->systemVar,cs_dst->nvariables,boolean);
2242  memcpy(cs_dst->systemVar,cs_src->systemVar,sizeof(boolean)*cs_dst->nvariables);
2243 
2244  NEW(cs_dst->notDummyVar,cs_dst->nvariables,boolean);
2245  memcpy(cs_dst->notDummyVar,cs_src->notDummyVar,sizeof(boolean)*cs_dst->nvariables);
2246 
2247  NEW(cs_dst->varType,cs_dst->nvariables,unsigned int);
2248  memcpy(cs_dst->varType,cs_src->varType,sizeof(unsigned int)*cs_dst->nvariables);
2249  }
2250 
2251  if (cs_dst->searchMode==MINIMIZATION_SEARCH)
2252  CopyEquation(&(cs_dst->eqMin),&(cs_src->eqMin));
2253 
2254 
2255  /*The simplified system*/
2256  NEW(cs_dst->orig2s,1,Tmapping);
2257  CopyMapping(cs_dst->orig2s,cs_src->orig2s);
2258 
2259  cs_dst->simp_empty=cs_src->simp_empty;
2260 
2261  CopyEquations(&(cs_dst->simp_equations),&(cs_src->simp_equations));
2262  CopyVariables(&(cs_dst->simp_variables),&(cs_src->simp_variables));
2263 
2264  cs_dst->simp_nequations=cs_src->simp_nequations;
2265  cs_dst->simp_nvariables=cs_src->simp_nvariables;
2266  cs_dst->simp_nee=cs_src->simp_nee;
2267 
2268  if (cs_src->simp_tp!=NULL)
2269  {
2270  NEW(cs_dst->simp_tp,cs_dst->simp_nvariables,unsigned int);
2271  memcpy(cs_dst->simp_tp,cs_src->simp_tp,sizeof(unsigned int)*cs_dst->simp_nvariables);
2272  }
2273  else
2274  cs_dst->simp_tp=NULL;
2275 
2276  CopyJacobian(&(cs_dst->J),&(cs_src->J));
2277 
2278  /*Cached elements from the original system*/
2279  cs_dst->orig_nequations=cs_src->orig_nequations;
2280  cs_dst->orig_nvariables=cs_src->orig_nvariables;
2281 
2282  if (cs_dst->orig_nvariables==0)
2283  {
2284  cs_dst->orig_systemVar=NULL;
2285  cs_dst->orig_notDummyVar=NULL;
2286  cs_dst->orig_varNames=NULL;
2287  }
2288  else
2289  {
2290  NEW(cs_dst->orig_systemVar,cs_dst->orig_nvariables,boolean);
2291  memcpy(cs_dst->orig_systemVar,cs_src->orig_systemVar,cs_dst->orig_nvariables*sizeof(boolean));
2292 
2293  NEW(cs_dst->orig_notDummyVar,cs_dst->orig_nvariables,boolean);
2294  memcpy(cs_dst->orig_notDummyVar,cs_src->orig_notDummyVar,cs_dst->orig_nvariables*sizeof(boolean));
2295 
2296  NEW(cs_dst->orig_varNames,cs_dst->orig_nvariables,char *);
2297  GetVariableNames(cs_dst->orig_varNames,&(cs_dst->orig_variables));
2298  }
2299  }
2300  else
2301  {
2302  cs_dst->orig2sd=NULL;
2303  cs_dst->orig2s=NULL;
2304 
2305  cs_dst->nequations=0;
2306  cs_dst->nvariables=0;
2307  cs_dst->systemVar=NULL;
2308  cs_dst->notDummyVar=NULL;
2309  cs_dst->varType=NULL;
2310  cs_dst->simp_tp=NULL;
2311 
2312  cs_dst->orig_nequations=0;
2313  cs_dst->orig_nvariables=0;
2314  cs_dst->orig_systemVar=NULL;
2315  cs_dst->orig_notDummyVar=NULL;
2316  cs_dst->orig_varNames=NULL;
2317  }
2318 }
2319 
2321 {
2322  unsigned int nvar1,nvar2;
2323  unsigned int i;
2324 
2325  InitCuikSystem(cs);
2326 
2327  MergeConstants(&(cs1->constants),&(cs2->constants),&(cs->constants));
2328 
2329  cs->empty=((cs1->empty)&&(cs2->empty));
2330  cs->simp_empty=((cs1->simp_empty)&&(cs2->simp_empty));
2331 
2332  cs->updated=FALSE;
2333  cs->consistent=TRUE; /*not used untill updated==TRUE*/
2334 
2335  cs->scalar=((ScalarEquations(&(cs1->orig_equations)))&&
2336  (ScalarEquations(&(cs->orig_equations))));
2337 
2338  cs->orig2sd=NULL;
2339 
2340  cs->nequations=0;
2341  cs->nvariables=0;
2342  cs->systemVar=NULL;
2343  cs->notDummyVar=NULL;
2344  cs->varType=NULL;
2345 
2346  cs->orig_nequations=0;
2347  cs->orig_nvariables=0;
2348  cs->orig_systemVar=NULL;
2349  cs->orig_notDummyVar=NULL;
2350 
2351  /* One of the sets of variables is supposed to be a sub-set of the
2352  other.
2353  If nvar1 is larger than the larger set was in cs1 and nothing has
2354  to be done*/
2356  nvar1=NVariables(&(cs1->orig_variables));
2357  nvar2=NVariables(&(cs2->orig_variables));
2358  for(i=nvar1;i<nvar2;i++)
2359  AddVariable2CS(GetVariable(i,&(cs2->orig_variables)),cs);
2360 
2363 
2364  cs->searchMode=cs1->searchMode;
2365  CopyEquation(&(cs->orig_eqMin),&(cs1->orig_eqMin));
2366  AccumulateEquations(&(cs2->orig_eqMin),1,&(cs->orig_eqMin));
2367 }
2368 
2369 double EvaluateEqMin(void *b,void *cs)
2370 {
2371  double v=0;
2372 
2373  if (!((TCuikSystem *)cs)->updated)
2374  Error("The CuikSystem must be updated before using EvaluateEqMin");
2375 
2376  if (!((TCuikSystem *)cs)->consistent)
2377  Error("The CuikSystem must be consistent before using EvaluateEqMin");
2378 
2379  if (((TCuikSystem *)cs)->searchMode==MINIMIZATION_SEARCH)
2380  {
2381  #if (EQ_MIN_IN_CENTER)
2382  unsigned int nv;
2383  double *p;
2384  unsigned int i;
2385 
2386  nv=NVariables(&(((TCuikSystem *)cs)->variables));
2387 
2388  NEW(p,nv,double);
2389 
2390  for(i=0;i<nv;i++)
2391  p[i]=IntervalCenter(GetBoxInterval(i,(Tbox *)b));
2392 
2393  v=EvaluateWholeEquation(p,&(((TCuikSystem *)cs)->eqMin));
2394 
2395  free(p);
2396  #else
2397  Tinterval ie;
2398  double ct;
2399 
2400  ct=GetEquationValue(&(((TCuikSystem *)cs)->eqMin));
2401 
2402  EvaluateEquationInt(GetBoxIntervals((Tbox *)b),&ie,&(((TCuikSystem *)cs)->eqMin));
2403  IntervalOffset(&ie,-ct,&ie);
2404 
2405  /*v=LowerLimit(&ie);*/
2406  /*v=UpperLimit(&ie);*/
2407  v=IntervalCenter(&ie);
2408  #endif
2409  }
2410  return(v);
2411 }
2412 
2413 boolean CmpBoxesEquation(void *b1,void *b2,void *cs)
2414 {
2415  return(EvaluateEqMin(b1,cs)<EvaluateEqMin(b2,cs));
2416 }
2417 
2418 void SetCSSearchMode(unsigned int sm,Tequation *eqMin,TCuikSystem *cs)
2419 {
2420  switch(sm)
2421  {
2422  case DEPTH_FIRST_SEARCH:
2424  break;
2425  case BREADTH_FIRST_SEARCH:
2427  break;
2428  case MINIMIZATION_SEARCH:
2430  DeleteEquation(&(cs->orig_eqMin));
2431  CopyEquation(&(cs->orig_eqMin),eqMin);
2432  break;
2433  default:
2434  Error("Unkonwn search mode in SetCSSearchMode");
2435  }
2436 }
2437 
2438 void AddTerm2SearchCriterion(double w,unsigned int v,double val,TCuikSystem *cs)
2439 {
2440  if (w!=0.0)
2441  {
2442  Tequation eq;
2443  Tmonomial m;
2444 
2445  InitEquation(&eq);
2446  SetEquationCmp(EQU,&eq);
2447  InitMonomial(&m);
2448 
2449  AddCt2Monomial(w,&m);
2450  AddVariable2Monomial(NFUN,v,2,&m);
2451  AddMonomial(&m,&eq);
2452  ResetMonomial(&m);
2453 
2454  AddCt2Monomial(-w*2*val,&m);
2455  AddVariable2Monomial(NFUN,v,1,&m);
2456  AddMonomial(&m,&eq);
2457  ResetMonomial(&m);
2458 
2459  AddCt2Monomial(w*val*val,&m);
2460  AddMonomial(&m,&eq);
2461  ResetMonomial(&m);
2462 
2464  AccumulateEquations(&eq,1,&(cs->orig_eqMin));
2465  else
2467 
2468  DeleteMonomial(&m);
2469  DeleteEquation(&eq);
2470  }
2471 }
2472 
2473 /*
2474  Adds a equation to the CuikSystem.
2475  Equations as added by the user are not dummified. This helps
2476  when simplifying the system (dummify variables can prevent some
2477  simplifications or make them more "obscure")
2478  The dummificitaion is applied after the simplification of the
2479  system (see SimplifyCuikSystem)
2480 */
2482 {
2483  if (cs->updated)
2484  UnUpdateCuikSystem(cs);
2485 
2486  if (eq!=NULL)
2487  AddEquation(eq,&(cs->orig_equations));
2488 }
2489 
2491 {
2492  if (cs->updated)
2493  UnUpdateCuikSystem(cs);
2494 
2495  if (eq!=NULL)
2496  {
2497  cs->scalar=FALSE;
2498  if (!SimplifiedMEquation(eq))
2499  Error("Adding a non-simplified matrix equation to the cuiksystem");
2500  //SimplifyMEquation(eq);
2501  AddMatrixEquation(eq,&(cs->orig_equations));
2502  }
2503 }
2504 
2505 
2506 /*
2507  Checks if the system has a variable with the given name.
2508  If so, we return the variable identifier.
2509  If not, we create a new variable
2510  */
2512 {
2513  unsigned int id;
2514 
2515  if (cs->updated)
2516  UnUpdateCuikSystem(cs);
2517 
2519 
2520  if (id==NO_UINT)
2521  id=AddVariable(v,&(cs->orig_variables));
2522 
2523  return(id);
2524 }
2525 
2526 /*
2527  Returns a copy of the variables stored in the CuikSystem
2528 */
2530 {
2531  CopyVariables(vs,&(cs->orig_variables));
2532 }
2533 
2534 void GetCSVariableNames(char **varNames,TCuikSystem *cs)
2535 {
2536  /* We do not use the cached orig_varnames since the cs
2537  can be non-updated */
2538  GetVariableNames(varNames,&(cs->orig_variables));
2539 }
2540 
2541 /*
2542  Returns the number of variables in the CuikSystem
2543 */
2545 {
2546  return(NVariables(&(cs->orig_variables)));
2547 }
2548 
2550 {
2551  return(GetNumSystemVariables(&(cs->orig_variables))+
2553 }
2554 
2555 /*
2556  Returns the number of system variables in the CuikSystem
2557 */
2559 {
2561 }
2562 
2563 /*
2564  Gets a copy of variable with ID 'n'
2565 */
2566 void GetCSVariable(unsigned int n,Tvariable *v,TCuikSystem *cs)
2567 {
2569 }
2570 
2571 /*
2572  Sets a new range for variable with ID 'n'
2573 */
2574 void SetCSVariableRange(unsigned int n,Tinterval *r,TCuikSystem *cs)
2575 {
2576  if (cs->updated)
2577  UnUpdateCuikSystem(cs); /* Variables with ct range are simplified in
2578  the UpdateCuikSystem */
2579 
2581 }
2582 
2583 /*
2584  Gets the ID of the variable with the given name
2585 */
2586 unsigned int GetCSVariableID(char *name,TCuikSystem *cs)
2587 {
2588  return(GetVariableID(name,&(cs->orig_variables)));
2589 }
2590 
2591 char *GetCSVariableName(unsigned int id,TCuikSystem *cs)
2592 {
2593  return(VariableName(id,&(cs->orig_variables)));
2594 }
2595 
2597 {
2598  unsigned int n;
2599 
2600  if (!UpdateCuikSystem(p,cs))
2601  Error("Inconsistent input cuiksystem");
2602 
2603  n=GetVariableID(v,&(cs->simp_variables));
2604  return((n!=NO_UINT)&&
2605  ((IsSystemVariable(n,&(cs->simp_variables)))||
2606  (IsSecondaryVariable(n,&(cs->simp_variables)))));
2607 }
2608 
2609 /*
2610  Returns an array of booleans with sv[i]=TRUE if variable
2611  'i' is a system variable.
2612  Also returns the number of variables in the CuikSystem
2613 */
2614 unsigned int GetCSSystemVars(boolean **sv,TCuikSystem *cs)
2615 {
2616  unsigned int i,n;
2617 
2618  n=NVariables(&(cs->orig_variables));
2619  NEW(*sv,n,boolean);
2620 
2621  for(i=0;i<n;i++)
2622  (*sv)[i]=((IsSystemVariable(i,&(cs->orig_variables)))||
2623  (IsSecondaryVariable(i,&(cs->orig_variables))));
2624  return(n);
2625 }
2626 
2627 unsigned int GetCSVarTopology(unsigned int vID,TCuikSystem *cs)
2628 {
2629  return(GetVariableTopology(GetVariable(vID,&(cs->orig_variables))));
2630 }
2631 
2632 /*
2633  Returns a copy of the equations in the cuiksystem
2634 */
2636 {
2637  CopyEquations(eqs,&(cs->orig_equations));
2638 }
2639 
2640 /*
2641  Returns a copy of equation 'n' in the CuikSystem
2642 */
2643 void GetCSEquation(unsigned int n,Tequation *eq,TCuikSystem *cs)
2644 {
2645  if (!cs->scalar)
2646  Error("GetCSEquation only operates on scalar systems");
2647 
2648  CopyEquation(eq,GetEquation(n,&(cs->orig_equations)));
2649 }
2650 
2652 {
2653  return(PolynomialEquations(&(cs->orig_equations)));
2654 }
2655 
2657 {
2658  return(cs->scalar);
2659 }
2660 
2661 /*
2662  Returns the number of equations in the CuikSytem
2663 */
2665 {
2666  return(NEquations(&(cs->orig_equations)));
2667 }
2668 
2670 {
2671  InitJacobian(&(cs->orig_variables),&(cs->orig_equations),J);
2672 }
2673 
2675  unsigned int **t,TCuikSystem *cs)
2676 {
2677  if (!UpdateCuikSystem(p,cs))
2678  Error("Inconsistent input cuiksystem");
2679 
2680  return(GetVariablesTopology(t,&(cs->simp_variables)));
2681 }
2682 
2684 {
2685  if (!cs->scalar)
2686  GetCSJacobian(J,cs);
2687  else
2688  {
2689  if (!UpdateCuikSystem(p,cs))
2690  Error("Inconsistent input cuiksystem");
2691  InitJacobian(&(cs->simp_variables),&(cs->simp_equations),J);
2692  }
2693 }
2694 
2695 void AddJacobianEquationsInt(Tparameters *p,boolean *selectedVars,TJacobian *J,
2696  TCuikSystem *cs)
2697 {
2698  Tinterval one_one,zero_one;
2699  unsigned int *nvl;
2700  Tvariable v;
2701  char vname[100];
2702  Tequation eq,eq1;
2703  unsigned int i,j,nvars,neqs;
2704 
2705  if (!cs->scalar)
2706  Error("AddJacobianEquationsInt only operates on scalar systems");
2707 
2708  GetJacobianSize(&neqs,&nvars,J);
2709 
2710  NewInterval(-1.0,1.0,&one_one);
2711  NewInterval( 0.0,1.0,&zero_one);
2712 
2713  /* From this point we start to modify cs and thus, all
2714  queries to the original cs must be done on csOrig*/
2715 
2716  /*define the lambda variables*/
2717  NEW(nvl,neqs,unsigned int);
2718  for(i=0;i<neqs;i++)
2719  {
2720  sprintf(vname,"_lambda_%u",i);
2721  NewVariable(SYSTEM_VAR,vname,&v);
2722  SetVariableInterval((i==0?&zero_one:&one_one),&v);
2723 
2724  nvl[i]=AddVariable2CS(&v,cs);
2725 
2726  DeleteVariable(&v);
2727  }
2728 
2729  /*Lambdas must be normalized*/
2730  GenerateGeneralNormEquation(neqs,nvl,1.0,&eq);
2731  AddEquation2CS(p,&eq,cs);
2732  DeleteEquation(&eq);
2733 
2734  /*define the linear combination using the lambda*/
2735  for(j=0;j<nvars;j++)
2736  {
2737  if ((selectedVars==NULL)||(selectedVars[j]))
2738  {
2739  InitEquation(&eq);
2740  SetEquationCmp(EQU,&eq);
2742 
2743  for(i=0;i<neqs;i++)
2744  {
2745  CopyEquation(&eq1,GetJacobianEquation(i,j,J));
2746  VarScaleEquation(nvl[i],&eq1);
2747  AccumulateEquations(&eq1,1.0,&eq);
2748  DeleteEquation(&eq1);
2749  }
2750 
2751  AddEquation2CS(p,&eq,cs);
2752  DeleteEquation(&eq);
2753  }
2754  }
2755 }
2756 
2757 
2758 void AddJacobianEquations(Tparameters *p,boolean *selectedVars,TCuikSystem *cs)
2759 {
2760  unsigned int neqs,nvars,i,ns;
2761  TCuikSystem csOrig;
2762  TJacobian J;
2763 
2764  if (!UpdateCuikSystem(p,cs))
2765  Error("Inconsistent input cuiksystem");
2766 
2767  CopyCuikSystem(&csOrig,cs);
2768  GetCSJacobian(&J,&csOrig);
2769  GetJacobianSize(&neqs,&nvars,&J);
2770 
2771  if (selectedVars==NULL)
2772  ns=nvars;
2773  else
2774  {
2775  ns=0;
2776  for(i=0;i<nvars;i++)
2777  if (selectedVars[i]) ns++;
2778  }
2779 
2780  if ((neqs>0)&&(ns>0))
2781  AddJacobianEquationsInt(p,selectedVars,&J,cs);
2782 
2783  DeleteJacobian(&J);
2784  UnUpdateCuikSystem(cs);
2785  DeleteCuikSystem(&csOrig);
2786 }
2787 
2788 void AddSimplifiedJacobianEquations(Tparameters *p,boolean *selectedVars,TCuikSystem *cs)
2789 {
2790  unsigned int neqs,nvars,i,ns;
2791  boolean *simpSelectedVars;
2792  unsigned int simpID,nt;
2793  TLinearConstraint lc;
2794  TCuikSystem csOrig;
2795 
2796  if (!cs->scalar)
2797  Error("AddSimplifiedJacobianEquations only operates on scalar systems");
2798 
2799  if (!UpdateCuikSystem(p,cs))
2800  Error("Inconsistent input cuiksystem");
2801 
2802  CopyCuikSystem(&csOrig,cs);
2803 
2804  neqs=cs->simp_nee;
2805  nvars=cs->simp_nvariables;
2806 
2807  NEW(simpSelectedVars,nvars,boolean);
2808  if (selectedVars==NULL)
2809  ns=nvars;
2810  else
2811  {
2812  /*Translated the selected vars from original to simplified*/
2813  for(i=0;i<nvars;i++)
2814  simpSelectedVars[i]=TRUE;
2815 
2816  ns=0;
2817  for(i=0;i<cs->orig_nvariables;i++)
2818  {
2819  if (!selectedVars[i])
2820  {
2821  GetOriginalVarRelation(i,&lc,cs->orig2s);
2823  /* if nt is zero the variable is set to constant
2824  and must be not taken into account*/
2825  if (nt>0)
2826  {
2827  ns++;
2828  /*We use the first variable in the simplified
2829  system as equivalent to the original variable*/
2830  simpID=GetLinearConstraintVariable(0,&lc);
2831  simpSelectedVars[simpID]=FALSE;
2832  #if (_DEBUG>0)
2833  printf(" Var %s -> %s\n",cs->orig_varNames[i],
2834  VariableName(simpID,&(cs->simp_variables)));
2835  #endif
2836  }
2837  #if (_DEBUG>0)
2838  else
2839  printf(" Var %s -> constant (not considered)\n",cs->orig_varNames[i]);
2840  #endif
2842  }
2843  }
2844  }
2845 
2846  if ((neqs>0)&&(ns>0))
2847  {
2848  TJacobian J;
2849  Tequation *eq;
2850  unsigned int j;
2851 
2852  /* We have to translate the simplified Jacobian matrix to the original system */
2853  InitJacobian(&(csOrig.simp_variables),&(csOrig.simp_equations),&J);
2854 
2855  for(i=0;i<neqs;i++)
2856  {
2857  for(j=0;j<nvars;j++)
2858  {
2859  eq=GetJacobianEquation(i,j,&J);
2860  RewriteEquation2Orig(csOrig.orig2s,eq,eq);
2861  }
2862  }
2863 
2864  #if (_DEBUG>1)
2865  {
2866  char **varNames;
2867 
2868  NEW(varNames,csOrig.orig_nvariables,char*);
2869  GetVariableNames(varNames,&(csOrig.orig_variables));
2870  for(j=0;j<nvars;j++)
2871  {
2872  for(i=0;i<neqs;i++)
2873  {
2874  printf("J[%u,%u] ->",i,j);
2875  PrintEquation(stdout,varNames,GetJacobianEquation(i,j,&J));
2876  }
2877  }
2878  free(varNames);
2879  }
2880  #endif
2881 
2882  AddJacobianEquationsInt(p,simpSelectedVars,&J,cs);
2883 
2884  DeleteJacobian(&J);
2885  }
2886 
2887  UnUpdateCuikSystem(cs);
2888  DeleteCuikSystem(&csOrig);
2889 }
2890 
2891 /*
2892  Reduces a box as much as possible (according to parameters in 'p').
2893  If the set of parameters is undefined, we use a default set.
2894  This procedure is used in CuikPlan.
2895 
2896  This is basically another flavor of ReduceBox where parameters
2897  are given in one structure instead of as different inputs
2898 
2899  The input box and output boxes are refered to the original
2900  cuiksystem (the non-simplified one)
2901 
2902  The output can be EMPTY_BOX
2903  REDUCED_BOX
2904  REDUCED_BOX_WITH_SOLUTION
2905  ERROR_IN_PROCESS
2906 
2907 */
2908 unsigned int MaxReduction(Tparameters *p,unsigned int varMask,double *reduction,Tbox *b,TCuikSystem *cs)
2909 {
2910  unsigned int e;
2911 
2912  if (!cs->scalar)
2913  Error("MaxReduction only operates on scalar systems");
2914 
2915  if (UpdateCuikSystem(p,cs))
2916  {
2917  Tbox bsimp;
2918  double sIn,sOut;
2919 
2920  SimpleFromOriginal(b,&bsimp,cs->orig2sd);
2921  sIn=GetBoxSize(cs->systemVar,&bsimp);
2922  e=ReduceBox(p,(varMask==0?~DUMMY_VAR:varMask),&bsimp,cs);
2923  if (e==EMPTY_BOX)
2924  *reduction=0;
2925  else
2926  {
2927  sOut=GetBoxSize(cs->systemVar,&bsimp);
2928  *reduction=sOut/sIn;
2929  /* ERROR_IN_PROCESS, REDUCED_BOX, REDUCED_BOX_WITH_SOLUTION*/
2930  UpdateOriginalFromSimple(&bsimp,b,cs->orig2sd);
2931  }
2932  DeleteBox(&bsimp);
2933  }
2934  else
2935  e=EMPTY_BOX;
2936 
2937  return(e);
2938 }
2939 
2940 
2941 boolean SampleCuikSystem(Tparameters *p,char *fname,Tlist *sb,
2942  unsigned int nsamples,unsigned int ntries,
2943  unsigned int ndof,TCuikSystem *cs)
2944 {
2945  Tbox init_box;
2946  boolean haveSample;
2947 
2948  if (!cs->scalar)
2949  Error("SampleCuikSystem only operates on scalar systems");
2950 
2951  GenerateInitialBox(&init_box,cs);
2952 
2953  haveSample=SampleCuikSystemInBox(p,fname,sb,nsamples,ntries,ndof,&init_box,cs);
2954 
2955  DeleteBox(&init_box);
2956 
2957  return(haveSample);
2958 }
2959 
2960 boolean SampleCuikSystemInBox(Tparameters *p,char *fname,Tlist *sb,
2961  unsigned int nsamples,unsigned int ntries,
2962  unsigned int ndof,
2963  Tbox *init_box,TCuikSystem *cs)
2964 {
2965  Tbox orig_ranges,b,*bAux,*box;
2966  Tlist solutions;
2967  unsigned int k,d,ns,nv,ne,i,nt;
2968  double v;
2969 
2970  Tfilename fsols_box;
2971  Tfilename fsols_point;
2972 
2973  FILE *f_out_box;
2974  FILE *f_out_point;
2975 
2976  Titerator it;
2977  boolean *systemVars;
2978  Tinterval *r;
2979  unsigned int spt;
2980  double ms;
2981  double reduction;
2982 
2983  boolean ok;
2984  boolean end;
2985 
2986  /* free var = variable not in equations -> must be fixed always */
2987  unsigned int *freeVars;
2988  unsigned int nFreeVars;
2989 
2990  unsigned int nsols;
2991 
2992  if (!cs->scalar)
2993  Error("SampleCuikSystemInBox only operates on scalar systems");
2994 
2995  /*Get a copy of the ranges in 'cs' for all the variables*/
2996  GenerateInitialBox(&orig_ranges,cs);
2997 
2998  /*Set the ranges for the 'cs' variables to those in the given box (init_box)*/
2999  if (cs->updated)
3000  UnUpdateCuikSystem(cs);
3001 
3002  nv=NVariables(&(cs->orig_variables));
3003  for(k=0;k<nv;k++)
3005 
3006  NEW(freeVars,nv,unsigned int);
3007  nFreeVars=0;
3008  for(k=0;k<nv;k++)
3009  {
3010  if ((GetVariableTypeN(k,&(cs->orig_variables))==SYSTEM_VAR)&&
3011  (!UsedVarInEquations(k,&(cs->orig_equations))))
3012  {
3013  freeVars[nFreeVars]=k;
3014  nFreeVars++;
3015  }
3016  }
3017 
3018  /*Change the SPLIT_TYPE parameter to random split*/
3019  spt=(unsigned int)GetParameter(CT_SPLIT_TYPE,p);
3020  ChangeParameter(CT_SPLIT_TYPE,2,p); /* select split dimentions at random */
3021  /*Change the SMALL_SIGMA parameter to ensure it is small enough.
3022  SMALL_SIGMA controls the precission of the isolated solutions and, thus, of the
3023  sampled points. */
3026  /*keep the original nsols parameter*/
3027  nsols=(unsigned int)GetParameter(CT_N_SOLUTIONS,p);
3028 
3029  if (fname==NULL)
3030  {
3031  f_out_box=NULL;
3032  f_out_point=NULL;
3033  }
3034  else
3035  {
3036  /* Check if we can create the output files */
3037  CreateFileName(NULL,fname,"_samples",SOL_EXT,&fsols_box);
3038  f_out_box=fopen(GetFileFullName(&fsols_box),"w");
3039  if (!f_out_box)
3040  Error("Could not open box sample output file");
3041 
3042  CreateFileName(NULL,fname,NULL,LINKS_EXT,&fsols_point);
3043  f_out_point=fopen(GetFileFullName(&fsols_point),"w");
3044  if (!f_out_point)
3045  Error("Could not open box sample output file");
3046  }
3047 
3048  if (sb!=NULL)
3049  InitListOfBoxes(sb);
3050 
3051  nv=GetCSNumVariables(cs);
3052  ne=GetCSNumEquations(cs);
3053 
3054  if (ndof==NO_UINT)
3055  {
3056  /* try to get the number of dof from the system (vars-equations) */
3057  if (ne>=nv)
3058  Error("Over/Well-constrained system. Please indicate the number of dof in the call");
3059  ndof=nv-ne;
3060  }
3061 
3062  if (MaxReduction(p,~DUMMY_VAR,&reduction,init_box,cs)==EMPTY_BOX)
3063  end=TRUE;
3064  else
3065  end=FALSE;
3066 
3067  if (nFreeVars>ndof)
3068  Error("More free vars than degrees of freedom");
3069 
3070  ns=0; /*samples already determined*/
3071  nt=0; /*number of tries so far*/
3072  while(!end)
3073  {
3074  do {
3075  ok=TRUE;
3076 
3077  /*fix 'ndof' variables*/
3078  CopyBox(&b,init_box);
3079  for(k=0;k<nv;k++)
3080  SetCSVariableRange(k,GetBoxInterval(k,&b),cs);
3081 
3082  for(k=0;k<nFreeVars;k++)
3083  {
3084  r=GetBoxInterval(freeVars[k],&b);
3085  v=randomInInterval(r);
3086 
3087  #if (_DEBUG>1)
3088  printf("[CS]->fixing %u to %f\n",freeVars[k],v);
3089  #endif
3090 
3091  NewInterval(v,v,r);
3092 
3093  SetCSVariableRange(freeVars[k],r,cs);
3094  }
3095 
3096  for(k=nFreeVars;((ok)&&(k<ndof));k++)
3097  {
3098  /*fix 'ndof' variables*/
3099 
3100  /* This updates the system, including simplification */
3101  d=ComputeSplitDim(p,&b,cs);
3102 
3103  if (d==NO_UINT)
3104  ok=FALSE;
3105  else
3106  {
3107  r=GetBoxInterval(d,&b);
3108  v=randomInInterval(r);
3109 
3110  #if (_DEBUG>1)
3111  printf("[CS]->fixing %u to %f\n",d,v);
3112  #endif
3113 
3114  NewInterval(v,v,r);
3115 
3116  SetCSVariableRange(d,r,cs);
3117  }
3118  }
3119 
3120  #if (_DEBUG>1)
3121  if (ok)
3122  {printf("[CS]");PrintBox(stdout,&b);}
3123  else
3124  printf("Inconsistent system while sampling\n");
3125  #endif
3126 
3127  DeleteBox(&b);
3128 
3129  nt++;
3130  if ((ntries!=NO_UINT)&&(nt==ntries))
3131  end=TRUE;
3132  } while((!ok)&&(!end));
3133 
3134  if (!end)
3135  {
3136 
3137  InitListOfBoxes(&solutions);
3138 
3139  /*See if there are solutions for the fixed ranges*/
3140  ChangeParameter(CT_N_SOLUTIONS,nsamples-ns,p);
3141  SolveCuikSystem(p,FALSE,NULL,NULL,f_out_box,&solutions,cs);
3142 
3143  /*We have to print the samples*/
3144  GetCSSystemVars(&systemVars,cs);
3145  InitIterator(&it,&solutions);
3146  First(&it);
3147  while((!EndOfList(&it))&&(!end))
3148  {
3149  bAux=(Tbox *)GetCurrent(&it);
3150 
3151  if (f_out_point!=NULL)
3152  {
3153  /* Print the center of the box into the point file */
3154  /* We only print values for the system vars */
3155  for(i=0;i<nv;i++)
3156  {
3157  if (systemVars[i])
3158  fprintf(f_out_point," %.12g",IntervalCenter(GetBoxInterval(i,bAux)));
3159  }
3160  fprintf(f_out_point,"\n");
3161  }
3162  if (sb!=NULL)
3163  {
3164  NEW(box,1,Tbox);
3165 
3166  CopyBox(box,bAux);
3167 
3168  AddLastElement((void *)box,sb);
3169  }
3170 
3171  ns++;
3172  if (ns<nsamples)
3173  end=TRUE;
3174  Advance(&it);
3175  }
3176  free(systemVars);
3177 
3178  fflush(f_out_point);
3179  fflush(f_out_box);
3180 
3181  /*Remove all boxes from the list and the current box*/
3182  DeleteListOfBoxes(&solutions);
3183  }
3184  }
3185 
3186  if (f_out_box!=NULL)
3187  {
3188  DeleteFileName(&fsols_box);
3189  fclose(f_out_box);
3190  }
3191  if (f_out_point!=NULL)
3192  {
3193  fclose(f_out_point);
3194  DeleteFileName(&fsols_point);
3195  }
3196 
3197  /*Restore the user given split_type, min_sigma, and n_solutions*/
3201 
3202  /*Return the original ranges to the 'cs' variables*/
3203  if (cs->updated)
3204  UnUpdateCuikSystem(cs);
3205 
3206  nv=NVariables(&(cs->orig_variables));
3207  for(k=0;k<nv;k++)
3208  SetVariableInterval(GetBoxInterval(k,&orig_ranges),GetVariable(k,&(cs->orig_variables)));
3209 
3210  free(freeVars);
3211 
3212  DeleteBox(&orig_ranges);
3213 
3214  return(!end);
3215 }
3216 
3217 
3218 boolean IncrementalSampleCuikSystem(Tparameters *p,char *fname,Tlist *sb,
3219  boolean *fixVars,
3220  unsigned int nsamples,unsigned int ntries,
3221  unsigned int ndof,TCuikSystem *cs)
3222 {
3223  Tbox init_box;
3224  boolean haveSample;
3225 
3226  if (!cs->scalar)
3227  Error("IncrementalSampleCuikSystem only operates on scalar systems");
3228 
3229  GenerateInitialBox(&init_box,cs);
3230 
3231  haveSample=IncrementalSampleCuikSystemInBox(p,fname,sb,fixVars,nsamples,ntries,ndof,&init_box,cs);
3232 
3233  DeleteBox(&init_box);
3234 
3235  return(haveSample);
3236 }
3237 
3239  boolean *fixVars,
3240  unsigned int nsamples,unsigned int ntries,
3241  unsigned int ndof,
3242  Tbox *init_box,TCuikSystem *cs)
3243 {
3244  Tbox b,*bAux,*box;
3245  Tlist solutions;
3246  unsigned int k,d,ns,nv,ne,i,nt;
3247  double v;
3248 
3249  double smallSigma;
3250 
3251  Tfilename fsols_box;
3252  Tfilename fsols_point;
3253 
3254  FILE *f_out_box;
3255  FILE *f_out_point;
3256 
3257  Titerator it;
3258  boolean *systemVars;
3259  Tinterval *r,newr;
3260  Tbox *bl;
3261  double ms;
3262  double reduction;
3263 
3264  boolean ok;
3265  boolean end;
3266 
3267  boolean *fv;
3268  unsigned int *ind2fv;
3269  unsigned int nfv;
3270  unsigned int newtonIterations;
3271  double *sol;
3272  Tbox b_sol;
3273  Tbox orig_ranges;
3274  unsigned int level,*attempts,maxAttemptsPerLevel;
3275  unsigned int nsols;
3276 
3277  if (!cs->scalar)
3278  Error("IncrementalSampleCuikSystemInBox only operates on scalar systems");
3279 
3280  /*Set the ranges for the 'cs' variables to those in the given box (init_box)*/
3281 
3282  if (!UpdateCuikSystem(p,cs))
3283  Error("Inconsistent cuiksystem in IncrementalSampleCuikSystemInBox");
3284 
3285  BoxFromVariables(&orig_ranges,&(cs->orig_variables));
3286 
3287  nv=NVariables(&(cs->orig_variables));
3288 
3289  NEW(sol,nv,double);
3290 
3292  smallSigma=100*GetParameter(CT_EPSILON,p);
3293  ChangeParameter(CT_SMALL_SIGMA,smallSigma,p);
3294 
3295  nsols=(unsigned int)GetParameter(CT_N_SOLUTIONS,p);
3296 
3297  newtonIterations=(unsigned int)GetParameter(CT_MAX_NEWTON_ITERATIONS,p);
3298 
3299  /* Array of variables that can be fixed */
3300  NEW(fv,nv,boolean);
3301  NEW(ind2fv,nv,unsigned int);
3302  nfv=0;
3303  for(k=0;k<nv;k++)
3304  {
3305  if ((IntervalSize(GetBoxInterval(k,init_box))>smallSigma)&&
3306  (IsInSimple(k,cs->orig2sd)))
3307  {
3308  if (((fixVars!=NULL)&&(fixVars[k]))||
3309  ((fixVars==NULL)&&(IsSystemVariable(k,&(cs->orig_variables)))))
3310  {
3311  fv[k]=TRUE;
3312  ind2fv[nfv]=k;
3313  nfv++;
3314  }
3315  else
3316  fv[k]=FALSE;
3317  }
3318  else
3319  fv[k]=FALSE;
3320  }
3321 
3322  if (fname==NULL)
3323  {
3324  f_out_box=NULL;
3325  f_out_point=NULL;
3326  }
3327  else
3328  {
3329  /* Check if we can create the output files */
3330  CreateFileName(NULL,fname,"_samples",SOL_EXT,&fsols_box);
3331  f_out_box=fopen(GetFileFullName(&fsols_box),"w");
3332  if (!f_out_box)
3333  Error("Could not open box sample output file");
3334 
3335  CreateFileName(NULL,fname,NULL,LINKS_EXT,&fsols_point);
3336  f_out_point=fopen(GetFileFullName(&fsols_point),"w");
3337  if (!f_out_point)
3338  Error("Could not open box sample output file");
3339  }
3340 
3341  if (sb!=NULL)
3342  InitListOfBoxes(sb);
3343 
3344  nv=GetCSNumVariables(cs);
3345  ne=GetCSNumEquations(cs);
3346 
3347  if (ndof==NO_UINT)
3348  {
3349  /* try to get the number of dof from the system (vars-equations) */
3350  if (ne>=nv)
3351  Error("Over/Well-constrained system. Please indicate the number of dof in the call");
3352  ndof=nv-ne;
3353  }
3354  #if (_DEBUG>1)
3355  printf("Initial reduction (s:%g)\n",GetBoxSumSide(NULL,init_box));
3356  #endif
3357  if (MaxReduction(p,~DUMMY_VAR,&reduction,init_box,cs)==EMPTY_BOX)
3358  end=TRUE;
3359  else
3360  end=FALSE;
3361 
3362  #if (_DEBUG>1)
3363  if (!end)
3364  {
3365  printf("[CS] (s:%g)",GetBoxSumSide(NULL,init_box));
3366  PrintBox(stdout,init_box);
3367  }
3368  else
3369  printf("Inconsistent system while sampling\n");
3370  fflush(stdout);
3371  #endif
3372 
3373  NEW(attempts,ndof,unsigned int);
3374  NEW(bl,ndof,Tbox);
3375  maxAttemptsPerLevel=10;
3376 
3377  ns=0; /*samples already determined*/
3378  nt=0; /*time we tried to sample so far*/
3379  while((!end)&&(ns<nsamples))
3380  {
3381  ok=TRUE;
3382 
3383  CopyBox(&b,init_box);
3384 
3385  /*fix ndof variables*/
3386  level=0;
3387  for(k=0;k<ndof;k++)
3388  attempts[k]=0;
3389 
3390  while((level<ndof)&&(!end))
3391  {
3392  /*Keep a copy of the box BEFORE fixing the new variable at the current
3393  'level'. This is used when backtracking to recover the state before this
3394  assigment*/
3395  CopyBox(&(bl[level]),&b);
3396 
3397  do {
3398  /*select one of the still not fixed variables*/
3399  d=ind2fv[randomMax(nfv-1)];
3400  r=GetBoxInterval(d,&b);
3401  } while (IntervalSize(r)<=smallSigma);
3402 
3403  #if (_DEBUG>1)
3404  printf("[CS](level %u/%u bs:%f)",level,ndof,GetBoxSumSide(NULL,&b));
3405  #endif
3406 
3407  {
3408  Tinterval raux;
3409  double c,l;
3410 
3411  c=IntervalCenter(r);
3412  l=IntervalSize(r)/2;
3413 
3414  NewInterval(c-l,c+l,&raux);
3415 
3416  /*Now change the current box ,b*/
3417  v=randomInInterval(&raux);
3418  }
3419  //v=randomInInterval(r);
3420 
3421  NewInterval(v,v,&newr);
3422 
3423  SetBoxInterval(d,&newr,&b);
3424 
3425  /* Change the variable ranges to those in b */
3426  UnUpdateCuikSystem(cs);
3427  VariablesFromBox(&b,&(cs->orig_variables));
3428 
3429  attempts[level]++;
3430 
3431  #if (_DEBUG>1)
3432  printf("(attempts %u/%u)->fixing %u to %f \n",
3433  attempts[level],maxAttemptsPerLevel,d,v);
3434  #endif
3435 
3436  /* And reduce the systems as much as possible */
3437  ok=((UpdateCuikSystem(p,cs))&&
3438  (MaxReduction(p,~DUMMY_VAR,&reduction,&b,cs)!=EMPTY_BOX));
3439 
3440  if (ok)
3441  {
3442  if (newtonIterations>0)
3443  {
3444  if ((CuikNewtonInBox(p,&b,sol,&b_sol,cs)&(CONVERGED_IN_GLOBAL|CONVERGED_IN_BOX))>0)
3445  {
3446  fprintf(f_out_box,"Newton (%g)",ErrorInSolution(&b_sol,cs));
3447  PrintBoxSubset(f_out_box,cs->orig_notDummyVar,cs->orig_varNames,&b_sol);
3448  fflush(f_out_box);
3449  }
3450  DeleteBox(&b_sol);
3451  }
3452  level++;
3453 
3454  #if (_DEBUG>1)
3455  printf("[CS] (l:%u s:%g)",level,GetBoxSumSide(NULL,&b));
3456  PrintBox(stdout,&b);
3457  #endif
3458  }
3459  else
3460  {
3461  #if (_DEBUG>1)
3462  printf("Inconsistent system while sampling. Backtracking\n");
3463  #endif
3464 
3465  /*Undo the last assigment*/
3466  DeleteBox(&b);
3467  CopyBox(&b,&(bl[level]));
3468  DeleteBox(&(bl[level]));
3469 
3470  if (attempts[level]>maxAttemptsPerLevel)
3471  {
3472  attempts[level]=0;
3473  if (level>0)
3474  {
3475  level--;
3476  /*move to the state just before the assigment at leve-1*/
3477  DeleteBox(&b);
3478  CopyBox(&b,&(bl[level]));
3479  DeleteBox(&(bl[level]));
3480  }
3481  }
3482  }
3483  }
3484 
3485  if (ok)
3486  {
3487  /* Last degree of freedom is not fixed but we solve the cuiksystem in the box */
3488  InitListOfBoxes(&solutions);
3489 
3490  /*See if there are solutions for the fixed ranges*/
3491  ChangeParameter(CT_N_SOLUTIONS,nsamples-ns,p);
3492  SolveCuikSystem(p,FALSE,NULL,&b,f_out_box,&solutions,cs);
3493 
3494  /*We have to print the samples*/
3495  GetCSSystemVars(&systemVars,cs);
3496  InitIterator(&it,&solutions);
3497  First(&it);
3498  while((!EndOfList(&it))&&(ns<nsamples))
3499  {
3500  bAux=(Tbox *)GetCurrent(&it);
3501 
3502  if (f_out_point!=NULL)
3503  {
3504  /* Print the center of the box into the point file */
3505  /* We only print values for the system vars */
3506  for(i=0;i<nv;i++)
3507  {
3508  if (systemVars[i])
3509  fprintf(f_out_point," %.12g",IntervalCenter(GetBoxInterval(i,bAux)));
3510  }
3511  fprintf(f_out_point,"\n");
3512  }
3513  if (sb!=NULL)
3514  {
3515  NEW(box,1,Tbox);
3516 
3517  CopyBox(box,bAux);
3518 
3519  AddLastElement((void *)box,sb);
3520  }
3521 
3522  ns++;
3523  Advance(&it);
3524  }
3525  free(systemVars);
3526 
3527  fflush(f_out_point);
3528  fflush(f_out_box);
3529 
3530  /*Remove all boxes from the list and the current box*/
3531  DeleteListOfBoxes(&solutions);
3532  }
3533 
3534  DeleteBox(&b);
3535 
3536  nt++;
3537  if ((ntries!=NO_UINT)&&(nt==ntries))
3538  end=TRUE;
3539 
3540  }
3541 
3542  if (f_out_box!=NULL)
3543  {
3544  DeleteFileName(&fsols_box);
3545  fclose(f_out_box);
3546  }
3547  if (f_out_point!=NULL)
3548  {
3549  fclose(f_out_point);
3550  DeleteFileName(&fsols_point);
3551  }
3552 
3553  /*Restore the user given small_sigma and number of solutions*/
3556 
3557  free(attempts);
3558  free(bl);
3559 
3560  free(fv);
3561  free(ind2fv);
3562  free(sol);
3563 
3564  UnUpdateCuikSystem(cs);
3565  VariablesFromBox(&orig_ranges,&(cs->orig_variables));
3566  DeleteBox(&orig_ranges);
3567 
3568  return(!end);
3569 }
3570 
3571 unsigned int CuikNewtonSimp(Tparameters *p,double *x,TCuikSystem *cs)
3572 {
3573  boolean converged=FALSE;
3574  unsigned int out=DIVERGED;
3575 
3576  if (UpdateCuikSystem(p,cs))
3577  {
3578  if (cs->simp_nequations==0)
3579  {
3580  out=CONVERGED_IN_BOX;
3581  }
3582  else
3583  {
3584  TNewton newton;
3585  double *Jx_d;
3586  double *b_d;
3587 
3588  double epsilon,nullSingularValue;
3589  double errorVal;
3590  int err;
3591  unsigned int nStep,nStepMax;
3592 
3593  #if (_DEBUG>2)
3594  {
3595  unsigned int i;
3596 
3597  printf("Starting Newton Iteration form: [");
3598  for(i=0;i<cs->simp_nvariables;i++)
3599  printf("%f ",x[i]);
3600  printf("]\n");
3601  }
3602  #endif
3603  InitNewton(cs->simp_nee,cs->simp_nvariables,&newton);
3604  Jx_d=GetNewtonMatrixBuffer(&newton);
3605  b_d=GetNewtonRHBuffer(&newton);
3606 
3607  epsilon=GetParameter(CT_EPSILON,p);
3608  nullSingularValue=epsilon/100;
3609 
3610  nStepMax=(unsigned int)GetParameter(CT_MAX_NEWTON_ITERATIONS,p);
3611  if (nStepMax==0)
3612  Error("Parameter MAX_NEWTON_ITERATIONS must be larger than 0 to use CuikNewton");
3613 
3614  nStep=0;
3615 
3616  while((!converged)&&(nStep<nStepMax))
3617  {
3618  EvaluateJacobianInVector(x,cs->simp_nee,cs->simp_nvariables,Jx_d,&(cs->J));
3620 
3621  err=NewtonStep(nullSingularValue,x,&errorVal,&newton);
3622 
3623  if (err)
3624  nStep=nStepMax;
3625  else
3626  {
3627  #if (_DEBUG>2)
3628  {
3629  unsigned int i;
3630 
3631  printf(" Iteration: %u error: %g",nStep,errorVal);
3632  #if (_DEBUG>3)
3633  printf(" x: [\n");
3634  for(i=0;i<cs->simp_nvariables;i++)
3635  printf("%f ",x[i]);
3636  printf("]");
3637  #endif
3638  printf("\n");
3639  }
3640  #endif
3641 
3642  /* stopCriterion test */
3643  if(errorVal<epsilon)
3644  converged=TRUE;
3645 
3646  nStep++;
3647  }
3648  }
3649 
3650  if (cs->simp_tp!=NULL)
3651  ArrayPi2Pi(cs->simp_nvariables,cs->simp_tp,x);
3652 
3653  DeleteNewton(&newton);
3654 
3655  if ((converged)&&(ErrorInSimpCSEquations(p,x,cs)<epsilon))
3656  out=CONVERGED_IN_GLOBAL;
3657  else
3658  out=DIVERGED;
3659  }
3660  }
3661  else
3662  Error("Inconsistent input cuiksystem");
3663 
3664  return(out);
3665 }
3666 
3667 unsigned int CuikNewtonInBox(Tparameters *p,Tbox *bIn,double *sol,Tbox *b_sol,TCuikSystem *cs)
3668 {
3669  boolean converged=FALSE;
3670  unsigned int out=DIVERGED;
3671 
3672  if (UpdateCuikSystem(p,cs))
3673  {
3674  if (cs->simp_nequations==0)
3675  {
3676  unsigned int i;
3677  Tbox ib;
3678 
3679  /* The system of equations is empty -> just pic a random value for the variables */
3680  BoxFromVariables(&ib,&(cs->orig_variables));
3681  for(i=0;i<cs->orig_nvariables;i++)
3682  sol[i]=randomInInterval(GetBoxInterval(i,&ib));
3683  DeleteBox(&ib);
3684  InitBoxFromPoint(cs->orig_nvariables,sol,b_sol);
3685  out=CONVERGED_IN_BOX;
3686  }
3687  else
3688  {
3689  Tbox bsimp;
3690  Tinterval c,*r;
3691 
3692  unsigned int nv,ne,bs;
3693  unsigned int i,j,k,l;
3694 
3695  double *m,*h;
3696  double *x;
3697  double *Jx_d;
3698  double *b_d;
3699 
3700  /* Newton required variables */
3701  TNewton newton;
3702 
3703  double epsilon;
3704  double nullSingularValue;
3705  double errorVal;
3706  unsigned int nStep,nStepMax;
3707  int err;
3708  Tbox box;
3709 
3710  /*randomSet(2341);*/
3711 
3712  bs=GetBoxNIntervals(bIn);
3713  if (cs->orig_nvariables<bs)
3714  Error("Box size missmatch in CuikNewtonInBox");
3715  else
3716  {
3717  if (cs->orig_nvariables>bs)
3718  {
3719  if (cs->orig_nvariables==(bs+GetNumDummyVariables(&(cs->orig_variables))))
3720  {
3721  /* We are given a solution box without the dummy vars */
3722  /* We need to compute the values for the dummy vars from the system ones*/
3723 
3724  GenerateInitialBox(&box,cs); /*all variables to initial ranges*/
3725 
3726  SetBoxSubset(cs->orig_systemVar,bIn,&box); /* only dummies to initial ranges */
3727  RegenerateSolution(p,&box,cs); /* define dummies from system variables */
3728  }
3729  else
3730  Error("Box size missmatch in CuikNewtonInBox");
3731  }
3732  else
3733  CopyBox(&box,bIn);
3734  }
3735 
3736  nStepMax=(unsigned int)GetParameter(CT_MAX_NEWTON_ITERATIONS,p);
3737  if (nStepMax==0)
3738  Error("Parameter MAX_NEWTON_ITERATIONS must be larger than 0 to use CuikNewton");
3739 
3740  SimpleFromOriginal(&box,&bsimp,cs->orig2s);
3741 
3742  /* We will add some variables and equations to deal with ranges */
3743  #if (NEWTON_WITHIN_RANGES)
3744  nv=cs->simp_nvariables+cs->simp_nvariables;
3745  ne=cs->simp_nee+cs->simp_nvariables;
3746  #else
3747  nv=cs->simp_nvariables;
3748  ne=cs->simp_nee;
3749  #endif
3750 
3751  /* Space for the variable ranges */
3752  NEW(m,cs->simp_nvariables,double);
3753  NEW(h,cs->simp_nvariables,double);
3754  for(i=0;i<cs->simp_nvariables;i++)
3755  {
3756  r=GetBoxInterval(i,&bsimp);
3757  m[i]=IntervalCenter(r);
3758  h[i]=IntervalSize(r)/2.0;
3759  }
3760 
3761  /* Matrices and arrays used in the Newton iteration */
3762  NEW(x,nv,double);
3763  for(i=0;i<cs->simp_nvariables;i++)
3764  x[i]=randomInInterval(GetBoxInterval(i,&bsimp));
3765 
3766  for(i=cs->simp_nvariables,j=0;i<nv;i++,j++)
3767  x[i]=asin((x[j]-m[j])/h[j]);
3768 
3769  #if (_DEBUG>2)
3770  {
3771  printf("Starting Newton Iteration form: x=[");
3772  for(i=0;i<nv;i++)
3773  printf("%f ",x[i]);
3774  printf("]\n");
3775  }
3776  #endif
3777 
3778  InitNewton(ne,nv,&newton);
3779  Jx_d=GetNewtonMatrixBuffer(&newton);
3780 
3781  k=ne*nv;
3782  if (k!=(cs->simp_nee*cs->simp_nvariables))
3783  {
3784  /* If we have limit constraints just init the full
3785  buffer to zero. */
3786  for(i=0;i<k;i++)
3787  Jx_d[i]=0;
3788  }
3789 
3790  b_d=GetNewtonRHBuffer(&newton);
3791 
3792  epsilon=GetParameter(CT_EPSILON,p);
3793  nullSingularValue=epsilon/100;
3794 
3795  nStep=0;
3796 
3797  while((!converged)&&(nStep<nStepMax))
3798  {
3799  EvaluateJacobianInVector(x,ne,nv,Jx_d,&(cs->J));
3801 
3802  for(i=cs->simp_nvariables,j=0,l=cs->simp_nee;i<nv;i++,j++,l++)
3803  {
3804  /* 'j' is the index for the original variable 'v' */
3805  /* 'i' is the index for the new variable 'q' */
3806  NewtonSetMatrix(l,j,1,&newton);
3807  NewtonSetMatrix(l,i,-h[j]*cos(x[i]),&newton);
3808 
3809  NewtonSetRH(l,(x[j]-m[j]-h[j]*sin(x[i])),&newton);
3810  }
3811 
3812  errorVal=Norm(ne,b_d);
3813 
3814  //fprintf(stderr,"%u Newton Error: %g\n",nStep,errorVal);
3815 
3816  converged=(errorVal<epsilon);
3817 
3818  if (!converged)
3819  {
3820  err=NewtonStep(nullSingularValue,x,&errorVal,&newton);
3821 
3822  if (err)
3823  nStep=nStepMax;
3824  else
3825  {
3826  #if (_DEBUG>2)
3827  printf(" Iteration: %u error: %g",nStep,errorVal);
3828  #if (_DEBUG>3)
3829  printf(" x: [\n");
3830  for(i=0;i<cs->simp_nvariables;i++)
3831  printf("%f ",x[i]);
3832  printf("]");
3833  #endif
3834  printf("\n");
3835  #endif
3836 
3837  nStep++;
3838  }
3839  }
3840  }
3841  DeleteNewton(&newton);
3842 
3843  if (cs->simp_tp!=NULL)
3844  ArrayPi2Pi(cs->simp_nvariables,cs->simp_tp,x);
3845 
3846  for(i=0;i<cs->simp_nvariables;i++)
3847  {
3848  NewInterval(x[i],x[i],&c);
3849  SetBoxInterval(i,&c,&bsimp);
3850  }
3851 
3852  CopyBox(b_sol,&box);
3853  UpdateOriginalFromSimple(&bsimp,b_sol,cs->orig2s);
3854 
3855  /* Some dummy variables in original system might not be present in
3856  the simplified system but their value can be trivially deduced
3857  from the system variables*/
3858  RegenerateSolution(p,b_sol,cs);
3859 
3860  /*At this point b_sol must be a punctual box (fully defined) */
3861  for(i=0;i<cs->orig_nvariables;i++)
3862  sol[i]=IntervalCenter(GetBoxInterval(i,b_sol));
3863 
3864  if (!converged) //((!converged)||(ErrorInSolution(b_sol,cs)>epsilon))
3865  out=DIVERGED;
3866  else
3867  {
3868  Tbox init_box;
3869 
3870  GenerateInitialBox(&init_box,cs);
3871 
3872  if ((ErrorInInequalities(b_sol,cs)>epsilon)||
3873  (!PointInBox(cs->orig_notDummyVar,cs->orig_nvariables,sol,epsilon,&init_box)))
3875  else
3876  {
3877  if (PointInBox(cs->orig_systemVar,cs->orig_nvariables,sol,epsilon,&box))
3878  out=CONVERGED_IN_BOX;
3879  else
3880  out=CONVERGED_IN_GLOBAL;
3881  }
3882 
3883  DeleteBox(&init_box);
3884  }
3885 
3886  DeleteBox(&bsimp);
3887  free(x);
3888  free(m);
3889  free(h);
3890  DeleteBox(&box);
3891  }
3892  }
3893  else
3894  Error("Inconsistent input cuiksystem");
3895 
3896  return(out);
3897 }
3898 
3899 /*
3900  Find a solution of a cuiksystem set of equations
3901  using the Newton-Rhapson method. If the system is undetermined,
3902  the method is based on the Moor-Penrose generalised inverse.
3903  The method can diverge if the initial point is far from the solutions of the system.
3904  */
3905 boolean CuikNewton(Tparameters *p,double *sol,Tbox *b_sol,TCuikSystem *cs)
3906 {
3907  Tbox init_box;
3908  boolean converged;
3909 
3910  GenerateInitialBox(&init_box,cs);
3911 
3912  converged=((CuikNewtonInBox(p,&init_box,sol,b_sol,cs)&(CONVERGED_IN_BOX|CONVERGED_IN_GLOBAL))>0);
3913 
3914  DeleteBox(&init_box);
3915 
3916  return(converged);
3917 }
3918 
3919 /*
3920  Takes a cuik systema and returns all solutions.
3921  Solutions are stored in file 'f_out' (if defined) and in list 'sol'
3922  (if defined), or both.
3923 
3924  This procedure is used for single-cpu versions of Cuik.
3925  */
3927  boolean restart,char *fstate,Tbox *searchSpace,
3928  FILE *f_out,Tlist *sol,TCuikSystem *cs)
3929 {
3930  Theap boxes;
3931  #if (_DEBUG>1)
3932  unsigned int nbox=0;
3933  #endif
3934 
3935  Tbox b;
3936  unsigned int c; /* output of ReduceBox */
3937  unsigned int current_level; /* level of the current box in the search three */
3938  Tbox init_box;
3939  boolean done;
3940  unsigned int nb; /*boxes for recovery mode*/
3941  unsigned int statePeriod;
3942  unsigned int nsols;
3943 
3944  if (!cs->scalar)
3945  Error("SolveCuikSystem only operates on scalar systems");
3946 
3947  /************************************************************************************/
3948  /************************************************************************************/
3949  /************************************************************************************/
3950 
3951  if (UpdateCuikSystem(p,cs))
3952  {
3953  if (cs->simp_empty)
3954  {
3955  Tbox b;
3956 
3957  /* The system of equations is empty -> just pic a random value for the variables */
3958  InitBox(cs->simp_nvariables,NULL,&b);
3959  PostProcessBox(p,REDUCED_BOX_WITH_SOLUTION,f_out,sol,&boxes,&b,cs);
3960  DeleteBox(&b);
3961  }
3962  else
3963  {
3964  switch(cs->searchMode)
3965  {
3966  case DEPTH_FIRST_SEARCH:
3967  InitHeapOfBoxes(CmpBoxDepthFirst,NULL,&boxes);
3968  break;
3969  case BREADTH_FIRST_SEARCH:
3970  InitHeapOfBoxes(CmpBoxBreadthFirst,NULL,&boxes);
3971  break;
3972  case MINIMIZATION_SEARCH:
3973  InitHeapOfBoxes(CmpBoxesEquation,(void *)cs,&boxes);
3974  break;
3975  }
3976 
3977  if ((restart)&&(fstate!=NULL))
3978  {
3979  Tlist pboxes;
3980 
3981  LoadCSState(fstate,&pboxes,cs);
3982  AddList2Heap(&pboxes,&boxes);
3983  DeleteListOfBoxes(&pboxes);
3984  }
3985  else
3986  {
3987  if (searchSpace==NULL)
3988  {
3989  BoxFromVariables(&init_box,&(cs->variables));
3990  AddBox2HeapOfBoxes(&init_box,&boxes);
3991  DeleteBox(&init_box);
3992  }
3993  else
3994  {
3995  SimpleFromOriginal(searchSpace,&init_box,cs->orig2sd);
3996  AddBox2HeapOfBoxes(&init_box,&boxes);
3997  }
3998 
3999  InitStatistics(1,HeapOfBoxesVolume(cs->systemVar,&boxes),&(cs->st));
4000  }
4001 
4002  current_level=0;
4003  nb=0;
4004  statePeriod=(unsigned int)GetParameter(CT_STATE_PERIOD,p);
4005  nsols=(unsigned int)GetParameter(CT_N_SOLUTIONS,p);
4006 
4007  done=FALSE;
4008 
4009  /*the cuik solve process itself*/
4010  do {
4011  if(!HeapEmpty(&boxes))
4012  {
4013  /*Get a new box*/
4014 
4015  if ((statePeriod>0)&&(fstate!=NULL))
4016  {
4017  nb++;
4018  if (nb==statePeriod)
4019  {
4020  Tlist pboxes;
4021 
4022  Heap2List(&pboxes,&boxes);
4023  SaveCSState(fstate,&pboxes,cs);
4024  DeleteListOfBoxes(&pboxes);
4025  nb=0;
4026  }
4027  }
4028 
4029  NewBoxProcessed(&(cs->st));
4030 
4031  /*The box is extracted to avoid another process
4032  to get the same box*/
4033  ExtractMinElement(&b,&boxes);
4034 
4035  current_level=GetBoxLevel(&b);
4036 
4037  NewMaxLevel(current_level,&(cs->st));
4038 
4039  #if (_DEBUG>0)
4040  fprintf(stderr,"<%g %g>[%u]",
4041  GetBoxVolume(cs->systemVar,&b),
4042  GetBoxSize(cs->systemVar,&b),
4043  current_level);
4044  #endif
4045 
4046  #if (_DEBUG>1)
4047  printf("\n\n%u: Analizing box: ",nbox);
4048  nbox++;
4049  PrintBox(stdout,&b);
4050  printf(" Box level : %u\n",GetBoxLevel(&b));
4051  printf(" Box volume : %g\n",GetBoxVolume(cs->systemVar,&b));
4052  printf(" Box size : %g\n",GetBoxSize(cs->systemVar,&b));
4054  printf(" Box penalty : %g\n",EvaluateEqMin((void *)&b,(void *)cs));
4055  #endif
4056 
4057  /* Reduce the box by all but dummy vars */;
4058  c=ReduceBox(p,~DUMMY_VAR,&b,cs);
4059 
4060  #if (_DEBUG>0)
4061  fprintf(stderr," -> ");
4062  #endif
4063 
4064  PostProcessBox(p,c,f_out,sol,&boxes,&b,cs);
4065 
4066  if (nsols>0)
4067  done=(GetNSolutionBoxes(&(cs->st))>=nsols);
4068 
4069  DeleteBox(&b);
4070  }
4071  } while((!done)&&(!HeapEmpty(&boxes))); /*Until everything is completed*/
4072 
4073  DeleteHeap(&boxes); /*for sanity we delete the empty list*/
4074 
4075  PrintStatistics((f_out==NULL?stdout:f_out),&(cs->st));
4076  }
4077  }
4078 }
4079 
4080 
4081 #if (_USE_MPI)
4082 /*
4083  In multi-cpu versions we use a MPI version of CuikSolve
4084  One of the CPUs takes the role of scheduler in charge of managing the
4085  lists of boxes pending to be processed and classifing boxes out of
4086  reduction process as solutions of boxes to be split.
4087  The rest of CPUs take the charge of reducing boxes (i.e., of executing
4088  ReduceBox to the boxes received from the scheduler) and of returning
4089  the resulting boxes (with the associated result code) to the
4090  scheduler.
4091 
4092  The scheduler has to take into account that some child nodes can
4093  'die' while reducing boxes. So it has to keep track of when a box was
4094  sent to a child and prediodically check whether a timeout was been reached.
4095  If so, the corresponding processor is considered as 'dead' the
4096  box sent is recovered and added to the list of boxes to be processed
4097  (i.e., send to other child nodes still alive).
4098 
4099  The scheduler should considere the extreme case where all child process die
4100  and should take care of "awakening" dead processors.
4101  This is not implemented.... till now.
4102 */
4104  boolean restart,char *fstate,Tbox *searchSpace,
4105  FILE *f_out,TCuikSystem *cs)
4106 {
4107  Theap boxes; /*queued boxes*/
4108  Tbox b;
4109  unsigned int c;
4110  unsigned int nr;
4111  unsigned int current_level;
4112 
4113  /*Variables to set up the parallelism*/
4114  unsigned int in_process; /* number of boxes send to child processors and still
4115  in process */
4116  unsigned int available; /* number of child processors available */
4117  signed int np; /* number of processors used (np-1) is the number of childs (process
4118  0 is the scheduler) */
4119  MPI_Status status; /* Output status of a MPI command */
4120  MPI_Request *request; /* Communication request. When the scheduller sends out a box
4121  is starts a communication request */
4122  unsigned int buffer_size; /* Size of a buffer of doubles where we store the
4123  information of a box */
4124  double **buffers_out; /* boxes send to the child processors */
4125  double **buffers_in; /* boxes returned by the child processors */
4126  signed int i,completed; /* MPI_Testany returns completed=1 if any of the pending
4127  communications requests has been replied (then, 'i'
4128  is the number of the compleated request)*/
4129  boolean *lost_packet; /* true if a box was never returned */
4130  time_t *time_stamp; /* time when we send out a box */
4131  time_t deadline; /* deadline (in seconds) used to compute the maximum expected
4132  return time for a box*/
4133  Tbox init_box; /* initial search space */
4134  unsigned int nb;
4135  unsigned int statePeriod;
4136 
4137  unsigned int nsols;
4138  boolean done;
4139 
4140  if (!cs->scalar)
4141  Error("MPI_SolveCuikSystem only operates on scalar systems");
4142 
4143  /*lowest possible priority but not as low as that of the children processors*/
4144  /*setpriority(PRIO_PROCESS,0,PRIO_MAX-1); */
4145 
4146  MPI_Comm_size(MPI_COMM_WORLD,&np); /*total number of processes (take into
4147  account that this routine is process number 0)*/
4148 
4149  if (UpdateCuikSystem(p,cs))
4150  {
4151  switch(cs->searchMode)
4152  {
4153  case DEPTH_FIRST_SEARCH:
4154  InitHeapOfBoxes(CmpBoxDepthFirst,NULL,&boxes);
4155  break;
4156  case BREADTH_FIRST_SEARCH:
4157  InitHeapOfBoxes(CmpBoxBreadthFirst,NULL,&boxes);
4158  break;
4159  case MINIMIZATION_SEARCH:
4160  InitHeapOfBoxes(CmpBoxesEquation,(void *)cs,&boxes);
4161  break;
4162  }
4163 
4164  BoxFromVariables(&init_box,&(cs->variables));
4165 
4166  /* Reserve space for the boxes sent/received and associated information */
4167  /*All boxes are the same size as the initial one*/
4168  buffer_size=GetBoxBufferSize(&init_box);
4169 
4170  if ((restart)&&(fstate!=NULL))
4171  {
4172  Tlist pboxes;
4173 
4174  LoadCSState(fstate,&pboxes,cs);
4175  AddList2Heap(&pboxes,&boxes);
4176  DeleteListOfBoxes(&pboxes);
4177  }
4178  else
4179  {
4180  if (searchSpace==NULL)
4181  AddBox2HeapOfBoxes(&init_box,&boxes);
4182  else
4183  {
4184  Tbox lb;
4185 
4186  SimpleFromOriginal(searchSpace,&lb,cs->orig2sd);
4187  AddBox2HeapOfBoxes(&lb,&boxes);
4188  DeleteBox(&lb);
4189  }
4190 
4191  InitStatistics(np,HeapOfBoxesVolume(cs->systemVar,&boxes),&(cs->st));
4192  }
4193  DeleteBox(&init_box);
4194 
4195  current_level=0;
4196 
4197  NEW(buffers_out,np,double*);
4198  NEW(buffers_in,np,double*);
4199 
4200  NEW(request,np,MPI_Request);
4201  NEW(lost_packet,np,boolean);
4202  NEW(time_stamp,np,time_t);
4203 
4204  for(i=1;i<np;i++)
4205  {
4206  request[i]=MPI_REQUEST_NULL;
4207  NEW(buffers_out[i],buffer_size,double);
4208  NEW(buffers_in[i],buffer_size,double);
4209  lost_packet[i]=FALSE;
4210  }
4211 
4212  request[0]=MPI_REQUEST_NULL; /* Processor 0 is the scheduler, do not use it as a child */
4213  in_process=0;
4214  available=(np-1);
4215 
4216  nb=0;
4217  statePeriod=GetParameter(CT_STATE_PERIOD,p);
4218 
4219  nsols=(unsigned int)GetParameter(CT_N_SOLUTIONS,p);
4220 
4221  done=FALSE;
4222 
4223  /*the cuik solve process itself*/
4224  do {
4225 
4226  if ((!done)&& /*we don't have enough solutions*/
4227  (!HeapEmpty(&boxes))&& /* there is something to be processes */
4228  (available>0)) /* and someone ready to process it */
4229  {
4230  /* Search for a free child (at least one must exists) */
4231  i=1;
4232  while ((i<np)&&(request[i]!=MPI_REQUEST_NULL)) i++;
4233 
4234  if (i==np)
4235  Error("wrong number of available processors");
4236  else
4237  {
4238  if ((statePeriod>0)&&(fstate!=NULL))
4239  {
4240  nb++;
4241  if (nb==statePeriod)
4242  {
4243  Tbox btmp;
4244  unsigned int ctmp,ntmp;
4245  Tlist pboxes;
4246  unsigned int kk;
4247 
4248  Heap2List(&pboxes,&boxes);
4249  for(kk=1;kk<np;kk++)
4250  {
4251  if (request[kk]!=MPI_REQUEST_NULL)
4252  {
4253  InitBox(cs->nvariables,NULL,&btmp);
4254  Buffer2Box(&ctmp,&ntmp,buffers_out[kk],&btmp);
4255  AddFirstElement((void *)&btmp,&pboxes);
4256  }
4257  }
4258  SaveCSState(fstate,&pboxes,cs);
4259  DeleteListOfBoxes(&pboxes);
4260  nb=0;
4261  }
4262  }
4263 
4264  /*Get a new box*/
4265  ExtractMinElement(&b,&boxes);
4266 
4267  /* Transform the box into a buffer of doubles */
4268  Box2Buffer(REDUCED_BOX,0,buffers_out[i],&b);
4269 
4270  /* Send it to the iddle child processors */
4271  /* and start waiting for the reply (i.e., the reduced box) */
4272  if (MPI_Send(buffers_out[i],buffer_size,MPI_DOUBLE,i,20,MPI_COMM_WORLD)==MPI_SUCCESS)
4273  {
4274  /* The box is sent and now we start waiting for it */
4275  if (MPI_Irecv(buffers_in[i],buffer_size,MPI_DOUBLE,i,20,MPI_COMM_WORLD,&(request[i]))==MPI_SUCCESS)
4276  {
4277 
4278  current_level=GetBoxLevel(&b);
4279 
4280  NewBoxProcessed(&(cs->st));
4281  NewMaxLevel(current_level,&(cs->st));
4282 
4283  /* One more box is in process */
4284  in_process++;
4285 
4286  /* One less processor is available */
4287  available--;
4288 
4289  #if (_DEBUG>0)
4290  /* Inform about the box launchig */
4291  fprintf(stderr,"b<v:%g, s:%g, l:%u>-> %u\n",
4292  GetBoxVolume(cs->systemVar,&b),
4293  GetBoxSize(cs->systemVar,&b),
4294  current_level,i);
4295  /*
4296  fprintf(stderr," A:%d P:%d Q:%u\n",
4297  available,in_process,HeapSize(&boxes));
4298  */
4299  #endif
4300 
4301  /* Store the time at which we sent out the box */
4302  time_stamp[i]=time(NULL);
4303  /* and, mark the sent box as not lost */
4304  lost_packet[i]=FALSE;
4305 
4306  #if (_DEBUG>1)
4307  printf("Sending box to processors : %u (iddle -> %u box_in_list-> %u)\n",
4308  i,np-in_process-1,HeapSize(&boxes));
4309  printf(" Box: ");
4310  PrintBox(stdout,&b);
4311  #endif
4312 
4313  DeleteBox(&b);
4314  }
4315  else
4316  {
4317  /*The box will never be back: add it to the list and free the processor */
4318  #if (_DEBUG>0)
4319  fprintf(stderr,"RF-> %u\n",i);
4320  #endif
4321 
4322  AddBox2HeapOfBoxes(&b,&boxes);
4323  request[i]=MPI_REQUEST_NULL;
4324  }
4325  }
4326  else
4327  {
4328  /* We didn't manage to send out the box, just return it to the boxes
4329  to be processes*/
4330  #if (_DEBUG>0)
4331  fprintf(stderr,"SF-> %u\n",i);
4332  #endif
4333 
4334  AddBox2HeapOfBoxes(&b,&boxes);
4335  }
4336  }
4337  }
4338 
4339  /* see if we already got any reply from the slaves :-) */
4340  if ((MPI_Testany(np,request,&i,&completed,&status)==MPI_SUCCESS)&&
4341  (completed))
4342  {
4343  /* if 'completed' 'i' is the request completed */
4344  #if (_DEBUG>1)
4345  printf("Receiving box from processors : %u\n",i);
4346  #endif
4347 
4348  available++; /* either if lost or not, the processor becomes available from now on */
4349  request[i]=MPI_REQUEST_NULL;
4350 
4351  if (!lost_packet[i])
4352  {
4353  /* If so, reconstruct the output box and analyze it */
4354  InitBox(cs->nvariables,NULL,&b);
4355 
4356  Buffer2Box(&c,&nr,buffers_in[i],&b);
4357 
4358  AddNBoxReductions(nr,&(cs->st));
4359 
4360  in_process--; /* lost ones are already discounted when classified as lost
4361  (to be able to finish even if lost ones exits)*/
4362 
4363  #if (_DEBUG>0)
4364  fprintf(stderr," %u->b<v:%g s:%g l:%u>(t=%lu): ",
4365  i,
4366  GetBoxVolume(cs->systemVar,&b),
4367  GetBoxSize(cs->systemVar,&b),
4368  GetBoxLevel(&b),
4369  time(NULL)-time_stamp[i]);
4370  /*
4371  fprintf(stderr," A:%d P:%d Q:%u\n",
4372  available,in_process,HeapSize(&boxes));
4373  */
4374  #endif
4375 
4376  PostProcessBox(p,c,f_out,NULL,&boxes,&b,cs);
4377 
4378  if (nsols>0)
4379  done=(GetNSolutionBoxes(&(cs->st))>=nsols);
4380 
4381  DeleteBox(&b);
4382  }
4383  #if (_DEBUG>0)
4384  else
4385  {
4386  /* We receive a box that was thought to be lost. The only think we
4387  do is to mark the child as available again by setting request[i]=MPI_REQUEST_NULL;
4388  below */
4389  fprintf(stderr," lb[%u](t=%lu > %u)\n",i,
4390  time(NULL)-time_stamp[i],MPI_TREAT_BOX_TIMEOUT(cs));
4391  }
4392  #endif
4393  }
4394 
4395 
4396  /* See if one of the send boxes is lost.
4397  All boxes sent before the deadline are considered lost.
4398  Lost boxes are recovered and added to the list of pending boxes.
4399  The processor to which the box was assigned is discarted for future
4400  jobs.
4401  */
4402  deadline=time(NULL)-(time_t)MPI_TREAT_BOX_TIMEOUT(cs);
4403  for(i=1;i<np;i++)
4404  {
4405  if ((!lost_packet[i])&& /*not already lost*/
4406  (request[i]!=MPI_REQUEST_NULL)&& /*and we are waiting for a reply*/
4407  (time_stamp[i]<deadline)) /*for too long actually*/
4408  {
4409  /* we considere the sub-task as lost, we add the box for further process by
4410  another computer and we mark the current message as lost */
4411 
4412  lost_packet[i]=TRUE; /* mark the packet as lost */
4413 
4414  InitBox(cs->nvariables,NULL,&b); /*Recover it and put it again in the list of
4415  boxes to be processed*/
4416  Buffer2Box(&c,&nr,buffers_out[i],&b); /* Error code (c) and additional information
4417  (nr) are empty (this is a box whose process
4418  was not concluded)*/
4419 
4420  NewLostBox(&(cs->st));
4421  #if (_DEBUG>0)
4422  fprintf(stderr," l[%u] (elapsed %lu)\n",i,time(NULL)-time_stamp[i]);
4423  #endif
4424 
4425  /* a time out tipically means a error in the simplex */
4426  PostProcessBox(p,ERROR_IN_PROCESS,f_out,NULL,&boxes,&b,cs);
4427 
4428  DeleteBox(&b);
4429 
4430  in_process--; /* one less box is in process (Observe that the processors is
4431  not marked as available) */
4432  }
4433  }
4434 
4435  fflush(stderr);
4436 
4437  /* Repeat until we don't have enough solutions and there are still some boxes from
4438  where we can get solutions.
4439  If we already have all the desired solutions but some boxes are in_process, we
4440  wait untill their process finishes (so that all processors are ready to receive
4441  the kill signal)
4442  */
4443  } while(((!done)&&(!HeapEmpty(&boxes)))||(in_process>0));
4444 
4445  /* Send a magic packed to the child so that they kill themselves */
4446 
4447  for(i=1;i<np;i++)
4448  {
4449  if (request[i]!=MPI_REQUEST_NULL)
4450  {
4451  /*we have a pending communication with this processor -> cancel it*/
4452  MPI_Cancel(&(request[i]));
4453  }
4454  else
4455  {
4456  /*nothing pending with this processors, just order it to die*/
4457  buffers_out[i][0]=-1.0;
4458  MPI_Send(buffers_out[i],buffer_size,MPI_DOUBLE,i,20,MPI_COMM_WORLD);
4459  }
4460  }
4461 
4462  free(request);
4463  free(lost_packet);
4464  free(time_stamp);
4465 
4466  DeleteHeap(&boxes);
4467 
4468  if (f_out!=NULL)
4469  {
4470  fprintf(f_out,"\n\nCuik executed in %u (%u) child processors\n",available,np-1);
4471  PrintStatistics(f_out,&(cs->st));
4472  }
4473 
4474  DeleteStatistics(&(cs->st));
4475 
4476  for(i=1;i<np;i++)
4477  {
4478  free(buffers_in[i]);
4479  free(buffers_out[i]);
4480  }
4481 
4482  free(buffers_in);
4483  free(buffers_out);
4484  }
4485 }
4486 
4487 /*
4488  This is the process executed by each child process
4489 */
4491 {
4492  if (!cs->scalar)
4493  Error("MPI_TreatBox only operates on scalar systems");
4494 
4495  if (UpdateCuikSystem(p,cs))
4496  {
4497  boolean end;
4498  unsigned int n;
4499  double *buffer;
4500  unsigned int c,nr;
4501  Tbox b;
4502  MPI_Status status;
4503 
4504  InitBox(cs->nvariables,NULL,&b);
4505  n=GetBoxBufferSize(&b);
4506  NEW(buffer,n,double);
4507 
4508  /* A child can be in execution in computers not only devoted to cuik.
4509  We reduce the priority of the cuik process to allow other processes
4510  (those own by the owners of the machines) to run smoothly.
4511  This way we can enlarge the cluster of procesors to almost all
4512  computers in our Lab.
4513  In the future we have to develop programs to check which computers
4514  are on and to add them into the machine list (i.e., the members of
4515  the cluster)
4516  */
4517  /*setpriority(PRIO_PROCESS,0,PRIO_MAX); */
4518 
4519  end=FALSE; /* While the kill buffer is not received */
4520 
4521  while(!end)
4522  {
4523  /* Get new box from the master */
4524  if (MPI_Recv(buffer,n,MPI_DOUBLE,0,20,MPI_COMM_WORLD,&status)==MPI_SUCCESS)
4525  {
4526  /* The first double is used to send control messages between the
4527  scheduler and the slave. If negative, the slave dies*/
4528  if (buffer[0]<0)
4529  end=TRUE;
4530  else
4531  {
4532  /* If not the kill buffer, recover the box out of the buffer */
4533  Buffer2Box(&c,&nr,buffer,&b);
4534 
4535  ResetNBoxReductions(&(cs->st)); /* Let's count the box reductions for this
4536  box shrink only*/
4537 
4538  #if (_DEBUG>1)
4539  printf("Box from main processor:");
4540  PrintBox(stdout,&b);
4541  fflush(stdout);
4542  #endif
4543 
4544  /* And reduce it */
4545  c=ReduceBox(p,~DUMMY_VAR,&b,cs);
4546 
4547  #if (_DEBUG>1)
4548  printf("Box out of ReduceBox (%u):",c);
4549  PrintBox(stdout,&b);
4550  fflush(stdout);
4551  #endif
4552 
4553  /* Then pack the result into another buffer */
4554  Box2Buffer(c,GetNBoxReductions(&(cs->st)),buffer,&b);
4555 
4556  /* and send it back to the master */
4557 
4558  #if (_DEBUG>1)
4559  if (MPI_Send(buffer,n,MPI_DOUBLE,0,20,MPI_COMM_WORLD)!=MPI_SUCCESS)
4560  printf("Sending to master failed\n");
4561  #else
4562  MPI_Send(buffer,n,MPI_DOUBLE,0,20,MPI_COMM_WORLD);
4563  #endif
4564 
4565  }
4566  }
4567  #if (_DEBUG>1)
4568  else
4569  printf("Receive from master failed\n");
4570  #endif
4571  }
4572 
4573  DeleteBox(&b);
4574  free(buffer);
4575  }
4576 }
4577 #endif
4578 
4579 /*
4580  Returns a box defined by the ranges for the variables given by the user in the
4581  input file. This is the initial search space.
4582  */
4584 {
4585  BoxFromVariables(box,&(cs->orig_variables));
4586 }
4587 
4589 {
4590  if (!UpdateCuikSystem(p,cs))
4591  Error("Inconsistent system in GenerateSimpInitialBox");
4592 
4593  BoxFromVariables(box,&(cs->simp_variables));
4594 }
4595 
4597 {
4598  boolean ok;
4599 
4600  if (UpdateCuikSystem(p,cs))
4601  {
4602  unsigned int i;
4603  Tbox bInit;
4604  double epsilon,rho;
4605  Tinterval all;
4606 
4607  /* If we set the dummies to range [-INF,INF] they are not
4608  adjusted via crop. */
4609  NewInterval(-INF/2,INF/2,&all);
4610 
4611  epsilon=GetParameter(CT_EPSILON,p);
4612  rho=GetParameter(CT_RHO,p);
4613 
4614  InitBox(cs->orig_nvariables,NULL,&bInit);
4615 
4616  /* System and secondary vars already have value in the input box but
4617  dummies and cartesian variables probably not.
4618  These are set to inf range to avoid any limit
4619  in their value when regenerated (boxes commming from
4620  Newton can have extrange values in dummy variables. */
4621  for(i=0;i<cs->orig_nvariables;i++)
4622  {
4623  if ((IsDummyVariable(i,&(cs->orig_variables)))||
4624  (IsCartesianVariable(i,&(cs->orig_variables))))
4625  SetBoxInterval(i,&all,b);
4626  }
4627  DeleteBox(&bInit);
4628 
4629  ok=TRUE;
4630 
4631  /* First we define the value of the dummy variables from the
4632  system ones. They are related by dymmy equations and, given
4633  the system varibles, the dummy variables can be easily defined
4634  via crop. */
4635  for(i=0;((ok)&&(i<cs->orig_nequations));i++)
4636  {
4637  if (IsDummyEquation(i,&(cs->orig_equations)))
4638  ok=(CropEquation(i,DUMMY_VAR,epsilon,rho,b,&(cs->orig_variables),&(cs->orig_equations))!=EMPTY_BOX);
4639  }
4640 
4641  /*
4642  Now we try to define the cartesian variables from the system
4643  and dummy ones.
4644  The separating plane cartesian variables can not be properly bounded
4645  this way but the ones for body vertexes can (and they are the relevant
4646  ones)
4647  */
4648  for(i=0;((ok)&&(i<cs->orig_nequations));i++)
4649  {
4650  if (IsCoordEquation(i,&(cs->orig_equations)))
4651  ok=(CropEquation(i,CARTESIAN_VAR,epsilon,rho,b,&(cs->orig_variables),&(cs->orig_equations))!=EMPTY_BOX);
4652  }
4653  }
4654  else
4655  ok=FALSE;
4656 
4657  return(ok);
4658 }
4659 
4660 unsigned int RegenerateSolutionPoint(Tparameters *p,double *pt,double **rp,TCuikSystem *cs)
4661 {
4662  unsigned int i,k;
4663  Tbox b;
4664 
4665  if (!UpdateCuikSystem(p,cs))
4666  Error("Inconsistent system in RegenerateSolutionPoint");
4667 
4668  NEW(*rp,cs->orig_nvariables,double);
4669 
4670  k=0;
4671  for(i=0;i<cs->orig_nvariables;i++)
4672  {
4673  /*The input box only has values for the input variables*/
4674  if (cs->orig_systemVar[i])
4675  {
4676  (*rp)[i]=pt[k];
4677  k++;
4678  }
4679  else
4680  (*rp)[i]=0.0; /*Non system variables are set to 0. This value is not used */
4681  }
4682 
4683  InitBoxFromPoint(cs->orig_nvariables,*rp,&b);
4684 
4685  if (!RegenerateSolution(p,&b,cs))
4686  Error("Invalid solution point in RegenerateSolutionPoint");
4687 
4688  k=0;
4689  for(i=0;i<cs->orig_nvariables;i++)
4690  {
4691  /*The input box only has values for the input variables (system and secondary)*/
4692  if (!cs->orig_systemVar[i])
4693  (*rp)[i]=IntervalCenter(GetBoxInterval(i,&b));
4694  }
4695  DeleteBox(&b);
4696 
4697  return(cs->orig_nvariables);
4698 }
4699 
4700 
4702 {
4703  if (UpdateCuikSystem(p,cs))
4704  {
4705  BoxFromVariables(boxO,&(cs->orig_variables));
4706  UpdateOriginalFromSimple(boxS,boxO,cs->orig2sd);
4707  }
4708  else
4709  Error("Inconsistent input cuiksystem");
4710 }
4711 
4712 unsigned int RegenerateOriginalPoint(Tparameters *p,double *s,double **o,TCuikSystem *cs)
4713 {
4714  if (UpdateCuikSystem(p,cs))
4715  {
4718  }
4719  else
4720  Error("Inconsistent input cuiksystem");
4721 
4722  return(cs->orig_nvariables);
4723 }
4724 
4725 unsigned int GenerateSimplifiedPoint(Tparameters *p,double *o,double **s,TCuikSystem *cs)
4726 {
4727  if (UpdateCuikSystem(p,cs))
4729  else
4730  Error("Inconsistent input cuiksystem");
4731 
4732  return(cs->simp_nvariables);
4733 }
4734 
4735 /*
4736  Selects one dimension along which to split box 'b'. The dimension can
4737  be selected according to the size or according to the error that each
4738  variable induce in each equation.
4739  This is similar to ComputeSplitDimInt but the output is an index in
4740  the original system and not in the simplified one
4741  */
4743 {
4744  Tbox bsimp;
4745  unsigned int d;
4746 
4747  if (!cs->scalar)
4748  Error("ComputeSplitDim only operates on scalar systems");
4749 
4750  if (UpdateCuikSystem(p,cs))
4751  {
4752  SimpleFromOriginal(b,&bsimp,cs->orig2sd);
4753  d=GetVarIDInOriginal(ComputeSplitDimInt(p,&bsimp,cs),cs->orig2sd);
4754  DeleteBox(&bsimp);
4755  }
4756  else
4757  d=NO_UINT;
4758 
4759  return(d);
4760 }
4761 
4762 /*
4763  Returns TRUE if point stored in vector 'v' is included in box 'b'
4764  Dummy variables are not taken into account in this procedure.
4765 */
4767 {
4768  unsigned int i,k,nv;
4769  double *val;
4770  boolean in;
4771 
4772  nv=NVariables(&(cs->orig_variables));
4773 
4774  if (nv!=GetBoxNIntervals(b))
4775  Error("Wrong box size in PointInSystemBox");
4776 
4777  k=0;
4778  in=TRUE;
4779  for(i=0;((in)&&(i<nv));i++)
4780  {
4781  if (!IsDummyVariable(i,&(cs->orig_variables)))
4782  {
4783  val=(double *)GetVectorElement(k,v);
4784  if (val==NULL)
4785  Error("Vector with incorrect number of elements in PointInSystemBox");
4786  in=IsInside(*val,0.0,GetBoxInterval(i,b));
4787  k++;
4788  }
4789  }
4790  return(in);
4791 }
4792 
4793 void EvaluateCSEquations(double *p,double *r,TCuikSystem *cs)
4794 {
4796 }
4797 
4798 void EvaluateSimpCSEquations(Tparameters *pr,double *p,double *r,TCuikSystem *cs)
4799 {
4800  if (!UpdateCuikSystem(pr,cs))
4801  Error("Inconsistent input cuiksystem");
4802 
4804 }
4805 
4806 void EvaluateSubSetSimpCSEquations(Tparameters *pr,boolean *se,double *p,double *r,TCuikSystem *cs)
4807 {
4808  if (!UpdateCuikSystem(pr,cs))
4809  Error("Inconsistent input cuiksystem");
4810 
4812 }
4813 
4814 void EvaluateCSJacobian(double *p,double ***J,TCuikSystem *cs)
4815 {
4816  TJacobian Ja;
4817 
4818  InitJacobian(&(cs->orig_variables),&(cs->orig_equations),&Ja);
4819 
4821  EvaluateJacobian(p,*J,&Ja);
4822 
4823  DeleteJacobian(&Ja);
4824 }
4825 
4826 double ErrorInCSEquations(double *p,TCuikSystem *cs)
4827 {
4828  double *r,e;
4829  unsigned int neq;
4830 
4831  neq=NEqualityEquations(&(cs->orig_equations));
4832  if (neq>0)
4833  {
4834  NEW(r,neq,double);
4836  e=Norm(neq,r);
4837  free(r);
4838  }
4839  else
4840  e=0.0;
4841 
4842  return(e);
4843 }
4844 
4846 {
4847  double *r,e;
4848  unsigned int neq;
4849 
4850  if (!UpdateCuikSystem(pr,cs))
4851  Error("Inconsistent input cuiksystem");
4852 
4853  neq=NEqualityEquations(&(cs->simp_equations));
4854  if (neq>0)
4855  {
4856  NEW(r,neq,double);
4858  e=Norm(neq,r);
4859  free(r);
4860  }
4861  else
4862  e=0.0;
4863 
4864  return(e);
4865 }
4866 
4867 double EvaluateCSCost(Tparameters *p,boolean simp,double *s,void *cs)
4868 {
4869  double v;
4870 
4871  if (simp)
4872  v=EvaluateWholeEquation(s,&(((TCuikSystem *)cs)->eqMin));
4873  else
4874  v=EvaluateWholeEquation(s,&(((TCuikSystem *)cs)->orig_eqMin));
4875 
4876  return(v);
4877 }
4878 
4879 /*
4880  Takes the center of the box and computes the error of this
4881  point when replaced into the system of equations.
4882  We use this function to check the validity of the outputs
4883  of Cuik.
4884  */
4886 {
4887  double *p; /*central point of the box*/
4888  unsigned int i;
4889  double *r,maxError;
4890  //double e;
4891  unsigned int nv,neq;
4892 
4893  maxError=0;
4894 
4895  nv=NVariables(&(cs->orig_variables));
4896  neq=NEqualityEquations(&(cs->orig_equations));
4897  if ((nv>0)&&(neq>0))
4898  {
4899  if (nv!=GetBoxNIntervals(b))
4900  Error("Wrong box size in ErrorInSolution");
4901 
4902  NEW(p,nv,double);
4903 
4904  for(i=0;i<nv;i++)
4905  p[i]=IntervalCenter(GetBoxInterval(i,b));
4906 
4907  NEW(r,neq,double);
4909  maxError=Norm(neq,r);
4910  /*
4911  for(i=0;i<neq;i++)
4912  {
4913  e=fabs(r[i]);
4914  if (e>maxError)
4915  maxError=e;
4916  }
4917  */
4918  free(r);
4919  free(p);
4920  }
4921  return(maxError);
4922 }
4923 
4925 {
4926  unsigned int i,neq,nv;
4927  double maxError,*r,*p;
4928 
4929  nv=NVariables(&(cs->orig_variables));
4930 
4931  if (nv!=GetBoxNIntervals(b))
4932  Error("Wrong box size in ErrorInInequalities");
4933 
4935  if (neq>0)
4936  {
4937  NEW(p,nv,double);
4938  for(i=0;i<cs->orig_nvariables;i++)
4939  p[i]=IntervalCenter(GetBoxInterval(i,b));
4940 
4941  NEW(r,neq,double);
4943  maxError=MaxVector(neq,r);
4944  free(r);
4945  free(p);
4946  }
4947  else
4948  maxError=0.0;
4949 
4950  return(maxError);
4951 }
4952 
4953 boolean InequalitiesHoldOnPoint(double *p,TCuikSystem *cs)
4954 {
4955  unsigned int neq;
4956  boolean hold;
4957  double *r;
4958 
4960  if (neq>0)
4961  {
4962  NEW(r,neq,double);
4964  hold=(MaxVector(neq,r)==0.0);
4965  free(r);
4966  }
4967  else
4968  hold=TRUE;
4969 
4970  return(hold);
4971 }
4972 
4974 {
4975  unsigned int neq;
4976  boolean hold;
4977  double *r;
4978 
4979  if (!UpdateCuikSystem(pr,cs))
4980  Error("Inconsistent input cuiksystem");
4981 
4983  if (neq>0)
4984  {
4985  NEW(r,neq,double);
4987  hold=(MaxVector(neq,r)==0.0);
4988  free(r);
4989  }
4990  else
4991  hold=TRUE;
4992 
4993  return(hold);
4994 }
4995 
4997 {
4998  unsigned int neq;
4999  double e,*r;
5000 
5001  if (!UpdateCuikSystem(pr,cs))
5002  Error("Inconsistent input cuiksystem");
5003 
5005  if (neq>0)
5006  {
5007  NEW(r,neq,double);
5009  e=MaxVector(neq,r);
5010  free(r);
5011  }
5012  else
5013  e=0.0;
5014 
5015  return(e);
5016 }
5017 
5018 /*
5019  Prints the information stored in the cuiksystem
5020  It prints it in a form than can be parsed again.
5021 */
5023 {
5024  unsigned int nv;
5025  char **varNames;
5026 
5027  PrintVariables(f,&(cs->orig_variables));
5028 
5029  nv=NVariables(&(cs->orig_variables));
5030  NEW(varNames,nv,char *);
5031  GetVariableNames(varNames,&(cs->orig_variables));
5032  PrintEquations(f,varNames,&(cs->orig_equations));
5034  {
5035  fprintf(f,"\n[SEARCH]\n\n MIN ");
5036  PrintMonomials(f,varNames,&(cs->orig_eqMin));
5037  }
5038  free(varNames);
5039  fprintf(f,"\n\n");
5040 }
5041 
5042 /*
5043  Prints the information stored in the cuiksystem together with the simplified
5044  CuikSystem.
5045  It prints it in a form than can be parsed again.
5046 */
5048 {
5049  unsigned int nv;
5050  char **varNames;
5051 
5052  #if (_DEBUG>1)
5053  fprintf(f,"%%****************************************\n");
5054  fprintf(f,"%% Original system \n");
5055  PrintVariables(f,&(cs->orig_variables));
5056 
5057  nv=NVariables(&(cs->orig_variables));
5058  NEW(varNames,nv,char *);
5059  GetVariableNames(varNames,&(cs->orig_variables));
5060  PrintEquations(f,varNames,&(cs->orig_equations));
5062  {
5063  fprintf(f,"\n[SEARCH]\n\n MIN ");
5064  PrintMonomials(f,varNames,&(cs->orig_eqMin));
5065  }
5066  fprintf(f,"\n\n");
5067  free(varNames);
5068  #endif
5069 
5070  if (UpdateCuikSystem(p,cs))
5071  {
5072  fprintf(f,"%%****************************************\n");
5073  fprintf(f,"%% Simplified system \n");
5074  fprintf(f,"%% SIMPLIFICATION_LEVEL: %u\n",
5075  (unsigned int)(GetParameter(CT_SIMPLIFICATION_LEVEL,p)));
5076  fprintf(f,"%% Variable reduction %u -> %u\n",cs->orig_nvariables,
5077  cs->simp_nvariables);
5078  fprintf(f,"%% Num syst+secu variables in original : %u \n",
5081  fprintf(f,"%% Num syst+secu variables in simplified: %u \n",
5084  PrintMapping(f,cs->orig2sd);
5085  PrintVariables(f,&(cs->variables));
5086 
5087  nv=NVariables(&(cs->variables));
5088  NEW(varNames,nv,char *);
5089  GetVariableNames(varNames,&(cs->variables));
5090  PrintEquations(f,varNames,&(cs->equations));
5092  {
5093  fprintf(f,"\n[SEARCH]\n\n MIN ");
5094  PrintMonomials(f,varNames,&(cs->eqMin));
5095  }
5096  free(varNames);
5097  }
5098  else
5099  fprintf(f,"INCONSISTENT INPUT CUIK SYSTEM\n");
5100 }
5101 
5103 {
5104  if (UpdateCuikSystem(p,cs))
5105  SaveMapping(f,cs->orig2sd);
5106  else
5107  fprintf(f,"INCONSISTENT INPUT CUIK SYSTEM\n");
5108 }
5109 
5110 /*
5111  Destructor of the cuiksystem structure.
5112 */
5114 {
5115  DeleteConstants(&(cs->constants));
5116 
5119  DeleteStatistics(&(cs->st));
5120 
5121  DeleteEquation(&(cs->orig_eqMin));
5122 
5123  UnUpdateCuikSystem(cs);
5124 
5125  cs->empty=TRUE;
5126  cs->simp_empty=TRUE;
5127 }
unsigned int SimplexNRows(TSimplex *s)
Gets the number of rows (i.e., constraints) of the simplex structure.
Definition: simplex_clp.c:107
Definition of the boolean type.
boolean PointInBox(boolean *used, unsigned int n, double *v, double tol, Tbox *b)
Checks if a point is included in a(sub-) box.
Definition: box.c:332
void First(Titerator *i)
Moves an iterator to the first position of its associated list.
Definition: list.c:356
double MaxVector(unsigned int m, double *v)
Value of the maximum element of a vector.
double EvaluateEqMin(void *b, void *cs)
Evaluates the equation to minimize in a given box.
Definition: cuiksystem.c:2369
Set of variables of a cuiksystem.
Definition: variables.h:38
void PrintCuikSystemWithSimplification(Tparameters *p, FILE *f, TCuikSystem *cs)
Prints the simplified cuiksystem.
Definition: cuiksystem.c:5047
void Box2Buffer(unsigned int c, unsigned int n, double *buffer, Tbox *b)
Converts a box into an array of doubles.
Definition: box.c:81
void InitConstants(Tconstants *cts)
Initializes a constant set.
Definition: constants.c:21
void DeleteLinearConstraint(TLinearConstraint *lc)
Destructor.
#define CARTESIAN_VAR
One of the possible type of variables.
Definition: variable.h:62
void RewriteEquation2Simp(double epsilon, Tmapping *map, Tequation *eqOut, Tequation *eq)
Applies the simplification mapping to an equation.
Definition: equation.c:237
Tinterval * GetBoxInterval(unsigned int n, Tbox *b)
Returns a pointer to one of the intervals defining the box.
Definition: box.c:270
void SetCSVariableRange(unsigned int n, Tinterval *r, TCuikSystem *cs)
Gets the range of a variable from a cuiksystem.
Definition: cuiksystem.c:2574
double EvaluateWholeEquation(double *varValues, Tequation *eq)
Evaluates an equation in a given point.
Definition: equation.c:1615
#define ERROR_IN_PROCESS
One of the possible results of reducing a box.
Definition: box.h:32
void EvaluateSubSetSimpCSEquations(Tparameters *pr, boolean *se, double *p, double *r, TCuikSystem *cs)
Evaluates a subset of the simplified equation set on a point.
Definition: cuiksystem.c:4806
#define SYSTEM_EQ
One of the possible type of equations.
Definition: equation.h:146
Tequation * GetEquation(unsigned int n, Tequations *eqs)
Gets an equation from the set.
Definition: equations.c:1697
boolean IsCSPolynomial(TCuikSystem *cs)
Identifies polynomial cuiksystems.
Definition: cuiksystem.c:2651
Tmonomial * GetMonomial(unsigned int i, Tequation *eq)
Gets a monomial from an equation.
Definition: equation.c:1584
void InitJacobian(Tvariables *vs, Tequations *eqs, TJacobian *j)
Constructor.
Definition: jacobian.c:16
#define CT_EPSILON
Numerical tolerance.
Definition: parameters.h:194
unsigned int NEquations(Tequations *eqs)
Number of equations in the set.
Definition: equations.c:1075
#define FALSE
FALSE.
Definition: boolean.h:30
void SaveMapping(FILE *f, Tmapping *m)
Saves the mapping information into a file.
Definition: csmapping.c:291
boolean SampleCuikSystem(Tparameters *p, char *fname, Tlist *sb, unsigned int nsamples, unsigned int ntries, unsigned int ndof, TCuikSystem *cs)
Generates samples for a cuiksystem.
Definition: cuiksystem.c:2941
double GetTime(Tstatistics *t)
Curent time.
Definition: statistics.c:74
#define CONVERGED_OUTSIDE_GLOBAL
One of the possible outputs of the Newton iteration.
Definition: cuiksystem.h:150
boolean IsSecondaryVariable(unsigned int n, Tvariables *vs)
Identifies secondary variables in a set.
Definition: variables.c:108
unsigned int NVariables(Tvariables *vs)
Gets the number of variables in a set.
Definition: variables.c:69
unsigned int ComputeSplitDimInt(Tparameters *p, Tbox *b, TCuikSystem *cs)
Computes the optimal split dimension for a given box.
Definition: cuiksystem.c:1972
boolean SimplifiedMEquation(TMequation *me)
Identifies simplified equations.
Definition: mequation.c:294
unsigned int GetVariableTypeN(unsigned int n, Tvariables *vs)
Gets the type of a particular variable.
Definition: variables.c:123
void SetEquationType(unsigned int type, Tequation *eq)
Changes the type of the equation (SYSTEM_EQ, CARTESIAN_EQ, DUMMY_EQ, DERIVED_EQ). ...
Definition: equation.c:1013
unsigned int GetCSNumEquations(TCuikSystem *cs)
Gets the number of equations already in the cuiksystem.
Definition: cuiksystem.c:2664
boolean IsCoordEquation(unsigned int i, Tequations *eqs)
Identify coordenalization equations.
Definition: equations.c:1126
void GetVariableNames(char **varNames, Tvariables *vs)
Gets the name for all the variables in the set.
Definition: variables.c:234
#define NEW(_var, _n, _type)
Allocates memory space.
Definition: defines.h:385
double GetBoxSize(boolean *used, Tbox *b)
Computes the size of the box.
Definition: box.c:607
Data structure to hold the information about the name of a file.
Definition: filename.h:248
#define CUT_POINT
Point, relative to the size of the selected box side, where we split a box.
Definition: cuiksystem.h:88
void ArrayPi2Pi(unsigned int n, unsigned int *t, double *a)
Applies PI2PI to an array.
unsigned int ReduceBoxEquationWise(Tparameters *p, Tbox *b, TCuikSystem *cs)
Reduces a box considering one equation at a time.
Definition: cuiksystem.c:714
void NewLostBox(Tstatistics *t)
Increases the number of lost boxes.
Definition: statistics.c:173
void EvaluateEquationInt(Tinterval *varValues, Tinterval *i_out, Tequation *eq)
Interval evaluation of an equation.
Definition: equation.c:1649
A linear constraint with an associated error.
Tequations orig_equations
Definition: cuiksystem.h:236
boolean * systemVar
Definition: cuiksystem.h:206
void DeleteEquation(Tequation *eq)
Destructor.
Definition: equation.c:1748
unsigned int GetVarIDInOriginal(unsigned int v, Tmapping *m)
Gets the original identifier of a simplified variable.
Definition: csmapping.c:135
#define CT_SMALL_SIGMA
Box size threshold.
Definition: parameters.h:229
void EvaluateJacobian(double *v, double **m, TJacobian *j)
Evaluates the Jacobian.
Definition: jacobian.c:79
boolean IsCartesianVariable(unsigned int n, Tvariables *vs)
Identifies cartesian variables in a set.
Definition: variables.c:118
void GenerateInitialBox(Tbox *box, TCuikSystem *cs)
Gives the search space in the form of a box.
Definition: cuiksystem.c:4583
void * GetVectorElement(unsigned int i, Tvector *vector)
Returns a pointer to a vector element.
Definition: vector.c:269
void SaveStatistics(FILE *f, Tstatistics *t)
Saves the statistics to a file in binary format.
Definition: statistics.c:260
void AddMatrixEquation(TMequation *equation, Tequations *eqs)
Adds a matrix equation to the set.
Definition: equations.c:1665
#define CT_SPLIT_TYPE
Split type.
Definition: parameters.h:346
unsigned int randomMax(unsigned int m)
Returns a random integer in the range [0,m].
Definition: random.c:77
void SetBoxInterval(unsigned int n, Tinterval *is, Tbox *b)
Replaces a particular interval in a box.
Definition: box.c:259
void SaveCSState(char *fname, Tlist *lb, TCuikSystem *cs)
Saves internal the cuik solver state information to a file.
Definition: cuiksystem.c:2093
boolean IncrementalSampleCuikSystem(Tparameters *p, char *fname, Tlist *sb, boolean *fixVars, unsigned int nsamples, unsigned int ntries, unsigned int ndof, TCuikSystem *cs)
Generates samples for a cuiksystem.
Definition: cuiksystem.c:3218
void MergeEquations(Tequations *eqs1, Tequations *eqs)
Copy constructor.
Definition: equations.c:815
boolean * orig_notDummyVar
Definition: cuiksystem.h:248
Tvariables simp_variables
Definition: cuiksystem.h:223
Collection of methods to work on Theap of boxes.
void ResetNBoxReductions(Tstatistics *t)
Resets the number of reduced boxes.
Definition: statistics.c:196
#define SYSTEM_VAR
One of the possible type of variables.
Definition: variable.h:24
void InitStatistics(unsigned int np, double v, Tstatistics *t)
Constructor.
Definition: statistics.c:21
CBLAS_INLINE double Norm(unsigned int s, double *v)
Computes the norm of a vector.
unsigned int AddVariable2CS(Tvariable *v, TCuikSystem *cs)
Adds a variable to the system.
Definition: cuiksystem.c:2511
void AddBox2HeapOfBoxes(Tbox *b, Theap *h)
Adds a box to a heap of boxes.
Definition: box_heap.c:26
void LoadStatistics(FILE *f, Tstatistics *t)
Loads the statistics from a file in binary format.
Definition: statistics.c:266
void CopyEquation(Tequation *eq_dst, Tequation *eq_orig)
Copy constructor.
Definition: equation.c:216
void CopyEquations(Tequations *eqs_dst, Tequations *eqs_src)
Copy constructor.
Definition: equations.c:768
double GetBoxSumSide(boolean *used, Tbox *b)
Computes the sum of the sides of the box.
Definition: box.c:974
void NewEmptyBox(Tstatistics *t)
Increases the number of empty boxes.
Definition: statistics.c:147
char ** orig_varNames
Definition: cuiksystem.h:250
#define CT_LR2TM_S
Threshold to switch from linear relaxations to Taylor models for saddle equations.
Definition: parameters.h:398
#define DEPTH_FIRST_SEARCH
Depth first search.
Definition: cuiksystem.h:54
void DeleteJacobian(TJacobian *j)
Destructor.
Definition: jacobian.c:249
void NewBoxReduction(Tstatistics *t)
Increases the number of reduced boxes.
Definition: statistics.c:181
double HeapOfBoxesVolume(boolean *used, Theap *h)
Computes the volume of a heap.
Definition: box_heap.c:31
double ErrorInSimpCSEquations(Tparameters *pr, double *p, TCuikSystem *cs)
Evaluates the norm of the error in a point for the simplified equations.
Definition: cuiksystem.c:4845
#define EQU
In a Tequation, the equation relational operator is equal.
Definition: equation.h:201
void SetOriginalVarRelation(unsigned int nvo, TLinearConstraint *lc, Tmapping *m)
Set a relation for a variable in the original set.
Definition: csmapping.c:100
#define DUMMY_EQ
One of the possible type of equations.
Definition: equation.h:164
void InitLinearConstraint(TLinearConstraint *lc)
Constructor.
void SaveCuikSystemSimplification(Tparameters *p, FILE *f, TCuikSystem *cs)
Saves the simplification information associated with a cuiksystem.
Definition: cuiksystem.c:5102
void AddSimplifiedJacobianEquations(Tparameters *p, boolean *selectedVars, TCuikSystem *cs)
Adds linear a linear combination of the Jacobian to the system.
Definition: cuiksystem.c:2788
void UpdateSplitWeight(unsigned int ne, double *splitDim, Tbox *b, Tequations *eqs)
Computes the linearization error induced by the variables of a given equation.
Definition: equations.c:2430
void DeleteCuikSystem(TCuikSystem *cs)
Destructor.
Definition: cuiksystem.c:5113
boolean LinearConstraintIncludes(unsigned int ind, TLinearConstraint *lc)
Checks if a variable is included in a linear constraint.
unsigned int GetCSSystemVars(boolean **sv, TCuikSystem *cs)
Identifies the system variables.
Definition: cuiksystem.c:2614
void UnUpdateCuikSystem(TCuikSystem *cs)
Removes the cached information for the cuiksystem.
Definition: cuiksystem.c:1763
Tequation * GetJacobianEquation(unsigned int r, unsigned int c, TJacobian *j)
Returns one element of the Jacobian.
Definition: jacobian.c:57
char * GetCSVariableName(unsigned int id, TCuikSystem *cs)
Gets a variable name.
Definition: cuiksystem.c:2591
#define EMPTY_BOX
One of the possible results of reducing a box.
Definition: box.h:25
void AllocateJacobianEvaluation(double ***m, TJacobian *j)
Allocate space for the Jacobian evaluation.
Definition: jacobian.c:65
Tequations equations
Definition: cuiksystem.h:201
boolean IsSimplificable(unsigned int simpLevel, unsigned int nTerms, boolean *systemVars, Tbox *cb, unsigned int *v, TLinearConstraint *lc, Tequation *eq)
Identify equations than can trigger variable simplifications.
Definition: equation.c:840
unsigned int GetNBoxReductions(Tstatistics *t)
Gets the number of reduced boxes.
Definition: statistics.c:186
unsigned int ReduceBox(Tparameters *p, unsigned int varMask, Tbox *b, TCuikSystem *cs)
Reduces a box as much as possible.
Definition: cuiksystem.c:833
Definition of the Tfilename type and the associated functions.
#define BREADTH_FIRST_SEARCH
Breadth first search.
Definition: cuiksystem.h:63
void InitVariables(Tvariables *vs)
Constructor.
Definition: variables.c:24
unsigned int GetLinearConstraintVariable(unsigned int i, TLinearConstraint *lc)
Gets the a particular variable index.
void RegenerateOriginalBox(Tparameters *p, Tbox *boxS, Tbox *boxO, TCuikSystem *cs)
Generates a box in the original cuiksystem from a box of the simplified one.
Definition: cuiksystem.c:4701
unsigned int NInequalityEquations(Tequations *eqs)
Number of inequalities in the set.
Definition: equations.c:1100
void PrintLinearConstraint(FILE *f, boolean eq, char **varName, TLinearConstraint *lc)
Prints a linear constraint.
void SetVariableInterval(Tinterval *i, Tvariable *v)
Sets the new range for the variable.
Definition: variable.c:70
#define TRUE
TRUE.
Definition: boolean.h:21
Mapping between the sets of variables in two different cuiksystems.
Definition: csmapping.h:53
double * GetNewtonRHBuffer(TNewton *n)
Buffer to store the Newton right hand.
Definition: algebra.c:1131
boolean CmpBoxBreadthFirst(void *b1, void *b2, void *userData)
Determines which box to explore first in breadth first mode.
Definition: box.c:1110
unsigned int searchMode
Definition: cuiksystem.h:193
void InitBox(unsigned int dim, Tinterval *is, Tbox *b)
Initializes a box.
Definition: box.c:23
#define CT_SAFE_SIMPLEX
Trade off between speed and numerical stability when using the simplex.
Definition: parameters.h:358
void PointFromVariables(double **p, Tvariables *vs)
Creates a point from the center of the ranges of a set of variables.
Definition: variables.c:294
Tmapping * orig2sd
Definition: cuiksystem.h:199
void InitEquation(Tequation *eq)
Constructor.
Definition: equation.c:86
void Error(const char *s)
General error function.
Definition: error.c:80
void AddMatrixEquation2CS(Tparameters *p, TMequation *eq, TCuikSystem *cs)
Adds a matrix equation to the system.
Definition: cuiksystem.c:2490
#define NFUN
No trigonometric function for the variable.
Definition: variable_set.h:36
#define TOPOLOGY_R
One of the possible topologies.
Definition: defines.h:122
void SetBoxSubset(boolean *used, Tbox *bset, Tbox *b)
Changes a sub-set of ranges in a given box.
Definition: box.c:238
void RemoveVariable(unsigned int n, Tvariables *vs)
Removes a variable from a set.
Definition: variables.c:257
boolean AddEquation2Simplex(unsigned int ne, double lr2tm_q, double lr2tm_s, double epsilon, unsigned int safeSimplex, double rho, Tbox *b, Tvariables *vs, TSimplex *lp, Tequations *eqs)
Adds an equation to the simplex.
Definition: equations.c:2185
void GetJacobianSize(unsigned int *nr, unsigned int *nc, TJacobian *j)
Returns the size of the Jacobian.
Definition: jacobian.c:43
void CuikSystemMerge(Tparameters *p, TCuikSystem *cs1, TCuikSystem *cs2, TCuikSystem *cs)
Produces the union of two cuik systems.
Definition: cuiksystem.c:2320
void CSRemoveUnusedVars(Tparameters *p, TCuikSystem *cs)
Removes non-used variables.
Definition: cuiksystem.c:1260
unsigned int orig_nequations
Definition: cuiksystem.h:239
#define CT_SIGMA
Box size threshold.
Definition: parameters.h:236
boolean simp_empty
Definition: cuiksystem.h:220
#define MINIMIZATION_SEARCH
Search based on a minimum value of a given equation.
Definition: cuiksystem.h:76
boolean UpdateCuikSystem(Tparameters *p, TCuikSystem *cs)
Updates the cached information for the cuiksystem.
Definition: cuiksystem.c:1844
void GetCSVariable(unsigned int n, Tvariable *v, TCuikSystem *cs)
Gets the a variable from a cuiksystem.
Definition: cuiksystem.c:2566
void EvaluateCSEquations(double *p, double *r, TCuikSystem *cs)
Evaluates the equation set on a point.
Definition: cuiksystem.c:4793
#define CONVERGED_IN_BOX
One of the possible outputs of the Newton iteration.
Definition: cuiksystem.h:124
A simplex tableau structure.
Definition: simplex.h:73
unsigned int CropEquation(unsigned int ne, unsigned int varType, double epsilon, double rho, Tbox *b, Tvariables *vs, Tequations *eqs)
Equation-wise crop.
Definition: equations.c:1175
#define ZERO
Floating point operations giving a value below this constant (in absolute value) are considered 0...
Definition: defines.h:37
void AddTerm2LinearConstraint(unsigned int ind, double val, TLinearConstraint *lc)
Adds a scaled variable to the linear constraint.
unsigned int GetBoxBufferSize(Tbox *b)
Returns the size of a box when converted to an array of doubles.
Definition: box.c:76
Matrix equation.
Definition: mequation.h:38
void SetEquationValue(double v, Tequation *eq)
Changes the right-hand value of the equation.
Definition: equation.c:1026
void AddNBoxReductions(unsigned int nr, Tstatistics *t)
Increases the number of reduced boxes.
Definition: statistics.c:201
Collection of methods to work on Tlist of boxes.
void VariablesFromBox(Tbox *b, Tvariables *vs)
Define the range for the variables from a box.
Definition: variables.c:305
unsigned int GetSimpCSTopology(Tparameters *p, unsigned int **t, TCuikSystem *cs)
Topology of the variables in the simplified system.
Definition: cuiksystem.c:2674
void AddMonomial(Tmonomial *f, Tequation *eq)
Adds a new monomial to the equation.
Definition: equation.c:1356
void SimplexCreate(double epsilon, unsigned int ncols, TSimplex *s)
Constructor.
Definition: simplex_clp.c:22
void CopyConstants(Tconstants *cts_dst, Tconstants *cts_src)
Copies a set of constants.
Definition: constants.c:30
void DeleteMapping(Tmapping *m)
Destructor.
Definition: csmapping.c:310
boolean VarIncluded(unsigned int id, Tvariable_set *vs)
Checks if a variable index is included in a set of variable indexes.
Definition: variable_set.c:183
unsigned int GetVariableFunctionN(unsigned int n, Tvariable_set *vs)
Gets a variable function from a variable set.
Definition: variable_set.c:474
void AddTerm2SearchCriterion(double w, unsigned int v, double val, TCuikSystem *cs)
Adds penalty terms to the search criterion.
Definition: cuiksystem.c:2438
unsigned int GetBoxNIntervals(Tbox *b)
Box dimensionality.
Definition: box.c:992
void CopyCuikSystem(TCuikSystem *cs_dst, TCuikSystem *cs_src)
Copy constructor.
Definition: cuiksystem.c:2204
void MPI_TreatBox(Tparameters *p, TCuikSystem *cs)
Determines the solution set for a cuiksystem. Child process.
Definition: cuiksystem.c:4490
boolean IsSystemVariable(unsigned int n, Tvariables *vs)
Identifies system variables in a set.
Definition: variables.c:103
unsigned int GetBoxLevel(Tbox *b)
Returns the box level.
Definition: box.c:676
int NewtonStep(double nullSingularValue, double *x, double *dif, TNewton *n)
One step in a Newton iteration.
unsigned int GetVariableN(unsigned int n, Tvariable_set *vs)
Gets a variable identifier from a variable set.
Definition: variable_set.c:447
Error and warning functions.