COCONUT API: cint_evaluator.h Source File

00001 // Complex Interval Evaluator implementation -*- C++ -*-
00002 
00003 // $Id: cint_evaluator.h 347 2008-02-06 15:20:01Z schilly $
00004 // Copyright (C) 2001-2003 Hermann Schichl, Knut Petras
00005 //
00006 // This file is part of the COCONUT API.  This library
00007 // is free software; you can redistribute it and/or modify it under the
00008 // terms of the Library GNU General Public License as published by the
00009 // Free Software Foundation; either version 2, or (at your option)
00010 // any later version.
00011 
00012 // This library is distributed in the hope that it will be useful,
00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 // Library GNU General Public License for more details.
00016 
00017 // As a special exception, you may use this file as part of a free software
00018 // library without restriction.  Specifically, if other files instantiate
00019 // templates or use macros or inline functions from this file, or you compile
00020 // this file and link it with other files to produce an executable, this
00021 // file does not by itself cause the resulting executable to be covered by
00022 // the Library GNU General Public License.  This exception does not however
00023 // invalidate any other reasons why the executable file might be covered by
00024 // the Library GNU General Public License.
00025 
00028 #ifndef _CINT_EVALUATOR_H_
00029 #define _CINT_EVALUATOR_H_
00030 
00031 #include <coconut_config.h>
00032 #include <evaluator.h>
00033 #include <expression.h>
00034 #include <model.h>
00035 #include <eval_main.h>
00036 #include <linalg.h>
00037 #include <math.h>
00038 #include <cinterval.h>
00039 #include <api_exception.h>
00040 
00041 using namespace vgtl;
00042 
00043 namespace coco {
00044 
00045 typedef cinterval (*cinterval_evaluator)(const std::vector<cinterval>* __x,
00046                                  const variable_indicator& __v);
00047 
00048 struct cinterval_eval_type
00049 {
00050   const std::vector<cinterval>* x;
00051   std::vector<cinterval>* cache;
00052   const model* mod;
00053   void *p;
00054   cinterval d;
00055   int info;
00056   cinterval r;
00057   unsigned int n;
00058 };
00059 
00060 class cinterval_eval : public
00061   cached_forward_evaluator_base<cinterval_eval_type,expression_node,cinterval,
00062                                 expression_const_walker>
00063 {
00064 private:
00065   typedef cached_forward_evaluator_base<cinterval_eval_type,expression_node,
00066                                       cinterval,expression_const_walker> _Base;
00067 
00068 protected:
00069   bool is_cached(const node_data_type& __data)
00070   {
00071     if(__data.operator_type == EXPRINFO_LIN ||
00072        __data.operator_type == EXPRINFO_QUAD)
00073       return true;
00074     else
00075       return false;
00076 #if 0
00077     if(eval_data.cache && __data.n_parents > 1 && __data.n_children > 0)
00078       return true;
00079     else
00080       return false;
00081 #endif
00082   }
00083 
00084 private:
00085   /*  cinterval __power(double __coeff, cinterval __x, int __exp)
00086   {
00087     if(__exp == 0)
00088       return 1.;
00089     else
00090     {
00091       cinterval k = __coeff*__x;
00092       switch(__exp)
00093       {
00094         case 1:
00095           return k;
00096           break;
00097         case 2:
00098           return sqr(k);
00099           break;
00100         case -1:
00101           return 1./k;
00102           break;
00103         case -2:
00104           return 1./sqr(k);
00105           break;
00106         default:
00107           return ipow(k, __exp);
00108           break;
00109       }
00110     }
00111     return 0;           // we never come here!
00112   }
00113   */
00114 public:
00115   cinterval_eval(const std::vector<cinterval>& __x, const variable_indicator& __v,
00116                 const model& __m, std::vector<cinterval>* __c) : _Base()
00117   {
00118     eval_data.x = &__x;
00119     eval_data.cache = __c;
00120     eval_data.mod = &__m;
00121     eval_data.d = 0.;
00122     eval_data.p = NULL;
00123     v_ind = &__v;
00124     eval_data.n = 0;
00125   }
00126 
00127   cinterval_eval(const cinterval_eval& __v) : _Base(__v) {}
00128   
00129   ~cinterval_eval() {}
00130 
00131   expression_const_walker short_cut_to(const expression_node& __data)
00132   { return eval_data.mod->node(0); }    // return anything
00133 
00134   void initialize() { return; }
00135   
00136   int initialize(const expression_node& __data)
00137   {
00138     if(__data.ev && (*__data.ev)[FUNC_CRANGE])
00139     // is there a short-cut evaluator defined?
00140     {
00141       eval_data.r =
00142         (*(cinterval_evaluator)(*__data.ev)[FUNC_CRANGE])(eval_data.x,
00143                                                                  *v_ind);
00144       return 0;     // don't perform the remaining graph walk
00145     }
00146     else
00147     {
00148       switch(__data.operator_type)
00149       {
00150         case EXPRINFO_SUM:
00151         case EXPRINFO_PROD:
00152         case EXPRINFO_MAX:
00153         case EXPRINFO_MIN:
00154         case EXPRINFO_INVERT:
00155         case EXPRINFO_DIV:
00156           eval_data.r = __data.params.nd();
00157           break;
00158         case EXPRINFO_MONOME:
00159         case EXPRINFO_IN:
00160         case EXPRINFO_AND:
00161         case EXPRINFO_NOGOOD:
00162           eval_data.r = 1.;
00163           break;
00164         case EXPRINFO_ALLDIFF:
00165           eval_data.p = (void*) new std::vector<cinterval>;
00166           ((std::vector<cinterval>*)eval_data.p)->reserve(__data.n_children);
00167           // no break!
00168         case EXPRINFO_MEAN:
00169         case EXPRINFO_IF:
00170         case EXPRINFO_OR:
00171         case EXPRINFO_NOT:
00172         case EXPRINFO_COUNT:
00173         case EXPRINFO_SCPROD:
00174         case EXPRINFO_NORM:
00175         case EXPRINFO_LEVEL:
00176           eval_data.r = 0.;
00177           break;
00178         case EXPRINFO_DET:
00179         case EXPRINFO_PSD:
00180           // here construct square cinterval matrix
00181           break;
00182         case EXPRINFO_COND:
00183         case EXPRINFO_FEM:
00184         case EXPRINFO_MPROD:
00185           // here construct rectangular cinterval matrix
00186           break;
00187       }
00188       return 1;
00189     }
00190   }
00191 
00192   void calculate(const expression_node& __data)
00193   {
00194     if(__data.operator_type > 0)
00195     {
00196       eval_data.r = __data.i_evaluate(-1, __data.params.nn(), *eval_data.x,
00197                                       *v_ind, eval_data.r, cintval(0.0), NULL);
00198     }
00199     eval_data.r.intersectwith(__data.f_bounds);
00200   }
00201 
00202   void retrieve_from_cache(const expression_node& __data)
00203   {
00204     if(__data.operator_type == EXPRINFO_LIN)
00205     {
00206       const linalg::matrix<double>::Row& lrw(eval_data.mod->lin[__data.params.nn()]);
00207       linalg::matrix<double>::Row::const_iterator _x, _e;
00208 
00209       eval_data.r = cinterval(0.);
00210       for(_x = lrw.begin(); _x != _e; ++_x)
00211         eval_data.r += (*eval_data.x)[_x.index()] * *_x;
00212     }
00213     else if(__data.operator_type == EXPRINFO_QUAD)
00214     {
00215       throw nyi_exception("CI evaluator: QUAD");
00216     }
00217     else
00218       eval_data.r = (*eval_data.cache)[__data.node_num];
00219   }
00220 
00221   int update(const cinterval& __rval)
00222   {
00223     eval_data.r = __rval;
00224     return 0;
00225   }
00226 
00227   int update(const expression_node& __data, const cinterval& __rval)
00228   {
00229     cinterval __x;
00230     int ret = 0;
00231     if(__data.operator_type < 0)
00232     {
00233       switch(__data.operator_type)
00234       {
00235         case EXPRINFO_CONSTANT:
00236           eval_data.r = __data.params.nd();
00237           break;
00238         case EXPRINFO_VARIABLE:
00239           eval_data.r = (*eval_data.x)[__data.params.nn()];
00240           break;
00241         case EXPRINFO_SUM:
00242         case EXPRINFO_MEAN:
00243           eval_data.r += __data.coeffs[eval_data.n++]*__rval;
00244           break;
00245         case EXPRINFO_PROD:
00246           // the coefficients MUST be collected in __data.params.nd()
00247           eval_data.r *= __rval;
00248           if(eval_data.r == 0.) ret = -1;
00249           break;
00250         case EXPRINFO_MONOME:
00251           // the coefficients MUST be collected in __data.params.nd()
00252           eval_data.r *= __power(__data.coeffs[eval_data.n], __rval,
00253                                  __data.params.n()[eval_data.n]);
00254           if(eval_data.r == 0.) ret = -1;
00255           eval_data.n++;
00256           break;
00257         case EXPRINFO_MAX:
00258           /* CHANGE HERE */
00259           { cinterval h = __rval * __data.coeffs[eval_data.n++];
00260             eval_data.r = imax(eval_data.r, h);
00261           }
00262           break;
00263         case EXPRINFO_MIN:
00264           /* CHANGE HERE */
00265           { cinterval h = __rval * __data.coeffs[eval_data.n++];
00266             eval_data.r = imin(eval_data.r, h);
00267           }
00268           break;
00269         case EXPRINFO_SCPROD:
00270           { cinterval h = __data.coeffs[eval_data.n]*__rval;
00271             // if the number of children is odd, the list is implicitly
00272             // lengthened by one null element (i.e. implicitly shortened)
00273             if(eval_data.n & 1) // y_n
00274               eval_data.r += eval_data.d*h;
00275             else // x_n
00276               eval_data.d = h;
00277           }
00278           eval_data.n++;
00279           break;
00280         case EXPRINFO_NORM:
00281           if(__data.params.nd() != 2.)
00282           {
00283             throw nyi_exception("CI evaluator: non-euclidean NORM");
00284           }
00285           else
00286           {
00287             cinterval h = __data.coeffs[eval_data.n++]*__rval;
00288             eval_data.r += sqr(h);
00289           }
00290           if(eval_data.n == __data.n_children && __data.params.nd() == 2.)
00291             eval_data.r = sqrt(eval_data.r);
00292           break;
00293         case EXPRINFO_INVERT:
00294           eval_data.r /= __rval;
00295           break;
00296         case EXPRINFO_DIV:
00297           // this evaluator requires, that the second coefficient
00298           // is put into the first.
00299           if(eval_data.n == 0)
00300             eval_data.r *= __rval;
00301           else
00302             eval_data.r /= __rval;
00303           ++eval_data.n;
00304           break;
00305         case EXPRINFO_SQUARE:
00306           { cinterval h = __data.coeffs[0]*__rval+__data.params.nd();
00307             eval_data.r = sqr(h);
00308           }
00309           break;
00310         case EXPRINFO_INTPOWER:
00311           eval_data.r = __power(__data.coeffs[0], __rval, __data.params.nn());
00312           break;
00313         case EXPRINFO_SQROOT:
00314           eval_data.r = sqrt(__data.coeffs[0]*__rval+__data.params.nd());
00315           break;
00316         case EXPRINFO_ABS:
00317           eval_data.r = abs(__data.coeffs[0]*__rval+__data.params.nd());
00318           break;
00319         case EXPRINFO_POW:
00320           { cinterval h = __rval * __data.coeffs[eval_data.n];
00321             if(eval_data.n == 0)
00322               eval_data.r = h+__data.params.nd();
00323             else
00324               eval_data.r = pow(eval_data.r, h);
00325           }
00326           ++eval_data.n;
00327           break;
00328         case EXPRINFO_EXP:
00329           eval_data.r = exp(__rval*__data.coeffs[0]+__data.params.nd());
00330           break;
00331         case EXPRINFO_LOG:
00332           eval_data.r = log(__rval*__data.coeffs[0]+__data.params.nd());
00333           break;
00334         case EXPRINFO_SIN:
00335           eval_data.r = sin(__rval*__data.coeffs[0]+__data.params.nd());
00336           break;
00337         case EXPRINFO_COS:
00338           eval_data.r = cos(__rval*__data.coeffs[0]+__data.params.nd());
00339           break;
00340         case EXPRINFO_ATAN2:
00341           { cinterval h = __rval * __data.coeffs[eval_data.n];
00342             if(eval_data.n == 0)
00343               eval_data.r = h;
00344             else
00345               eval_data.r = atan2(eval_data.r, h);
00346             ++eval_data.n;
00347           }
00348           break;
00349         case EXPRINFO_GAUSS:
00350           { cinterval h = (__data.coeffs[0]*__rval-__data.params.d()[0])/
00351                                 __data.params.d()[1];
00352             eval_data.r = gauss(h);
00353           }
00354           break;
00355         case EXPRINFO_POLY:
00356           throw nyi_exception("CI evaluator: POLY");
00357           break;
00358         case EXPRINFO_LIN:
00359         case EXPRINFO_QUAD:
00360           // we will never be here
00361           break;
00362         case EXPRINFO_IN:
00363           {
00364             __x = __data.coeffs[eval_data.n]*__rval;
00365             const cinterval& i(__data.params.i()[eval_data.n]);
00366             if(i.disjoint(__x))     // it is not in
00367             {
00368               eval_data.r = -1.;
00369               ret = -1;
00370             }
00371             else if(eval_data.r != cinterval(-1.))
00372             {
00373               if(!i.superset(__x))  // in and out are possible
00374               {
00375                 if(i.sup() == __x.inf() || i.inf() == __x.sup())
00376                 // intersection is only on the border
00377                   eval_data.r = cinterval(-1.,0.);
00378                 else if(eval_data.r.sup() > 0)
00379                   eval_data.r = cinterval(-1.,1.);
00380               }
00381               else if(i.inf() == __x.inf() || i.sup() == __x.sup())
00382               // they are not disjoint, and i is a superset of __x
00383               // but a border is possible
00384               {
00385                 if(eval_data.r.inf() == 1.)
00386                 {
00387                   if(__x.is_thin())       // only border
00388                     eval_data.r = 0.;
00389                   else                    // interior and border possible
00390                     eval_data.r = cinterval(0.,1.);
00391                 }
00392                 else if(__x.is_thin() && eval_data.r.inf() == 0.)
00393                   eval_data.r = 0.;
00394               }
00395               // if __x is in the interior of i eval_data.r does not change
00396             }
00397           }
00398           eval_data.n++;
00399           break;
00400         case EXPRINFO_IF:
00401           __x = __rval * __data.coeffs[eval_data.n];
00402           if(eval_data.n == 0) // this is the condition
00403           {
00404             const cinterval& i(__data.params.ni());
00405             if(i.superset(__x))      // always true
00406               eval_data.info = 0;
00407             else if(i.disjoint(__x)) // never true
00408             {
00409               eval_data.info = 1;
00410               ret = 1;               // skip if part
00411             }
00412             else                     // both is possible
00413               eval_data.info = 2;
00414           }
00415           else if(eval_data.n == 1)  // this is the if part
00416           {
00417             eval_data.r = __x;
00418             if(eval_data.info == 0)  // it was always true
00419               ret = -1;              // don't continue walking
00420           }
00421           else // this is the else part
00422           {
00423             if(eval_data.info == 1) // only the else part
00424               eval_data.r = __x;
00425             else                      // both parts
00426               eval_data.r.hullwith(__x);
00427           }
00428           eval_data.n += ret+1;
00429           break;
00430         case EXPRINFO_AND:
00431           { __x = __data.coeffs[eval_data.n]*__rval;
00432             const cinterval& i(__data.params.i()[eval_data.n]);
00433             if(eval_data.r != 0.)
00434             {
00435               if(i.disjoint(__x))
00436               {
00437                 eval_data.r = 0;
00438                 ret = -1;         // result cannot become more false
00439               }
00440               else if(!i.superset(__x))
00441                 eval_data.r = cinterval(0.,1.);
00442               // else this does not change the value
00443             }
00444           }
00445           eval_data.n++;
00446           break;
00447         case EXPRINFO_OR:
00448           { __x = __data.coeffs[eval_data.n]*__rval;
00449             const cinterval& i(__data.params.i()[eval_data.n]);
00450             if(eval_data.r != 1.)
00451             {
00452               if(i.superset(__x))
00453               {
00454                 eval_data.r = 1;
00455                 ret = -1;         // result cannot become more true
00456               }
00457               else if(!i.disjoint(__x))
00458                 eval_data.r = cinterval(0.,1.);
00459             }
00460           }
00461           eval_data.n++;
00462           break;
00463         case EXPRINFO_NOT:
00464           { __x = __data.coeffs[0]*__rval;
00465             const cinterval& i(__data.params.ni());
00466             if(i.superset(__x))
00467               eval_data.r = 0;
00468             else if(i.disjoint(__x))
00469               eval_data.r = 1;
00470             else
00471               eval_data.r = cinterval(0.,1.);
00472           }
00473           break;
00474         case EXPRINFO_IMPLIES:
00475           { const cinterval& i(__data.params.i()[eval_data.n]);
00476             __x = __rval * __data.coeffs[eval_data.n];
00477             if(eval_data.n == 0)
00478             {
00479               if(i.disjoint(__x))
00480               {
00481                 eval_data.r = 1.;
00482                 ret = -1;        // ex falso quodlibet
00483               }
00484               else if(!i.superset(__x))
00485                 eval_data.r = cinterval(0.,1.);
00486               else
00487                 eval_data.r = cinterval(0.);
00488             }
00489             else // we only get here if first clause might be true
00490             {
00491               if(i.superset(__x))
00492                 eval_data.r = 1.;
00493               else if(!i.disjoint(__x))
00494                 eval_data.r = cinterval(0.,1.);
00495               // else the value does not change
00496             }
00497             ++eval_data.n;
00498           }
00499           break;
00500         case EXPRINFO_COUNT:
00501           { __x = __data.coeffs[eval_data.n]*__rval;
00502             const cinterval& i(__data.params.i()[eval_data.n]);
00503             if(i.superset(__x))
00504               eval_data.r += 1;
00505             else if(!i.disjoint(__x))
00506               eval_data.r += cinterval(0.,1.);
00507           }
00508           eval_data.n++;
00509           break;
00510         case EXPRINFO_ALLDIFF:
00511           { __x = __data.coeffs[eval_data.n]*__rval;
00512             for(std::vector<cinterval>::const_iterator _b =
00513                         ((std::vector<cinterval>*)eval_data.p)->begin();
00514                 _b != ((std::vector<cinterval>*)eval_data.p)->end(); ++_b)
00515             {
00516               cinterval __h = abs(__x-*_b);
00517               if(__h.sup() <= __data.params.nd())  // they are always equal
00518               {
00519                 eval_data.r = 0;
00520                 ret = -1;               // don't continue
00521                 break;
00522               }
00523               else if(__h.inf() <= __data.params.nd()) // they are not always different
00524                 eval_data.r = cinterval(0.,1.);
00525             }
00526             if(ret != -1)
00527               ((std::vector<cinterval>*) eval_data.p)->push_back(__x);
00528           }
00529           eval_data.n++;
00530           if(eval_data.n == __data.n_children || ret == -1)
00531             delete (std::vector<cinterval>*) eval_data.p;
00532           break;
00533         case EXPRINFO_HISTOGRAM:
00534           throw nyi_exception("CI evaluator: HISTOGRAM");
00535           break;
00536         case EXPRINFO_LEVEL:
00537           { int lo = (int)eval_data.r.inf();
00538             int hi = (int)eval_data.r.sup();
00539             __x = __data.coeffs[eval_data.n]*__rval;
00540             cinterval _h;
00541 
00542             // first adjust upper bound
00543             if(hi != INT_MAX)
00544             {
00545               while(hi < __data.params.im().nrows())
00546               { 
00547                 _h = __data.params.im()[hi][eval_data.n];
00548                 if(_h.superset(__x))
00549                   break;
00550                 hi++;
00551               }
00552               if(hi == __data.params.im().nrows())
00553                 hi = INT_MAX;
00554             }
00555             // then adjust lower bound
00556             if(lo != INT_MAX)
00557             {
00558               while(lo < __data.params.im().nrows())
00559               { 
00560                 _h = __data.params.im()[lo][eval_data.n];
00561                 if(!_h.disjoint(__x))
00562                   break;
00563                 lo++;
00564               }
00565               if(lo == __data.params.im().nrows())
00566                 lo = INT_MAX;
00567             }
00568             if(hi == INT_MAX && lo == INT_MAX)
00569             {
00570               eval_data.r = INT_MAX;
00571               ret = -1;
00572             }
00573             else
00574               eval_data.r = cinterval(lo+0.,hi+0.);
00575           }
00576           eval_data.n++;
00577           break;
00578         case EXPRINFO_NEIGHBOR:
00579           // this is suboptimal. It does not take into account
00580           // all the possibilities to proove that the result actually
00581           // EQUALS 1
00582           if(eval_data.n == 0)
00583             eval_data.r = __data.coeffs[0]*__rval;
00584           else
00585           {
00586             cinterval h = eval_data.r;
00587             eval_data.r = 0;
00588             __x = __data.coeffs[1]*__rval;
00589             for(unsigned int i = 0; i < __data.params.n().size(); i+=2)
00590             {
00591               if(h == __data.params.n()[i]+0. &&
00592                  __x == __data.params.n()[i+1]+0.)
00593               {
00594                 eval_data.r = 1;
00595                 break;
00596               }
00597               else if(h.contains(__data.params.n()[i]+0.) &&
00598                       __x.contains(__data.params.n()[i+1]+0.))
00599               {
00600                 eval_data.r = cinterval(0.,1.);
00601                 break;
00602               }
00603             }
00604           }
00605           eval_data.n++;
00606           break;
00607         case EXPRINFO_NOGOOD:
00608           __x = __data.coeffs[eval_data.n]*__rval;
00609           if(eval_data.r != 0.)
00610           { 
00611             if(!__x.contains(__data.params.n()[eval_data.n]+0.))
00612             {
00613               eval_data.r = 0.;
00614               ret = -1;
00615             }
00616             else if(__x != __data.params.n()[eval_data.n]+0.)
00617               eval_data.r = cinterval(0.,1.);
00618           }
00619           eval_data.n++;
00620           break;
00621         case EXPRINFO_EXPECTATION:
00622           throw nyi_exception("CI evaluator: EXPECTATION");
00623           break;
00624         case EXPRINFO_INTEGRAL:
00625           throw nyi_exception("CI evaluator: INTEGRAL");
00626           break;
00627         case EXPRINFO_LOOKUP:
00628         case EXPRINFO_PWLIN:
00629         case EXPRINFO_SPLINE:
00630         case EXPRINFO_PWCONSTLC:
00631         case EXPRINFO_PWCONSTRC:
00632           throw nyi_exception("CI evaluator: Table Operations");
00633           break;
00634         case EXPRINFO_DET:
00635         case EXPRINFO_COND:
00636         case EXPRINFO_PSD:
00637         case EXPRINFO_MPROD:
00638         case EXPRINFO_FEM:
00639           throw nyi_exception("CI evaluator: Matrix Operations");
00640           break;
00641         case EXPRINFO_RE:
00642         case EXPRINFO_IM:
00643         case EXPRINFO_ARG:
00644         case EXPRINFO_CPLXCONJ:
00645           throw nyi_exception("CI evaluator: Complex Operations");
00646           break;
00647         case EXPRINFO_CMPROD:
00648         case EXPRINFO_CGFEM:
00649           throw nyi_exception("CI evaluator: Const Matrix Operations");
00650           break;
00651         default:
00652           throw api_exception(apiee_evaluator,
00653               std::string("CI evaluator: unknown function type ")+
00654               convert_to_str(__data.operator_type));
00655           break;
00656       }
00657     }
00658     else if(__data.operator_type > 0)
00659       // update the function arguments
00660       eval_data.r = __data.i_evaluate(eval_data.n++, __data.params.nn(),
00661                         *eval_data.x, *v_ind, eval_data.r, __rval, NULL);
00662     // use the cache only if it is worthwhile
00663     //    if(eval_data.cache && __data.n_parents > 1 && __data.n_children > 0)
00664     if(eval_data.cache)
00665       (*eval_data.cache)[__data.node_num] = eval_data.r;
00666     return ret;
00667   }
00668 
00669   cinterval calculate_value(bool eval_all)
00670   {
00671     return eval_data.r;
00672   }
00673 };
00674 
00675 } // namespace coco
00676 
00677 #endif /* _CINT_EVALUATOR_H_ */