1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 * -------------------------------------------------------------------------- * 3 * Lepton * 4 * -------------------------------------------------------------------------- * 5 * This is part of the Lepton expression parser originating from * 6 * Simbios, the NIH National Center for Physics-Based Simulation of * 7 * Biological Structures at Stanford, funded under the NIH Roadmap for * 8 * Medical Research, grant U54 GM072970. See https://simtk.org. * 9 * * 10 * Portions copyright (c) 2013-2016 Stanford University and the Authors. * 11 * Authors: Peter Eastman * 12 * Contributors: * 13 * * 14 * Permission is hereby granted, free of charge, to any person obtaining a * 15 * copy of this software and associated documentation files (the "Software"), * 16 * to deal in the Software without restriction, including without limitation * 17 * the rights to use, copy, modify, merge, publish, distribute, sublicense, * 18 * and/or sell copies of the Software, and to permit persons to whom the * 19 * Software is furnished to do so, subject to the following conditions: * 20 * * 21 * The above copyright notice and this permission notice shall be included in * 22 * all copies or substantial portions of the Software. * 23 * * 24 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * 25 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * 26 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * 27 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, * 28 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR * 29 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE * 30 * USE OR OTHER DEALINGS IN THE SOFTWARE. * 31 * -------------------------------------------------------------------------- * 32 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ 33 #ifndef __PLUMED_lepton_Operation_h 34 #define __PLUMED_lepton_Operation_h 35 36 /* -------------------------------------------------------------------------- * 37 * lepton * 38 * -------------------------------------------------------------------------- * 39 * This is part of the lepton expression parser originating from * 40 * Simbios, the NIH National Center for Physics-Based Simulation of * 41 * Biological Structures at Stanford, funded under the NIH Roadmap for * 42 * Medical Research, grant U54 GM072970. See https://simtk.org. * 43 * * 44 * Portions copyright (c) 2009-2015 Stanford University and the Authors. * 45 * Authors: Peter Eastman * 46 * Contributors: * 47 * * 48 * Permission is hereby granted, free of charge, to any person obtaining a * 49 * copy of this software and associated documentation files (the "Software"), * 50 * to deal in the Software without restriction, including without limitation * 51 * the rights to use, copy, modify, merge, publish, distribute, sublicense, * 52 * and/or sell copies of the Software, and to permit persons to whom the * 53 * Software is furnished to do so, subject to the following conditions: * 54 * * 55 * The above copyright notice and this permission notice shall be included in * 56 * all copies or substantial portions of the Software. * 57 * * 58 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * 59 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * 60 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * 61 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, * 62 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR * 63 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE * 64 * USE OR OTHER DEALINGS IN THE SOFTWARE. * 65 * -------------------------------------------------------------------------- */ 66 67 #include "windowsIncludes.h" 68 #include "CustomFunction.h" 69 #include "Exception.h" 70 #include <cmath> 71 #include <map> 72 #include <string> 73 #include <vector> 74 #include <sstream> 75 #include <algorithm> 76 #include <limits> 77 78 namespace PLMD { 79 namespace lepton { 80 81 class ExpressionTreeNode; 82 83 /** 84 * An Operation represents a single step in the evaluation of an expression, such as a function, 85 * an operator, or a constant value. Each Operation takes some number of values as arguments 86 * and produces a single value. 87 * 88 * This is an abstract class with subclasses for specific operations. 89 */ 90 91 class LEPTON_EXPORT Operation { 92 public: ~Operation()93 virtual ~Operation() { 94 } 95 /** 96 * This enumeration lists all Operation subclasses. This is provided so that switch statements 97 * can be used when processing or analyzing parsed expressions. 98 */ 99 enum Id {CONSTANT, VARIABLE, CUSTOM, ADD, SUBTRACT, MULTIPLY, DIVIDE, POWER, NEGATE, SQRT, EXP, LOG, 100 SIN, COS, SEC, CSC, TAN, COT, ASIN, ACOS, ATAN, SINH, COSH, TANH, ERF, ERFC, STEP, DELTA, NANDELTA, SQUARE, CUBE, RECIPROCAL, 101 ADD_CONSTANT, MULTIPLY_CONSTANT, POWER_CONSTANT, MIN, MAX, ABS, FLOOR, CEIL, SELECT, 102 ACOT, ASEC, ACSC, COTH, SECH, CSCH, ASINH, ACOSH, ATANH, ACOTH, ASECH, ACSCH, ATAN2}; 103 /** 104 * Get the name of this Operation. 105 */ 106 virtual std::string getName() const = 0; 107 /** 108 * Get this Operation's ID. 109 */ 110 virtual Id getId() const = 0; 111 /** 112 * Get the number of arguments this operation expects. 113 */ 114 virtual int getNumArguments() const = 0; 115 /** 116 * Create a clone of this Operation. 117 */ 118 virtual Operation* clone() const = 0; 119 /** 120 * Perform the computation represented by this operation. 121 * 122 * @param args the array of arguments 123 * @param variables a map containing the values of all variables 124 * @return the result of performing the computation. 125 */ 126 virtual double evaluate(double* args, const std::map<std::string, double>& variables) const = 0; 127 /** 128 * Return an ExpressionTreeNode which represents the analytic derivative of this Operation with respect to a variable. 129 * 130 * @param children the child nodes 131 * @param childDerivs the derivatives of the child nodes with respect to the variable 132 * @param variable the variable with respect to which the derivate should be taken 133 */ 134 virtual ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const = 0; 135 /** 136 * Get whether this operation should be displayed with infix notation. 137 */ isInfixOperator()138 virtual bool isInfixOperator() const { 139 return false; 140 } 141 /** 142 * Get whether this is a symmetric binary operation, such that exchanging its arguments 143 * does not affect the result. 144 */ isSymmetric()145 virtual bool isSymmetric() const { 146 return false; 147 } 148 virtual bool operator!=(const Operation& op) const { 149 return op.getId() != getId(); 150 } 151 virtual bool operator==(const Operation& op) const { 152 return !(*this != op); 153 } 154 class Constant; 155 class Variable; 156 class Custom; 157 class Add; 158 class Subtract; 159 class Multiply; 160 class Divide; 161 class Power; 162 class Negate; 163 class Sqrt; 164 class Exp; 165 class Log; 166 class Sin; 167 class Cos; 168 class Sec; 169 class Csc; 170 class Tan; 171 class Cot; 172 class Asin; 173 class Acos; 174 class Atan; 175 class Sinh; 176 class Cosh; 177 class Tanh; 178 class Erf; 179 class Erfc; 180 class Step; 181 class Delta; 182 class Nandelta; 183 class Square; 184 class Cube; 185 class Reciprocal; 186 class AddConstant; 187 class MultiplyConstant; 188 class PowerConstant; 189 class Min; 190 class Max; 191 class Abs; 192 class Floor; 193 class Ceil; 194 class Select; 195 class Acot; 196 class Asec; 197 class Acsc; 198 class Coth; 199 class Sech; 200 class Csch; 201 class Asinh; 202 class Acosh; 203 class Atanh; 204 class Acoth; 205 class Asech; 206 class Acsch; 207 class Atan2; 208 }; 209 210 class LEPTON_EXPORT Operation::Constant : public Operation { 211 public: Constant(double value)212 Constant(double value) : value(value) { 213 } getName()214 std::string getName() const { 215 std::stringstream name; 216 name << value; 217 return name.str(); 218 } getId()219 Id getId() const { 220 return CONSTANT; 221 } getNumArguments()222 int getNumArguments() const { 223 return 0; 224 } clone()225 Operation* clone() const { 226 return new Constant(value); 227 } evaluate(double * args,const std::map<std::string,double> & variables)228 double evaluate(double* args, const std::map<std::string, double>& variables) const { 229 return value; 230 } 231 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; getValue()232 double getValue() const { 233 return value; 234 } 235 bool operator!=(const Operation& op) const { 236 const Constant* o = dynamic_cast<const Constant*>(&op); 237 return (o == NULL || o->value != value); 238 } 239 private: 240 double value; 241 }; 242 243 class LEPTON_EXPORT Operation::Variable : public Operation { 244 public: Variable(const std::string & name)245 Variable(const std::string& name) : name(name) { 246 } getName()247 std::string getName() const { 248 return name; 249 } getId()250 Id getId() const { 251 return VARIABLE; 252 } getNumArguments()253 int getNumArguments() const { 254 return 0; 255 } clone()256 Operation* clone() const { 257 return new Variable(name); 258 } evaluate(double * args,const std::map<std::string,double> & variables)259 double evaluate(double* args, const std::map<std::string, double>& variables) const { 260 std::map<std::string, double>::const_iterator iter = variables.find(name); 261 if (iter == variables.end()) 262 throw Exception("No value specified for variable "+name); 263 return iter->second; 264 } 265 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 266 bool operator!=(const Operation& op) const { 267 const Variable* o = dynamic_cast<const Variable*>(&op); 268 return (o == NULL || o->name != name); 269 } 270 private: 271 std::string name; 272 }; 273 274 class LEPTON_EXPORT Operation::Custom : public Operation { 275 public: Custom(const std::string & name,CustomFunction * function)276 Custom(const std::string& name, CustomFunction* function) : name(name), function(function), isDerivative(false), derivOrder(function->getNumArguments(), 0) { 277 } Custom(const Custom & base,int derivIndex)278 Custom(const Custom& base, int derivIndex) : name(base.name), function(base.function->clone()), isDerivative(true), derivOrder(base.derivOrder) { 279 derivOrder[derivIndex]++; 280 } ~Custom()281 ~Custom() { 282 delete function; 283 } getName()284 std::string getName() const { 285 return name; 286 } getId()287 Id getId() const { 288 return CUSTOM; 289 } getNumArguments()290 int getNumArguments() const { 291 return function->getNumArguments(); 292 } clone()293 Operation* clone() const { 294 Custom* clone = new Custom(name, function->clone()); 295 clone->isDerivative = isDerivative; 296 clone->derivOrder = derivOrder; 297 return clone; 298 } evaluate(double * args,const std::map<std::string,double> & variables)299 double evaluate(double* args, const std::map<std::string, double>& variables) const { 300 if (isDerivative) 301 return function->evaluateDerivative(args, &derivOrder[0]); 302 return function->evaluate(args); 303 } 304 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; getDerivOrder()305 const std::vector<int>& getDerivOrder() const { 306 return derivOrder; 307 } 308 bool operator!=(const Operation& op) const { 309 const Custom* o = dynamic_cast<const Custom*>(&op); 310 return (o == NULL || o->name != name || o->isDerivative != isDerivative || o->derivOrder != derivOrder); 311 } 312 private: 313 std::string name; 314 CustomFunction* function; 315 bool isDerivative; 316 std::vector<int> derivOrder; 317 }; 318 319 class LEPTON_EXPORT Operation::Add : public Operation { 320 public: Add()321 Add() { 322 } getName()323 std::string getName() const { 324 return "+"; 325 } getId()326 Id getId() const { 327 return ADD; 328 } getNumArguments()329 int getNumArguments() const { 330 return 2; 331 } clone()332 Operation* clone() const { 333 return new Add(); 334 } evaluate(double * args,const std::map<std::string,double> & variables)335 double evaluate(double* args, const std::map<std::string, double>& variables) const { 336 return args[0]+args[1]; 337 } 338 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; isInfixOperator()339 bool isInfixOperator() const { 340 return true; 341 } isSymmetric()342 bool isSymmetric() const { 343 return true; 344 } 345 }; 346 347 class LEPTON_EXPORT Operation::Subtract : public Operation { 348 public: Subtract()349 Subtract() { 350 } getName()351 std::string getName() const { 352 return "-"; 353 } getId()354 Id getId() const { 355 return SUBTRACT; 356 } getNumArguments()357 int getNumArguments() const { 358 return 2; 359 } clone()360 Operation* clone() const { 361 return new Subtract(); 362 } evaluate(double * args,const std::map<std::string,double> & variables)363 double evaluate(double* args, const std::map<std::string, double>& variables) const { 364 return args[0]-args[1]; 365 } 366 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; isInfixOperator()367 bool isInfixOperator() const { 368 return true; 369 } 370 }; 371 372 class LEPTON_EXPORT Operation::Multiply : public Operation { 373 public: Multiply()374 Multiply() { 375 } getName()376 std::string getName() const { 377 return "*"; 378 } getId()379 Id getId() const { 380 return MULTIPLY; 381 } getNumArguments()382 int getNumArguments() const { 383 return 2; 384 } clone()385 Operation* clone() const { 386 return new Multiply(); 387 } evaluate(double * args,const std::map<std::string,double> & variables)388 double evaluate(double* args, const std::map<std::string, double>& variables) const { 389 return args[0]*args[1]; 390 } 391 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; isInfixOperator()392 bool isInfixOperator() const { 393 return true; 394 } isSymmetric()395 bool isSymmetric() const { 396 return true; 397 } 398 }; 399 400 class LEPTON_EXPORT Operation::Divide : public Operation { 401 public: Divide()402 Divide() { 403 } getName()404 std::string getName() const { 405 return "/"; 406 } getId()407 Id getId() const { 408 return DIVIDE; 409 } getNumArguments()410 int getNumArguments() const { 411 return 2; 412 } clone()413 Operation* clone() const { 414 return new Divide(); 415 } evaluate(double * args,const std::map<std::string,double> & variables)416 double evaluate(double* args, const std::map<std::string, double>& variables) const { 417 return args[0]/args[1]; 418 } 419 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; isInfixOperator()420 bool isInfixOperator() const { 421 return true; 422 } 423 }; 424 425 class LEPTON_EXPORT Operation::Power : public Operation { 426 public: Power()427 Power() { 428 } getName()429 std::string getName() const { 430 return "^"; 431 } getId()432 Id getId() const { 433 return POWER; 434 } getNumArguments()435 int getNumArguments() const { 436 return 2; 437 } clone()438 Operation* clone() const { 439 return new Power(); 440 } evaluate(double * args,const std::map<std::string,double> & variables)441 double evaluate(double* args, const std::map<std::string, double>& variables) const { 442 return std::pow(args[0], args[1]); 443 } 444 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; isInfixOperator()445 bool isInfixOperator() const { 446 return true; 447 } 448 }; 449 450 class LEPTON_EXPORT Operation::Negate : public Operation { 451 public: Negate()452 Negate() { 453 } getName()454 std::string getName() const { 455 return "-"; 456 } getId()457 Id getId() const { 458 return NEGATE; 459 } getNumArguments()460 int getNumArguments() const { 461 return 1; 462 } clone()463 Operation* clone() const { 464 return new Negate(); 465 } evaluate(double * args,const std::map<std::string,double> & variables)466 double evaluate(double* args, const std::map<std::string, double>& variables) const { 467 return -args[0]; 468 } 469 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 470 }; 471 472 class LEPTON_EXPORT Operation::Sqrt : public Operation { 473 public: Sqrt()474 Sqrt() { 475 } getName()476 std::string getName() const { 477 return "sqrt"; 478 } getId()479 Id getId() const { 480 return SQRT; 481 } getNumArguments()482 int getNumArguments() const { 483 return 1; 484 } clone()485 Operation* clone() const { 486 return new Sqrt(); 487 } evaluate(double * args,const std::map<std::string,double> & variables)488 double evaluate(double* args, const std::map<std::string, double>& variables) const { 489 return std::sqrt(args[0]); 490 } 491 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 492 }; 493 494 class LEPTON_EXPORT Operation::Exp : public Operation { 495 public: Exp()496 Exp() { 497 } getName()498 std::string getName() const { 499 return "exp"; 500 } getId()501 Id getId() const { 502 return EXP; 503 } getNumArguments()504 int getNumArguments() const { 505 return 1; 506 } clone()507 Operation* clone() const { 508 return new Exp(); 509 } evaluate(double * args,const std::map<std::string,double> & variables)510 double evaluate(double* args, const std::map<std::string, double>& variables) const { 511 return std::exp(args[0]); 512 } 513 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 514 }; 515 516 class LEPTON_EXPORT Operation::Log : public Operation { 517 public: Log()518 Log() { 519 } getName()520 std::string getName() const { 521 return "log"; 522 } getId()523 Id getId() const { 524 return LOG; 525 } getNumArguments()526 int getNumArguments() const { 527 return 1; 528 } clone()529 Operation* clone() const { 530 return new Log(); 531 } evaluate(double * args,const std::map<std::string,double> & variables)532 double evaluate(double* args, const std::map<std::string, double>& variables) const { 533 return std::log(args[0]); 534 } 535 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 536 }; 537 538 class LEPTON_EXPORT Operation::Sin : public Operation { 539 public: Sin()540 Sin() { 541 } getName()542 std::string getName() const { 543 return "sin"; 544 } getId()545 Id getId() const { 546 return SIN; 547 } getNumArguments()548 int getNumArguments() const { 549 return 1; 550 } clone()551 Operation* clone() const { 552 return new Sin(); 553 } evaluate(double * args,const std::map<std::string,double> & variables)554 double evaluate(double* args, const std::map<std::string, double>& variables) const { 555 return std::sin(args[0]); 556 } 557 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 558 }; 559 560 class LEPTON_EXPORT Operation::Cos : public Operation { 561 public: Cos()562 Cos() { 563 } getName()564 std::string getName() const { 565 return "cos"; 566 } getId()567 Id getId() const { 568 return COS; 569 } getNumArguments()570 int getNumArguments() const { 571 return 1; 572 } clone()573 Operation* clone() const { 574 return new Cos(); 575 } evaluate(double * args,const std::map<std::string,double> & variables)576 double evaluate(double* args, const std::map<std::string, double>& variables) const { 577 return std::cos(args[0]); 578 } 579 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 580 }; 581 582 class LEPTON_EXPORT Operation::Sec : public Operation { 583 public: Sec()584 Sec() { 585 } getName()586 std::string getName() const { 587 return "sec"; 588 } getId()589 Id getId() const { 590 return SEC; 591 } getNumArguments()592 int getNumArguments() const { 593 return 1; 594 } clone()595 Operation* clone() const { 596 return new Sec(); 597 } evaluate(double * args,const std::map<std::string,double> & variables)598 double evaluate(double* args, const std::map<std::string, double>& variables) const { 599 return 1.0/std::cos(args[0]); 600 } 601 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 602 }; 603 604 class LEPTON_EXPORT Operation::Csc : public Operation { 605 public: Csc()606 Csc() { 607 } getName()608 std::string getName() const { 609 return "csc"; 610 } getId()611 Id getId() const { 612 return CSC; 613 } getNumArguments()614 int getNumArguments() const { 615 return 1; 616 } clone()617 Operation* clone() const { 618 return new Csc(); 619 } evaluate(double * args,const std::map<std::string,double> & variables)620 double evaluate(double* args, const std::map<std::string, double>& variables) const { 621 return 1.0/std::sin(args[0]); 622 } 623 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 624 }; 625 626 class LEPTON_EXPORT Operation::Tan : public Operation { 627 public: Tan()628 Tan() { 629 } getName()630 std::string getName() const { 631 return "tan"; 632 } getId()633 Id getId() const { 634 return TAN; 635 } getNumArguments()636 int getNumArguments() const { 637 return 1; 638 } clone()639 Operation* clone() const { 640 return new Tan(); 641 } evaluate(double * args,const std::map<std::string,double> & variables)642 double evaluate(double* args, const std::map<std::string, double>& variables) const { 643 return std::tan(args[0]); 644 } 645 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 646 }; 647 648 class LEPTON_EXPORT Operation::Cot : public Operation { 649 public: Cot()650 Cot() { 651 } getName()652 std::string getName() const { 653 return "cot"; 654 } getId()655 Id getId() const { 656 return COT; 657 } getNumArguments()658 int getNumArguments() const { 659 return 1; 660 } clone()661 Operation* clone() const { 662 return new Cot(); 663 } evaluate(double * args,const std::map<std::string,double> & variables)664 double evaluate(double* args, const std::map<std::string, double>& variables) const { 665 return 1.0/std::tan(args[0]); 666 } 667 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 668 }; 669 670 class LEPTON_EXPORT Operation::Asin : public Operation { 671 public: Asin()672 Asin() { 673 } getName()674 std::string getName() const { 675 return "asin"; 676 } getId()677 Id getId() const { 678 return ASIN; 679 } getNumArguments()680 int getNumArguments() const { 681 return 1; 682 } clone()683 Operation* clone() const { 684 return new Asin(); 685 } evaluate(double * args,const std::map<std::string,double> & variables)686 double evaluate(double* args, const std::map<std::string, double>& variables) const { 687 return std::asin(args[0]); 688 } 689 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 690 }; 691 692 class LEPTON_EXPORT Operation::Acos : public Operation { 693 public: Acos()694 Acos() { 695 } getName()696 std::string getName() const { 697 return "acos"; 698 } getId()699 Id getId() const { 700 return ACOS; 701 } getNumArguments()702 int getNumArguments() const { 703 return 1; 704 } clone()705 Operation* clone() const { 706 return new Acos(); 707 } evaluate(double * args,const std::map<std::string,double> & variables)708 double evaluate(double* args, const std::map<std::string, double>& variables) const { 709 return std::acos(args[0]); 710 } 711 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 712 }; 713 714 class LEPTON_EXPORT Operation::Atan : public Operation { 715 public: Atan()716 Atan() { 717 } getName()718 std::string getName() const { 719 return "atan"; 720 } getId()721 Id getId() const { 722 return ATAN; 723 } getNumArguments()724 int getNumArguments() const { 725 return 1; 726 } clone()727 Operation* clone() const { 728 return new Atan(); 729 } evaluate(double * args,const std::map<std::string,double> & variables)730 double evaluate(double* args, const std::map<std::string, double>& variables) const { 731 return std::atan(args[0]); 732 } 733 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 734 }; 735 736 class LEPTON_EXPORT Operation::Sinh : public Operation { 737 public: Sinh()738 Sinh() { 739 } getName()740 std::string getName() const { 741 return "sinh"; 742 } getId()743 Id getId() const { 744 return SINH; 745 } getNumArguments()746 int getNumArguments() const { 747 return 1; 748 } clone()749 Operation* clone() const { 750 return new Sinh(); 751 } evaluate(double * args,const std::map<std::string,double> & variables)752 double evaluate(double* args, const std::map<std::string, double>& variables) const { 753 return std::sinh(args[0]); 754 } 755 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 756 }; 757 758 class LEPTON_EXPORT Operation::Cosh : public Operation { 759 public: Cosh()760 Cosh() { 761 } getName()762 std::string getName() const { 763 return "cosh"; 764 } getId()765 Id getId() const { 766 return COSH; 767 } getNumArguments()768 int getNumArguments() const { 769 return 1; 770 } clone()771 Operation* clone() const { 772 return new Cosh(); 773 } evaluate(double * args,const std::map<std::string,double> & variables)774 double evaluate(double* args, const std::map<std::string, double>& variables) const { 775 return std::cosh(args[0]); 776 } 777 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 778 }; 779 780 class LEPTON_EXPORT Operation::Tanh : public Operation { 781 public: Tanh()782 Tanh() { 783 } getName()784 std::string getName() const { 785 return "tanh"; 786 } getId()787 Id getId() const { 788 return TANH; 789 } getNumArguments()790 int getNumArguments() const { 791 return 1; 792 } clone()793 Operation* clone() const { 794 return new Tanh(); 795 } evaluate(double * args,const std::map<std::string,double> & variables)796 double evaluate(double* args, const std::map<std::string, double>& variables) const { 797 return std::tanh(args[0]); 798 } 799 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 800 }; 801 802 class LEPTON_EXPORT Operation::Erf : public Operation { 803 public: Erf()804 Erf() { 805 } getName()806 std::string getName() const { 807 return "erf"; 808 } getId()809 Id getId() const { 810 return ERF; 811 } getNumArguments()812 int getNumArguments() const { 813 return 1; 814 } clone()815 Operation* clone() const { 816 return new Erf(); 817 } 818 double evaluate(double* args, const std::map<std::string, double>& variables) const; 819 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 820 }; 821 822 class LEPTON_EXPORT Operation::Erfc : public Operation { 823 public: Erfc()824 Erfc() { 825 } getName()826 std::string getName() const { 827 return "erfc"; 828 } getId()829 Id getId() const { 830 return ERFC; 831 } getNumArguments()832 int getNumArguments() const { 833 return 1; 834 } clone()835 Operation* clone() const { 836 return new Erfc(); 837 } 838 double evaluate(double* args, const std::map<std::string, double>& variables) const; 839 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 840 }; 841 842 class LEPTON_EXPORT Operation::Step : public Operation { 843 public: Step()844 Step() { 845 } getName()846 std::string getName() const { 847 return "step"; 848 } getId()849 Id getId() const { 850 return STEP; 851 } getNumArguments()852 int getNumArguments() const { 853 return 1; 854 } clone()855 Operation* clone() const { 856 return new Step(); 857 } evaluate(double * args,const std::map<std::string,double> & variables)858 double evaluate(double* args, const std::map<std::string, double>& variables) const { 859 return (args[0] >= 0.0 ? 1.0 : 0.0); 860 } 861 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 862 }; 863 864 class LEPTON_EXPORT Operation::Delta : public Operation { 865 public: Delta()866 Delta() { 867 } getName()868 std::string getName() const { 869 return "delta"; 870 } getId()871 Id getId() const { 872 return DELTA; 873 } getNumArguments()874 int getNumArguments() const { 875 return 1; 876 } clone()877 Operation* clone() const { 878 return new Delta(); 879 } evaluate(double * args,const std::map<std::string,double> & variables)880 double evaluate(double* args, const std::map<std::string, double>& variables) const { 881 return (args[0] == 0.0 ? 1.0/0.0 : 0.0); 882 } 883 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 884 }; 885 886 class LEPTON_EXPORT Operation::Nandelta : public Operation { 887 public: Nandelta()888 Nandelta() { 889 } getName()890 std::string getName() const { 891 return "nandelta"; 892 } getId()893 Id getId() const { 894 return NANDELTA; 895 } getNumArguments()896 int getNumArguments() const { 897 return 1; 898 } clone()899 Operation* clone() const { 900 return new Nandelta(); 901 } evaluate(double * args,const std::map<std::string,double> & variables)902 double evaluate(double* args, const std::map<std::string, double>& variables) const { 903 return (args[0] == 0.0 ? std::numeric_limits<double>::quiet_NaN() : 0.0); 904 } 905 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 906 }; 907 908 class LEPTON_EXPORT Operation::Square : public Operation { 909 public: Square()910 Square() { 911 } getName()912 std::string getName() const { 913 return "square"; 914 } getId()915 Id getId() const { 916 return SQUARE; 917 } getNumArguments()918 int getNumArguments() const { 919 return 1; 920 } clone()921 Operation* clone() const { 922 return new Square(); 923 } evaluate(double * args,const std::map<std::string,double> & variables)924 double evaluate(double* args, const std::map<std::string, double>& variables) const { 925 return args[0]*args[0]; 926 } 927 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 928 }; 929 930 class LEPTON_EXPORT Operation::Cube : public Operation { 931 public: Cube()932 Cube() { 933 } getName()934 std::string getName() const { 935 return "cube"; 936 } getId()937 Id getId() const { 938 return CUBE; 939 } getNumArguments()940 int getNumArguments() const { 941 return 1; 942 } clone()943 Operation* clone() const { 944 return new Cube(); 945 } evaluate(double * args,const std::map<std::string,double> & variables)946 double evaluate(double* args, const std::map<std::string, double>& variables) const { 947 return args[0]*args[0]*args[0]; 948 } 949 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 950 }; 951 952 class LEPTON_EXPORT Operation::Reciprocal : public Operation { 953 public: Reciprocal()954 Reciprocal() { 955 } getName()956 std::string getName() const { 957 return "recip"; 958 } getId()959 Id getId() const { 960 return RECIPROCAL; 961 } getNumArguments()962 int getNumArguments() const { 963 return 1; 964 } clone()965 Operation* clone() const { 966 return new Reciprocal(); 967 } evaluate(double * args,const std::map<std::string,double> & variables)968 double evaluate(double* args, const std::map<std::string, double>& variables) const { 969 return 1.0/args[0]; 970 } 971 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 972 }; 973 974 class LEPTON_EXPORT Operation::AddConstant : public Operation { 975 public: AddConstant(double value)976 AddConstant(double value) : value(value) { 977 } getName()978 std::string getName() const { 979 std::stringstream name; 980 name << value << "+"; 981 return name.str(); 982 } getId()983 Id getId() const { 984 return ADD_CONSTANT; 985 } getNumArguments()986 int getNumArguments() const { 987 return 1; 988 } clone()989 Operation* clone() const { 990 return new AddConstant(value); 991 } evaluate(double * args,const std::map<std::string,double> & variables)992 double evaluate(double* args, const std::map<std::string, double>& variables) const { 993 return args[0]+value; 994 } 995 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; getValue()996 double getValue() const { 997 return value; 998 } 999 bool operator!=(const Operation& op) const { 1000 const AddConstant* o = dynamic_cast<const AddConstant*>(&op); 1001 return (o == NULL || o->value != value); 1002 } 1003 private: 1004 double value; 1005 }; 1006 1007 class LEPTON_EXPORT Operation::MultiplyConstant : public Operation { 1008 public: MultiplyConstant(double value)1009 MultiplyConstant(double value) : value(value) { 1010 } getName()1011 std::string getName() const { 1012 std::stringstream name; 1013 name << value << "*"; 1014 return name.str(); 1015 } getId()1016 Id getId() const { 1017 return MULTIPLY_CONSTANT; 1018 } getNumArguments()1019 int getNumArguments() const { 1020 return 1; 1021 } clone()1022 Operation* clone() const { 1023 return new MultiplyConstant(value); 1024 } evaluate(double * args,const std::map<std::string,double> & variables)1025 double evaluate(double* args, const std::map<std::string, double>& variables) const { 1026 return args[0]*value; 1027 } 1028 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; getValue()1029 double getValue() const { 1030 return value; 1031 } 1032 bool operator!=(const Operation& op) const { 1033 const MultiplyConstant* o = dynamic_cast<const MultiplyConstant*>(&op); 1034 return (o == NULL || o->value != value); 1035 } 1036 private: 1037 double value; 1038 }; 1039 1040 class LEPTON_EXPORT Operation::PowerConstant : public Operation { 1041 public: PowerConstant(double value)1042 PowerConstant(double value) : value(value) { 1043 intValue = (int) value; 1044 isIntPower = (intValue == value); 1045 } getName()1046 std::string getName() const { 1047 std::stringstream name; 1048 name << "^" << value; 1049 return name.str(); 1050 } getId()1051 Id getId() const { 1052 return POWER_CONSTANT; 1053 } getNumArguments()1054 int getNumArguments() const { 1055 return 1; 1056 } clone()1057 Operation* clone() const { 1058 return new PowerConstant(value); 1059 } evaluate(double * args,const std::map<std::string,double> & variables)1060 double evaluate(double* args, const std::map<std::string, double>& variables) const { 1061 if (isIntPower) { 1062 // Integer powers can be computed much more quickly by repeated multiplication. 1063 1064 int exponent = intValue; 1065 double base = args[0]; 1066 if (exponent < 0) { 1067 exponent = -exponent; 1068 base = 1.0/base; 1069 } 1070 double result = 1.0; 1071 while (exponent != 0) { 1072 if ((exponent&1) == 1) 1073 result *= base; 1074 base *= base; 1075 exponent = exponent>>1; 1076 } 1077 return result; 1078 } 1079 else 1080 return std::pow(args[0], value); 1081 } 1082 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; getValue()1083 double getValue() const { 1084 return value; 1085 } 1086 bool operator!=(const Operation& op) const { 1087 const PowerConstant* o = dynamic_cast<const PowerConstant*>(&op); 1088 return (o == NULL || o->value != value); 1089 } isInfixOperator()1090 bool isInfixOperator() const { 1091 return true; 1092 } 1093 private: 1094 double value; 1095 int intValue; 1096 bool isIntPower; 1097 }; 1098 1099 class LEPTON_EXPORT Operation::Min : public Operation { 1100 public: Min()1101 Min() { 1102 } getName()1103 std::string getName() const { 1104 return "min"; 1105 } getId()1106 Id getId() const { 1107 return MIN; 1108 } getNumArguments()1109 int getNumArguments() const { 1110 return 2; 1111 } clone()1112 Operation* clone() const { 1113 return new Min(); 1114 } evaluate(double * args,const std::map<std::string,double> & variables)1115 double evaluate(double* args, const std::map<std::string, double>& variables) const { 1116 // parens around (std::min) are workaround for horrible microsoft max/min macro trouble 1117 return (std::min)(args[0], args[1]); 1118 } 1119 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 1120 }; 1121 1122 class LEPTON_EXPORT Operation::Max : public Operation { 1123 public: Max()1124 Max() { 1125 } getName()1126 std::string getName() const { 1127 return "max"; 1128 } getId()1129 Id getId() const { 1130 return MAX; 1131 } getNumArguments()1132 int getNumArguments() const { 1133 return 2; 1134 } clone()1135 Operation* clone() const { 1136 return new Max(); 1137 } evaluate(double * args,const std::map<std::string,double> & variables)1138 double evaluate(double* args, const std::map<std::string, double>& variables) const { 1139 // parens around (std::min) are workaround for horrible microsoft max/min macro trouble 1140 return (std::max)(args[0], args[1]); 1141 } 1142 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 1143 }; 1144 1145 class LEPTON_EXPORT Operation::Abs : public Operation { 1146 public: Abs()1147 Abs() { 1148 } getName()1149 std::string getName() const { 1150 return "abs"; 1151 } getId()1152 Id getId() const { 1153 return ABS; 1154 } getNumArguments()1155 int getNumArguments() const { 1156 return 1; 1157 } clone()1158 Operation* clone() const { 1159 return new Abs(); 1160 } evaluate(double * args,const std::map<std::string,double> & variables)1161 double evaluate(double* args, const std::map<std::string, double>& variables) const { 1162 return std::abs(args[0]); 1163 } 1164 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 1165 }; 1166 1167 class LEPTON_EXPORT Operation::Floor : public Operation { 1168 public: 1169 Floor()1170 Floor() { 1171 } getName()1172 std::string getName() const { 1173 return "floor"; 1174 } getId()1175 Id getId() const { 1176 return FLOOR; 1177 } getNumArguments()1178 int getNumArguments() const { 1179 return 1; 1180 } clone()1181 Operation* clone() const { 1182 return new Floor(); 1183 } evaluate(double * args,const std::map<std::string,double> & variables)1184 double evaluate(double* args, const std::map<std::string, double>& variables) const { 1185 return std::floor(args[0]); 1186 } 1187 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 1188 }; 1189 1190 class LEPTON_EXPORT Operation::Ceil : public Operation { 1191 public: Ceil()1192 Ceil() { 1193 } getName()1194 std::string getName() const { 1195 return "ceil"; 1196 } getId()1197 Id getId() const { 1198 return CEIL; 1199 } getNumArguments()1200 int getNumArguments() const { 1201 return 1; 1202 } clone()1203 Operation* clone() const { 1204 return new Ceil(); 1205 } evaluate(double * args,const std::map<std::string,double> & variables)1206 double evaluate(double* args, const std::map<std::string, double>& variables) const { 1207 return std::ceil(args[0]); 1208 } 1209 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 1210 }; 1211 1212 class LEPTON_EXPORT Operation::Select : public Operation { 1213 public: Select()1214 Select() { 1215 } getName()1216 std::string getName() const { 1217 return "select"; 1218 } getId()1219 Id getId() const { 1220 return SELECT; 1221 } getNumArguments()1222 int getNumArguments() const { 1223 return 3; 1224 } clone()1225 Operation* clone() const { 1226 return new Select(); 1227 } evaluate(double * args,const std::map<std::string,double> & variables)1228 double evaluate(double* args, const std::map<std::string, double>& variables) const { 1229 return (args[0] != 0.0 ? args[1] : args[2]); 1230 } 1231 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; 1232 }; 1233 1234 #define LEPTON_CLASS_OPERATION(Name,name,NAME,nargs,impl) \ 1235 class LEPTON_EXPORT Operation::Name : public Operation { \ 1236 public: \ 1237 Name() { \ 1238 } \ 1239 std::string getName() const { \ 1240 return #name; \ 1241 } \ 1242 Id getId() const { \ 1243 return NAME; \ 1244 } \ 1245 int getNumArguments() const { \ 1246 return nargs; \ 1247 } \ 1248 Operation* clone() const { \ 1249 return new Name(); \ 1250 } \ 1251 double evaluate(double* args, const std::map<std::string, double>& variables) const { \ 1252 return impl; \ 1253 } \ 1254 ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const; \ 1255 } 1256 1257 LEPTON_CLASS_OPERATION(Acot,acot,ACOT,1,std::atan(1.0/args[0])); 1258 LEPTON_CLASS_OPERATION(Asec,asec,ASEC,1,std::acos(1.0/args[0])); 1259 LEPTON_CLASS_OPERATION(Acsc,acsc,ACSC,1,std::asin(1.0/args[0])); 1260 LEPTON_CLASS_OPERATION(Coth,coth,ACOT,1,1.0/std::tanh(args[0])); 1261 LEPTON_CLASS_OPERATION(Sech,sech,SECH,1,1.0/std::cosh(args[0])); 1262 LEPTON_CLASS_OPERATION(Csch,csch,CSCH,1,1.0/std::sinh(args[0])); 1263 1264 LEPTON_CLASS_OPERATION(Asinh,asinh,ASINH,1,std::asinh(args[0])); 1265 LEPTON_CLASS_OPERATION(Acosh,acosh,ACOSH,1,std::acosh(args[0])); 1266 LEPTON_CLASS_OPERATION(Atanh,atanh,ATANH,1,std::atanh(args[0])); 1267 1268 LEPTON_CLASS_OPERATION(Acoth,acoth,ACOTH,1,0.5*std::log((args[0]+1.0)/(args[0]-1.0))); 1269 LEPTON_CLASS_OPERATION(Asech,asech,ASECH,1,std::log(std::sqrt(1.0/args[0]-1.0)*std::sqrt(1.0/args[0]+1.0)+1.0/args[0])); 1270 LEPTON_CLASS_OPERATION(Acsch,acsch,ACSCH,1,std::log(1.0/args[0]+std::sqrt(1.0/(args[0]*args[0])+1.0))); 1271 1272 LEPTON_CLASS_OPERATION(Atan2,atan2,ATAN2,2,std::atan2(args[0],args[1])); 1273 1274 } // namespace lepton 1275 } // namespace PLMD 1276 1277 #endif /*LEPTON_OPERATION_H_*/ 1278