Changeset 3669 for trunk/test_more
 Timestamp:
 Mar 8, 2015 1:43:45 PM (5 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/test_more/optimize.cpp
r3642 r3669 4 4 5 5 CppAD is distributed under multiple licenses. This distribution is under 6 the terms of the 6 the terms of the 7 7 Eclipse Public License Version 1.0. 8 8 … … 17 17 18 18 namespace { 19 // include conditional skip optimization 20 bool conditional_skip; 21 19 22 // accuracy for almost equal checks 20 23 double eps = 10. * std::numeric_limits<double>::epsilon(); … … 30 33 using CppAD::AD; 31 34 using CppAD::vector; 32 33 // independent variable vector 35 36 // independent variable vector 34 37 vector< AD<double> > ax(2), ay(1); 35 38 ax[0] = 1.0; 36 39 ax[1] = 2.0; 37 40 Independent(ax); 38 41 39 42 // first conditional expression 40 43 AD<double> ac1 = CondExpLe(ax[0], ax[1], 2.0 * ax[0], 3.0 * ax[1] ); 41 44 42 45 // second conditional expression 43 46 AD<double> ac2 = CondExpGe(ax[0], ax[1], 4.0 * ax[0], 5.0 * ax[1] ); 44 47 45 48 // third conditional expression 46 49 AD<double> ac3 = CondExpLt(ax[0], ax[1], 6.0 * ac1, 7.0 * ac2 ); 47 50 48 51 // create function object f : ax > ay 49 52 ay[0] = ac3; 50 53 CppAD::ADFun<double> f(ax, ay); 51 54 52 55 // now optimize the operation sequence 53 f.optimize(); 54 56 if( conditional_skip ) 57 f.optimize(); 58 else 59 f.optimize("no_conditional_skip"); 60 55 61 // now zero order forward 56 62 vector<double> x(2), y(1); … … 71 77 c2 = 4.0 * x[0]; 72 78 else c2 = 5.0 * x[1]; 73 79 74 80 // third conditional expression 75 81 double c3; … … 77 83 c3 = 6.0 * c1; 78 84 else c3 = 7.0 * c2; 79 85 80 86 ok &= y[0] == c3; 81 87 } … … 86 92 // the operands in the logical comparison because of the CondExp 87 93 // sparsity pattern. 88 void j_algo( 94 void j_algo( 89 95 const CppAD::vector< CppAD::AD<double> >& ax , 90 96 CppAD::vector< CppAD::AD<double> >& ay ) 91 97 { ay[0] = CondExpGt(ax[0], ax[1], ax[2], ax[3]); } 92 98 93 99 bool atomic_cond_exp_sparsity(void) 94 100 { bool ok = true; 95 101 using CppAD::AD; 96 102 using CppAD::vector; 97 103 98 104 // Create a checkpoint version of the function g 99 105 vector< AD<double> > au(4), av(1); … … 101 107 au[i] = AD<double>(i); 102 108 CppAD::checkpoint<double> j_check("j_check", j_algo, au, av); 103 104 // independent variable vector 109 110 // independent variable vector 105 111 vector< AD<double> > ax(2), ay(1); 106 112 ax[0] = 1.; 107 113 ax[1] = 1.; 108 114 Independent(ax); 109 115 110 116 // call atomic function that does not get used 111 for(size_t i = 0; i < 4; i++) 117 for(size_t i = 0; i < 4; i++) 112 118 au[i] = ax[0] + AD<double>(i + 1) * ax[1]; 113 119 j_check(au, ay); 114 120 115 121 // create function object f : ax > ay 116 122 CppAD::ADFun<double> f(ax, ay); 117 123 118 124 119 125 // now optimize the operation sequence 120 126 j_check.option( atomic_sparsity_option ); 121 f.optimize(); 122 127 if( conditional_skip ) 128 f.optimize(); 129 else 130 f.optimize("no_conditional_skip"); 131 123 132 // check result where true case is used; i.e., au[0] > au[1] 124 133 vector<double> x(2), y(1); … … 127 136 y = f.Forward(0, x); 128 137 ok &= y[0] == x[0] + double(3) * x[1]; 129 130 138 139 131 140 // check result where false case is used; i.e., au[0] <= au[1] 132 141 x[0] = 1.; … … 134 143 y = f.Forward(0, x); 135 144 ok &= y[0] == x[0] + double(4) * x[1]; 136 145 137 146 return ok; 138 147 } … … 160 169 CppAD::checkpoint<double> h_check("h_check", h_algo, ax, ah); 161 170 162 // independent variable vector 171 // independent variable vector 163 172 Independent(ax); 164 173 … … 168 177 169 178 // conditional expression 170 ay[0] = CondExpLt(ax[0], ax[1], ag[0], ah[0]); 171 179 ay[0] = CondExpLt(ax[0], ax[1], ag[0], ah[0]); 180 172 181 // create function object f : ax > ay 173 182 CppAD::ADFun<double> f; 174 183 f.Dependent(ax, ay); 175 184 176 185 // use zero order to evaluate when condition is true 177 186 CppAD::vector<double> x(2), dx(2); … … 188 197 189 198 // now optimize the operation sequence 190 f.optimize(); 199 if( conditional_skip ) 200 f.optimize(); 201 else 202 f.optimize("no_conditional_skip"); 191 203 192 204 // optimized zero order forward when condition is false … … 210 222 ok &= dx[0] == 1.; 211 223 ok &= dx[1] == 1.; 212 224 213 225 return ok; 214 226 } 215 227 //  216 228 // Test of optimizing out arguments to an atomic function 217 void g_algo( 229 void g_algo( 218 230 const CppAD::vector< CppAD::AD<double> >& ax , 219 231 CppAD::vector< CppAD::AD<double> >& ay ) … … 224 236 using CppAD::AD; 225 237 using CppAD::vector; 226 238 227 239 // Create a checkpoint version of the function g 228 240 vector< AD<double> > ax(2), ay(2), az(1); … … 230 242 ax[1] = 1.; 231 243 CppAD::checkpoint<double> g_check("g_check", g_algo, ax, ay); 232 233 // independent variable vector 244 245 // independent variable vector 234 246 Independent(ax); 235 247 236 248 // call atomic function that does not get used 237 249 g_check(ax, ay); 238 250 239 251 // conditional expression 240 az[0] = CondExpLt(ax[0], ax[1], ax[0] + ax[1], ax[0]  ax[1]); 241 252 az[0] = CondExpLt(ax[0], ax[1], ax[0] + ax[1], ax[0]  ax[1]); 253 242 254 // create function object f : ax > az 243 255 CppAD::ADFun<double> f(ax, az); … … 246 258 // (include ay[0] and ay[1]) 247 259 size_t n_before = f.size_var(); 248 260 249 261 // now optimize the operation sequence 250 262 g_check.option( atomic_sparsity_option ); 251 f.optimize(); 252 253 // number of variables after optimization 263 if( conditional_skip ) 264 f.optimize(); 265 else 266 f.optimize("no_conditional_skip"); 267 268 // number of variables after optimization 254 269 // (does not include ay[0] and ay[1]) 255 270 size_t n_after = f.size_var(); 256 271 ok &= n_after + 2 == n_before; 257 272 258 273 // check optimization works ok 259 274 vector<double> x(2), z(1); … … 262 277 z = f.Forward(0, x); 263 278 ok &= z[0] == x[0]  x[1]; 264 279 265 280 return ok; 266 281 } … … 284 299 g_check(au, aw); 285 300 286 // now create f(x) = x_0  x_1 301 // now create f(x) = x_0  x_1 287 302 ay[0] = aw[0]; 288 303 CppAD::ADFun<double> f(ax, ay); … … 290 305 // number of variables before optimization 291 306 size_t n_before = f.size_var(); 292 307 293 308 // now optimize f so that the calculation of au[1] is removed 294 309 g_check.option( atomic_sparsity_option ); 295 f.optimize(); 310 if( conditional_skip ) 311 f.optimize(); 312 else 313 f.optimize("no_conditional_skip"); 296 314 297 315 // check difference in number of variables … … 304 322 x[1] = 6.0; 305 323 y = f.Forward(0, x); 306 ok &= (y[0] == x[0] + x[1]); 324 ok &= (y[0] == x[0] + x[1]); 307 325 308 326 return ok; … … 324 342 // unary operator where operand is arg[0] 325 343 // (note that sin corresponds to two tape variables) 326 not_used = CppAD::abs(x[0]); 344 not_used = CppAD::abs(x[0]); 327 345 y[0] = sin(x[0]); 328 346 original += 3; … … 375 393 size_t opt; 376 394 size_t i, j; 377 395 378 396 // domain space vector 379 397 size_t n = 6; 380 398 CppAD::vector< AD<double> > X(n); 381 399 for(j = 0; j < n; j++) 382 X[j] = 1. / double(j + 1); 383 400 X[j] = 1. / double(j + 1); 401 384 402 // declare independent variables and start tape recording 385 403 CppAD::Independent(X); 386 387 // range space vector 404 405 // range space vector 388 406 size_t m = n; 389 407 CppAD::vector< AD<double> > Y(m); 390 408 depend_fun(X, Y, original, opt); 391 409 392 410 // create f: X > Y and stop tape recording 393 411 CppAD::ADFun<double> F; 394 F.Dependent(X, Y); 395 412 F.Dependent(X, Y); 413 396 414 CppAD::vector<double> x(n), y(m), check(m); 397 415 for(j = 0; j < n; j++) … … 401 419 for(i = 0; i < m; i++) 402 420 ok &= NearEqual(y[i], check[i], eps, eps); 403 421 404 422 // Check size before optimization 405 423 ok &= F.size_var() == original; 406 424 407 425 // Optimize the operation sequence 408 F.optimize(); 409 426 if( conditional_skip ) 427 F.optimize(); 428 else 429 F.optimize("no_conditional_skip"); 430 410 431 // Check size after optimization 411 432 ok &= F.size_var() == opt; 412 433 413 434 // check result now 414 435 // (should have already been checked if NDEBUG not defined) … … 416 437 for(i = 0; i < m; i++) 417 438 ok &= NearEqual(y[i], check[i], eps, eps); 418 439 419 440 return ok; 420 441 } … … 430 451 CppAD::vector< AD<double> > X(n); 431 452 for(j = 0; j < n; j++) 432 X[j] = double(j); 453 X[j] = double(j); 433 454 434 455 // range space vector … … 442 463 for(j = 0; j < n; j++) 443 464 V[j] = 0; 444 465 445 466 // declare independent variables and start tape recording 446 467 CppAD::Independent(X); … … 461 482 Y[i] = U[I]; 462 483 } 463 484 464 485 // create f: X > Y and stop tape recording 465 // Y[ X[0] ] = X[1] and other components of Y are zero. 486 // Y[ X[0] ] = X[1] and other components of Y are zero. 466 487 CppAD::ADFun<double> F; 467 F.Dependent(X, Y); 488 F.Dependent(X, Y); 468 489 469 490 // Check number of VecAD vectors plus number of VecAD elements 470 ok &= (F.size_VecAD() == 2 + n + m); 471 491 ok &= (F.size_VecAD() == 2 + n + m); 492 472 493 CppAD::vector<double> x(n), y(m); 473 494 for(j = 0; j < n; j++) … … 481 502 } 482 503 483 F.optimize(); 504 if( conditional_skip ) 505 F.optimize(); 506 else 507 F.optimize("no_conditional_skip"); 484 508 485 509 // Check number of VecAD vectors plus number of VecAD elements 486 ok &= (F.size_VecAD() == 1 + m); 510 ok &= (F.size_VecAD() == 1 + m); 487 511 y = F.Forward(0, x); 488 512 for(i = 0; i < m; i++) … … 491 515 else ok &= (y[i] == x[1]); 492 516 } 493 517 494 518 return ok; 495 519 } … … 502 526 size_t n = 3; 503 527 size_t j; 504 528 505 529 vector< AD<double> > X(n), Y(n); 506 vector<double> x(n), y(n); 507 530 vector<double> x(n), y(n); 531 508 532 for(j = 0; j < n; j++) 509 533 X[j] = x[j] = double(j+2); 510 534 511 535 CppAD::Independent(X); 512 536 513 537 Y[0] = pow(X[0], 2.0); 514 538 Y[1] = pow(2.0, X[1]); 515 539 Y[2] = pow(X[0], X[1]); 516 540 517 541 CppAD::ADFun<double> F(X, Y); 518 F.optimize(); 542 if( conditional_skip ) 543 F.optimize(); 544 else 545 F.optimize("no_conditional_skip"); 519 546 y = F.Forward(0, x); 520 547 … … 526 553 527 554 // check reverse mode derivative 528 vector<double> w(n), dw(n); 555 vector<double> w(n), dw(n); 529 556 w[0] = 0.; 530 557 w[1] = 0.; … … 540 567 check = 0.; 541 568 ok &= NearEqual( dw[2], check, eps, eps ); 542 569 543 570 return ok; 544 571 } … … 564 591 vector<double> y_original = F.Forward(0, x); 565 592 size_t size_original = F.size_var(); 566 F.optimize(); 593 if( conditional_skip ) 594 F.optimize(); 595 else 596 F.optimize("no_conditional_skip"); 567 597 ok &= F.size_var() + 5 == size_original; 568 598 vector<double> y = F.Forward(0, x); … … 582 612 583 613 // unary operator where operand is arg[0] and one result 584 Scalar a1 = CppAD::exp(x[0]); 614 Scalar a1 = CppAD::exp(x[0]); 585 615 original += 1; 586 616 opt += 1; … … 651 681 size_t opt; 652 682 size_t i, j; 653 683 654 684 // domain space vector 655 685 size_t n = 7; 656 686 CppAD::vector< AD<double> > X(n); 657 687 for(j = 0; j < n; j++) 658 X[j] = 1. / double(j + 1); 659 688 X[j] = 1. / double(j + 1); 689 660 690 // declare independent variables and start tape recording 661 691 CppAD::Independent(X); 662 663 // range space vector 692 693 // range space vector 664 694 size_t m = n; 665 695 CppAD::vector< AD<double> > Y(m); 666 696 duplicate_fun(X, Y, original, opt); 667 697 668 698 // create f: X > Y and stop tape recording 669 699 CppAD::ADFun<double> F; 670 F.Dependent(X, Y); 671 700 F.Dependent(X, Y); 701 672 702 CppAD::vector<double> x(n), y(m), check(m); 673 703 for(j = 0; j < n; j++) … … 677 707 for(i = 0; i < m; i++) 678 708 ok &= NearEqual(y[i], check[i], eps, eps); 679 709 680 710 // Check size before optimization 681 711 ok &= F.size_var() == (n + 1 + original); 682 712 683 713 // Optimize the operation sequence 684 F.optimize(); 685 714 if( conditional_skip ) 715 F.optimize(); 716 else 717 F.optimize("no_conditional_skip"); 718 686 719 // Check size after optimization 687 720 ok &= F.size_var() == (n + 1 + opt); 688 721 689 722 // check result now 690 723 // (should have already been checked if NDEBUG not defined) … … 692 725 for(i = 0; i < m; i++) 693 726 ok &= NearEqual(y[i], check[i], eps, eps); 694 727 695 728 return ok; 696 729 } … … 707 740 CppAD::vector< AD<double> > X(n); 708 741 for(j = 0; j < n; j++) 709 X[j] = double(j + 2); 742 X[j] = double(j + 2); 710 743 711 744 // range space vector … … 727 760 // create a duplicate that can only be dectected using new 728 761 // argument indices 729 AD<double> B2 = A2 / 2.; 730 731 // Make a new variable for result 762 AD<double> B2 = A2 / 2.; 763 764 // Make a new variable for result 732 765 // and make it depend on all the variables 733 766 Y[0] = B1 + B2; … … 735 768 // create f: X > Y and stop tape recording 736 769 CppAD::ADFun<double> F; 737 F.Dependent(X, Y); 770 F.Dependent(X, Y); 738 771 739 772 // check number of variables in original function 740 ok &= (F.size_var() == 1 + n + m + 4 ); 741 773 ok &= (F.size_var() == 1 + n + m + 4 ); 774 742 775 CppAD::vector<double> x(n), y(m); 743 776 for(j = 0; j < n; j++) … … 748 781 ok &= ( y[i] == Value( Y[i] ) ); 749 782 750 F.optimize(); 783 if( conditional_skip ) 784 F.optimize(); 785 else 786 F.optimize("no_conditional_skip"); 751 787 752 788 // check number of variables in optimized version 753 ok &= (F.size_var() == 1 + n + m + 2 ); 789 ok &= (F.size_var() == 1 + n + m + 2 ); 754 790 755 791 y = F.Forward(0, x); … … 771 807 CppAD::vector< AD<double> > X(n); 772 808 for(j = 0; j < n; j++) 773 X[j] = double(j + 2); 809 X[j] = double(j + 2); 774 810 775 811 // range space vector … … 791 827 // create a duplicate that can only be dectected using new 792 828 // argument indices 793 AD<double> B2 = 2. * A2; 794 795 // Make a new variable for result 829 AD<double> B2 = 2. * A2; 830 831 // Make a new variable for result 796 832 // and make it depend on all the variables 797 833 Y[0] = B1 + B2; … … 799 835 // create f: X > Y and stop tape recording 800 836 CppAD::ADFun<double> F; 801 F.Dependent(X, Y); 837 F.Dependent(X, Y); 802 838 803 839 // check number of variables in original function 804 ok &= (F.size_var() == 1 + n + m + 4 ); 805 840 ok &= (F.size_var() == 1 + n + m + 4 ); 841 806 842 CppAD::vector<double> x(n), y(m); 807 843 for(j = 0; j < n; j++) … … 812 848 ok &= ( y[i] == Value( Y[i] ) ); 813 849 814 F.optimize(); 850 if( conditional_skip ) 851 F.optimize(); 852 else 853 F.optimize("no_conditional_skip"); 815 854 816 855 // check number of variables in optimized version 817 ok &= (F.size_var() == 1 + n + m + 2 ); 856 ok &= (F.size_var() == 1 + n + m + 2 ); 818 857 819 858 y = F.Forward(0, x); … … 843 882 CppAD::Independent(X); 844 883 845 // check a huge number of same operation with different operands 884 // check a huge number of same operation with different operands 846 885 size_t n_operations = std::min( 847 886 size_t(CPPAD_HASH_TABLE_SIZE) + 5, … … 854 893 // create f: X > Y and stop tape recording 855 894 CppAD::ADFun<double> F; 856 F.Dependent(X, Y); 895 F.Dependent(X, Y); 857 896 858 897 // check number of variables in original function 859 ok &= (F.size_var() == 1 + n + n_operations ); 860 898 ok &= (F.size_var() == 1 + n + n_operations ); 899 861 900 CppAD::vector<double> x(n), y(m); 862 901 x[0] = 1.; … … 865 904 ok &= ( y[0] == Value( Y[0] ) ); 866 905 867 F.optimize(); 906 if( conditional_skip ) 907 F.optimize(); 908 else 909 F.optimize("no_conditional_skip"); 868 910 869 911 // check same number of variables in optimized version 870 ok &= (F.size_var() == 1 + n + n_operations ); 912 ok &= (F.size_var() == 1 + n + n_operations ); 871 913 872 914 y = F.Forward(0, x); … … 887 929 CppAD::vector< AD<double> > X(n); 888 930 for(j = 0; j < n; j++) 889 X[j] = double(j + 2); 931 X[j] = double(j + 2); 890 932 891 933 size_t n_original = 1 + n; … … 915 957 916 958 CppAD::ADFun<double> F; 917 F.Dependent(X, Y); 959 F.Dependent(X, Y); 918 960 919 961 // check number of variables in original function 920 ok &= (F.size_var() == n_original ); 921 962 ok &= (F.size_var() == n_original ); 963 922 964 CppAD::vector<double> x(n), y(m); 923 965 for(j = 0; j < n; j++) … … 928 970 ok &= ( y[i] == Value( Y[i] ) ); 929 971 930 F.optimize(); 972 if( conditional_skip ) 973 F.optimize(); 974 else 975 F.optimize("no_conditional_skip"); 931 976 932 977 // check number of variables in optimized version 933 ok &= (F.size_var() == n_optimize ); 978 ok &= (F.size_var() == n_optimize ); 934 979 935 980 y = F.Forward(0, x); … … 942 987 bool forward_csum(void) 943 988 { bool ok = true; 944 989 945 990 using namespace CppAD; 946 947 // independent variable vector 991 992 // independent variable vector 948 993 CppAD::vector< AD<double> > X(2); 949 X[0] = 0.; 994 X[0] = 0.; 950 995 X[1] = 1.; 951 996 Independent(X); 952 997 953 998 // compute sum of elements in X 954 999 CppAD::vector< AD<double> > Y(1); 955 1000 Y[0] = X[0] + X[0] + X[1]; 956 1001 957 1002 // create function object F : X > Y 958 1003 ADFun<double> F(X, Y); 959 1004 960 1005 // now optimize the operation sequence 961 F.optimize(); 962 963 // use zero order to evaluate F[ (3, 4) ] 1006 if( conditional_skip ) 1007 F.optimize(); 1008 else 1009 F.optimize("no_conditional_skip"); 1010 1011 // use zero order to evaluate F[ (3, 4) ] 964 1012 CppAD::vector<double> x0( F.Domain() ); 965 1013 CppAD::vector<double> y0( F.Range() ); … … 968 1016 y0 = F.Forward(0, x0); 969 1017 ok &= NearEqual(y0[0] , x0[0]+x0[0]+x0[1], eps, eps); 970 1018 971 1019 // evaluate derivative of F in X[0] direction 972 1020 CppAD::vector<double> x1( F.Domain() ); … … 976 1024 y1 = F.Forward(1, x1); 977 1025 ok &= NearEqual(y1[0] , x1[0]+x1[0]+x1[1], eps, eps); 978 1026 979 1027 // evaluate second derivative of F in X[0] direction 980 1028 CppAD::vector<double> x2( F.Domain() ); … … 985 1033 double F_00 = 2. * y2[0]; 986 1034 ok &= NearEqual(F_00, 0., eps, eps); 987 1035 988 1036 return ok; 989 1037 } … … 991 1039 bool reverse_csum(void) 992 1040 { bool ok = true; 993 1041 994 1042 using namespace CppAD; 995 996 // independent variable vector 1043 1044 // independent variable vector 997 1045 CppAD::vector< AD<double> > X(2); 998 X[0] = 0.; 1046 X[0] = 0.; 999 1047 X[1] = 1.; 1000 1048 Independent(X); 1001 1049 1002 1050 // compute sum of elements in X 1003 1051 CppAD::vector< AD<double> > Y(1); 1004 1052 Y[0] = X[0]  (X[0]  X[1]); 1005 1053 1006 1054 // create function object F : X > Y 1007 1055 ADFun<double> F(X, Y); 1008 1056 1009 1057 // now optimize the operation sequence 1010 F.optimize(); 1011 1012 // use zero order to evaluate F[ (3, 4) ] 1058 if( conditional_skip ) 1059 F.optimize(); 1060 else 1061 F.optimize("no_conditional_skip"); 1062 1063 // use zero order to evaluate F[ (3, 4) ] 1013 1064 CppAD::vector<double> x0( F.Domain() ); 1014 1065 CppAD::vector<double> y0( F.Range() ); … … 1017 1068 y0 = F.Forward(0, x0); 1018 1069 ok &= NearEqual(y0[0] , x0[0]x0[0]+x0[1], eps, eps); 1019 1020 // evaluate derivative of F 1070 1071 // evaluate derivative of F 1021 1072 CppAD::vector<double> dF( F.Domain() ); 1022 1073 CppAD::vector<double> w( F.Range() ); … … 1025 1076 ok &= NearEqual(dF[0] , 0., eps, eps); 1026 1077 ok &= NearEqual(dF[1] , 1., eps, eps); 1027 1078 1028 1079 return ok; 1029 1080 } … … 1031 1082 { bool ok = true; 1032 1083 using namespace CppAD; 1033 1084 1034 1085 // dimension of the domain space 1035 size_t n = 3; 1036 1086 size_t n = 3; 1087 1037 1088 // dimension of the range space 1038 1089 size_t m = 3; 1039 1040 // independent variable vector 1090 1091 // independent variable vector 1041 1092 CppAD::vector< AD<double> > X(n); 1042 X[0] = 2.; 1093 X[0] = 2.; 1043 1094 X[1] = 3.; 1044 1095 X[2] = 4.; 1045 1096 Independent(X); 1046 1097 1047 1098 // dependent variable vector 1048 1099 CppAD::vector< AD<double> > Y(m); 1049 1100 1050 1101 // check results vector 1051 1102 CppAD::vector< bool > Check(m * n); 1052 1103 1053 1104 // initialize index into Y 1054 1105 size_t index = 0; 1055 1056 // Y[0] 1106 1107 // Y[0] 1057 1108 Y[index] = X[0] + X[1] + 5.; 1058 1109 Check[index * n + 0] = true; … … 1060 1111 Check[index * n + 2] = false; 1061 1112 index++; 1062 1063 // Y[1] 1113 1114 // Y[1] 1064 1115 Y[index] = Y[0]  (X[1] + X[2]); 1065 1116 Check[index * n + 0] = true; … … 1068 1119 index++; 1069 1120 1070 // Y[2] 1121 // Y[2] 1071 1122 // 2DO: There is a subtitle issue that has to do with using reverse 1072 1123 // jacobian sparsity patterns during the optimization process. … … 1078 1129 Check[index * n + 2] = true; 1079 1130 index++; 1080 1131 1081 1132 // check final index 1082 1133 assert( index == m ); 1083 1134 1084 1135 // create function object F : X > Y 1085 1136 ADFun<double> F(X, Y); 1086 F.optimize(); 1087 1137 if( conditional_skip ) 1138 F.optimize(); 1139 else 1140 F.optimize("no_conditional_skip"); 1141 1088 1142 //  1089 // dependency matrix for the identity function 1143 // dependency matrix for the identity function 1090 1144 CppAD::vector< std::set<size_t> > Sx(n); 1091 1145 size_t i, j; … … 1094 1148 Sx[i].insert(i); 1095 1149 } 1096 1150 1097 1151 // evaluate the dependency matrix for F(x) 1098 1152 CppAD::vector< std::set<size_t> > Sy(m); 1099 1153 Sy = F.ForSparseJac(n, Sx); 1100 1154 1101 1155 // check values 1102 1156 bool found; … … 1106 1160 ok &= (found == Check[i * n + j]); 1107 1161 } 1108 } 1109 1162 } 1163 1110 1164 return ok; 1111 1165 } … … 1113 1167 { bool ok = true; 1114 1168 using namespace CppAD; 1115 1169 1116 1170 // dimension of the domain space 1117 size_t n = 3; 1118 1171 size_t n = 3; 1172 1119 1173 // dimension of the range space 1120 1174 size_t m = 3; 1121 1122 // independent variable vector 1175 1176 // independent variable vector 1123 1177 CppAD::vector< AD<double> > X(n); 1124 X[0] = 2.; 1178 X[0] = 2.; 1125 1179 X[1] = 3.; 1126 1180 X[2] = 4.; 1127 1181 Independent(X); 1128 1182 1129 1183 // dependent variable vector 1130 1184 CppAD::vector< AD<double> > Y(m); 1131 1185 1132 1186 // check results vector 1133 1187 CppAD::vector< bool > Check(m * n); 1134 1188 1135 1189 // initialize index into Y 1136 1190 size_t index = 0; 1137 1138 // Y[0] 1191 1192 // Y[0] 1139 1193 Y[index] = X[0] + X[1] + 5.; 1140 1194 Check[index * n + 0] = true; … … 1142 1196 Check[index * n + 2] = false; 1143 1197 index++; 1144 1145 // Y[1] 1198 1199 // Y[1] 1146 1200 Y[index] = Y[0]  (X[1] + X[2]); 1147 1201 Check[index * n + 0] = true; … … 1150 1204 index++; 1151 1205 1152 // Y[2] 1206 // Y[2] 1153 1207 Y[index] = CondExpLe(X[0], X[1], X[1]+X[1], X[2]X[2]); 1154 1208 Check[index * n + 0] = false; … … 1156 1210 Check[index * n + 2] = true; 1157 1211 index++; 1158 1212 1159 1213 // check final index 1160 1214 assert( index == m ); 1161 1215 1162 1216 // create function object F : X > Y 1163 1217 ADFun<double> F(X, Y); 1164 F.optimize(); 1165 1218 if( conditional_skip ) 1219 F.optimize(); 1220 else 1221 F.optimize("no_conditional_skip"); 1222 1166 1223 //  1167 1224 // dependency matrix for the identity function … … 1172 1229 Py[ i * m + j ] = (i == j); 1173 1230 } 1174 1231 1175 1232 // evaluate the dependency matrix for F(x) 1176 1233 CppAD::vector< bool > Px(m * n); 1177 1234 Px = F.RevSparseJac(m, Py); 1178 1235 1179 1236 // check values 1180 1237 for(i = 0; i < m; i++) 1181 1238 { for(j = 0; j < n; j++) 1182 1239 ok &= (Px[i * n + j] == Check[i * n + j]); 1183 } 1184 1240 } 1241 1185 1242 return ok; 1186 1243 } … … 1189 1246 using CppAD::AD; 1190 1247 size_t i, j; 1191 1248 1192 1249 size_t n = 3; 1193 CppAD::vector< AD<double> > X(n); 1250 CppAD::vector< AD<double> > X(n); 1194 1251 X[0] = 1.; 1195 1252 X[1] = 2.; 1196 1253 X[2] = 3.; 1197 1254 CppAD::Independent(X); 1198 1255 1199 1256 size_t m = 1; 1200 1257 CppAD::vector< AD<double> > Y(m); 1201 Y[0] = CondExpGe( X[0], X[1], 1202 X[0] + (2. + X[1] + 3.) * X[1], 1258 Y[0] = CondExpGe( X[0], X[1], 1259 X[0] + (2. + X[1] + 3.) * X[1], 1203 1260 X[0] + (2. + X[2] + 3.) * X[1] 1204 1261 ); 1205 1262 1206 1263 CppAD::vector<bool> check(n * n); 1207 1264 check[0 * n + 0] = false; // partial w.r.t. x[0], x[0] … … 1216 1273 check[2 * n + 1] = true; // x[2], x[1] 1217 1274 check[2 * n + 2] = false; // x[2], x[2] 1218 1275 1219 1276 // create function object F : X > Y 1220 1277 CppAD::ADFun<double> F(X, Y); 1221 F.optimize(); 1222 1278 if( conditional_skip ) 1279 F.optimize(); 1280 else 1281 F.optimize("no_conditional_skip"); 1282 1223 1283 // sparsity pattern for the identity function U(x) = x 1224 1284 CppAD::vector<bool> Px(n * n); … … 1226 1286 for(j = 0; j < n; j++) 1227 1287 Px[ i * n + j ] = (i == j); 1228 1288 1229 1289 // compute sparsity pattern for Jacobian of F(U(x)) 1230 1290 CppAD::vector<bool> P_jac(m * n); 1231 1291 P_jac = F.ForSparseJac(n, Px); 1232 1233 // compute sparsity pattern for Hessian of F_k ( U(x) ) 1292 1293 // compute sparsity pattern for Hessian of F_k ( U(x) ) 1234 1294 CppAD::vector<bool> Py(m); 1235 1295 CppAD::vector<bool> Pxx(n * n); … … 1239 1299 for(i = 0; i < n * n; i++) 1240 1300 ok &= (Pxx[i] == check[i]); 1241 1301 1242 1302 return ok; 1243 1303 } … … 1248 1308 1249 1309 AD<double> zero(0.), one(1.), two(2.), three(3.); 1250 1310 1251 1311 size_t n = 4; 1252 CppAD::vector< AD<double> > X(n); 1312 CppAD::vector< AD<double> > X(n); 1253 1313 X[0] = zero; 1254 1314 X[1] = one; … … 1256 1316 X[3] = three; 1257 1317 CppAD::Independent(X); 1258 1318 1259 1319 size_t m = 4; 1260 1320 CppAD::vector< AD<double> > Y(m); … … 1265 1325 1266 1326 CppAD::ADFun<double> f(X, Y); 1267 f.optimize(); 1327 if( conditional_skip ) 1328 f.optimize(); 1329 else 1330 f.optimize("no_conditional_skip"); 1268 1331 1269 1332 CppAD::vector<double> x(n), y(m); … … 1284 1347 return ok; 1285 1348 } 1286 // check that CondExp properly handels expressions that get 1349 // check that CondExp properly handels expressions that get 1287 1350 // removed during opitmization 1288 1351 bool cond_exp_removed(void) … … 1290 1353 using CppAD::AD; 1291 1354 AD<double> zero(0.); 1292 1355 1293 1356 size_t n = 1; 1294 CppAD::vector< AD<double> > X(n); 1357 CppAD::vector< AD<double> > X(n); 1295 1358 X[0] = 1.0; 1296 1359 CppAD::Independent(X); 1297 1360 1298 1361 size_t m = 1; 1299 1362 CppAD::vector< AD<double> > Y(m); … … 1304 1367 1305 1368 CppAD::ADFun<double> f(X, Y); 1306 f.optimize(); 1369 if( conditional_skip ) 1370 f.optimize(); 1371 else 1372 f.optimize("no_conditional_skip"); 1307 1373 1308 1374 CppAD::vector<double> x(n), y(m), w(m), dw(n); … … 1336 1402 bool old_atomic_forward( 1337 1403 size_t id , 1338 size_t k , 1404 size_t k , 1339 1405 size_t n , 1340 1406 size_t m , 1341 1407 const CppAD::vector<bool>& vx , 1342 1408 CppAD::vector<bool>& vy , 1343 const CppAD::vector<double>& tx , 1409 const CppAD::vector<double>& tx , 1344 1410 CppAD::vector<double>& ty ) 1345 1411 { assert(n == 3 && m == 2); 1346 if( k > 0 ) 1412 if( k > 0 ) 1347 1413 return false; 1348 1414 … … 1352 1418 // y[1] = x[1] + x[2] 1353 1419 ty[1] = tx[1] + tx[2]; 1354 1420 1355 1421 if( vy.size() > 0 ) 1356 1422 { vy[0] = (vx[0]  vx[1]); 1357 1423 vy[1] = (vx[1]  vx[2]); 1358 1424 } 1359 return true; 1425 return true; 1360 1426 } 1361 1427 1362 1428 bool old_atomic_reverse( 1363 1429 size_t id , 1364 size_t k , 1365 size_t n , 1366 size_t m , 1367 const CppAD::vector<double>& tx , 1430 size_t k , 1431 size_t n , 1432 size_t m , 1433 const CppAD::vector<double>& tx , 1368 1434 const CppAD::vector<double>& ty , 1369 1435 CppAD::vector<double>& px , … … 1398 1464 my_union(r[2], r[2], s[1]); 1399 1465 1400 return true; 1466 return true; 1401 1467 } 1402 1468 … … 1421 1487 old_atomic_for_jac_sparse , 1422 1488 old_atomic_rev_jac_sparse , 1423 old_atomic_rev_hes_sparse 1489 old_atomic_rev_hes_sparse 1424 1490 ) 1425 1491 … … 1472 1538 ax[j] = 1. / 3.; 1473 1539 CppAD::Independent(ax); 1474 1540 1475 1541 // dependent variable vector 1476 1542 size_t m = 1; … … 1483 1549 } 1484 1550 CppAD::ADFun<double> f(ax, ay); 1485 1551 1486 1552 // Used to fail assert in optimize that forward mode results 1487 1553 // are identically equal 1488 f.optimize(); 1489 1554 if( conditional_skip ) 1555 f.optimize(); 1556 else 1557 f.optimize("no_conditional_skip"); 1558 1490 1559 return ok; 1491 1560 } … … 1497 1566 { bool ok = true; 1498 1567 using CppAD::vector; 1499 1568 1500 1569 vector< CppAD::AD<double> > ax(1), ay(1); 1501 ax[0] = 0.0; 1570 ax[0] = 0.0; 1502 1571 CppAD::Independent(ax); 1503 ay[0] = floor(ax[0]) + floor(ax[0]); 1572 ay[0] = floor(ax[0]) + floor(ax[0]); 1504 1573 CppAD::ADFun<double> f(ax, ay); 1505 1574 1506 1575 size_t size_before = f.size_var(); 1507 f.optimize(); 1576 if( conditional_skip ) 1577 f.optimize(); 1578 else 1579 f.optimize("no_conditional_skip"); 1508 1580 size_t size_after = f.size_var(); 1509 1581 ok &= size_after + 1 == size_before; 1510 1582 1511 1583 vector<double> x(1), y(1); 1512 1584 x[0] = 2.2; 1513 1585 y = f.Forward(0, x); 1514 1586 ok &= y[0] == 6.0; 1515 1587 1516 1588 return ok; 1517 1589 } 1518 1590 //  1519 void i_algo( 1591 void i_algo( 1520 1592 const CppAD::vector< CppAD::AD<double> >& ax , 1521 1593 CppAD::vector< CppAD::AD<double> >& ay ) … … 1527 1599 using CppAD::AD; 1528 1600 using CppAD::vector; 1529 1601 1530 1602 // Create a checkpoint version of the function i_algo 1531 1603 vector< AD<double> > au(1), av(1), aw(1); 1532 1604 au[0] = 1.0; 1533 1605 CppAD::checkpoint<double> i_check("i_check", i_algo, au, av); 1534 1535 // independent variable vector 1606 1607 // independent variable vector 1536 1608 vector< AD<double> > ax(2), ay(1); 1537 1609 ax[0] = 1.0; 1538 1610 ax[1] = 2.0; 1539 1611 Independent(ax); 1540 1612 1541 1613 // call atomic function that does not get used 1542 1614 au[0] = ax[0]; … … 1556 1628 x[1] = 2.0; 1557 1629 y_before = f.Forward(0, x); 1558 f.optimize(); 1630 if( conditional_skip ) 1631 f.optimize(); 1632 else 1633 f.optimize("no_conditional_skip"); 1559 1634 y_after = f.Forward(0, x); 1560 1635 1561 1636 ok &= y_before[0] == y_after[0]; 1562 1637 1563 1638 return ok; 1564 1639 } … … 1570 1645 using CppAD::AD; 1571 1646 using CppAD::vector; 1572 1647 1573 1648 // Create a checkpoint version of the function i_algo 1574 1649 vector< AD<double> > au(1), av(1), aw(1); 1575 1650 au[0] = 1.0; 1576 1651 CppAD::checkpoint<double> i_check("i_check", i_algo, au, av); 1577 1652 1578 1653 vector< AD<double> > ax(2), ay(1); 1579 AD<double> zero = 0.0; 1654 AD<double> zero = 0.0; 1580 1655 ax[0] = 1.0; 1581 1656 ax[1] = 1.0; … … 1597 1672 dx[1] = 2.0; 1598 1673 dy_before = f.Forward(1, dx); 1599 f.optimize(); 1674 if( conditional_skip ) 1675 f.optimize(); 1676 else 1677 f.optimize("no_conditional_skip"); 1600 1678 y_after = f.Forward(0, x); 1601 1679 dy_after = f.Forward(1, dx); … … 1620 1698 using CppAD::vector; 1621 1699 using CppAD::AD; 1622 1700 1623 1701 vector< AD<double> > ax(n), ay(1); 1624 1702 for(size_t j = 0; j < n; j++) … … 1627 1705 ay[0] = my_max(ax) + my_max(ax); 1628 1706 CppAD::ADFun<double> f(ax, ay); 1629 1630 f.optimize(); 1631 1707 1708 if( conditional_skip ) 1709 f.optimize(); 1710 else 1711 f.optimize("no_conditional_skip"); 1712 1632 1713 vector<double> x(n), w(1), dx(n); 1633 1714 for(size_t j = 0;j < n; j++) … … 1652 1733 1653 1734 // f(x) = x[0] + x[0] if x[0] >= 3 1654 // = x[0] + x[1] otherwise 1735 // = x[0] + x[1] otherwise 1655 1736 vector< AD<double> > ax(2), ay(3); 1656 1737 ax[0] = 1.0; … … 1666 1747 ay[2] = CppAD::CondExpGe(ax[0], three, exp(ax[0]), exp(ax[0]) ); 1667 1748 CppAD::ADFun<double> f(ax, ay); 1668 f.optimize(); 1749 if( conditional_skip ) 1750 f.optimize(); 1751 else 1752 f.optimize("no_conditional_skip"); 1669 1753 1670 1754 // check case where x[0] >= 3 … … 1692 1776 { bool ok = true; 1693 1777 atomic_sparsity_option = CppAD::atomic_base<double>::bool_sparsity_enum; 1778 conditional_skip = true; 1779 1780 // atomic sparsity loop 1694 1781 for(size_t i = 0; i < 2; i++) 1695 { // check conditional expression sparsity pattern 1696 // (used to optimize calls to atomic functions). 1782 { // check conditional expression sparsity pattern 1783 // (used to optimize calls to atomic functions). 1697 1784 ok &= atomic_cond_exp_sparsity(); 1698 1785 // check optimizing out entire atomic function … … 1701 1788 ok &= atomic_no_used(); 1702 1789 ok &= atomic_arguments(); 1703 atomic_sparsity_option = 1790 atomic_sparsity_option = 1704 1791 CppAD::atomic_base<double>::set_sparsity_enum; 1705 1792 } 1706 // check nested conditional expressions 1707 ok &= nested_cond_exp(); 1708 // check reverse dependency analysis optimization 1709 ok &= depend_one(); 1710 ok &= depend_two(); 1711 ok &= depend_three(); 1712 ok &= depend_four(); 1713 // check removal of duplicate expressions 1714 ok &= duplicate_one(); 1715 ok &= duplicate_two(); 1716 ok &= duplicate_three(); 1717 ok &= duplicate_four(); 1718 // convert sequence of additions to cummulative summation 1719 ok &= cummulative_sum(); 1720 ok &= forward_csum(); 1721 ok &= reverse_csum(); 1722 // sparsity patterns 1723 ok &= forward_sparse_jacobian(); 1724 ok &= reverse_sparse_jacobian(); 1725 ok &= reverse_sparse_hessian(); 1726 // check that CondExp properly detects dependencies 1727 ok &= cond_exp_depend(); 1728 // check that it properly handles expressions that have been removed 1729 ok &= cond_exp_removed(); 1730 // check old_atomic functions 1731 ok &= old_atomic_test(); 1732 // case where results are not identically equal 1733 ok &= not_identically_equal(); 1734 // case where a discrete function is used 1735 ok &= discrete_function(); 1736 // check conditional skip of an atomic function 1737 ok &= cond_exp_skip_atomic(); 1738 // check conditional dependence through atomic function 1739 ok &= cond_exp_atomic_dependence(); 1740 // check reverse mode conditional skipping 1741 ok &= cond_exp_reverse(); 1742 // check case where an expresion needed by both true and false case 1743 ok &= cond_exp_both_true_and_false(); 1793 1794 // conditional skip loop 1795 for(size_t i = 0; i < 2; i++) 1796 { // check nested conditional expressions 1797 ok &= nested_cond_exp(); 1798 // check reverse dependency analysis optimization 1799 ok &= depend_one(); 1800 ok &= depend_two(); 1801 ok &= depend_three(); 1802 ok &= depend_four(); 1803 // check removal of duplicate expressions 1804 ok &= duplicate_one(); 1805 ok &= duplicate_two(); 1806 ok &= duplicate_three(); 1807 ok &= duplicate_four(); 1808 // convert sequence of additions to cummulative summation 1809 ok &= cummulative_sum(); 1810 ok &= forward_csum(); 1811 ok &= reverse_csum(); 1812 // sparsity patterns 1813 ok &= forward_sparse_jacobian(); 1814 ok &= reverse_sparse_jacobian(); 1815 ok &= reverse_sparse_hessian(); 1816 // check that CondExp properly detects dependencies 1817 ok &= cond_exp_depend(); 1818 // check that it properly handles expressions that have been removed 1819 ok &= cond_exp_removed(); 1820 // check old_atomic functions 1821 ok &= old_atomic_test(); 1822 // case where results are not identically equal 1823 ok &= not_identically_equal(); 1824 // case where a discrete function is used 1825 ok &= discrete_function(); 1826 // check conditional skip of an atomic function 1827 ok &= cond_exp_skip_atomic(); 1828 // check conditional dependence through atomic function 1829 ok &= cond_exp_atomic_dependence(); 1830 // check reverse mode conditional skipping 1831 ok &= cond_exp_reverse(); 1832 // check case where an expresion needed by both true and false case 1833 ok &= cond_exp_both_true_and_false(); 1834 // 1835 conditional_skip = false; 1836 } 1744 1837 // 1745 1838 CppAD::user_atomic<double>::clear();
Note: See TracChangeset
for help on using the changeset viewer.