1 /** 2 * 3 * Copyright (c) 2005-2021 by Pierre-Henri WUILLEMIN(_at_LIP6) & Christophe GONZALES(_at_AMU) 4 * info_at_agrum_dot_org 5 * 6 * This library is free software: you can redistribute it and/or modify 7 * it under the terms of the GNU Lesser General Public License as published by 8 * the Free Software Foundation, either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * This library is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU Lesser General Public License for more details. 15 * 16 * You should have received a copy of the GNU Lesser General Public License 17 * along with this library. If not, see <http://www.gnu.org/licenses/>. 18 * 19 */ 20 21 22 #include <agrum/agrum.h> 23 24 #include <agrum/tools/core/utils_string.h> 25 #include <agrum/CN/credalNet.h> 26 27 namespace gum { 28 namespace credal { 29 30 template < typename GUM_SCALAR > CredalNet()31 CredalNet< GUM_SCALAR >::CredalNet() { 32 _initParams_(); 33 34 _src_bn_ = BayesNet< GUM_SCALAR >(); 35 _src_bn_min_ = BayesNet< GUM_SCALAR >(); 36 _src_bn_max_ = BayesNet< GUM_SCALAR >(); 37 38 GUM_CONSTRUCTOR(CredalNet); 39 } 40 41 template < typename GUM_SCALAR > addVariable(const std::string & name,const Size & card)42 NodeId CredalNet< GUM_SCALAR >::addVariable(const std::string& name, const Size& card) { 43 LabelizedVariable var(name, "node " + name, card); 44 45 NodeId a = _src_bn_.add(var); 46 NodeId b = _src_bn_min_.add(var); 47 NodeId c = _src_bn_max_.add(var); 48 49 if (a != b || a != c /*|| b != c*/) 50 GUM_ERROR(OperationNotAllowed, 51 "addVariable : not the same id over all networks : " << a << ", " << b << ", " 52 << c); 53 54 return a; 55 } 56 57 template < typename GUM_SCALAR > addArc(const NodeId & tail,const NodeId & head)58 void CredalNet< GUM_SCALAR >::addArc(const NodeId& tail, const NodeId& head) { 59 _src_bn_.addArc(tail, head); 60 _src_bn_min_.addArc(tail, head); 61 _src_bn_max_.addArc(tail, head); 62 } 63 64 template < typename GUM_SCALAR > setCPTs(const NodeId & id,const std::vector<std::vector<std::vector<GUM_SCALAR>>> & cpt)65 void CredalNet< GUM_SCALAR >::setCPTs( 66 const NodeId& id, 67 const std::vector< std::vector< std::vector< GUM_SCALAR > > >& cpt) { 68 const Potential< GUM_SCALAR >* const potential(&_src_bn_.cpt(id)); 69 70 auto var_dSize = _src_bn_.variable(id).domainSize(); 71 auto entry_size = potential->domainSize() / var_dSize; 72 73 if (cpt.size() != entry_size) 74 GUM_ERROR(SizeError, 75 "setCPTs : entry sizes of cpts does not match for node id : " 76 << id << " : " << cpt.size() << " != " << entry_size); 77 78 for (const auto& cset: cpt) { 79 if (cset.size() == 0) 80 GUM_ERROR(SizeError, 81 "setCPTs : vertices in credal set does not match for node id : " 82 << id << " with 0 vertices"); 83 84 for (const auto& vertex: cset) { 85 if (vertex.size() != var_dSize) 86 GUM_ERROR(SizeError, 87 "setCPTs : variable modalities in cpts does " 88 "not match for node id : " 89 << id << " with vertex " << vertex << " : " << vertex.size() 90 << " != " << var_dSize); 91 92 GUM_SCALAR sum = 0; 93 94 for (const auto& prob: vertex) { 95 sum += prob; 96 } 97 98 if (std::fabs(sum - 1) > 1e-6) 99 GUM_ERROR(CPTError, 100 "setCPTs : a vertex coordinates does not " 101 "sum to one for node id : " 102 << id << " with vertex " << vertex); 103 } 104 } 105 106 _credalNet_src_cpt_.insert(id, cpt); 107 } 108 109 template < typename GUM_SCALAR > setCPT(const NodeId & id,Size & entry,const std::vector<std::vector<GUM_SCALAR>> & cpt)110 void CredalNet< GUM_SCALAR >::setCPT(const NodeId& id, 111 Size& entry, 112 const std::vector< std::vector< GUM_SCALAR > >& cpt) { 113 const Potential< GUM_SCALAR >* const potential(&_src_bn_.cpt(id)); 114 115 auto var_dSize = _src_bn_.variable(id).domainSize(); 116 auto entry_size = potential->domainSize() / var_dSize; 117 118 if (entry >= entry_size) 119 GUM_ERROR(SizeError, 120 "setCPT : entry is greater or equal than entry size " 121 "(entries start at 0 up to entry_size - 1) : " 122 << entry << " >= " << entry_size); 123 124 if (cpt.size() == 0) GUM_ERROR(SizeError, "setCPT : empty credal set for entry : " << entry) 125 126 for (const auto& vertex: cpt) { 127 if (vertex.size() != var_dSize) 128 GUM_ERROR(SizeError, 129 "setCPT : variable modalities in cpts does not " 130 "match for node id : " 131 << id << " with vertex " << vertex << " at entry " << entry << " : " 132 << vertex.size() << " != " << var_dSize); 133 134 GUM_SCALAR sum = 0; 135 136 for (const auto& prob: vertex) { 137 sum += prob; 138 } 139 140 if (std::fabs(sum - 1) > 1e-6) 141 GUM_ERROR(CPTError, 142 "setCPT : a vertex coordinates does not sum to one for node id : " 143 << id << " at entry " << entry << " with vertex " << vertex); 144 } 145 146 // !! auto does NOT use adress (if available) unless explicitly asked !! 147 auto& node_cpt = _credalNet_src_cpt_.getWithDefault( 148 id, 149 std::vector< std::vector< std::vector< GUM_SCALAR > > >(entry_size)); 150 151 if (node_cpt[entry].size() != 0) 152 GUM_ERROR(DuplicateElement, 153 "setCPT : vertices of entry id " << entry 154 << " already set to : " << node_cpt[entry] 155 << ", cannot insert : " << cpt); 156 157 node_cpt[entry] = cpt; 158 159 /// __credalNet_src_cpt.set ( id, node_cpt ); 160 } 161 162 template < typename GUM_SCALAR > setCPT(const NodeId & id,Instantiation ins,const std::vector<std::vector<GUM_SCALAR>> & cpt)163 void CredalNet< GUM_SCALAR >::setCPT(const NodeId& id, 164 Instantiation ins, 165 const std::vector< std::vector< GUM_SCALAR > >& cpt) { 166 const Potential< GUM_SCALAR >* const potential(&_src_bn_.cpt(id)); 167 168 auto var_dSize = _src_bn_.variable(id).domainSize(); 169 auto entry_size = potential->domainSize() / var_dSize; 170 171 // to be sure of entry index reorder ins according to the bayes net 172 // potentials 173 // ( of the credal net ) 174 // it WONT throw an error if the sequences are not equal not because of 175 // order 176 // but content, so we double check (before & after order correction) 177 // beware of slaves & master 178 Instantiation ref(potential); 179 ref.forgetMaster(); 180 181 ins.forgetMaster(); 182 183 const auto& vseq = ref.variablesSequence(); 184 185 if (ins.variablesSequence() != vseq) { 186 ins.reorder(ref); 187 188 if (ins.variablesSequence() != vseq) 189 GUM_ERROR(OperationNotAllowed, 190 "setCPT : instantiation : " 191 << ins << " is not valid for node id " << id 192 << " which accepts instantiations such as (order is not " 193 "important) : " 194 << ref); 195 } 196 197 Idx entry = 0, jump = 1; 198 199 for (Idx i = 0, end = ins.nbrDim(); i < end; i++) { 200 if (_src_bn_.nodeId(ins.variable(i)) == id) continue; 201 202 entry += ins.val(i) * jump; 203 204 jump *= ins.variable(i).domainSize(); 205 } 206 207 if (entry >= entry_size) 208 GUM_ERROR(SizeError, 209 "setCPT : entry is greater or equal than entry size " 210 "(entries start at 0 up to entry_size - 1) : " 211 << entry << " >= " << entry_size); 212 213 if (cpt.size() == 0) GUM_ERROR(SizeError, "setCPT : empty credal set for entry : " << entry) 214 215 for (const auto& vertex: cpt) { 216 if (vertex.size() != var_dSize) 217 GUM_ERROR(SizeError, 218 "setCPT : variable modalities in cpts does not " 219 "match for node id : " 220 << id << " with vertex " << vertex << " at entry " << entry << " : " 221 << vertex.size() << " != " << var_dSize); 222 223 GUM_SCALAR sum = 0; 224 225 for (const auto& prob: vertex) { 226 sum += prob; 227 } 228 229 if (std::fabs(sum - 1) > 1e-6) 230 GUM_ERROR(CPTError, 231 "setCPT : a vertex coordinates does not sum to one for node id : " 232 << id << " at entry " << entry << " with vertex " << vertex); 233 } 234 235 auto& node_cpt = _credalNet_src_cpt_.getWithDefault( 236 id, 237 std::vector< std::vector< std::vector< GUM_SCALAR > > >(entry_size)); 238 239 if (node_cpt[entry].size() != 0) 240 GUM_ERROR(DuplicateElement, 241 "setCPT : vertices of entry : " << ins << " id " << entry 242 << " already set to : " << node_cpt[entry] 243 << ", cannot insert : " << cpt); 244 245 node_cpt[entry] = cpt; 246 247 /// __credalNet_src_cpt.set ( id, node_cpt ); 248 } 249 250 template < typename GUM_SCALAR > fillConstraints(const NodeId & id,const std::vector<GUM_SCALAR> & lower,const std::vector<GUM_SCALAR> & upper)251 void CredalNet< GUM_SCALAR >::fillConstraints(const NodeId& id, 252 const std::vector< GUM_SCALAR >& lower, 253 const std::vector< GUM_SCALAR >& upper) { 254 try { 255 _src_bn_min_.cpt(id).fillWith(lower); 256 _src_bn_max_.cpt(id).fillWith(upper); 257 } catch (const SizeError&) { 258 GUM_ERROR(SizeError, 259 "fillConstraints : sizes does not match in fillWith for node id : " << id); 260 } 261 } 262 263 template < typename GUM_SCALAR > fillConstraint(const NodeId & id,const Idx & entry,const std::vector<GUM_SCALAR> & lower,const std::vector<GUM_SCALAR> & upper)264 void CredalNet< GUM_SCALAR >::fillConstraint(const NodeId& id, 265 const Idx& entry, 266 const std::vector< GUM_SCALAR >& lower, 267 const std::vector< GUM_SCALAR >& upper) { 268 Potential< GUM_SCALAR >* const potential_min( 269 const_cast< Potential< GUM_SCALAR >* const >(&_src_bn_min_.cpt(id))); 270 Potential< GUM_SCALAR >* const potential_max( 271 const_cast< Potential< GUM_SCALAR >* const >(&_src_bn_max_.cpt(id))); 272 273 auto var_dSize = _src_bn_.variable(id).domainSize(); 274 275 if (lower.size() != var_dSize || upper.size() != var_dSize) 276 GUM_ERROR(SizeError, 277 "setCPT : variable modalities in cpts does not match for node id : " 278 << id << " with sizes of constraints : ( " << lower.size() << " || " 279 << upper.size() << " ) != " << var_dSize); 280 281 auto entry_size = potential_min->domainSize() / var_dSize; 282 283 if (entry >= entry_size) 284 GUM_ERROR(SizeError, 285 "setCPT : entry is greater or equal than entry size " 286 "(entries start at 0 up to entry_size - 1) : " 287 << entry << " >= " << entry_size); 288 289 Instantiation min(potential_min); 290 Instantiation max(potential_max); 291 min.setFirst(); 292 max.setFirst(); 293 294 Idx pos = 0; 295 296 while (pos != entry) { 297 ++min; 298 ++max; 299 ++pos; 300 } 301 302 for (Size i = 0; i < var_dSize; i++) { 303 potential_min->set(min, lower[i]); 304 potential_max->set(max, upper[i]); 305 ++min; 306 ++max; 307 } 308 } 309 310 template < typename GUM_SCALAR > fillConstraint(const NodeId & id,Instantiation ins,const std::vector<GUM_SCALAR> & lower,const std::vector<GUM_SCALAR> & upper)311 void CredalNet< GUM_SCALAR >::fillConstraint(const NodeId& id, 312 Instantiation ins, 313 const std::vector< GUM_SCALAR >& lower, 314 const std::vector< GUM_SCALAR >& upper) { 315 const Potential< GUM_SCALAR >* const potential(&_src_bn_.cpt(id)); 316 /* 317 auto var_dSize = _src_bn_.variable ( id ).domainSize(); 318 auto entry_size = potential->domainSize() / var_dSize; 319 */ 320 // to be sure of entry index reorder ins according to the bayes net 321 // potentials 322 // ( of the credal net ) 323 // it WONT throw an error if the sequences are not equal not because of 324 // order 325 // but content, so we double check (before & after order correction) 326 // beware of slaves & master 327 Instantiation ref(potential); 328 ref.forgetMaster(); 329 330 ins.forgetMaster(); 331 332 const auto& vseq = ref.variablesSequence(); 333 334 if (ins.variablesSequence() != vseq) { 335 ins.reorder(ref); 336 337 if (ins.variablesSequence() != vseq) 338 GUM_ERROR(OperationNotAllowed, 339 "setCPT : instantiation : " 340 << ins << " is not valid for node id " << id 341 << " which accepts instantiations such as (order is not " 342 "important) : " 343 << ref); 344 } 345 346 Idx entry = 0, jump = 1; 347 348 for (Idx i = 0, end = ins.nbrDim(); i < end; i++) { 349 if (_src_bn_.nodeId(ins.variable(i)) == id) continue; 350 351 entry += ins.val(i) * jump; 352 353 jump *= ins.variable(i).domainSize(); 354 } 355 356 /* 357 if ( entry >= entry_size ) 358 GUM_ERROR ( SizeError, "setCPT : entry is greater or equal than entry 359 size 360 (entries start at 0 up to entry_size - 1) : " << entry << " >= " << 361 entry_size 362 ); 363 364 if ( lower.size() != var_dSize || upper.size() != var_dSize ) 365 GUM_ERROR ( SizeError, "setCPT : variable modalities in cpts does not 366 match 367 for node id : " << id << " with sizes of constraints : ( "<< lower.size() 368 << " 369 || " << upper.size() << " ) != " << var_dSize ); 370 */ 371 fillConstraint(id, entry, lower, upper); 372 } 373 374 //////////////////////////////////////////////// 375 /// bnet accessors / shortcuts 376 377 template < typename GUM_SCALAR > instantiation(const NodeId & id)378 INLINE Instantiation CredalNet< GUM_SCALAR >::instantiation(const NodeId& id) { 379 return Instantiation(_src_bn_.cpt(id)); 380 } 381 382 template < typename GUM_SCALAR > domainSize(const NodeId & id)383 INLINE Size CredalNet< GUM_SCALAR >::domainSize(const NodeId& id) { 384 return _src_bn_.variable(id).domainSize(); 385 } 386 387 /////////////////////////////////////////////// 388 389 template < typename GUM_SCALAR > CredalNet(const std::string & src_min_num,const std::string & src_max_den)390 CredalNet< GUM_SCALAR >::CredalNet(const std::string& src_min_num, 391 const std::string& src_max_den) { 392 _initParams_(); 393 _initCNNets_(src_min_num, src_max_den); 394 395 GUM_CONSTRUCTOR(CredalNet); 396 } 397 398 template < typename GUM_SCALAR > CredalNet(const BayesNet<GUM_SCALAR> & src_min_num,const BayesNet<GUM_SCALAR> & src_max_den)399 CredalNet< GUM_SCALAR >::CredalNet(const BayesNet< GUM_SCALAR >& src_min_num, 400 const BayesNet< GUM_SCALAR >& src_max_den) { 401 _initParams_(); 402 _initCNNets_(src_min_num, src_max_den); 403 404 GUM_CONSTRUCTOR(CredalNet); 405 } 406 407 template < typename GUM_SCALAR > ~CredalNet()408 CredalNet< GUM_SCALAR >::~CredalNet() { 409 if (_current_bn_ != nullptr) delete _current_bn_; 410 411 if (_credalNet_current_cpt_ != nullptr) delete _credalNet_current_cpt_; 412 413 if (_current_nodeType_ != nullptr) delete _current_nodeType_; 414 415 GUM_DESTRUCTOR(CredalNet); 416 } 417 418 // from BNs with numerators & denominators or cpts & denominators to credal 419 template < typename GUM_SCALAR > bnToCredal(const GUM_SCALAR beta,const bool oneNet,const bool keepZeroes)420 void CredalNet< GUM_SCALAR >::bnToCredal(const GUM_SCALAR beta, 421 const bool oneNet, 422 const bool keepZeroes) { 423 GUM_SCALAR epsi_min = 1.; 424 GUM_SCALAR epsi_max = 0.; 425 GUM_SCALAR epsi_moy = 0.; 426 GUM_SCALAR epsi_den = 0.; 427 428 for (auto node: src_bn().nodes()) { 429 const Potential< GUM_SCALAR >* const potential(&_src_bn_.cpt(node)); 430 431 Potential< GUM_SCALAR >* const potential_min( 432 const_cast< Potential< GUM_SCALAR >* const >(&_src_bn_min_.cpt(node))); 433 Potential< GUM_SCALAR >* const potential_max( 434 const_cast< Potential< GUM_SCALAR >* const >(&_src_bn_max_.cpt(node))); 435 436 Size var_dSize = _src_bn_.variable(node).domainSize(); 437 Size entry_size = potential->domainSize() / var_dSize; 438 439 Instantiation ins(potential); 440 Instantiation ins_min(potential_min); 441 Instantiation ins_max(potential_max); 442 443 ins.setFirst(); 444 ins_min.setFirst(); 445 ins_max.setFirst(); 446 447 std::vector< GUM_SCALAR > vertex(var_dSize); 448 449 for (Size entry = 0; entry < entry_size; entry++) { 450 GUM_SCALAR den; 451 452 if (oneNet) 453 den = 0; 454 else 455 den = potential_max->get(ins_max); 456 457 Size nbm = 0; 458 459 for (Size modality = 0; modality < var_dSize; modality++) { 460 vertex[modality] = potential->get(ins); 461 462 if (oneNet) { 463 den += vertex[modality]; 464 465 if (vertex[modality] < 1 && vertex[modality] > 0) 466 GUM_ERROR(OperationNotAllowed, 467 "bnToCredal : the BayesNet contains " 468 "probabilities and not event counts " 469 "although user precised oneNet = " 470 << oneNet); 471 } 472 473 if (vertex[modality] > 0) nbm++; 474 475 ++ins; 476 } 477 478 /// check sum is 1 if not oneNet (we are not using counts) 479 if (!oneNet) { 480 GUM_SCALAR sum = 0; 481 482 for (auto modality = vertex.cbegin(), theEnd = vertex.cend(); modality != theEnd; 483 ++modality) { 484 sum += *modality; 485 } 486 487 if (std::fabs(1. - sum) > _epsRedund_) { 488 GUM_ERROR(CPTError, 489 _src_bn_.variable(node).name() 490 << "(" << _epsRedund_ << ") does not sum to one for" 491 << " " << entry << std::endl 492 << vertex << std::endl 493 << ins << std::endl); 494 } 495 } 496 497 /// end check sum is 1 498 499 GUM_SCALAR epsilon; 500 501 if (beta == 0) 502 epsilon = 0; 503 else if (den == 0 || beta == 1) 504 epsilon = GUM_SCALAR(1.0); 505 else 506 epsilon = GUM_SCALAR(std::pow(beta, std::log1p(den))); 507 508 epsi_moy += epsilon; 509 epsi_den += 1; 510 511 if (epsilon > epsi_max) epsi_max = epsilon; 512 513 if (epsilon < epsi_min) epsi_min = epsilon; 514 515 GUM_SCALAR min, max; 516 517 for (Size modality = 0; modality < var_dSize; modality++) { 518 if ((vertex[modality] > 0 && nbm > 1) || !keepZeroes) { 519 min = GUM_SCALAR((1. - epsilon) * vertex[modality]); 520 521 if (oneNet) min = GUM_SCALAR(min * 1.0 / den); 522 523 max = GUM_SCALAR(min + epsilon); 524 } else { // if ( ( vertex[modality] == 0 && keepZeroes ) || ( 525 // vertex[modality] > 0 && nbm <= 1 ) || ( vertex[modality] == 0 526 // && nbm <= 1 ) ) { 527 min = vertex[modality]; 528 529 if (oneNet) min = GUM_SCALAR(min * 1.0 / den); 530 531 max = min; 532 } 533 534 potential_min->set(ins_min, min); 535 potential_max->set(ins_max, max); 536 537 ++ins_min; 538 ++ins_max; 539 } // end of : for each modality 540 541 } // end of : for each entry 542 543 } // end of : for each variable 544 545 _epsilonMin_ = epsi_min; 546 _epsilonMax_ = epsi_max; 547 _epsilonMoy_ = (GUM_SCALAR)epsi_moy / (GUM_SCALAR)epsi_den; 548 549 _intervalToCredal_(); 550 } 551 552 template < typename GUM_SCALAR > lagrangeNormalization()553 void CredalNet< GUM_SCALAR >::lagrangeNormalization() { 554 for (auto node: _src_bn_.nodes()) { 555 const Potential< GUM_SCALAR >* const potential(&_src_bn_.cpt(node)); 556 557 auto var_dSize = _src_bn_.variable(node).domainSize(); 558 auto entry_size = potential->domainSize() / var_dSize; 559 560 Instantiation ins(potential); 561 562 ins.setFirst(); 563 564 std::vector< GUM_SCALAR > vertex(var_dSize); 565 566 for (Size entry = 0; entry < entry_size; entry++) { 567 GUM_SCALAR den = 0; 568 bool zeroes = false; 569 Instantiation ins_prev = ins; 570 571 for (Size modality = 0; modality < var_dSize; modality++) { 572 vertex[modality] = potential->get(ins); 573 574 if (vertex[modality] < 1 && vertex[modality] > 0) 575 GUM_ERROR(OperationNotAllowed, 576 "lagrangeNormalization : the BayesNet " 577 "contains probabilities and not event " 578 "counts."); 579 580 den += vertex[modality]; 581 582 if (!zeroes && vertex[modality] == 0) { zeroes = true; } 583 584 ++ins; 585 } 586 587 if (zeroes) { 588 ins = ins_prev; 589 590 for (Size modality = 0; modality < var_dSize; modality++) { 591 potential->set(ins, potential->get(ins) + 1); 592 ++ins; 593 } 594 } 595 596 } // end of : for each entry 597 598 } // end of : for each variable 599 } 600 601 template < typename GUM_SCALAR > idmLearning(const Idx s,const bool keepZeroes)602 void CredalNet< GUM_SCALAR >::idmLearning(const Idx s, const bool keepZeroes) { 603 for (auto node: _src_bn_.nodes()) { 604 const Potential< GUM_SCALAR >* const potential(&_src_bn_.cpt(node)); 605 606 Potential< GUM_SCALAR >* const potential_min( 607 const_cast< Potential< GUM_SCALAR >* const >(&_src_bn_min_.cpt(node))); 608 Potential< GUM_SCALAR >* const potential_max( 609 const_cast< Potential< GUM_SCALAR >* const >(&_src_bn_max_.cpt(node))); 610 611 Size var_dSize = _src_bn_.variable(node).domainSize(); 612 Size entry_size = potential->domainSize() / var_dSize; 613 614 Instantiation ins(potential); 615 Instantiation ins_min(potential_min); 616 Instantiation ins_max(potential_max); 617 618 ins.setFirst(); 619 ins_min.setFirst(); 620 ins_max.setFirst(); 621 622 std::vector< GUM_SCALAR > vertex(var_dSize); 623 624 for (Size entry = 0; entry < entry_size; entry++) { 625 GUM_SCALAR den = 0; 626 Size nbm = 0; 627 628 for (Size modality = 0; modality < var_dSize; modality++) { 629 vertex[modality] = potential->get(ins); 630 631 if (vertex[modality] < 1 && vertex[modality] > 0) 632 GUM_ERROR(OperationNotAllowed, 633 "idmLearning : the BayesNet contains " 634 "probabilities and not event counts."); 635 636 den += vertex[modality]; 637 638 if (vertex[modality] > 0) nbm++; 639 640 ++ins; 641 } 642 643 if (nbm > 1 || !keepZeroes) den += s; 644 645 GUM_SCALAR min, max; 646 647 for (Size modality = 0; modality < var_dSize; modality++) { 648 min = vertex[modality]; 649 max = min; 650 651 if ((vertex[modality] > 0 && nbm > 1) || !keepZeroes) { max += s; } 652 653 min = GUM_SCALAR(min * 1.0 / den); 654 max = GUM_SCALAR(max * 1.0 / den); 655 656 potential_min->set(ins_min, min); 657 potential_max->set(ins_max, max); 658 659 ++ins_min; 660 ++ins_max; 661 } // end of : for each modality 662 663 } // end of : for each entry 664 665 } // end of : for each variable 666 667 _epsilonMin_ = GUM_SCALAR(s); 668 _epsilonMax_ = GUM_SCALAR(s); 669 _epsilonMoy_ = GUM_SCALAR(s); 670 _intervalToCredal_(); 671 } 672 673 /* no need for lrs : (max ... min ... max) vertices from bnToCredal() */ 674 template < typename GUM_SCALAR > _intervalToCredal_()675 void CredalNet< GUM_SCALAR >::_intervalToCredal_() { 676 if (!_credalNet_src_cpt_.empty()) _credalNet_src_cpt_.clear(); 677 678 _credalNet_src_cpt_.resize(_src_bn_.size()); 679 680 for (auto node: _src_bn_.nodes()) { 681 const Potential< GUM_SCALAR >* const potential_min(&_src_bn_min_.cpt(node)); 682 const Potential< GUM_SCALAR >* const potential_max(&_src_bn_max_.cpt(node)); 683 684 Size var_dSize = _src_bn_.variable(node).domainSize(); 685 Size entry_size = potential_min->domainSize() / var_dSize; 686 687 std::vector< std::vector< std::vector< GUM_SCALAR > > > var_cpt(entry_size); 688 689 Instantiation ins_min(potential_min); 690 Instantiation ins_max(potential_max); 691 692 ins_min.setFirst(); 693 ins_max.setFirst(); 694 695 std::vector< GUM_SCALAR > lower(var_dSize); 696 std::vector< GUM_SCALAR > upper(var_dSize); 697 698 for (Size entry = 0; entry < entry_size; entry++) { 699 for (Size modality = 0; modality < var_dSize; modality++, ++ins_min, ++ins_max) { 700 lower[modality] = potential_min->get(ins_min); 701 upper[modality] = potential_max->get(ins_max); 702 } 703 704 bool all_equals = true; 705 std::vector< std::vector< GUM_SCALAR > > vertices; 706 707 for (Size modality = 0; modality < var_dSize; modality++) { 708 if (std::fabs(upper[modality] - lower[modality]) < 1e-6) continue; 709 710 all_equals = false; 711 std::vector< GUM_SCALAR > vertex(var_dSize); 712 vertex[modality] = upper[modality]; 713 714 for (Size mod = 0; mod < var_dSize; mod++) { 715 if (modality != mod) vertex[mod] = lower[mod]; 716 } 717 718 GUM_SCALAR total = 0; 719 720 auto vsize = vertex.size(); 721 722 for (Size i = 0; i < vsize; i++) 723 total += vertex[i]; 724 725 if (std::fabs(total - 1.) > 1e-6) 726 GUM_ERROR(CPTError, 727 _src_bn_.variable(node).name() 728 << " does not sum to one for " << entry << std::endl 729 << vertex << std::endl); 730 731 vertices.push_back(vertex); 732 } 733 734 if (all_equals) { 735 std::vector< GUM_SCALAR > vertex(var_dSize); 736 737 for (Size modality = 0; modality < var_dSize; modality++) 738 vertex[modality] = lower[modality]; 739 740 GUM_SCALAR total = 0.; 741 742 auto vsize = vertex.size(); 743 744 for (Size i = 0; i < vsize; i++) 745 total += vertex[i]; 746 747 if (std::fabs(total - 1.) > 1e-6) 748 GUM_ERROR(CPTError, 749 _src_bn_.variable(node).name() 750 << " does not sum to one for " << entry << std::endl 751 << vertex << std::endl); 752 753 vertices.push_back(vertex); 754 } 755 756 var_cpt[entry] = vertices; 757 } 758 759 _credalNet_src_cpt_.insert(node, var_cpt); 760 761 } // end of : for each variable (node) 762 763 // get precise/credal/vacuous status of each variable 764 _sort_varType_(); 765 _separatelySpecified_ = true; 766 } 767 768 /* uses lrsWrapper */ 769 template < typename GUM_SCALAR > intervalToCredal()770 void CredalNet< GUM_SCALAR >::intervalToCredal() { 771 if (!_credalNet_src_cpt_.empty()) _credalNet_src_cpt_.clear(); 772 773 _credalNet_src_cpt_.resize(_src_bn_.size()); 774 775 LRSWrapper< GUM_SCALAR > lrsWrapper; 776 777 for (auto node: _src_bn_.nodes()) { 778 const Potential< GUM_SCALAR >* const potential_min(&_src_bn_min_.cpt(node)); 779 const Potential< GUM_SCALAR >* const potential_max(&_src_bn_max_.cpt(node)); 780 781 Size var_dSize = _src_bn_.variable(node).domainSize(); 782 Size entry_size = potential_min->domainSize() / var_dSize; 783 784 std::vector< std::vector< std::vector< GUM_SCALAR > > > var_cpt(entry_size); 785 786 Instantiation ins_min(potential_min); 787 Instantiation ins_max(potential_max); 788 789 ins_min.setFirst(); 790 ins_max.setFirst(); 791 792 lrsWrapper.setUpH(var_dSize); 793 794 for (Size entry = 0; entry < entry_size; entry++) { 795 for (Size modality = 0; modality < var_dSize; modality++) { 796 if (potential_min->get(ins_min) > potential_max->get(ins_max)) { 797 GUM_ERROR(CPTError, 798 "For variable " 799 << _src_bn_.variable(node).name() << " (at " << ins_min 800 << "), the min is greater than the max : " << potential_min->get(ins_min) 801 << ">" << potential_max->get(ins_max) << "."); 802 } 803 lrsWrapper.fillH(potential_min->get(ins_min), potential_max->get(ins_max), modality); 804 ++ins_min; 805 ++ins_max; 806 } 807 808 lrsWrapper.H2V(); 809 var_cpt[entry] = lrsWrapper.getOutput(); 810 lrsWrapper.nextHInput(); 811 } 812 813 _credalNet_src_cpt_.insert(node, var_cpt); 814 815 } // end of : for each variable (node) 816 817 // get precise/credal/vacuous status of each variable 818 _sort_varType_(); 819 _separatelySpecified_ = true; 820 } 821 822 /* call lrs */ 823 template < typename GUM_SCALAR > intervalToCredalWithFiles()824 void CredalNet< GUM_SCALAR >::intervalToCredalWithFiles() { 825 if (!_credalNet_src_cpt_.empty()) _credalNet_src_cpt_.clear(); 826 827 _credalNet_src_cpt_.resize(_src_bn_.size()); 828 829 for (auto node: _src_bn_.nodes()) { 830 const Potential< GUM_SCALAR >* const potential_min(&_src_bn_min_.cpt(node)); 831 const Potential< GUM_SCALAR >* const potential_max(&_src_bn_max_.cpt(node)); 832 833 auto var_dSize = _src_bn_.variable(node).domainSize(); 834 auto entry_size = potential_min->domainSize() / var_dSize; 835 836 std::vector< std::vector< std::vector< GUM_SCALAR > > > var_cpt(entry_size); 837 838 Instantiation ins_min(potential_min); 839 Instantiation ins_max(potential_max); 840 841 ins_min.setFirst(); 842 ins_max.setFirst(); 843 844 // use iterator 845 for (Size entry = 0; entry < entry_size; entry++) { 846 std::vector< std::vector< GUM_SCALAR > > vertices; 847 std::vector< GUM_SCALAR > vertex(var_dSize); // if not interval 848 849 std::vector< std::vector< GUM_SCALAR > > inequalities( 850 var_dSize * 2, 851 std::vector< GUM_SCALAR >(var_dSize + 1, 0)); 852 853 std::vector< GUM_SCALAR > sum_ineq1(var_dSize + 1, -1); 854 std::vector< GUM_SCALAR > sum_ineq2(var_dSize + 1, 1); 855 sum_ineq1[0] = 1; 856 sum_ineq2[0] = -1; 857 858 bool isInterval = false; 859 860 for (Size modality = 0; modality < var_dSize; modality++) { 861 inequalities[modality * 2][0] = -potential_min->get(ins_min); 862 inequalities[modality * 2 + 1][0] = potential_max->get(ins_max); 863 inequalities[modality * 2][modality + 1] = 1; 864 inequalities[modality * 2 + 1][modality + 1] = -1; 865 866 vertex[modality] = inequalities[modality * 2 + 1][0]; 867 868 if (!isInterval 869 && (-inequalities[modality * 2][0] != inequalities[modality * 2 + 1][0])) 870 isInterval = true; 871 872 ++ins_min; 873 ++ins_max; 874 } 875 876 inequalities.push_back(sum_ineq1); 877 inequalities.push_back(sum_ineq2); 878 879 if (!isInterval) { 880 vertices.push_back(vertex); 881 } else { 882 try { 883 _H2Vlrs_(inequalities, vertices); 884 // __H2Vcdd ( inequalities, vertices ); 885 } catch (const std::exception& err) { 886 std::cout << err.what() << std::endl; 887 throw; 888 } 889 890 } // end of : is interval 891 892 if (entry == 0 && vertices.size() >= 2) { 893 auto tmp = vertices[0]; 894 vertices[0] = vertices[1]; 895 vertices[1] = tmp; 896 } 897 898 var_cpt[entry] = vertices; 899 900 } // end of : for each entry 901 902 _credalNet_src_cpt_.insert(node, var_cpt); 903 // std::cout << _src_bn_.variable(node_idIt).name() << std::endl; 904 // std::cout << var_cpt << std::endl; 905 906 } // end of : for each variable (node) 907 908 // get precise/credal/vacuous status of each variable 909 _sort_varType_(); 910 _separatelySpecified_ = true; 911 } 912 913 /** 914 * to call after bnToCredal( GUM_SCALAR beta ) 915 * save a BN with lower probabilities and a BN with upper ones 916 */ 917 template < typename GUM_SCALAR > saveBNsMinMax(const std::string & min_path,const std::string & max_path)918 void CredalNet< GUM_SCALAR >::saveBNsMinMax(const std::string& min_path, 919 const std::string& max_path) const { 920 BIFWriter< GUM_SCALAR > writer; 921 922 std::string minfilename = min_path; //"min.bif"; 923 std::string maxfilename = max_path; //"max.bif"; 924 std::ofstream min_file(minfilename.c_str(), std::ios::out | std::ios::trunc); 925 std::ofstream max_file(maxfilename.c_str(), std::ios::out | std::ios::trunc); 926 927 if (!min_file.good()) 928 GUM_ERROR(IOError, "bnToCredal() : could not open stream : min_file : " << minfilename); 929 930 if (!max_file.good()) { 931 min_file.close(); 932 GUM_ERROR(IOError, "bnToCredal() : could not open stream : min_file : " << maxfilename); 933 } 934 935 try { 936 writer.write(min_file, _src_bn_min_); 937 writer.write(max_file, _src_bn_max_); 938 } catch (Exception& err) { 939 GUM_SHOWERROR(err); 940 min_file.close(); 941 max_file.close(); 942 throw(err); 943 } 944 945 min_file.close(); 946 max_file.close(); 947 } 948 949 template < typename GUM_SCALAR > approximatedBinarization()950 void CredalNet< GUM_SCALAR >::approximatedBinarization() { 951 // don't forget to delete the old one ( _current_), if necessary at the end 952 auto bin_bn = new BayesNet< GUM_SCALAR >(); 953 954 // __bnCopy ( * _bin_bn_ ); 955 // delete old one too 956 auto credalNet_bin_cpt 957 = new NodeProperty< std::vector< std::vector< std::vector< GUM_SCALAR > > > >(); 958 959 // delete old one too 960 auto bin_nodeType = new NodeProperty< NodeType >(); 961 962 const BayesNet< GUM_SCALAR >* current_bn; 963 // const NodeProperty< nodeType > * _current_nodeType_; 964 const NodeProperty< std::vector< std::vector< std::vector< GUM_SCALAR > > > >* 965 credalNet_current_cpt; 966 967 if (this->_current_bn_ == nullptr) 968 current_bn = &this->_src_bn_; 969 else 970 current_bn = this->_current_bn_; 971 972 if (this->_credalNet_current_cpt_ == nullptr) 973 credalNet_current_cpt = &this->_credalNet_src_cpt_; 974 else 975 credalNet_current_cpt = this->_credalNet_current_cpt_; 976 977 /*if ( this-> _current_nodeType_ == nullptr ) 978 _current_nodeType_ = & this-> _nodeType_; 979 else 980 _current_nodeType_ = this-> _current_nodeType_;*/ 981 982 if (!_var_bits_.empty()) _var_bits_.clear(); 983 984 bin_bn->beginTopologyTransformation(); 985 986 for (auto node: current_bn->nodes()) { 987 auto var_dSize = current_bn->variable(node).domainSize(); 988 989 if (var_dSize != 2) { 990 unsigned long b, c; 991 superiorPow((unsigned long)var_dSize, b, c); 992 Size nb_bits{Size(b)}; 993 994 std::string bit_name; 995 std::vector< NodeId > bits(nb_bits); 996 997 for (Size bit = 0; bit < nb_bits; bit++) { 998 bit_name = current_bn->variable(node).name() + "-b"; 999 std::stringstream ss; 1000 ss << bit; 1001 bit_name += ss.str(); 1002 1003 LabelizedVariable var_bit(bit_name, "node " + bit_name, 2); 1004 NodeId iD = bin_bn->add(var_bit); 1005 1006 bits[bit] = iD; 1007 } // end of : for each bit 1008 1009 _var_bits_.insert(node, bits); 1010 1011 } // end of : if variable is not binary 1012 else { 1013 std::string bit_name = current_bn->variable(node).name(); 1014 LabelizedVariable var_bit(bit_name, "node " + bit_name, 2); 1015 NodeId iD = bin_bn->add(var_bit); 1016 1017 _var_bits_.insert(node, std::vector< NodeId >(1, iD)); 1018 } 1019 1020 } // end of : for each original variable 1021 1022 for (auto node: current_bn->nodes()) { 1023 NodeSet parents = current_bn->parents(node); 1024 1025 if (!parents.empty()) { 1026 for (auto par: current_bn->parents(node)) { 1027 for (Size parent_bit = 0, spbits = Size(_var_bits_[par].size()); parent_bit < spbits; 1028 parent_bit++) 1029 for (Size var_bit = 0, mbits = Size(_var_bits_[node].size()); var_bit < mbits; 1030 var_bit++) 1031 bin_bn->addArc(_var_bits_[par][parent_bit], _var_bits_[node][var_bit]); 1032 } 1033 } 1034 1035 // arcs with one's bits 1036 auto bitsize = _var_bits_[node].size(); 1037 1038 for (Size bit_c = 1; bit_c < bitsize; bit_c++) 1039 for (Size bit_p = 0; bit_p < bit_c; bit_p++) 1040 bin_bn->addArc(_var_bits_[node][bit_p], _var_bits_[node][bit_c]); 1041 1042 } // end of : for each original variable 1043 1044 bin_bn->endTopologyTransformation(); 1045 1046 // binarization of cpts 1047 1048 auto varsize = current_bn->size(); 1049 1050 for (Size var = 0; var < varsize; var++) { 1051 auto bitsize = _var_bits_[var].size(); 1052 1053 for (Size i = 0; i < bitsize; i++) { 1054 Potential< GUM_SCALAR > const* potential(&bin_bn->cpt(_var_bits_[var][i])); 1055 Instantiation ins(potential); 1056 ins.setFirst(); 1057 1058 auto entry_size = potential->domainSize() / 2; 1059 std::vector< std::vector< std::vector< GUM_SCALAR > > > var_cpt(entry_size); 1060 1061 Size old_conf = 0; 1062 1063 for (Size conf = 0; conf < entry_size; conf++) { 1064 std::vector< std::vector< GUM_SCALAR > > pvar_cpt; 1065 auto verticessize = (*credalNet_current_cpt)[var][old_conf].size(); 1066 1067 for (Size old_distri = 0; old_distri < verticessize; old_distri++) { 1068 const std::vector< GUM_SCALAR >& vertex 1069 = (*credalNet_current_cpt)[var][old_conf][old_distri]; 1070 auto vertexsize = vertex.size(); 1071 1072 std::vector< Idx > incc(vertexsize, 0); 1073 1074 for (Size preced = 0; preced < i; preced++) { 1075 auto bit_pos = ins.pos(bin_bn->variable(_var_bits_[var][preced])); 1076 auto val = ins.val(bit_pos); 1077 1078 Size pas = Size(int2Pow((unsigned long)preced)); 1079 Size elem; 1080 1081 if (val == 0) 1082 elem = 0; 1083 else 1084 elem = pas; 1085 1086 while (elem < vertexsize) { 1087 incc[elem]++; 1088 elem++; 1089 1090 if (elem % pas == 0) elem += pas; 1091 } 1092 } 1093 1094 Size pas = Size(int2Pow((unsigned long)i)); 1095 1096 std::vector< GUM_SCALAR > distri(2, 0); 1097 int pos = 1; 1098 1099 for (Size elem = 0; elem < vertexsize; elem++) { 1100 if (elem % pas == 0) pos = -pos; 1101 1102 if (incc[elem] == i) 1103 (pos < 0) ? (distri[0] += vertex[elem]) : (distri[1] += vertex[elem]); 1104 } 1105 1106 if (i > 0) { 1107 GUM_SCALAR den = distri[0] + distri[1]; 1108 1109 if (den == 0) { 1110 distri[0] = 0; 1111 distri[1] = 0; 1112 } else { 1113 distri[0] /= den; 1114 distri[1] /= den; 1115 } 1116 } 1117 1118 pvar_cpt.push_back(distri); 1119 1120 } // end of old distris 1121 1122 // get min/max approx, 2 vertices 1123 std::vector< std::vector< GUM_SCALAR > > vertices(2, std::vector< GUM_SCALAR >(2, 1)); 1124 vertices[1][1] = 0; 1125 1126 auto new_verticessize = pvar_cpt.size(); 1127 1128 for (Size v = 0; v < new_verticessize; v++) { 1129 if (pvar_cpt[v][1] < vertices[0][1]) vertices[0][1] = pvar_cpt[v][1]; 1130 1131 if (pvar_cpt[v][1] > vertices[1][1]) vertices[1][1] = pvar_cpt[v][1]; 1132 } 1133 1134 vertices[0][0] = 1 - vertices[0][1]; 1135 vertices[1][0] = 1 - vertices[1][1]; 1136 1137 pvar_cpt = vertices; 1138 1139 var_cpt[conf] = pvar_cpt; 1140 1141 ++ins; 1142 ++ins; 1143 1144 old_conf++; 1145 1146 if (old_conf == (*credalNet_current_cpt)[var].size()) old_conf = 0; 1147 1148 } // end of new parent conf 1149 1150 credalNet_bin_cpt->insert(_var_bits_[var][i], var_cpt); 1151 1152 } // end of bit i 1153 1154 } // end of old variable 1155 1156 bin_bn->beginTopologyTransformation(); 1157 1158 /* indicatrices variables */ 1159 auto old_varsize = _var_bits_.size(); 1160 1161 for (Size i = 0; i < old_varsize; i++) { 1162 auto bitsize = _var_bits_[i].size(); 1163 1164 // binary variable 1165 if (bitsize == 1) continue; 1166 1167 auto old_card = _src_bn_.variable(i).domainSize(); 1168 1169 for (Size mod = 0; mod < old_card; mod++) { 1170 std::stringstream ss; 1171 ss << _src_bn_.variable(i).name(); 1172 ss << "-v"; 1173 ss << mod; 1174 1175 LabelizedVariable var(ss.str(), "node " + ss.str(), 2); 1176 const NodeId indic = bin_bn->add(var); 1177 1178 // arcs from one's bits 1179 for (Size bit = 0; bit < bitsize; bit++) 1180 bin_bn->addArc(_var_bits_[i][bit], indic); 1181 1182 // cpt 1183 Size num = Size(int2Pow(long(bitsize))); 1184 1185 std::vector< std::vector< std::vector< GUM_SCALAR > > > icpt(num); 1186 1187 for (Size entry = 0; entry < num; entry++) { 1188 std::vector< std::vector< GUM_SCALAR > > vertices(1, std::vector< GUM_SCALAR >(2, 0)); 1189 1190 if (mod == entry) 1191 vertices[0][1] = 1; 1192 else 1193 vertices[0][0] = 1; 1194 1195 icpt[entry] = vertices; 1196 } 1197 1198 credalNet_bin_cpt->insert(indic, icpt); 1199 1200 bin_nodeType->insert(indic, NodeType::Indic); 1201 } // end of each modality, i.e. as many indicatrice 1202 } 1203 1204 bin_bn->endTopologyTransformation(); 1205 1206 if (this->_current_bn_ != nullptr) delete this->_current_bn_; 1207 1208 this->_current_bn_ = bin_bn; 1209 1210 if (this->_credalNet_current_cpt_ != nullptr) delete this->_credalNet_current_cpt_; 1211 1212 this->_credalNet_current_cpt_ = credalNet_bin_cpt; 1213 1214 if (this->_current_nodeType_ != nullptr) delete this->_current_nodeType_; 1215 1216 this->_current_nodeType_ = bin_nodeType; 1217 1218 _sort_varType_(); // will fill _bin_nodeType_ except for NodeType::Indic 1219 // variables 1220 1221 computeBinaryCPTMinMax(); 1222 } 1223 1224 template < typename GUM_SCALAR > 1225 const NodeProperty< std::vector< std::vector< std::vector< GUM_SCALAR > > > >& credalNet_currentCpt()1226 CredalNet< GUM_SCALAR >::credalNet_currentCpt() const { 1227 if (_credalNet_current_cpt_ != nullptr) return *_credalNet_current_cpt_; 1228 1229 return _credalNet_src_cpt_; 1230 } 1231 1232 template < typename GUM_SCALAR > 1233 const NodeProperty< std::vector< std::vector< std::vector< GUM_SCALAR > > > >& credalNet_srcCpt()1234 CredalNet< GUM_SCALAR >::credalNet_srcCpt() const { 1235 return _credalNet_src_cpt_; 1236 } 1237 1238 template < typename GUM_SCALAR > 1239 typename CredalNet< GUM_SCALAR >::NodeType currentNodeType(const NodeId & id)1240 CredalNet< GUM_SCALAR >::currentNodeType(const NodeId& id) const { 1241 if (_current_nodeType_ != nullptr) return (*(_current_nodeType_))[id]; 1242 1243 return _original_nodeType_[id]; 1244 } 1245 1246 template < typename GUM_SCALAR > 1247 typename CredalNet< GUM_SCALAR >::NodeType nodeType(const NodeId & id)1248 CredalNet< GUM_SCALAR >::nodeType(const NodeId& id) const { 1249 return _original_nodeType_[id]; 1250 } 1251 1252 template < typename GUM_SCALAR > isSeparatelySpecified()1253 const bool CredalNet< GUM_SCALAR >::isSeparatelySpecified() const { 1254 return _separatelySpecified_; 1255 } 1256 1257 template < typename GUM_SCALAR > hasComputedBinaryCPTMinMax()1258 const bool CredalNet< GUM_SCALAR >::hasComputedBinaryCPTMinMax() const { 1259 return _hasComputedBinaryCPTMinMax_; 1260 } 1261 1262 // only if CN is binary !! 1263 template < typename GUM_SCALAR > computeBinaryCPTMinMax()1264 void CredalNet< GUM_SCALAR >::computeBinaryCPTMinMax() { 1265 _binCptMin_.resize(current_bn().size()); 1266 _binCptMax_.resize(current_bn().size()); 1267 1268 for (auto node: current_bn().nodes()) { 1269 auto pConf = credalNet_currentCpt()[node].size(); 1270 std::vector< GUM_SCALAR > min(pConf); 1271 std::vector< GUM_SCALAR > max(pConf); 1272 1273 for (Size pconf = 0; pconf < pConf; pconf++) { 1274 GUM_SCALAR v1, v2; 1275 v1 = credalNet_currentCpt()[node][pconf][0][1]; 1276 1277 if (credalNet_currentCpt()[node][pconf].size() > 1) 1278 v2 = credalNet_currentCpt()[node][pconf][1][1]; 1279 else 1280 v2 = v1; 1281 1282 GUM_SCALAR delta = v1 - v2; 1283 min[pconf] = (delta >= 0) ? v2 : v1; 1284 max[pconf] = (delta >= 0) ? v1 : v2; 1285 } 1286 1287 _binCptMin_[node] = min; 1288 _binCptMax_[node] = max; 1289 } 1290 1291 _hasComputedBinaryCPTMinMax_ = true; 1292 } 1293 1294 template < typename GUM_SCALAR > 1295 const std::vector< std::vector< GUM_SCALAR > >& get_binaryCPT_min()1296 CredalNet< GUM_SCALAR >::get_binaryCPT_min() const { 1297 return _binCptMin_; 1298 } 1299 1300 template < typename GUM_SCALAR > 1301 const std::vector< std::vector< GUM_SCALAR > >& get_binaryCPT_max()1302 CredalNet< GUM_SCALAR >::get_binaryCPT_max() const { 1303 return _binCptMax_; 1304 } 1305 1306 template < typename GUM_SCALAR > epsilonMin()1307 const GUM_SCALAR& CredalNet< GUM_SCALAR >::epsilonMin() const { 1308 return _epsilonMin_; 1309 } 1310 1311 template < typename GUM_SCALAR > epsilonMax()1312 const GUM_SCALAR& CredalNet< GUM_SCALAR >::epsilonMax() const { 1313 return _epsilonMax_; 1314 } 1315 1316 template < typename GUM_SCALAR > epsilonMean()1317 const GUM_SCALAR& CredalNet< GUM_SCALAR >::epsilonMean() const { 1318 return _epsilonMoy_; 1319 } 1320 1321 template < typename GUM_SCALAR > toString()1322 std::string CredalNet< GUM_SCALAR >::toString() const { 1323 std::stringstream output; 1324 const BayesNet< GUM_SCALAR >* _current_bn_; 1325 const NodeProperty< std::vector< std::vector< std::vector< GUM_SCALAR > > > >* 1326 _credalNet_current_cpt_; 1327 1328 if (this->_current_bn_ == nullptr) 1329 _current_bn_ = &this->_src_bn_; 1330 else 1331 _current_bn_ = this->_current_bn_; 1332 1333 if (this->_credalNet_current_cpt_ == nullptr) 1334 _credalNet_current_cpt_ = &this->_credalNet_src_cpt_; 1335 else 1336 _credalNet_current_cpt_ = this->_credalNet_current_cpt_; 1337 1338 for (auto node: _current_bn_->nodes()) { 1339 const Potential< GUM_SCALAR >* potential(&_current_bn_->cpt(node)); 1340 auto pconfs = potential->domainSize() / _current_bn_->variable(node).domainSize(); 1341 1342 output << "\n" << _current_bn_->variable(node) << "\n"; 1343 1344 Instantiation ins(potential); 1345 ins.forgetMaster(); 1346 ins.erase(_current_bn_->variable(node)); 1347 ins.setFirst(); 1348 1349 for (Size pconf = 0; pconf < pconfs; pconf++) { 1350 output << ins << " : "; 1351 output << (*_credalNet_current_cpt_)[node][pconf] << "\n"; 1352 1353 if (pconf < pconfs - 1) ++ins; 1354 } 1355 } 1356 1357 output << "\n"; 1358 1359 return output.str(); 1360 } 1361 1362 template < typename GUM_SCALAR > current_bn()1363 const BayesNet< GUM_SCALAR >& CredalNet< GUM_SCALAR >::current_bn() const { 1364 if (_current_bn_ != nullptr) return *_current_bn_; 1365 1366 return _src_bn_; 1367 } 1368 1369 template < typename GUM_SCALAR > src_bn()1370 const BayesNet< GUM_SCALAR >& CredalNet< GUM_SCALAR >::src_bn() const { 1371 return _src_bn_; 1372 } 1373 1374 /////////// protected stuff ////////// 1375 1376 /////////// private stuff //////////// 1377 1378 template < typename GUM_SCALAR > _initParams_()1379 void CredalNet< GUM_SCALAR >::_initParams_() { 1380 _epsilonMin_ = 0; 1381 _epsilonMax_ = 0; 1382 _epsilonMoy_ = 0; 1383 1384 _epsRedund_ = GUM_SCALAR(1e-6); 1385 1386 // farey algorithm 1387 _epsF_ = GUM_SCALAR(1e-6); 1388 _denMax_ = GUM_SCALAR(1e6); // beware LRSWrapper 1389 1390 // continued fractions, beware LRSWrapper 1391 // decimal paces ( _epsC_ * _precisionC_ == 1) 1392 _precisionC_ = GUM_SCALAR(1e6); 1393 _deltaC_ = 5; 1394 1395 // old custom algorithm 1396 _precision_ = GUM_SCALAR(1e6); // beware LRSWrapper 1397 1398 _current_bn_ = nullptr; 1399 _credalNet_current_cpt_ = nullptr; 1400 _current_nodeType_ = nullptr; 1401 1402 _hasComputedBinaryCPTMinMax_ = false; 1403 } 1404 1405 template < typename GUM_SCALAR > _initCNNets_(const std::string & src_min_num,const std::string & src_max_den)1406 void CredalNet< GUM_SCALAR >::_initCNNets_(const std::string& src_min_num, 1407 const std::string& src_max_den) { 1408 BIFReader< GUM_SCALAR > reader(&_src_bn_, src_min_num); 1409 std::string other; 1410 1411 if (src_max_den.compare("") != 0) 1412 other = src_max_den; 1413 else 1414 other = src_min_num; 1415 1416 BIFReader< GUM_SCALAR > reader_min(&_src_bn_min_, src_min_num); 1417 BIFReader< GUM_SCALAR > reader_max(&_src_bn_max_, other); 1418 1419 reader.proceed(); 1420 reader_min.proceed(); 1421 reader_max.proceed(); 1422 } 1423 1424 template < typename GUM_SCALAR > _initCNNets_(const BayesNet<GUM_SCALAR> & src_min_num,const BayesNet<GUM_SCALAR> & src_max_den)1425 void CredalNet< GUM_SCALAR >::_initCNNets_(const BayesNet< GUM_SCALAR >& src_min_num, 1426 const BayesNet< GUM_SCALAR >& src_max_den) { 1427 _src_bn_ = src_min_num; 1428 _src_bn_min_ = src_min_num; 1429 1430 if (src_max_den.size() > 0) 1431 _src_bn_max_ = src_max_den; 1432 else 1433 _src_bn_max_ = src_min_num; 1434 } 1435 1436 template < typename GUM_SCALAR > _find_dNode_card_(const std::vector<std::vector<std::vector<GUM_SCALAR>>> & var_cpt)1437 int CredalNet< GUM_SCALAR >::_find_dNode_card_( 1438 const std::vector< std::vector< std::vector< GUM_SCALAR > > >& var_cpt) const { 1439 Size vertices_size = 0; 1440 1441 for (auto entry = var_cpt.cbegin(), theEnd = var_cpt.cend(); entry != theEnd; ++entry) { 1442 if (entry->size() > vertices_size) vertices_size = Size(entry->size()); 1443 } 1444 1445 return int(vertices_size); 1446 } 1447 1448 template < typename GUM_SCALAR > _bnCopy_(BayesNet<GUM_SCALAR> & dest)1449 void CredalNet< GUM_SCALAR >::_bnCopy_(BayesNet< GUM_SCALAR >& dest) { 1450 const BayesNet< GUM_SCALAR >* _current_bn_; 1451 1452 if (this->_current_bn_ == nullptr) 1453 _current_bn_ = &this->_src_bn_; 1454 else 1455 _current_bn_ = this->_current_bn_; 1456 1457 for (auto node: _current_bn_->nodes()) 1458 dest.add(_current_bn_->variable(node)); 1459 1460 dest.beginTopologyTransformation(); 1461 1462 for (auto node: _current_bn_->nodes()) { 1463 for (auto parent_idIt: _current_bn_->cpt(node).variablesSequence()) { 1464 if (_current_bn_->nodeId(*parent_idIt) != node) 1465 dest.addArc(_current_bn_->nodeId(*parent_idIt), node); 1466 } // end of : for each parent in order of appearence 1467 } // end of : for each variable 1468 1469 dest.endTopologyTransformation(); 1470 } 1471 1472 /* 1473 // cdd can use real values, not just rationals / integers 1474 template< typename GUM_SCALAR > 1475 void CredalNet< GUM_SCALAR >:: _H2Vcdd_ ( const std::vector< std::vector< 1476 GUM_SCALAR > > & h_rep, std::vector< std::vector< GUM_SCALAR > > & v_rep ) 1477 const { 1478 dd_set_global_constants(); 1479 1480 dd_MatrixPtr M, G; 1481 dd_PolyhedraPtr poly; 1482 dd_ErrorType err; 1483 1484 unsigned int rows = h_rep.size(); 1485 unsigned int cols = 0; 1486 if( h_rep.size() > 0 ) 1487 cols = h_rep[0].size(); 1488 1489 M = dd_CreateMatrix( rows, cols); 1490 1491 for ( unsigned int row = 0; row < rows; row++ ) 1492 for ( unsigned int col = 0; col < cols; col++ ) 1493 dd_set_d( M->matrix[row][col], h_rep[row][col] ); 1494 1495 M->representation = dd_Inequality; 1496 1497 poly = dd_DDMatrix2Poly(M, &err); 1498 G = dd_CopyGenerators(poly); 1499 1500 rows = G->rowsize; 1501 cols = G->colsize; 1502 1503 v_rep.clear(); 1504 for ( unsigned int row = 0; row < rows; row++ ) { 1505 std::vector< GUM_SCALAR > aRow(cols - 1); 1506 1507 if ( *G->matrix[row][0] != 1 ) 1508 GUM_ERROR(OperationNotAllowed, " __H2Vcdd : not reading a vertex") 1509 1510 for ( unsigned int col = 0; col < cols - 1; col++ ) 1511 aRow[col] = *G->matrix[row][ col + 1 ]; 1512 1513 v_rep.push_back(aRow); 1514 } 1515 1516 dd_FreeMatrix(M); 1517 dd_FreeMatrix(G); 1518 dd_FreePolyhedra(poly); 1519 1520 dd_free_global_constants(); 1521 } 1522 */ 1523 1524 template < typename GUM_SCALAR > _H2Vlrs_(const std::vector<std::vector<GUM_SCALAR>> & h_rep,std::vector<std::vector<GUM_SCALAR>> & v_rep)1525 void CredalNet< GUM_SCALAR >::_H2Vlrs_(const std::vector< std::vector< GUM_SCALAR > >& h_rep, 1526 std::vector< std::vector< GUM_SCALAR > >& v_rep) const { 1527 // write H rep file 1528 int64_t num, den; 1529 1530 std::string sinefile = getUniqueFileName(); // generate unique file name, we 1531 // need to add .ine or .ext for lrs 1532 // to know which input it is (Hrep 1533 // to Vrep or Vrep to Hrep) 1534 sinefile += ".ine"; 1535 1536 std::ofstream h_file(sinefile.c_str(), std::ios::out | std::ios::trunc); 1537 1538 if (!h_file.good()) 1539 GUM_ERROR(IOError, " __H2Vlrs : could not open lrs input file : " << sinefile) 1540 1541 h_file << "H - representation\n"; 1542 h_file << "begin\n"; 1543 h_file << h_rep.size() << ' ' << h_rep[0].size() << " rational\n"; 1544 1545 for (auto it = h_rep.cbegin(), theEnd = h_rep.cend(); it != theEnd; ++it) { 1546 for (auto it2 = it->cbegin(), theEnd2 = it->cend(); it2 != theEnd2; ++it2) { 1547 // get integer fraction from decimal value 1548 // smallest numerator & denominator is farley, also 1549 // best precision 1550 Rational< GUM_SCALAR >::farey(num, 1551 den, 1552 ((*it2 > 0) ? *it2 : -*it2), 1553 int64_t(_denMax_), 1554 _epsF_); 1555 1556 h_file << ((*it2 > 0) ? num : -num) << '/' << den << ' '; 1557 } 1558 1559 h_file << '\n'; 1560 } 1561 1562 h_file << "end\n"; 1563 h_file.close(); 1564 1565 // call lrs 1566 // lrs arguments 1567 char* args[3]; 1568 1569 std::string soft_name = "lrs"; 1570 std::string extfile(sinefile); 1571 extfile += ".ext"; 1572 1573 args[0] = new char[soft_name.size()]; 1574 args[1] = new char[sinefile.size()]; 1575 args[2] = new char[extfile.size()]; 1576 1577 strcpy(args[0], soft_name.c_str()); 1578 strcpy(args[1], sinefile.c_str()); 1579 strcpy(args[2], extfile.c_str()); 1580 1581 // standard cout to null (avoid lrs flooding) 1582 int old_cout, new_cout; 1583 fflush(stdout); 1584 old_cout = dup(1); 1585 1586 new_cout = open("/dev/null", O_WRONLY); 1587 dup2(new_cout, 1); 1588 close(new_cout); 1589 1590 lrs_main(3, args); 1591 1592 // restore standard cout 1593 fflush(stdout); 1594 dup2(old_cout, 1); 1595 close(old_cout); 1596 1597 delete[] args[2]; 1598 delete[] args[1]; 1599 delete[] args[0]; 1600 1601 // read V rep file 1602 std::ifstream v_file(extfile.c_str() /*extfilename.c_str()*/, std::ios::in); 1603 1604 if (!v_file.good()) GUM_ERROR(IOError, " __H2Vlrs : could not open lrs ouput file : ") 1605 1606 std::string line, tmp; 1607 char * cstr, *p; 1608 GUM_SCALAR probability; 1609 1610 std::string::size_type pos; 1611 bool keep_going = true; 1612 // int vertices; 1613 1614 std::vector< GUM_SCALAR > vertex; 1615 1616 v_file.ignore(256, 'l'); 1617 1618 while (v_file.good() && keep_going) { 1619 getline(v_file, line); 1620 1621 if (line.size() == 0) 1622 continue; 1623 else if (line.compare("end") == 0) { 1624 keep_going = false; 1625 // this is to get vertices number : 1626 /*getline ( v_file, line ); 1627 std::string::size_type pos, end_pos; 1628 pos = line.find ( "vertices = " ); 1629 end_pos = line.find ( "rays", pos + 9 ); 1630 vertices = atoi ( line.substr ( pos + 9, end_pos - pos - 9 ).c_str() 1631 );*/ 1632 break; 1633 } else if (line[1] != '1') { 1634 GUM_ERROR(IOError, 1635 " __H2Vlrs : reading something other than a vertex from " 1636 "lrs output file : "); 1637 } 1638 1639 line = line.substr(2); 1640 cstr = new char[line.size() + 1]; 1641 strcpy(cstr, line.c_str()); 1642 1643 p = strtok(cstr, " "); 1644 1645 while (p != nullptr) { 1646 tmp = p; 1647 1648 if (tmp.compare("1") == 0 || tmp.compare("0") == 0) 1649 probability = GUM_SCALAR(atof(tmp.c_str())); 1650 else { 1651 pos = tmp.find("/"); 1652 probability = GUM_SCALAR(atof(tmp.substr(0, pos).c_str()) 1653 / atof(tmp.substr(pos + 1, tmp.size()).c_str())); 1654 } 1655 1656 vertex.push_back(probability); 1657 p = strtok(nullptr, " "); 1658 } // end of : for all tokens 1659 1660 delete[] p; 1661 delete[] cstr; 1662 1663 bool is_redund = false; 1664 1665 #pragma omp parallel 1666 { 1667 int this_thread = getThreadNumber(); 1668 int num_threads = getNumberOfRunningThreads(); 1669 1670 auto begin_pos = (this_thread + 0) * v_rep.size() / num_threads; 1671 auto end_pos = (this_thread + 1) * v_rep.size() / num_threads; 1672 1673 for (auto p = begin_pos; p < end_pos; p++) { 1674 #pragma omp flush(is_redund) 1675 1676 if (is_redund) break; 1677 1678 bool thread_redund = true; 1679 1680 auto vsize = vertex.size(); 1681 1682 for (Size modality = 0; modality < vsize; modality++) { 1683 if (std::fabs(vertex[modality] - v_rep[p][modality]) > _epsRedund_) { 1684 thread_redund = false; 1685 break; 1686 } 1687 } 1688 1689 if (thread_redund) { 1690 is_redund = true; 1691 #pragma omp flush(is_redund) 1692 } 1693 } // end of : each thread for 1694 } // end of : parallel 1695 1696 if (!is_redund) v_rep.push_back(vertex); 1697 1698 vertex.clear(); 1699 1700 } // end of : file 1701 1702 v_file.close(); 1703 1704 if (std::remove(sinefile.c_str()) != 0) GUM_ERROR(IOError, "error removing : " + sinefile) 1705 1706 if (std::remove(extfile.c_str()) != 0) GUM_ERROR(IOError, "error removing : " + extfile) 1707 } 1708 1709 template < typename GUM_SCALAR > _sort_varType_()1710 void CredalNet< GUM_SCALAR >::_sort_varType_() { 1711 NodeProperty< NodeType >* _current_nodeType_; 1712 const NodeProperty< std::vector< std::vector< std::vector< GUM_SCALAR > > > >* 1713 _credalNet_current_cpt_; 1714 1715 const BayesNet< GUM_SCALAR >* _current_bn_; 1716 1717 if (this->_current_bn_ == nullptr) 1718 _current_bn_ = &_src_bn_; 1719 else 1720 _current_bn_ = this->_current_bn_; 1721 1722 if (this->_credalNet_current_cpt_ == nullptr) 1723 _credalNet_current_cpt_ = &_credalNet_src_cpt_; 1724 else 1725 _credalNet_current_cpt_ = this->_credalNet_current_cpt_; 1726 1727 if (this->_current_nodeType_ == nullptr) 1728 _current_nodeType_ = &_original_nodeType_; 1729 else 1730 _current_nodeType_ = this->_current_nodeType_; 1731 1732 /*if ( ! _current_nodeType_->empty() ) 1733 _current_nodeType_->clear();*/ 1734 1735 for (auto node: _current_bn_->nodes()) { 1736 // indicatrices are already present 1737 if (_current_nodeType_->exists(node)) continue; 1738 1739 bool precise = true, vacuous = true; 1740 1741 for (auto entry = (*_credalNet_current_cpt_)[node].cbegin(), 1742 theEnd2 = (*_credalNet_current_cpt_)[node].cend(); 1743 entry != theEnd2; 1744 ++entry) { 1745 auto vertices = entry->size(); 1746 auto var_dSize = (*entry)[0].size(); 1747 1748 if (precise && vertices > 1) precise = false; 1749 1750 if (vacuous && vertices == var_dSize) { 1751 std::vector< bool > elem(var_dSize, false); 1752 1753 for (auto vertex = entry->cbegin(), vEnd = entry->cend(); vertex != vEnd; ++vertex) { 1754 for (auto probability = vertex->cbegin(), pEnd = vertex->cend(); probability != pEnd; 1755 ++probability) { 1756 if (*probability == 1) { 1757 elem[probability - vertex->begin()] = true; 1758 break; 1759 } 1760 } // end of : for each modality 1761 1762 break; // not vacuous 1763 } // end of : for each vertex 1764 1765 for (auto /*std::vector< bool >::const_iterator*/ probability = elem.cbegin(); 1766 probability != elem.cend(); 1767 ++probability) 1768 if (*probability == false) vacuous = false; 1769 1770 } // end of : if vertices == dSize 1771 else 1772 vacuous = false; 1773 1774 if (vacuous == false && precise == false) { 1775 _current_nodeType_->insert(node, NodeType::Credal); 1776 break; 1777 } 1778 1779 } // end of : for each parents entry 1780 1781 if (vacuous) 1782 _current_nodeType_->insert(node, NodeType::Vacuous); 1783 else if (precise) 1784 _current_nodeType_->insert(node, NodeType::Precise); 1785 1786 } // end of : for each variable 1787 } 1788 1789 } // namespace credal 1790 } // namespace gum 1791