00001 #include "linear_constraint.h"
00002
00003 #include "defines.h"
00004 #include "interval.h"
00005
00006 #include <string.h>
00007 #include <math.h>
00008
00017 void InitLinearConstraint(TLinearConstraint *lc)
00018 {
00019 lc->max=INIT_NUM_TERMS_LC;
00020 NEW(lc->ind,lc->max,unsigned int);
00021 NEW(lc->val,lc->max,double);
00022 ResetLinearConstraint(lc);
00023 }
00024
00025 void ResetLinearConstraint(TLinearConstraint *lc)
00026 {
00027 lc->n=0;
00028 NewInterval(0,0,&(lc->error));
00029 }
00030
00031 void CopyLinearConstraint(TLinearConstraint *lc_dst,TLinearConstraint *lc_src)
00032 {
00033 lc_dst->n=lc_src->n;
00034 lc_dst->max=lc_src->max;
00035 NEW(lc_dst->ind,lc_src->max,unsigned int);
00036 NEW(lc_dst->val,lc_src->max,double);
00037
00038 memcpy(lc_dst->ind,lc_src->ind,lc_dst->n*sizeof(unsigned int));
00039 memcpy(lc_dst->val,lc_src->val,lc_dst->n*sizeof(double));
00040
00041 CopyInterval(&(lc_dst->error),&(lc_src->error));
00042 }
00043
00044 unsigned int GetNumTermsInLinearConstraint(TLinearConstraint *lc)
00045 {
00046 return(lc->n);
00047 }
00048
00049 double *GetLinearConstraintCoefficients(TLinearConstraint *lc)
00050 {
00051 return(lc->val);
00052 }
00053
00054 double GetLinearConstraintCoefficient(unsigned int i,TLinearConstraint *lc)
00055 {
00056 if (i<lc->n)
00057 return(lc->val[i]);
00058 else
00059 {
00060 Error("Index out of range in GetLinearConstraintCoefficient");
00061 return(0);
00062 }
00063 }
00064
00065 unsigned int *GetLinearConstraintVariables(TLinearConstraint *lc)
00066 {
00067 return(lc->ind);
00068 }
00069
00070 unsigned int GetLinearConstraintVariable(unsigned int i,TLinearConstraint *lc)
00071 {
00072 if (i<lc->n)
00073 return(lc->ind[i]);
00074 else
00075 {
00076 Error("Index out of range in GetLinearConstraintCoefficient");
00077 return(NO_UINT);
00078 }
00079 }
00080
00081 void GetLinearConstraintError(Tinterval *error,TLinearConstraint *lc)
00082 {
00083 CopyInterval(error,&(lc->error));
00084 }
00085
00086 void SetLinearConstraintError(Tinterval *error,TLinearConstraint *lc)
00087 {
00088 CopyInterval(&(lc->error),error);
00089 }
00090
00091 void AddCt2LinearConstraint(double ct,TLinearConstraint *lc)
00092 {
00093 IntervalOffset(&(lc->error),-ct,&(lc->error));
00094 }
00095
00096 void AddTerm2LinearConstraint(unsigned int ind,double val,TLinearConstraint *lc)
00097 {
00098 if (val!=0.0)
00099 {
00100 boolean found;
00101 unsigned int i;
00102
00103 found=FALSE;
00104 i=0;
00105 while((!found)&&(i<lc->n))
00106 {
00107 found=(lc->ind[i]==ind);
00108 if (!found) i++;
00109 }
00110 if (!found)
00111 {
00112 if (lc->n==lc->max)
00113 {
00114 MEM_DUP(lc->ind,lc->max,unsigned int);
00115 MEM_EXPAND(lc->val,lc->max,double);
00116 }
00117 lc->ind[lc->n]=ind;
00118 lc->val[lc->n]=val;
00119 lc->n++;
00120 }
00121 else
00122 lc->val[i]+=val;
00123 }
00124 }
00125
00126 double RemoveTermFromLinearConstraint(unsigned int ind,TLinearConstraint *lc)
00127 {
00128 unsigned int i,j;
00129 boolean found;
00130 double ct;
00131
00132 found=FALSE;
00133 i=0;
00134 while((!found)&&(i<lc->n))
00135 {
00136 found=(lc->ind[i]==ind);
00137 if (!found) i++;
00138 }
00139
00140 if (found)
00141 {
00142 ct=lc->val[i];
00143 for(j=i+1;j<lc->n;j++)
00144 {
00145 lc->ind[j-1]=lc->ind[j];
00146 lc->val[j-1]=lc->val[j];
00147 }
00148 lc->n--;
00149 }
00150 else
00151 ct=0.0;
00152
00153 return(ct);
00154 }
00155
00156 boolean LinearConstraintIncludes(unsigned int ind,TLinearConstraint *lc)
00157 {
00158 unsigned int i;
00159 boolean found;
00160
00161 found=FALSE;
00162 i=0;
00163 while((!found)&&(i<lc->n))
00164 {
00165 found=(lc->ind[i]==ind);
00166 if (!found) i++;
00167 }
00168
00169 return(found);
00170 }
00171
00172 void InvertLinearConstraint(TLinearConstraint *lc)
00173 {
00174 unsigned int i;
00175
00176 for(i=0;i<lc->n;i++)
00177 lc->val[i]=-lc->val[i];
00178
00179 IntervalInvert(&(lc->error),&(lc->error));
00180 }
00181
00182 void ScaleLinearConstraint(double a,TLinearConstraint *lc)
00183 {
00184 unsigned int i;
00185 TLinearConstraint lc2;
00186
00187 InitLinearConstraint(&lc2);
00188
00189 for(i=0;i<lc->n;i++)
00190 AddTerm2LinearConstraint(lc->ind[i],lc->val[i]*a,&lc2);
00191
00192 IntervalScale(&(lc->error),a,&(lc2.error));
00193
00194 DeleteLinearConstraint(lc);
00195 CopyLinearConstraint(lc,&lc2);
00196 DeleteLinearConstraint(&lc2);
00197 }
00198
00199 void AddLinearConstraints(TLinearConstraint *lc1,TLinearConstraint *lc2)
00200 {
00201 unsigned int i;
00202
00203 IntervalAdd(&(lc1->error),&(lc2->error),&(lc2->error));
00204
00205 for(i=0;i<lc1->n;i++)
00206 AddTerm2LinearConstraint(lc1->ind[i],lc1->val[i],lc2);
00207 }
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222 void CleanLinearConstraint(double epsilon,Tinterval *is,TLinearConstraint *lc)
00223 {
00224 TLinearConstraint lcInt;
00225 Tinterval error,r;
00226 unsigned int i,k;
00227
00228 InitLinearConstraint(&lcInt);
00229
00230 CopyInterval(&error,&(lc->error));
00231
00232 for(i=0;i<lc->n;i++)
00233 {
00234 k=lc->ind[i];
00235
00236 if ((fabs(lc->val[i])<=epsilon)||
00237 (IntervalSize(&(is[k]))<=epsilon))
00238 {
00239 IntervalScale(&(is[k]),-lc->val[i],&r);
00240 IntervalAdd(&r,&error,&error);
00241 }
00242 else
00243 AddTerm2LinearConstraint(k,lc->val[i],&lcInt);
00244 }
00245
00246 CopyInterval(&(lcInt.error),&error);
00247
00248 DeleteLinearConstraint(lc);
00249 CopyLinearConstraint(lc,&lcInt);
00250 DeleteLinearConstraint(&lcInt);
00251 }
00252
00253 boolean SimplifyLinearConstraint(boolean *full,Tinterval *is,TLinearConstraint *lc)
00254 {
00255 boolean simplified;
00256
00257
00258
00259
00260
00261 *full=TRUE;
00262 simplified=FALSE;
00263
00264 if (lc->n==1)
00265 {
00266 Tinterval a,one;
00267 unsigned int k;
00268
00269 simplified=TRUE;
00270
00271 NewInterval(1,1,&one);
00272
00273
00274
00275 NewInterval(lc->val[0]-ZERO,lc->val[0]-ZERO,&a);
00276 IntervalDivision(&one,&a,&a);
00277 IntervalProduct(&a,&(lc->error),&a);
00278
00279 k=lc->ind[0];
00280
00281 *full=Intersection(&a,&(is[k]),&(is[k]));
00282 }
00283
00284 return(simplified);
00285 }
00286
00287 unsigned int CropLinearConstraint(double epsilon,Tbox *b,
00288 TLinearConstraint *lc)
00289 {
00290 if ((lc->n==0)||(IntervalSize(&(lc->error))>=INF))
00291 return(NOT_REDUCED_BOX);
00292 else
00293 {
00294 Tinterval out,new_range,range,a,ct,one,zero;
00295 unsigned int j,k,x_j;
00296 unsigned int status;
00297 unsigned int m;
00298 Tinterval *is;
00299 boolean haveInf;
00300
00301 m=GetBoxNIntervals(b);
00302 is=GetBoxIntervals(b);
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314 status=NOT_REDUCED_BOX;
00315
00316 for(j=0;((j<lc->n)&&(status!=EMPTY_BOX));j++)
00317 {
00318 x_j=lc->ind[j];
00319
00320
00321 if (IntervalSize(&(is[x_j]))>=epsilon)
00322 {
00323 CopyInterval(&out,&(lc->error));
00324
00325
00326
00327
00328 NewInterval(-ZERO,ZERO,&zero);
00329 IntervalAdd(&out,&zero,&out);
00330
00331 haveInf=FALSE;
00332
00333
00334 for(k=0;((!haveInf)&&(k<lc->n));k++)
00335 {
00336 if (k!=j)
00337 {
00338 if (IntervalSize(&(is[lc->ind[k]]))<INF)
00339 {
00340
00341
00342
00343 NewInterval(-lc->val[k]-ZERO,-lc->val[k]+ZERO,&ct);
00344 IntervalProduct(&(is[lc->ind[k]]),&ct,&a);
00345 IntervalAdd(&out,&a,&out);
00346 }
00347 else
00348 haveInf=TRUE;
00349 }
00350 }
00351
00352 if (!haveInf)
00353 {
00354 if (fabs(lc->val[j])<=10*epsilon)
00355 {
00356
00357
00358
00359 if ((!IsInside(0.0,&out))&&
00360 (fabs(LowerLimit(&out))>10*epsilon)&&
00361 (fabs(UpperLimit(&out))>10*epsilon))
00362 status=EMPTY_BOX;
00363 }
00364 else
00365 {
00366 NewInterval(1,1,&one);
00367
00368
00369
00370 NewInterval(lc->val[j]-ZERO,lc->val[j]+ZERO,&ct);
00371 IntervalDivision(&one,&ct,&ct);
00372 IntervalProduct(&ct,&out,&range);
00373
00374 if (Intersection(&(is[x_j]),&range,&new_range))
00375 {
00376 double s1,s2;
00377
00378
00379 s1=IntervalSize(&(is[x_j]));
00380 s2=IntervalSize(&new_range);
00381 if (s2<(s1-ZERO))
00382 {
00383 CopyInterval(&(is[x_j]),&new_range);
00384 status=REDUCED_BOX;
00385 }
00386 }
00387 else
00388 status=EMPTY_BOX;
00389 }
00390 }
00391 }
00392 }
00393
00394 return(status);
00395 }
00396 }
00397
00398
00399
00400
00401
00402
00403 boolean CmpLinearConstraints(double *scaleOne2Two,TLinearConstraint *lc1,TLinearConstraint *lc2)
00404 {
00405 boolean equal;
00406
00407 equal=FALSE;
00408
00409 if (lc1->n==lc2->n)
00410 {
00411 unsigned int j,l;
00412 double n1,n2,n12,c;
00413
00414 n12=0;
00415 equal=TRUE;
00416 j=0;
00417 while((equal)&&(j<lc1->n))
00418 {
00419
00420 equal=FALSE;
00421 l=0;
00422 while((!equal)&&(l<lc2->n))
00423 {
00424 equal=(lc1->ind[j]==lc2->ind[l]);
00425 if (equal)
00426 n12+=(lc1->val[j]*lc2->val[l]);
00427 else
00428 l++;
00429 }
00430 j++;
00431 }
00432
00433 if(equal)
00434 {
00435 n1=0.0;
00436 n2=0.0;
00437 for(j=0;j<lc1->n;j++)
00438 {
00439 n1+=(lc1->val[j]*lc1->val[j]);
00440 n2+=(lc2->val[j]*lc2->val[j]);
00441 }
00442 n1=sqrt(n1);
00443 n2=sqrt(n2);
00444
00445
00446
00447
00448 c=n12/(n1*n2);
00449 equal=((fabs(c-1.0)<=ZERO)||(fabs(c+1.0)<=ZERO));
00450
00451 if (equal)
00452 {
00453
00454
00455
00456
00457 *scaleOne2Two=(c<0?-1.0:1.0)*n2/n1;
00458
00459
00460 }
00461 }
00462 }
00463 return(equal);
00464 }
00465
00466 double EvaluateLinearConstraint(double *varValues,TLinearConstraint *lc)
00467 {
00468 double v;
00469 unsigned int kk;
00470
00471 v=0;
00472 for(kk=0;kk<lc->n;kk++)
00473 v+=(lc->val[kk]*varValues[lc->ind[kk]]);
00474
00475 return(v);
00476 }
00477
00478 void EvaluateLinearConstraintInt(Tinterval *varValues,Tinterval *i_out,TLinearConstraint *lc)
00479 {
00480 unsigned int kk;
00481 Tinterval a,r;
00482
00483 NewInterval(0,0,i_out);
00484 for(kk=0;kk<lc->n;kk++)
00485 {
00486
00487
00488
00489 NewInterval(lc->val[kk]-ZERO,lc->val[kk]+ZERO,&a);
00490 IntervalProduct(&(varValues[lc->ind[kk]]),&a,&r);
00491 IntervalAdd(&r,i_out,i_out);
00492 }
00493 }
00494
00495 void PrintLinearConstraint(FILE *f,boolean eq,char **varName,TLinearConstraint *lc)
00496 {
00497 unsigned int kk;
00498 Tinterval e;
00499
00500 for(kk=0;kk<lc->n;kk++)
00501 {
00502 if (lc->val[kk]==1)
00503 {
00504 if (kk>0)
00505 fprintf(f,"+");
00506 }
00507 else
00508 {
00509 if (lc->val[kk]==-1)
00510 fprintf(f,"-");
00511 else
00512 fprintf(f,"%+.12g*",lc->val[kk]);
00513 }
00514
00515 if (varName==NULL)
00516 fprintf(f,"v_%u",lc->ind[kk]);
00517 else
00518 fprintf(f,"%s",varName[lc->ind[kk]]);
00519
00520 }
00521 if (eq)
00522 {
00523 fprintf(f," = ");
00524 PrintInterval(f,&(lc->error));
00525 fprintf(f," (s:%g)",IntervalSize(&(lc->error)));
00526 }
00527 else
00528 {
00529 if (IntervalSize(&(lc->error))<ZERO)
00530 {
00531 double c;
00532
00533 c=IntervalCenter(&(lc->error));
00534 if ((c!=0.0)||(lc->n==0))
00535 fprintf(f," %+.12g",-c);
00536 }
00537 else
00538 {
00539 IntervalInvert(&(lc->error),&e);
00540 fprintf(f," + ");
00541 PrintInterval(f,&e);
00542 }
00543 }
00544 fprintf(f,"\n");
00545 }
00546
00547 void SaveLinearConstraint(FILE *f,TLinearConstraint *lc)
00548 {
00549 unsigned int i;
00550
00551 fprintf(f,"%u %u ",lc->n,lc->max);
00552
00553 for(i=0;i<lc->n;i++)
00554 fprintf(f,"%u %f ",lc->ind[i],lc->val[i]);
00555
00556 fprintf(f," %f %f ",LowerLimit(&(lc->error)),UpperLimit(&(lc->error)));
00557 }
00558
00559 void LoadLinearConstraint(FILE *f,TLinearConstraint *lc)
00560 {
00561 unsigned int i;
00562 double l,u;
00563
00564 fscanf(f,"%u %u",&(lc->n),&(lc->max));
00565 NEW(lc->ind,lc->max,unsigned int);
00566 NEW(lc->val,lc->max,double);
00567
00568 for(i=0;i<lc->n;i++)
00569 {
00570 fscanf(f,"%u",&(lc->ind[i]));
00571 fscanf(f,"%lf",&(lc->val[i]));
00572 }
00573 fscanf(f,"%lf",&l);
00574 fscanf(f,"%lf",&u);
00575 NewInterval(l,u,&(lc->error));
00576 }
00577
00578
00579 void DeleteLinearConstraint(TLinearConstraint *lc)
00580 {
00581 free(lc->ind);
00582 free(lc->val);
00583 DeleteInterval(&(lc->error));
00584 }
00585
00586