From ad3450294edb6a564d4660e4afcab115b3449f9a Mon Sep 17 00:00:00 2001 From: Carlo Camilloni Date: Thu, 25 Jun 2026 14:02:04 +0200 Subject: [PATCH 1/6] TTSKETCH: strong optimisation of the contraction calculation --- .../ttsketch/rt-ttmetad-2/COLVAR.reference | 24 + regtest/ttsketch/rt-ttmetad-2/HILLS | 4 + regtest/ttsketch/rt-ttmetad-2/HILLS.reference | 4 + regtest/ttsketch/rt-ttmetad-2/Makefile | 1 + regtest/ttsketch/rt-ttmetad-2/config | 4 + .../ttsketch/rt-ttmetad-2/forces.reference | 504 +++++++++++++++++ regtest/ttsketch/rt-ttmetad-2/plumed.dat | 24 + regtest/ttsketch/rt-ttmetad-2/traj.gro | 525 ++++++++++++++++++ regtest/ttsketch/rt-ttmetad-2/ttsketch.h5 | Bin 0 -> 66008 bytes src/ttsketch/TTHelper.cpp | 165 ++++-- src/ttsketch/TTHelper.h | 7 +- src/ttsketch/TTMetaD.cpp | 10 +- 12 files changed, 1237 insertions(+), 35 deletions(-) create mode 100644 regtest/ttsketch/rt-ttmetad-2/COLVAR.reference create mode 100644 regtest/ttsketch/rt-ttmetad-2/HILLS create mode 100644 regtest/ttsketch/rt-ttmetad-2/HILLS.reference create mode 100644 regtest/ttsketch/rt-ttmetad-2/Makefile create mode 100644 regtest/ttsketch/rt-ttmetad-2/config create mode 100644 regtest/ttsketch/rt-ttmetad-2/forces.reference create mode 100644 regtest/ttsketch/rt-ttmetad-2/plumed.dat create mode 100644 regtest/ttsketch/rt-ttmetad-2/traj.gro create mode 100644 regtest/ttsketch/rt-ttmetad-2/ttsketch.h5 diff --git a/regtest/ttsketch/rt-ttmetad-2/COLVAR.reference b/regtest/ttsketch/rt-ttmetad-2/COLVAR.reference new file mode 100644 index 0000000000..44a2f81dc0 --- /dev/null +++ b/regtest/ttsketch/rt-ttmetad-2/COLVAR.reference @@ -0,0 +1,24 @@ +#! FIELDS time phi psi tt.bias +#! SET min_phi -pi +#! SET max_phi pi + 0.000000 -1.2379 1.9470 5.0125 + 1.000000 -1.4839 1.9381 4.7386 + 2.000000 -1.3243 1.9663 5.7173 + 3.000000 -1.3340 1.9885 6.6124 + 4.000000 -1.4613 1.8901 6.9762 + 5.000000 -1.2202 2.0306 7.0226 + 6.000000 -1.3883 1.8776 8.4608 + 7.000000 -1.5481 1.9488 8.0675 + 8.000000 -1.8429 1.9025 4.8013 + 9.000000 -2.2424 1.9813 2.9942 + 10.000000 -1.1482 1.9684 7.7028 + 11.000000 -1.7580 1.9807 6.3092 + 12.000000 -1.3186 1.9486 7.0545 + 13.000000 -2.9911 1.9547 -0.0205 + 14.000000 -1.4112 2.0415 7.6081 + 15.000000 -2.5995 1.9166 2.0867 + 16.000000 -1.4608 1.9433 8.5441 + 17.000000 -1.3791 2.0625 8.9533 + 18.000000 -1.6771 1.9433 8.4947 + 19.000000 -1.5241 1.9642 10.3355 + 20.000000 -1.1997 1.9296 9.0332 diff --git a/regtest/ttsketch/rt-ttmetad-2/HILLS b/regtest/ttsketch/rt-ttmetad-2/HILLS new file mode 100644 index 0000000000..8841ba1482 --- /dev/null +++ b/regtest/ttsketch/rt-ttmetad-2/HILLS @@ -0,0 +1,4 @@ +#! FIELDS time phi psi sigma_phi sigma_psi height biasf +#! SET min_phi -pi +#! SET max_phi pi + 20.000000 -1.199652 1.929639 0.200000 0.200000 1.119533 10.000000 diff --git a/regtest/ttsketch/rt-ttmetad-2/HILLS.reference b/regtest/ttsketch/rt-ttmetad-2/HILLS.reference new file mode 100644 index 0000000000..7a38fbfc1e --- /dev/null +++ b/regtest/ttsketch/rt-ttmetad-2/HILLS.reference @@ -0,0 +1,4 @@ +#! FIELDS time phi psi sigma_phi sigma_psi height biasf +#! SET min_phi -pi +#! SET max_phi pi + 20.000000 -1.199652 1.929639 0.200000 0.200000 0.956750 10.000000 diff --git a/regtest/ttsketch/rt-ttmetad-2/Makefile b/regtest/ttsketch/rt-ttmetad-2/Makefile new file mode 100644 index 0000000000..3703b27cea --- /dev/null +++ b/regtest/ttsketch/rt-ttmetad-2/Makefile @@ -0,0 +1 @@ +include ../../scripts/test.make diff --git a/regtest/ttsketch/rt-ttmetad-2/config b/regtest/ttsketch/rt-ttmetad-2/config new file mode 100644 index 0000000000..6c0c53f78a --- /dev/null +++ b/regtest/ttsketch/rt-ttmetad-2/config @@ -0,0 +1,4 @@ +plumed_needs=libitensor +plumed_modules=ttsketch +type=driver +arg="--plumed plumed.dat --trajectory-stride 500 --timestep 0.002 --igro traj.gro --dump-forces forces --dump-forces-fmt %.4lf" diff --git a/regtest/ttsketch/rt-ttmetad-2/forces.reference b/regtest/ttsketch/rt-ttmetad-2/forces.reference new file mode 100644 index 0000000000..6de8227a4f --- /dev/null +++ b/regtest/ttsketch/rt-ttmetad-2/forces.reference @@ -0,0 +1,504 @@ +22 +0.3392 0.8024 -1.1416 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 4.2049 -8.4969 9.1777 +X 0.0000 0.0000 0.0000 +X -9.0418 15.8493 -12.9996 +X 0.0000 0.0000 0.0000 +X 10.6386 -13.4845 -0.2516 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -5.8017 6.1321 4.0735 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +-0.6416 -0.9488 1.5904 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -5.1548 11.5325 -14.6415 +X 0.0000 0.0000 0.0000 +X 12.4431 -24.0232 23.6159 +X 0.0000 0.0000 0.0000 +X -19.1651 19.4689 -2.9547 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 11.8769 -6.9783 -6.0197 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +0.3718 0.6267 -0.9985 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 4.7791 -11.5022 12.4824 +X 0.0000 0.0000 0.0000 +X -11.5151 25.4026 -23.7572 +X 0.0000 0.0000 0.0000 +X 21.3906 -20.0400 9.2219 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -14.6545 6.1396 2.0528 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +0.0484 -0.5185 0.4702 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 3.4750 -9.0333 9.8413 +X 0.0000 0.0000 0.0000 +X -10.3429 27.6455 -31.3154 +X 0.0000 0.0000 0.0000 +X 29.6775 -19.5861 27.6992 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -22.8096 0.9739 -6.2250 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +-2.7135 -1.0574 3.7709 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -15.6121 33.3615 -35.9836 +X 0.0000 0.0000 0.0000 +X 40.1842 -83.3356 79.0423 +X 0.0000 0.0000 0.0000 +X -81.7563 63.2698 -31.7764 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 57.1841 -13.2957 -11.2823 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +2.6284 2.6799 -5.3083 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 34.3505 -66.4796 59.7741 +X 0.0000 0.0000 0.0000 +X -82.5730 153.4644 -122.1257 +X 0.0000 0.0000 0.0000 +X 132.8211 -121.1657 52.2670 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -84.5987 34.1809 10.0847 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +-0.6941 2.8385 -2.1445 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -4.9084 14.8476 -14.1293 +X 0.0000 0.0000 0.0000 +X 19.9532 -66.1387 72.2110 +X 0.0000 0.0000 0.0000 +X -80.1786 45.0847 -75.4930 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 65.1338 6.2064 17.4113 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +-2.5925 -11.9624 14.5549 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -36.0889 82.9798 -93.0792 +X 0.0000 0.0000 0.0000 +X 62.8688 -138.2918 102.1982 +X 0.0000 0.0000 0.0000 +X -57.1003 122.6003 62.6054 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 30.3204 -67.2883 -71.7245 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +-1.5792 -6.6808 8.2601 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -18.0016 43.1472 -61.9284 +X 0.0000 0.0000 0.0000 +X 35.5299 -88.0116 85.8696 +X 0.0000 0.0000 0.0000 +X -43.1716 78.5709 21.8681 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 25.6433 -33.7064 -45.8093 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +0.6009 -3.8073 3.2064 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -21.0973 20.0305 -36.9089 +X 0.0000 0.0000 0.0000 +X 26.4570 -32.0899 41.0576 +X 0.0000 0.0000 0.0000 +X 7.6563 25.9075 25.6145 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -13.0161 -13.8482 -29.7632 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +2.3645 12.4411 -14.8056 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 68.6988 -143.6330 59.3899 +X 0.0000 0.0000 0.0000 +X -121.3598 261.3663 -59.8884 +X 0.0000 0.0000 0.0000 +X 89.6163 -199.4288 -86.4637 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -36.9553 81.6955 86.9622 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +0.9392 -2.1684 1.2292 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -10.5652 18.7203 -9.5495 +X 0.0000 0.0000 0.0000 +X 10.4496 -23.7792 5.2535 +X 0.0000 0.0000 0.0000 +X 12.0832 11.1432 18.5406 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -11.9676 -6.0843 -14.2446 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +-0.6482 0.4757 0.1726 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 1.7857 -4.5650 -0.8338 +X 0.0000 0.0000 0.0000 +X 0.6643 2.8198 2.6318 +X 0.0000 0.0000 0.0000 +X -8.6810 0.2521 -2.3954 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 6.2310 1.4931 0.5974 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +0.2321 -0.3647 0.1326 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -22.0060 14.4412 -19.3899 +X 0.0000 0.0000 0.0000 +X 25.9506 -17.9393 22.4898 +X 0.0000 0.0000 0.0000 +X 13.8806 -4.4235 13.4003 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -17.8251 7.9216 -16.5003 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +1.0972 -0.6129 -0.4843 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -2.4509 8.9772 0.2333 +X 0.0000 0.0000 0.0000 +X -10.8155 27.6448 -9.3088 +X 0.0000 0.0000 0.0000 +X 53.6576 -43.0775 7.4468 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -40.3912 6.4555 1.6287 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +0.8216 -0.5816 -0.2399 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -21.3939 30.2228 -16.1955 +X 0.0000 0.0000 0.0000 +X 25.6523 -45.3761 15.1324 +X 0.0000 0.0000 0.0000 +X 10.7215 6.0247 24.0897 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -14.9799 9.1286 -23.0267 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +0.1359 -1.2743 1.1384 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -3.5802 21.1283 -3.0993 +X 0.0000 0.0000 0.0000 +X 7.7594 -50.8804 3.3247 +X 0.0000 0.0000 0.0000 +X -19.9121 38.6250 17.5872 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 15.7329 -8.8730 -17.8125 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +1.7760 -2.0384 0.2624 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 3.5634 -11.4141 -2.7927 +X 0.0000 0.0000 0.0000 +X -27.7695 81.5986 11.9611 +X 0.0000 0.0000 0.0000 +X 85.0568 -79.8432 -35.3935 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -60.8506 9.6586 26.2252 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +0.8205 -0.3530 -0.4675 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -18.2913 73.7745 -1.7838 +X 0.0000 0.0000 0.0000 +X 25.9443 -131.2828 -19.4490 +X 0.0000 0.0000 0.0000 +X -8.6308 59.2111 80.5624 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.9778 -1.7028 -59.3296 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +0.5974 0.2830 -0.8803 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -7.6297 37.2833 5.1969 +X 0.0000 0.0000 0.0000 +X 8.8891 -65.9177 -19.7034 +X 0.0000 0.0000 0.0000 +X 1.3687 30.2035 44.6085 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X -2.6281 -1.5690 -30.1020 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +22 +-2.2878 -2.5878 4.8756 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 16.0775 -118.6758 -48.1327 +X 0.0000 0.0000 0.0000 +X -6.4971 152.6978 107.1400 +X 0.0000 0.0000 0.0000 +X -67.8256 -34.4113 -136.9059 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 58.2452 0.3892 77.8986 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 +X 0.0000 0.0000 0.0000 diff --git a/regtest/ttsketch/rt-ttmetad-2/plumed.dat b/regtest/ttsketch/rt-ttmetad-2/plumed.dat new file mode 100644 index 0000000000..1d02ca1479 --- /dev/null +++ b/regtest/ttsketch/rt-ttmetad-2/plumed.dat @@ -0,0 +1,24 @@ +phi: TORSION ATOMS=5,7,9,15 NOPBC +psi: ANGLE ATOMS=7,9,15 + +tt: TTMETAD ... + RESTART=YES + ARG=phi,psi + SIGMA=0.20,0.20 + HEIGHT=1.20 + PACE=500 + TEMP=300.0 + FILE=HILLS + FMT=%12.6f + BIASFACTOR=10 + SKETCH_RANK=2 + SKETCH_INITRANK=4 + SKETCH_PACE=5000 + INTERVAL_MIN=-3.14159,0.0 + INTERVAL_MAX=3.14159,3.14159 + SKETCH_NBASIS=5 + SKETCH_ALPHA=1.0 + DETERMINISTIC +... + +PRINT STRIDE=500 FILE=COLVAR ARG=phi,psi,tt.bias FMT=%8.4f diff --git a/regtest/ttsketch/rt-ttmetad-2/traj.gro b/regtest/ttsketch/rt-ttmetad-2/traj.gro new file mode 100644 index 0000000000..c1c387a014 --- /dev/null +++ b/regtest/ttsketch/rt-ttmetad-2/traj.gro @@ -0,0 +1,525 @@ +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 0.00000 + 22 + 1ACE HH31 1 1.474 1.585 1.200 + 1ACE CH3 2 1.483 1.508 1.277 + 1ACE HH32 3 1.476 1.561 1.372 + 1ACE HH33 4 1.578 1.455 1.278 + 1ACE C 5 1.353 1.428 1.279 + 1ACE O 6 1.263 1.449 1.357 + 2ALA N 7 1.343 1.328 1.191 + 2ALA H 8 1.415 1.321 1.120 + 2ALA CA 9 1.233 1.239 1.159 + 2ALA HA 10 1.144 1.302 1.155 + 2ALA CB 11 1.244 1.182 1.013 + 2ALA HB1 12 1.341 1.136 0.992 + 2ALA HB2 13 1.159 1.117 0.994 + 2ALA HB3 14 1.242 1.265 0.942 + 2ALA C 15 1.207 1.140 1.271 + 2ALA O 16 1.214 1.017 1.241 + 3NME N 17 1.191 1.177 1.398 + 3NME H 18 1.192 1.275 1.421 + 3NME CH3 19 1.189 1.086 1.518 + 3NME HH31 20 1.170 0.983 1.487 + 3NME HH32 21 1.283 1.087 1.574 + 3NME HH33 22 1.108 1.127 1.578 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 1.00000 + 22 + 1ACE HH31 1 1.480 1.571 1.214 + 1ACE CH3 2 1.481 1.493 1.289 + 1ACE HH32 3 1.502 1.528 1.390 + 1ACE HH33 4 1.551 1.417 1.255 + 1ACE C 5 1.344 1.432 1.275 + 1ACE O 6 1.250 1.462 1.345 + 2ALA N 7 1.342 1.327 1.193 + 2ALA H 8 1.430 1.313 1.144 + 2ALA CA 9 1.233 1.244 1.166 + 2ALA HA 10 1.144 1.307 1.173 + 2ALA CB 11 1.240 1.189 1.017 + 2ALA HB1 12 1.327 1.124 1.000 + 2ALA HB2 13 1.150 1.128 1.005 + 2ALA HB3 14 1.251 1.267 0.941 + 2ALA C 15 1.221 1.133 1.271 + 2ALA O 16 1.217 1.015 1.238 + 3NME N 17 1.204 1.174 1.395 + 3NME H 18 1.200 1.275 1.398 + 3NME CH3 19 1.188 1.089 1.516 + 3NME HH31 20 1.083 1.086 1.543 + 3NME HH32 21 1.233 0.990 1.511 + 3NME HH33 22 1.241 1.141 1.596 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 2.00000 + 22 + 1ACE HH31 1 1.532 1.520 1.209 + 1ACE CH3 2 1.478 1.493 1.300 + 1ACE HH32 3 1.465 1.586 1.356 + 1ACE HH33 4 1.548 1.426 1.350 + 1ACE C 5 1.352 1.423 1.279 + 1ACE O 6 1.252 1.461 1.340 + 2ALA N 7 1.351 1.326 1.190 + 2ALA H 8 1.442 1.293 1.160 + 2ALA CA 9 1.232 1.244 1.160 + 2ALA HA 10 1.146 1.310 1.151 + 2ALA CB 11 1.241 1.190 1.016 + 2ALA HB1 12 1.333 1.132 1.008 + 2ALA HB2 13 1.160 1.123 0.986 + 2ALA HB3 14 1.242 1.280 0.955 + 2ALA C 15 1.203 1.138 1.270 + 2ALA O 16 1.161 1.021 1.240 + 3NME N 17 1.230 1.171 1.396 + 3NME H 18 1.257 1.266 1.417 + 3NME CH3 19 1.217 1.090 1.512 + 3NME HH31 20 1.144 1.011 1.493 + 3NME HH32 21 1.307 1.029 1.526 + 3NME HH33 22 1.212 1.146 1.605 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 3.00000 + 22 + 1ACE HH31 1 1.439 1.582 1.175 + 1ACE CH3 2 1.474 1.516 1.254 + 1ACE HH32 3 1.480 1.585 1.338 + 1ACE HH33 4 1.569 1.465 1.242 + 1ACE C 5 1.364 1.419 1.280 + 1ACE O 6 1.277 1.446 1.367 + 2ALA N 7 1.358 1.323 1.194 + 2ALA H 8 1.443 1.313 1.140 + 2ALA CA 9 1.235 1.243 1.164 + 2ALA HA 10 1.150 1.310 1.170 + 2ALA CB 11 1.240 1.197 1.019 + 2ALA HB1 12 1.316 1.119 1.016 + 2ALA HB2 13 1.145 1.157 0.982 + 2ALA HB3 14 1.279 1.276 0.955 + 2ALA C 15 1.201 1.137 1.272 + 2ALA O 16 1.172 1.021 1.232 + 3NME N 17 1.218 1.166 1.402 + 3NME H 18 1.240 1.259 1.434 + 3NME CH3 19 1.186 1.086 1.518 + 3NME HH31 20 1.225 0.984 1.527 + 3NME HH32 21 1.193 1.134 1.616 + 3NME HH33 22 1.081 1.058 1.509 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 4.00000 + 22 + 1ACE HH31 1 1.549 1.508 1.196 + 1ACE CH3 2 1.500 1.486 1.290 + 1ACE HH32 3 1.487 1.571 1.357 + 1ACE HH33 4 1.563 1.415 1.343 + 1ACE C 5 1.362 1.425 1.270 + 1ACE O 6 1.265 1.465 1.340 + 2ALA N 7 1.349 1.324 1.182 + 2ALA H 8 1.432 1.287 1.138 + 2ALA CA 9 1.221 1.249 1.168 + 2ALA HA 10 1.138 1.318 1.184 + 2ALA CB 11 1.201 1.194 1.025 + 2ALA HB1 12 1.276 1.117 1.005 + 2ALA HB2 13 1.096 1.165 1.014 + 2ALA HB3 14 1.229 1.265 0.947 + 2ALA C 15 1.217 1.141 1.275 + 2ALA O 16 1.234 1.024 1.243 + 3NME N 17 1.183 1.174 1.400 + 3NME H 18 1.184 1.274 1.412 + 3NME CH3 19 1.187 1.078 1.509 + 3NME HH31 20 1.248 0.990 1.490 + 3NME HH32 21 1.220 1.120 1.604 + 3NME HH33 22 1.088 1.035 1.527 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 5.00000 + 22 + 1ACE HH31 1 1.449 1.585 1.168 + 1ACE CH3 2 1.479 1.518 1.248 + 1ACE HH32 3 1.523 1.577 1.328 + 1ACE HH33 4 1.565 1.461 1.213 + 1ACE C 5 1.364 1.422 1.284 + 1ACE O 6 1.305 1.438 1.389 + 2ALA N 7 1.347 1.326 1.187 + 2ALA H 8 1.423 1.328 1.122 + 2ALA CA 9 1.226 1.241 1.162 + 2ALA HA 10 1.139 1.308 1.162 + 2ALA CB 11 1.236 1.193 1.023 + 2ALA HB1 12 1.314 1.117 1.012 + 2ALA HB2 13 1.137 1.167 0.986 + 2ALA HB3 14 1.273 1.278 0.966 + 2ALA C 15 1.195 1.133 1.268 + 2ALA O 16 1.173 1.016 1.239 + 3NME N 17 1.204 1.175 1.393 + 3NME H 18 1.211 1.275 1.403 + 3NME CH3 19 1.188 1.090 1.513 + 3NME HH31 20 1.089 1.044 1.509 + 3NME HH32 21 1.267 1.014 1.509 + 3NME HH33 22 1.189 1.145 1.607 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 6.00000 + 22 + 1ACE HH31 1 1.517 1.511 1.181 + 1ACE CH3 2 1.490 1.488 1.284 + 1ACE HH32 3 1.482 1.582 1.339 + 1ACE HH33 4 1.569 1.421 1.316 + 1ACE C 5 1.359 1.414 1.282 + 1ACE O 6 1.272 1.447 1.358 + 2ALA N 7 1.351 1.320 1.186 + 2ALA H 8 1.434 1.297 1.133 + 2ALA CA 9 1.220 1.251 1.159 + 2ALA HA 10 1.139 1.323 1.167 + 2ALA CB 11 1.220 1.194 1.018 + 2ALA HB1 12 1.298 1.120 1.003 + 2ALA HB2 13 1.120 1.158 0.994 + 2ALA HB3 14 1.224 1.286 0.960 + 2ALA C 15 1.201 1.139 1.270 + 2ALA O 16 1.190 1.022 1.239 + 3NME N 17 1.203 1.178 1.393 + 3NME H 18 1.207 1.277 1.409 + 3NME CH3 19 1.211 1.102 1.515 + 3NME HH31 20 1.111 1.064 1.534 + 3NME HH32 21 1.275 1.017 1.492 + 3NME HH33 22 1.264 1.151 1.597 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 7.00000 + 22 + 1ACE HH31 1 1.483 1.590 1.185 + 1ACE CH3 2 1.505 1.501 1.245 + 1ACE HH32 3 1.538 1.533 1.344 + 1ACE HH33 4 1.580 1.430 1.209 + 1ACE C 5 1.379 1.418 1.267 + 1ACE O 6 1.298 1.443 1.349 + 2ALA N 7 1.360 1.320 1.187 + 2ALA H 8 1.426 1.297 1.114 + 2ALA CA 9 1.224 1.253 1.180 + 2ALA HA 10 1.147 1.326 1.205 + 2ALA CB 11 1.174 1.215 1.037 + 2ALA HB1 12 1.245 1.152 0.983 + 2ALA HB2 13 1.084 1.154 1.032 + 2ALA HB3 14 1.153 1.311 0.992 + 2ALA C 15 1.212 1.141 1.280 + 2ALA O 16 1.219 1.022 1.244 + 3NME N 17 1.191 1.176 1.408 + 3NME H 18 1.207 1.275 1.423 + 3NME CH3 19 1.162 1.068 1.509 + 3NME HH31 20 1.229 0.982 1.503 + 3NME HH32 21 1.162 1.109 1.610 + 3NME HH33 22 1.056 1.044 1.499 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 8.00000 + 22 + 1ACE HH31 1 1.576 1.425 1.169 + 1ACE CH3 2 1.523 1.455 1.260 + 1ACE HH32 3 1.547 1.556 1.294 + 1ACE HH33 4 1.566 1.393 1.338 + 1ACE C 5 1.385 1.422 1.255 + 1ACE O 6 1.308 1.499 1.305 + 2ALA N 7 1.346 1.311 1.189 + 2ALA H 8 1.419 1.248 1.159 + 2ALA CA 9 1.210 1.260 1.193 + 2ALA HA 10 1.139 1.326 1.242 + 2ALA CB 11 1.152 1.252 1.051 + 2ALA HB1 12 1.230 1.210 0.987 + 2ALA HB2 13 1.066 1.185 1.052 + 2ALA HB3 14 1.127 1.354 1.024 + 2ALA C 15 1.206 1.136 1.282 + 2ALA O 16 1.195 1.023 1.229 + 3NME N 17 1.210 1.152 1.420 + 3NME H 18 1.221 1.243 1.463 + 3NME CH3 19 1.185 1.048 1.518 + 3NME HH31 20 1.195 0.948 1.475 + 3NME HH32 21 1.261 1.070 1.593 + 3NME HH33 22 1.088 1.064 1.565 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 9.00000 + 22 + 1ACE HH31 1 1.515 1.474 1.040 + 1ACE CH3 2 1.535 1.461 1.147 + 1ACE HH32 3 1.561 1.560 1.184 + 1ACE HH33 4 1.612 1.386 1.165 + 1ACE C 5 1.406 1.430 1.217 + 1ACE O 6 1.361 1.502 1.307 + 2ALA N 7 1.345 1.316 1.190 + 2ALA H 8 1.384 1.254 1.121 + 2ALA CA 9 1.217 1.277 1.242 + 2ALA HA 10 1.187 1.348 1.319 + 2ALA CB 11 1.111 1.278 1.134 + 2ALA HB1 12 1.129 1.198 1.062 + 2ALA HB2 13 1.018 1.269 1.189 + 2ALA HB3 14 1.121 1.372 1.079 + 2ALA C 15 1.217 1.133 1.309 + 2ALA O 16 1.293 1.044 1.265 + 3NME N 17 1.132 1.116 1.408 + 3NME H 18 1.076 1.195 1.437 + 3NME CH3 19 1.112 1.003 1.490 + 3NME HH31 20 1.156 0.910 1.456 + 3NME HH32 21 1.153 1.026 1.588 + 3NME HH33 22 1.005 0.985 1.500 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 10.00000 + 22 + 1ACE HH31 1 1.536 1.485 1.135 + 1ACE CH3 2 1.521 1.460 1.240 + 1ACE HH32 3 1.528 1.557 1.289 + 1ACE HH33 4 1.602 1.397 1.276 + 1ACE C 5 1.385 1.401 1.261 + 1ACE O 6 1.341 1.416 1.373 + 2ALA N 7 1.322 1.332 1.167 + 2ALA H 8 1.364 1.316 1.076 + 2ALA CA 9 1.190 1.273 1.177 + 2ALA HA 10 1.125 1.352 1.216 + 2ALA CB 11 1.134 1.230 1.038 + 2ALA HB1 12 1.167 1.128 1.016 + 2ALA HB2 13 1.026 1.238 1.046 + 2ALA HB3 14 1.174 1.287 0.953 + 2ALA C 15 1.183 1.157 1.283 + 2ALA O 16 1.111 1.061 1.259 + 3NME N 17 1.264 1.157 1.396 + 3NME H 18 1.320 1.241 1.405 + 3NME CH3 19 1.264 1.070 1.510 + 3NME HH31 20 1.163 1.057 1.548 + 3NME HH32 21 1.311 0.974 1.487 + 3NME HH33 22 1.326 1.109 1.592 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 11.00000 + 22 + 1ACE HH31 1 1.607 1.359 1.178 + 1ACE CH3 2 1.535 1.442 1.176 + 1ACE HH32 3 1.521 1.470 1.072 + 1ACE HH33 4 1.569 1.533 1.224 + 1ACE C 5 1.402 1.412 1.242 + 1ACE O 6 1.368 1.454 1.350 + 2ALA N 7 1.319 1.332 1.177 + 2ALA H 8 1.363 1.279 1.103 + 2ALA CA 9 1.194 1.289 1.231 + 2ALA HA 10 1.153 1.355 1.308 + 2ALA CB 11 1.092 1.300 1.113 + 2ALA HB1 12 1.137 1.243 1.032 + 2ALA HB2 13 1.001 1.243 1.131 + 2ALA HB3 14 1.070 1.405 1.095 + 2ALA C 15 1.198 1.143 1.290 + 2ALA O 16 1.237 1.057 1.213 + 3NME N 17 1.171 1.126 1.422 + 3NME H 18 1.164 1.213 1.472 + 3NME CH3 19 1.200 1.009 1.496 + 3NME HH31 20 1.118 0.938 1.493 + 3NME HH32 21 1.284 0.953 1.454 + 3NME HH33 22 1.230 1.031 1.598 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 12.00000 + 22 + 1ACE HH31 1 1.491 1.569 1.138 + 1ACE CH3 2 1.511 1.462 1.134 + 1ACE HH32 3 1.614 1.449 1.166 + 1ACE HH33 4 1.495 1.444 1.028 + 1ACE C 5 1.416 1.394 1.229 + 1ACE O 6 1.445 1.392 1.345 + 2ALA N 7 1.299 1.357 1.181 + 2ALA H 8 1.278 1.358 1.082 + 2ALA CA 9 1.196 1.301 1.267 + 2ALA HA 10 1.185 1.364 1.355 + 2ALA CB 11 1.063 1.319 1.182 + 2ALA HB1 12 1.057 1.248 1.100 + 2ALA HB2 13 0.976 1.288 1.240 + 2ALA HB3 14 1.037 1.423 1.161 + 2ALA C 15 1.226 1.161 1.304 + 2ALA O 16 1.317 1.094 1.258 + 3NME N 17 1.135 1.110 1.387 + 3NME H 18 1.064 1.174 1.419 + 3NME CH3 19 1.121 0.972 1.429 + 3NME HH31 20 1.163 0.900 1.358 + 3NME HH32 21 1.177 0.954 1.520 + 3NME HH33 22 1.016 0.942 1.428 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 13.00000 + 22 + 1ACE HH31 1 1.528 1.537 1.077 + 1ACE CH3 2 1.543 1.448 1.137 + 1ACE HH32 3 1.634 1.463 1.195 + 1ACE HH33 4 1.560 1.362 1.072 + 1ACE C 5 1.421 1.428 1.229 + 1ACE O 6 1.396 1.520 1.301 + 2ALA N 7 1.355 1.314 1.219 + 2ALA H 8 1.385 1.241 1.156 + 2ALA CA 9 1.249 1.276 1.311 + 2ALA HA 10 1.297 1.288 1.409 + 2ALA CB 11 1.122 1.359 1.284 + 2ALA HB1 12 1.072 1.325 1.193 + 2ALA HB2 13 1.053 1.352 1.367 + 2ALA HB3 14 1.153 1.463 1.275 + 2ALA C 15 1.204 1.131 1.290 + 2ALA O 16 1.225 1.070 1.185 + 3NME N 17 1.131 1.077 1.382 + 3NME H 18 1.104 1.132 1.463 + 3NME CH3 19 1.083 0.939 1.380 + 3NME HH31 20 1.012 0.925 1.298 + 3NME HH32 21 1.168 0.871 1.372 + 3NME HH33 22 1.036 0.915 1.475 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 14.00000 + 22 + 1ACE HH31 1 1.543 1.468 1.123 + 1ACE CH3 2 1.548 1.443 1.229 + 1ACE HH32 3 1.566 1.540 1.276 + 1ACE HH33 4 1.634 1.382 1.257 + 1ACE C 5 1.420 1.385 1.285 + 1ACE O 6 1.412 1.375 1.404 + 2ALA N 7 1.324 1.361 1.200 + 2ALA H 8 1.345 1.354 1.101 + 2ALA CA 9 1.189 1.323 1.244 + 2ALA HA 10 1.163 1.394 1.323 + 2ALA CB 11 1.094 1.344 1.125 + 2ALA HB1 12 1.098 1.270 1.045 + 2ALA HB2 13 0.995 1.352 1.170 + 2ALA HB3 14 1.120 1.440 1.079 + 2ALA C 15 1.166 1.169 1.284 + 2ALA O 16 1.055 1.145 1.330 + 3NME N 17 1.263 1.077 1.279 + 3NME H 18 1.356 1.102 1.247 + 3NME CH3 19 1.246 0.942 1.335 + 3NME HH31 20 1.142 0.911 1.345 + 3NME HH32 21 1.293 0.868 1.270 + 3NME HH33 22 1.292 0.934 1.434 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 15.00000 + 22 + 1ACE HH31 1 1.508 1.417 1.043 + 1ACE CH3 2 1.539 1.448 1.143 + 1ACE HH32 3 1.565 1.553 1.146 + 1ACE HH33 4 1.630 1.395 1.169 + 1ACE C 5 1.438 1.410 1.253 + 1ACE O 6 1.429 1.471 1.356 + 2ALA N 7 1.348 1.327 1.217 + 2ALA H 8 1.362 1.285 1.126 + 2ALA CA 9 1.233 1.289 1.298 + 2ALA HA 10 1.269 1.295 1.401 + 2ALA CB 11 1.107 1.379 1.286 + 2ALA HB1 12 1.057 1.360 1.191 + 2ALA HB2 13 1.034 1.369 1.366 + 2ALA HB3 14 1.141 1.482 1.282 + 2ALA C 15 1.192 1.141 1.266 + 2ALA O 16 1.220 1.089 1.156 + 3NME N 17 1.131 1.075 1.364 + 3NME H 18 1.117 1.126 1.450 + 3NME CH3 19 1.090 0.931 1.376 + 3NME HH31 20 1.156 0.880 1.446 + 3NME HH32 21 0.996 0.929 1.431 + 3NME HH33 22 1.090 0.874 1.283 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 16.00000 + 22 + 1ACE HH31 1 1.587 1.386 1.166 + 1ACE CH3 2 1.554 1.425 1.262 + 1ACE HH32 3 1.559 1.534 1.266 + 1ACE HH33 4 1.627 1.390 1.335 + 1ACE C 5 1.417 1.372 1.300 + 1ACE O 6 1.390 1.364 1.421 + 2ALA N 7 1.334 1.344 1.205 + 2ALA H 8 1.365 1.367 1.111 + 2ALA CA 9 1.190 1.324 1.235 + 2ALA HA 10 1.160 1.386 1.319 + 2ALA CB 11 1.105 1.363 1.112 + 2ALA HB1 12 1.129 1.297 1.029 + 2ALA HB2 13 0.999 1.351 1.131 + 2ALA HB3 14 1.125 1.467 1.088 + 2ALA C 15 1.162 1.180 1.282 + 2ALA O 16 1.054 1.128 1.258 + 3NME N 17 1.261 1.106 1.318 + 3NME H 18 1.349 1.152 1.336 + 3NME CH3 19 1.250 0.960 1.324 + 3NME HH31 20 1.264 0.908 1.229 + 3NME HH32 21 1.322 0.926 1.399 + 3NME HH33 22 1.153 0.923 1.357 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 17.00000 + 22 + 1ACE HH31 1 1.614 1.382 1.173 + 1ACE CH3 2 1.559 1.432 1.253 + 1ACE HH32 3 1.547 1.532 1.213 + 1ACE HH33 4 1.625 1.444 1.339 + 1ACE C 5 1.426 1.380 1.288 + 1ACE O 6 1.411 1.348 1.406 + 2ALA N 7 1.325 1.370 1.200 + 2ALA H 8 1.341 1.388 1.102 + 2ALA CA 9 1.189 1.318 1.239 + 2ALA HA 10 1.170 1.358 1.338 + 2ALA CB 11 1.084 1.388 1.158 + 2ALA HB1 12 1.113 1.369 1.055 + 2ALA HB2 13 0.994 1.336 1.189 + 2ALA HB3 14 1.075 1.493 1.186 + 2ALA C 15 1.170 1.163 1.252 + 2ALA O 16 1.097 1.103 1.171 + 3NME N 17 1.235 1.099 1.349 + 3NME H 18 1.305 1.153 1.399 + 3NME CH3 19 1.220 0.964 1.393 + 3NME HH31 20 1.315 0.912 1.406 + 3NME HH32 21 1.168 0.953 1.489 + 3NME HH33 22 1.178 0.897 1.318 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 18.00000 + 22 + 1ACE HH31 1 1.555 1.419 1.126 + 1ACE CH3 2 1.554 1.431 1.235 + 1ACE HH32 3 1.554 1.535 1.267 + 1ACE HH33 4 1.638 1.381 1.283 + 1ACE C 5 1.423 1.378 1.292 + 1ACE O 6 1.416 1.347 1.409 + 2ALA N 7 1.322 1.351 1.211 + 2ALA H 8 1.350 1.370 1.115 + 2ALA CA 9 1.185 1.318 1.251 + 2ALA HA 10 1.171 1.358 1.352 + 2ALA CB 11 1.091 1.400 1.149 + 2ALA HB1 12 1.059 1.346 1.060 + 2ALA HB2 13 1.008 1.443 1.205 + 2ALA HB3 14 1.154 1.480 1.109 + 2ALA C 15 1.163 1.166 1.255 + 2ALA O 16 1.060 1.108 1.214 + 3NME N 17 1.255 1.102 1.328 + 3NME H 18 1.335 1.159 1.354 + 3NME CH3 19 1.246 0.964 1.369 + 3NME HH31 20 1.154 0.913 1.344 + 3NME HH32 21 1.323 0.899 1.327 + 3NME HH33 22 1.251 0.953 1.477 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 19.00000 + 22 + 1ACE HH31 1 1.559 1.326 1.133 + 1ACE CH3 2 1.552 1.413 1.200 + 1ACE HH32 3 1.538 1.504 1.142 + 1ACE HH33 4 1.631 1.428 1.273 + 1ACE C 5 1.426 1.375 1.280 + 1ACE O 6 1.428 1.345 1.398 + 2ALA N 7 1.319 1.363 1.209 + 2ALA H 8 1.322 1.382 1.110 + 2ALA CA 9 1.191 1.329 1.265 + 2ALA HA 10 1.186 1.364 1.368 + 2ALA CB 11 1.075 1.391 1.178 + 2ALA HB1 12 1.098 1.383 1.072 + 2ALA HB2 13 0.980 1.348 1.211 + 2ALA HB3 14 1.066 1.499 1.192 + 2ALA C 15 1.169 1.174 1.275 + 2ALA O 16 1.086 1.118 1.203 + 3NME N 17 1.233 1.105 1.363 + 3NME H 18 1.312 1.150 1.405 + 3NME CH3 19 1.242 0.956 1.350 + 3NME HH31 20 1.219 0.899 1.440 + 3NME HH32 21 1.165 0.918 1.284 + 3NME HH33 22 1.333 0.921 1.301 + 10.00000 10.00000 10.00000 +Generated by trjconv : Gromacs Runs One Microsecond At Cannonball Speeds t= 20.00000 + 22 + 1ACE HH31 1 1.622 1.459 1.287 + 1ACE CH3 2 1.546 1.434 1.214 + 1ACE HH32 3 1.578 1.365 1.135 + 1ACE HH33 4 1.509 1.531 1.180 + 1ACE C 5 1.430 1.358 1.283 + 1ACE O 6 1.444 1.310 1.394 + 2ALA N 7 1.315 1.370 1.215 + 2ALA H 8 1.323 1.427 1.132 + 2ALA CA 9 1.178 1.334 1.258 + 2ALA HA 10 1.163 1.376 1.357 + 2ALA CB 11 1.079 1.393 1.158 + 2ALA HB1 12 1.096 1.348 1.060 + 2ALA HB2 13 0.979 1.367 1.192 + 2ALA HB3 14 1.087 1.501 1.165 + 2ALA C 15 1.163 1.177 1.270 + 2ALA O 16 1.073 1.114 1.216 + 3NME N 17 1.249 1.110 1.347 + 3NME H 18 1.324 1.165 1.388 + 3NME CH3 19 1.242 0.964 1.360 + 3NME HH31 20 1.325 0.920 1.416 + 3NME HH32 21 1.151 0.938 1.414 + 3NME HH33 22 1.221 0.913 1.266 + 10.00000 10.00000 10.00000 diff --git a/regtest/ttsketch/rt-ttmetad-2/ttsketch.h5 b/regtest/ttsketch/rt-ttmetad-2/ttsketch.h5 new file mode 100644 index 0000000000000000000000000000000000000000..bc600f9723ab4d8f53d85069babd5886f5d3733b GIT binary patch literal 66008 zcmeHQ3v?9K8J-QITfn+t#S~wI4_qy9mPk1%R(E-*OGR7-+gMZ*A(Chs0*M;z(QZ{_ zQG1LpnuCXmt&dc##8OYw)??y>c6}Y=qZ%#NSc?s`*wj`;E2T60|KARKGdq*V4(#mQ z&B^|A@BQ!0+;9H-KknRncTIHmnFo&;H$u@re!p^n5=b55YY%_Ae7;l`(kL(I11eRN z1gBey=)kAwe~MDf`J?&#BHx1Z&p1QTIB|YZuBLPkr=!ffL#d>NjzIK`88N+$W!zB6 zQuU}5qryTfaN!M5S=W}=d_dGHa{^+`DMU7Mnw4y zar+aj*8oRNW$n_2Dn(gQQ~fne394%IsP0U8C) z+H*@GG)kc-@6&ikY5!MMuBfkGR;vWG2Sm>(lHx->gC0od;tajM{z3ULO)u%k;(7GD z(({Xoic+k8w#37H3fX)`&m(zGZvoXktL0F*wLC3iHWoAboSl z?|gp2`oLB5yP&3kZ<>)0YEPRoaej3}rLyI-m3u~vnt2UQkgY+?G;MnJHD8* zanBFdf4;1;V?<>2^zEwx!=4Ubbn%DLTOJ(U@#$^19{1T}qdGoSHm=+K$Xo4`-#TUM zx>>W@Uq1DT=NG-VqJ84g-GLv)mbE`P>W=8918!?yv$nqe&Dc%h8{@(M;T{HjPup*9 zPkyUp%R)BmQ_$Rs9Ebo5Nq;*1Fn%c4zp*iju=;#GVVUZ(n}o?c=wGXZ`)Y zO%IQGHhf>*3zx`VH z@Y)mRh)Bz{L_i`S5s(N-1cr`)Psz1k>36kXnR`3}F2}=s{G*%w>Nty_-S~$(ju_Ma z@X$q9p20;3h~poKpAsxbXi9bMl6w8PzG2ylilvq5_%;HD0n&8bF2sXYal0M^B5awn zpt81p*$ORIA8g!i-hZl&poS#4J$C&^If@t&pK>)xJc5AOe^hRF@~fEuVEj<)N0p7* zP&m<&v^g3~R0AJFf~(=z#bIUR`DUa}w=poc`f4878K{P0T} z62N(m{ygVbHUJ!>-u@m(EOxbjL-MVUp3*0JG*RGQ$LS`O@?p39lQa23oQ{#1NRO)a zvYsdI!u^Kd41DmO%s-ym9OtIn*}3Lpi809hl6VUN(av#qQd-qa9IBU4dI`NIQYYG*m22`H|?ycJc&yZJR{rjSNbK43&U zQc`c7^n2+0en%|;y*xavm%EQkH8<-Q^iJ>bSy@sAgq=1|h4oT4ZWRR+a9~a$SAi3$ z9P&`g>9Ui|>+P#t=w%4%(BB&dT`gC+$0H>1(1h7|uUzVKJj~3Clf)wkh<*(0Q=?!HviU2KLE?lsQ?x%;evCfV0d{La#nM!i>9`-G8RGu=F%ehm z6LI_`cCz(@g}!K+XwB!fQhU2i^PKc~SNga~_IXcHmM5Ld`1GCk)NB;y;Su_|I6h(Y z+mTxT@+kNbvR%ooQmHYd7L@M*RisD`5nKXR{(a;h~GDyV%-IxcUq@f^C5nZre^tP z=T*6=hhjVeR}F^;BXEFux?;LO}FmV4q3_d*D?3F4`X}KPJE9_vNnECt~~_m~Pz< zpf9SvY0c-gQhU4I-}rqc%Tr;*@8Z0H!N>3Bb}tbjNlUQ6+1kCC#0Wo|t4ZP!1Vpik~qd)YWA9cuH*PpFT)2#gdnKjEQ8m3LXb|~-R z527}4!+gA6(#QKBu=2s;NoAN;uCtC-!|CztsUwxCdw4-zdK_{Xt-lx3LxcPH+zVl|Gb`?ef*DNe1s>pq`aoZyDr^`FKC#3=o4 z2Fznb>^~~EKf-zya5%m0u8aFC`iF8ZU0D_Md*n*W2VWJZ6IbI-y+tvuB;sW?uB4AF z1Hb%TclRz231n)fzVq0c^74U;>+U)(pxXrQJ)HKJQ9kTeJhi);@qAi1om@(MK%;7Y zUN$@DH-+~qzz6g8;*49IuV`oI0^cRZV9xJQp3n^lh;}Z-Hyva-LK~c#w-@4mj8SYW zKQX`O?QJapy==_eE6X-O%>n_vQ*{v}#GFDeD;JuVBUj~v#*2|-PNAm~mryzEMk}W~ z%gyWUt6b=1I1+e+g8Ll!ew{m9Emyh6;}eO8Cd}#ggudM_$HRR7LbQTNwwA(pn>sieE&OJm>jB7VU&3L$lh-Eci!H?EKh|oZ%-VbFvgQ|v+gcZMRFHs@w2mgoz=t$znH5@;t>Qy zyNCF_hsVVLgH!Q)FYip|{`h@t0oXkwelNZPPu`ruKMK~Emm?QmKnLq3POL>Tk>p0z1XYlTiBYV3y(BIoQ+SPKEdpsscJT$>J;`hle$HP2+_b;d3 z8|O}J=f^bDHg4&Czqpzt9zj6#V_=_J1$*FB{4Uxb_r~wjT&+*Ub$1EYBiXt^fD5`) zvzW4Rc`x7IZud8SKb_^NFyeP{-oW7FcXPX^G?J_(SmbQ&-b!MGpY1(K^>&Lzy9a&N zBFx)zYTcbUPt52Ko?Um>eGQ4E{$4HJE&Gmt#~cw{aot^l^@$@FE8vsrdUqh}?#^N> z2T1s(sv!U2F@Z-6fBpyLPjORk! z?nag)z<_l_D{AP6Olp{FzZbu7iZ%yUrzIN2`!T#XLw%8Ay-Vp0R8q#*vl_?Tdq3d( z`?TTQ`^NZr+nQ_~ln*`ClCb96U#0frm;-$uTkLAN$~_*c%keNDhi&1074N6u7BL@h z-;}O}{Bx4M$Gmsr3EjiF1pPS8$OpBjO_?~qTHDa+DP_;jZR+JcI6<}wHClza)wNeB zjqmtk%EmoSzxZP8wpWjcd~xh$Qx>e=(*E&7vFS5T9^WzHqhmf;x2d?}_>UiGJMiiE z+wb}0&o3-qd`hsVTzBN+_PJ|!-F4flYr|{C&wO#_^!vjPT^M{D z_b}jl3Qn4Syzh=LcH6hQgj&swSY5p6(6NVY4p04P=W{zdUI|~l=*1)NYE20oTi#e0-tynaW}iLjqi|WruD}JeH-|&_ zJb%mUH@qG`?%m*A5owv02uK7Z0uljY5RBms6qDaWSnQ zIr2O9!L@WH^p6;)>jeGY^kgPq8>hRT!kv1HVti4=%WA%eK4T0x=kI}ga$1%Q;m7kd3`{{%VKrb8X@mhY0CvQ%nciMhtUXEPo z<=6(~m{XOrllPmKBNuwA>jBDvy;dHiwAupB{X?Ob;YbX9Kb=xn%T@02h)6s%(KgoO zMO}`E`TMYD?#FTN#&&*8@L^(<<6R7x$B5|1z&=F;d*IZ1JkkDG`7!yupH9rx`b1oh zSNwD9jsSg8_lPy0*Glc}c7NC7m9ab(#(F$)e8Lz{%FTW{O+20vXMwY`d+OuF2*22C za_aT=igpk2dnwBi@Np`BkMKCN`{VbD0elKgI z{u}3RZ0E=H{Dv6ek8?FiJc5Ae$G|=*EJuLBsrX&AKkkj+f9Ps`BF67s!t;k5fI?q{ zeuo@Mbp-6~c7NmdH7rks5xNwr>A>~6)hhC#}p?~AO3h=?YyXafg@;C!vXXm=Q zh%w0gl6VUN(awdq`Muo#fQD1^_I7)D)uhr!)`Qky7|3K{(YT-UPeHJHz>Hzf&Y&EZdc1y z?(rCYD?L+47Nz>W7i-$rkG@~m%RL?kyBrVm`3r$=f-}yj+W9dX-X})*m0V2{k02oW zF|bc7_o=4z671Od_Rg}B|VJT3+poQmJu{!Lei@rA7H zUVg{#6AHlY8S#7fHavNA3je5gyLmZs5oIMJ@(UG`7& zdiyHZ-tG-_-Q8qY%T@02m?rVigxOekce=~*Fpu9mxgW>5mf{{hwy+=&1H literal 0 HcmV?d00001 diff --git a/src/ttsketch/TTHelper.cpp b/src/ttsketch/TTHelper.cpp index 8d51193ea5..2a6dafc3b7 100644 --- a/src/ttsketch/TTHelper.cpp +++ b/src/ttsketch/TTHelper.cpp @@ -35,52 +35,155 @@ MPS ttSumRead(const std::string& filename, unsigned count) { return tt; } -// Evaluate TT by contracting each core with its 1D basis vector phi(x_i): -// result = sum_{i_1,...,i_d} G_1[i_1] * ... * G_d[i_d] * phi_1(x_1,i_1) * ... * phi_d(x_d,i_d) -// The contraction is performed left-to-right so intermediate results stay rank-1 scalars. -double ttEval(const MPS& tt, const std::vector& basis, const std::vector& elements, bool conv) { +// Precomputed basis evaluations, reusable across ttEval and ttGrad +struct BasisCache { + std::vector phi; // regular basis vectors + std::vector dphi; // derivative basis vectors (empty if not needed) +}; + +BasisCache buildBasisCache(const MPS& tt, + const std::vector& basis, + const std::vector& elements, + bool conv, + bool with_grad = false) { int d = length(tt); auto s = siteInds(tt); - std::vector basis_evals(d); - for(int i = 1; i <= d; ++i) { - basis_evals[i - 1] = ITensor(s(i)); - for(int j = 1; j <= dim(s(i)); ++j) { - basis_evals[i - 1].set(s(i) = j, basis[i - 1](elements[i - 1], j, conv)); + BasisCache cache; + cache.phi.resize(d); + if (with_grad) { + cache.dphi.resize(d); + } + + // This loop is embarrassingly parallel and cheap; OMP overhead likely not worth it + // unless d is large (d > ~32). Profile before adding #pragma omp parallel for here. + for (int i = 0; i < d; ++i) { + cache.phi[i] = ITensor(s(i + 1)); + if (with_grad) { + cache.dphi[i] = ITensor(s(i + 1)); + } + for (int j = 1; j <= dim(s(i + 1)); ++j) { + cache.phi[i].set(s(i + 1) = j, basis[i](elements[i], j, conv)); + if (with_grad) { + cache.dphi[i].set(s(i + 1) = j, basis[i].grad(elements[i], j, conv)); + } } } - auto result = tt(1) * basis_evals[0]; - for(int i = 2; i <= d; ++i) { - result *= tt(i) * basis_evals[i - 1]; + return cache; +} + +// Evaluate TT using a prebuilt BasisCache (avoids recomputing phi when called +// immediately before ttGrad). +double ttEvalCached(const MPS& tt, const BasisCache& cache) { + int d = length(tt); + auto result = tt(1) * cache.phi[0]; + for (int i = 2; i <= d; ++i) { + result *= tt(i) * cache.phi[i - 1]; } return elt(result); } -// Gradient of ttEval w.r.t. elements. For each dimension k, one TT contraction is performed -// with the k-th basis vector replaced by its derivative d phi_k/dx_k (chain rule). -// All basis evaluations are precomputed to avoid redundant work across the d contractions. -std::vector ttGrad(const MPS& tt, const std::vector& basis, const std::vector& elements, bool conv) { +// Original ttEval kept for standalone use (no redundant work when called alone) +double ttEval(const MPS& tt, + const std::vector& basis, + const std::vector& elements, + bool conv) { + auto cache = buildBasisCache(tt, basis, elements, conv, /*with_grad=*/false); + return ttEvalCached(tt, cache); +} + +// Gradient using prefix-suffix sweep: O(d) contractions instead of O(d^2). +// +// For a TT of length d, define: +// L[k] = (G_1 * phi_1) * ... * (G_k * phi_k) [left partial product up to k] +// R[k] = (G_k * phi_k) * ... * (G_d * phi_d) [right partial product from k] +// +// Then grad[k] = L[k-1] * (G_k * dphi_k) * R[k+1] +// +// We precompute all L[] (left-to-right pass) and all R[] (right-to-left pass), +// then combine them with a single derivative insertion per dimension. +std::vector ttGradCached(const MPS& tt, const BasisCache& cache) { int d = length(tt); - auto s = siteInds(tt); std::vector grad(d, 0.0); - std::vector basis_evals(d), basisd_evals(d); - for(int i = 1; i <= d; ++i) { - basis_evals[i - 1] = basisd_evals[i - 1] = ITensor(s(i)); - for(int j = 1; j <= dim(s(i)); ++j) { - basis_evals[i - 1].set(s(i) = j, basis[i - 1](elements[i - 1], j, conv)); - basisd_evals[i - 1].set(s(i) = j, basis[i - 1].grad(elements[i - 1], j, conv)); - } + + // --- Left prefix products: L[i] = contraction of sites 1..i with phi --- + // L[0] is a scalar 1 (identity for the left boundary) + // L[i] has the bond index between site i and i+1 dangling (or is scalar for i==d) + std::vector L(d + 1), R(d + 2); + + L[0] = ITensor(1.0); // scalar identity + for (int i = 1; i <= d; ++i) { + L[i] = L[i - 1] * (tt(i) * cache.phi[i - 1]); } - for(int k = 1; k <= d; ++k) { - // replace dimension k's basis vector with its derivative, keep all others unchanged - auto result = tt(1) * (k == 1 ? basisd_evals[0] : basis_evals[0]); - for(int i = 2; i <= d; ++i) { - result *= tt(i) * (k == i ? basisd_evals[i - 1] : basis_evals[i - 1]); - } - grad[k - 1] = elt(result); + + // --- Right suffix products: R[i] = contraction of sites i..d with phi --- + R[d + 1] = ITensor(1.0); + for (int i = d; i >= 1; --i) { + R[i] = (tt(i) * cache.phi[i - 1]) * R[i + 1]; + } + + // --- Combine: grad[k] = L[k-1] * (G_k * dphi_k) * R[k+1] --- + // Each combination is O(r^2 * n) instead of the full O(d * r^2 * n) + for (int k = 1; k <= d; ++k) { + grad[k - 1] = elt(L[k - 1] * (tt(k) * cache.dphi[k - 1]) * R[k + 1]); } + return grad; } +std::pair> + ttEvalAndGrad(const itensor::MPS& tt, + const std::vector& basis, + const std::vector& elements, +bool conv) { + int d = length(tt); + auto s = siteInds(tt); + + // Cache cores once — tt(i) does index lookup every call + std::vector cores(d); + for (int i = 0; i < d; ++i) { + cores[i] = tt(i + 1); + } + + // Build phi and dphi + std::vector phi(d), dphi(d); + for (int i = 0; i < d; ++i) { + phi[i] = dphi[i] = ITensor(s(i + 1)); + for (int j = 1; j <= dim(s(i + 1)); ++j) { + phi[i].set(s(i + 1) = j, basis[i](elements[i], j, conv)); + dphi[i].set(s(i + 1) = j, basis[i].grad(elements[i], j, conv)); + } + } + + // Prefix-suffix sweep + std::vector L(d + 1), R(d + 1); + L[0] = ITensor(1.0); + for (int i = 0; i < d; ++i) { + L[i + 1] = L[i] * (cores[i] * phi[i]); + } + R[d] = ITensor(1.0); + for (int i = d - 1; i >= 0; --i) { + R[i] = (cores[i] * phi[i]) * R[i + 1]; + } + + double val = elt(L[d]); + + std::vector grad(d); + for (int k = 0; k < d; ++k) { + grad[k] = elt(L[k] * (cores[k] * dphi[k]) * R[k + 1]); + } + + return {val, grad}; +} + +// Standalone gradient (for places that don't need the value) +std::vector ttGrad(const MPS& tt, + const std::vector& basis, + const std::vector& elements, + bool conv) { + auto cache = buildBasisCache(tt, basis, elements, conv, /*with_grad=*/true); + return ttGradCached(tt, cache); +} + // Compute covariance matrix, marginal means, and partition function of the TT distribution. // Precomputes ITensors for the three moment integrals int0/int1/int2 per dimension, // then evaluates expectations as TT contractions with these integral vectors. diff --git a/src/ttsketch/TTHelper.h b/src/ttsketch/TTHelper.h index a05007e1fb..e2c1781890 100644 --- a/src/ttsketch/TTHelper.h +++ b/src/ttsketch/TTHelper.h @@ -5,6 +5,7 @@ #ifdef __PLUMED_HAS_LIBITENSOR #ifdef __PLUMED_HAS_LIBHDF5 #include +#include namespace PLMD { namespace ttsketch { @@ -34,7 +35,11 @@ std::tuple, std::vector, double> covMat(const itensor::MP // Fill `grid` (bins x bins) with the 2D marginal density of `tt` for dimensions // pos1 and pos2, obtained by integrating out all other dimensions analytically. void marginal2d(const itensor::MPS& tt, const std::vector& basis, int pos1, int pos2, Matrix& grid, bool conv); - +std::pair> + ttEvalAndGrad(const itensor::MPS& tt, + const std::vector& basis, + const std::vector& elements, + bool conv); } } #endif // __PLUMED_HAS_LIBHDF5 diff --git a/src/ttsketch/TTMetaD.cpp b/src/ttsketch/TTMetaD.cpp index 6238f80133..4b6c7eedc3 100644 --- a/src/ttsketch/TTMetaD.cpp +++ b/src/ttsketch/TTMetaD.cpp @@ -618,9 +618,13 @@ double TTMetaD::getBias(const std::vector& cv) { } double TTMetaD::getBiasAndDerivatives(const std::vector& cv, std::vector& der) { - double bias = length(this->vb_) == 0 ? 0.0 : ttEval(this->vb_, this->sketch_basis_, cv, this->sketch_conv_); - if(length(this->vb_) != 0) { - der = ttGrad(this->vb_, this->sketch_basis_, cv, this->sketch_conv_); + double bias = 0.0; + if (length(this->vb_) != 0) { + auto [tt_bias, tt_grad] = ttEvalAndGrad(this->vb_, this->sketch_basis_, cv, this->sketch_conv_); + bias = tt_bias; + der = std::move(tt_grad); + } else { + der.assign(this->d_, 0.0); } unsigned nt = OpenMP::getNumThreads(); if(this->hills_.size() < 2 * nt || nt == 1) { From a5d327557080901ad721f2e80fb3e01eb7200a34 Mon Sep 17 00:00:00 2001 From: Carlo Camilloni Date: Thu, 25 Jun 2026 18:53:56 +0200 Subject: [PATCH 2/6] code cleaning and openMP --- src/ttsketch/TTHelper.cpp | 90 ++------------------------------------- src/ttsketch/TTHelper.h | 3 -- 2 files changed, 3 insertions(+), 90 deletions(-) diff --git a/src/ttsketch/TTHelper.cpp b/src/ttsketch/TTHelper.cpp index 2a6dafc3b7..ff5c75490a 100644 --- a/src/ttsketch/TTHelper.cpp +++ b/src/ttsketch/TTHelper.cpp @@ -4,6 +4,7 @@ See the COPYRIGHT file distributed with this software for license details. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ #include "TTHelper.h" +#include "tools/OpenMP.h" #ifdef __PLUMED_HAS_LIBITENSOR #ifdef __PLUMED_HAS_LIBHDF5 @@ -35,53 +36,6 @@ MPS ttSumRead(const std::string& filename, unsigned count) { return tt; } -// Precomputed basis evaluations, reusable across ttEval and ttGrad -struct BasisCache { - std::vector phi; // regular basis vectors - std::vector dphi; // derivative basis vectors (empty if not needed) -}; - -BasisCache buildBasisCache(const MPS& tt, - const std::vector& basis, - const std::vector& elements, - bool conv, - bool with_grad = false) { - int d = length(tt); - auto s = siteInds(tt); - BasisCache cache; - cache.phi.resize(d); - if (with_grad) { - cache.dphi.resize(d); - } - - // This loop is embarrassingly parallel and cheap; OMP overhead likely not worth it - // unless d is large (d > ~32). Profile before adding #pragma omp parallel for here. - for (int i = 0; i < d; ++i) { - cache.phi[i] = ITensor(s(i + 1)); - if (with_grad) { - cache.dphi[i] = ITensor(s(i + 1)); - } - for (int j = 1; j <= dim(s(i + 1)); ++j) { - cache.phi[i].set(s(i + 1) = j, basis[i](elements[i], j, conv)); - if (with_grad) { - cache.dphi[i].set(s(i + 1) = j, basis[i].grad(elements[i], j, conv)); - } - } - } - return cache; -} - -// Evaluate TT using a prebuilt BasisCache (avoids recomputing phi when called -// immediately before ttGrad). -double ttEvalCached(const MPS& tt, const BasisCache& cache) { - int d = length(tt); - auto result = tt(1) * cache.phi[0]; - for (int i = 2; i <= d; ++i) { - result *= tt(i) * cache.phi[i - 1]; - } - return elt(result); -} - // Original ttEval kept for standalone use (no redundant work when called alone) double ttEval(const MPS& tt, const std::vector& basis, @@ -99,37 +53,6 @@ double ttEval(const MPS& tt, // // Then grad[k] = L[k-1] * (G_k * dphi_k) * R[k+1] // -// We precompute all L[] (left-to-right pass) and all R[] (right-to-left pass), -// then combine them with a single derivative insertion per dimension. -std::vector ttGradCached(const MPS& tt, const BasisCache& cache) { - int d = length(tt); - std::vector grad(d, 0.0); - - // --- Left prefix products: L[i] = contraction of sites 1..i with phi --- - // L[0] is a scalar 1 (identity for the left boundary) - // L[i] has the bond index between site i and i+1 dangling (or is scalar for i==d) - std::vector L(d + 1), R(d + 2); - - L[0] = ITensor(1.0); // scalar identity - for (int i = 1; i <= d; ++i) { - L[i] = L[i - 1] * (tt(i) * cache.phi[i - 1]); - } - - // --- Right suffix products: R[i] = contraction of sites i..d with phi --- - R[d + 1] = ITensor(1.0); - for (int i = d; i >= 1; --i) { - R[i] = (tt(i) * cache.phi[i - 1]) * R[i + 1]; - } - - // --- Combine: grad[k] = L[k-1] * (G_k * dphi_k) * R[k+1] --- - // Each combination is O(r^2 * n) instead of the full O(d * r^2 * n) - for (int k = 1; k <= d; ++k) { - grad[k - 1] = elt(L[k - 1] * (tt(k) * cache.dphi[k - 1]) * R[k + 1]); - } - - return grad; -} - std::pair> ttEvalAndGrad(const itensor::MPS& tt, const std::vector& basis, @@ -137,6 +60,7 @@ std::pair> bool conv) { int d = length(tt); auto s = siteInds(tt); + const unsigned nt = OpenMP::getNumThreads(); // Cache cores once — tt(i) does index lookup every call std::vector cores(d); @@ -146,6 +70,7 @@ bool conv) { // Build phi and dphi std::vector phi(d), dphi(d); + #pragma omp parallel for num_threads(nt) schedule(static) for (int i = 0; i < d; ++i) { phi[i] = dphi[i] = ITensor(s(i + 1)); for (int j = 1; j <= dim(s(i + 1)); ++j) { @@ -175,15 +100,6 @@ bool conv) { return {val, grad}; } -// Standalone gradient (for places that don't need the value) -std::vector ttGrad(const MPS& tt, - const std::vector& basis, - const std::vector& elements, - bool conv) { - auto cache = buildBasisCache(tt, basis, elements, conv, /*with_grad=*/true); - return ttGradCached(tt, cache); -} - // Compute covariance matrix, marginal means, and partition function of the TT distribution. // Precomputes ITensors for the three moment integrals int0/int1/int2 per dimension, // then evaluates expectations as TT contractions with these integral vectors. diff --git a/src/ttsketch/TTHelper.h b/src/ttsketch/TTHelper.h index e2c1781890..04265d0efb 100644 --- a/src/ttsketch/TTHelper.h +++ b/src/ttsketch/TTHelper.h @@ -22,9 +22,6 @@ itensor::MPS ttSumRead(const std::string& filename, unsigned count); // Evaluate the TT bias at point `elements` by contracting each core with its // 1D basis vector. If conv=true, basis functions are evaluated in convolution mode. double ttEval(const itensor::MPS& tt, const std::vector& basis, const std::vector& elements, bool conv); -// Gradient of ttEval w.r.t. `elements`. Computed by replacing the k-th dimension's -// basis vector with its derivative while keeping all others unchanged (chain rule). -std::vector ttGrad(const itensor::MPS& tt, const std::vector& basis, const std::vector& elements, bool conv); // Compute the covariance matrix, mean vector, and partition function of the TT distribution. // Normalizes tt by Z = int tt(x) dx, then returns: // sigma(k,l) = E[x_k * x_l] - E[x_k]*E[x_l] (covariance matrix) From 327b7ae8666737e36c2d977a77976e0917995a30 Mon Sep 17 00:00:00 2001 From: Carlo Camilloni Date: Thu, 25 Jun 2026 19:08:00 +0200 Subject: [PATCH 3/6] fixed after code cleaning --- src/ttsketch/TTHelper.cpp | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/src/ttsketch/TTHelper.cpp b/src/ttsketch/TTHelper.cpp index ff5c75490a..ad9ac5ac71 100644 --- a/src/ttsketch/TTHelper.cpp +++ b/src/ttsketch/TTHelper.cpp @@ -41,8 +41,27 @@ double ttEval(const MPS& tt, const std::vector& basis, const std::vector& elements, bool conv) { - auto cache = buildBasisCache(tt, basis, elements, conv, /*with_grad=*/false); - return ttEvalCached(tt, cache); + const int d = length(tt); + auto s = siteInds(tt); + + std::vector cores(d), phi(d); + + for (int i = 0; i < d; ++i) { + cores[i] = tt(i + 1); + + phi[i] = ITensor(s(i + 1)); + for (int j = 1; j <= dim(s(i + 1)); ++j) { + phi[i].set(s(i + 1) = j, + basis[i](elements[i], j, conv)); + } + } + + auto result = cores[0] * phi[0]; + for (int i = 1; i < d; ++i) { + result *= cores[i] * phi[i]; + } + + return elt(result); } // Gradient using prefix-suffix sweep: O(d) contractions instead of O(d^2). From 371bf768b9d89c3559888baaeb83e663f7c6a909 Mon Sep 17 00:00:00 2001 From: Carlo Camilloni Date: Wed, 1 Jul 2026 14:30:52 +0200 Subject: [PATCH 4/6] more openMP: this squeeze some more performances with many CVs and openmp --- src/ttsketch/TTHelper.cpp | 226 +++++++++++++++++--------------------- 1 file changed, 101 insertions(+), 125 deletions(-) diff --git a/src/ttsketch/TTHelper.cpp b/src/ttsketch/TTHelper.cpp index ad9ac5ac71..409736d255 100644 --- a/src/ttsketch/TTHelper.cpp +++ b/src/ttsketch/TTHelper.cpp @@ -21,8 +21,7 @@ void ttWrite(const std::string& filename, const MPS& tt, unsigned count) { MPS ttRead(const std::string& filename, unsigned count) { auto f = h5_open(filename, 'r'); - auto tt = h5_read(f, "tt_" + std::to_string(count - 1)); - return tt; + return h5_read(f, "tt_" + std::to_string(count - 1)); } void ttSumWrite(const std::string& filename, const MPS& tt, unsigned count) { @@ -32,46 +31,17 @@ void ttSumWrite(const std::string& filename, const MPS& tt, unsigned count) { MPS ttSumRead(const std::string& filename, unsigned count) { auto f = h5_open(filename, 'r'); - auto tt = h5_read(f, "vb_" + std::to_string(count - 1)); - return tt; + return h5_read(f, "vb_" + std::to_string(count - 1)); } -// Original ttEval kept for standalone use (no redundant work when called alone) -double ttEval(const MPS& tt, - const std::vector& basis, - const std::vector& elements, - bool conv) { - const int d = length(tt); - auto s = siteInds(tt); - - std::vector cores(d), phi(d); - - for (int i = 0; i < d; ++i) { - cores[i] = tt(i + 1); - - phi[i] = ITensor(s(i + 1)); - for (int j = 1; j <= dim(s(i + 1)); ++j) { - phi[i].set(s(i + 1) = j, - basis[i](elements[i], j, conv)); - } - } - - auto result = cores[0] * phi[0]; - for (int i = 1; i < d; ++i) { - result *= cores[i] * phi[i]; - } - - return elt(result); -} - -// Gradient using prefix-suffix sweep: O(d) contractions instead of O(d^2). -// -// For a TT of length d, define: -// L[k] = (G_1 * phi_1) * ... * (G_k * phi_k) [left partial product up to k] -// R[k] = (G_k * phi_k) * ... * (G_d * phi_d) [right partial product from k] +// Evaluate TT and gradient together using a prefix-suffix sweep: O(d) contractions. // -// Then grad[k] = L[k-1] * (G_k * dphi_k) * R[k+1] +// For a TT of length d: +// L[i] = (G_1 * phi_1) * ... * (G_i * phi_i) [left partial product up to i] +// R[i] = (G_i * phi_i) * ... * (G_d * phi_d) [right partial product from i] // +// val = elt(L[d]) +// grad[k] = L[k-1] * (G_k * dphi_k) * R[k+1] std::pair> ttEvalAndGrad(const itensor::MPS& tt, const std::vector& basis, @@ -81,16 +51,11 @@ bool conv) { auto s = siteInds(tt); const unsigned nt = OpenMP::getNumThreads(); - // Cache cores once — tt(i) does index lookup every call - std::vector cores(d); - for (int i = 0; i < d; ++i) { - cores[i] = tt(i + 1); - } - - // Build phi and dphi - std::vector phi(d), dphi(d); + // Cache cores and build phi/dphi — all independent across dimensions + std::vector cores(d), phi(d), dphi(d); #pragma omp parallel for num_threads(nt) schedule(static) for (int i = 0; i < d; ++i) { + cores[i] = tt(i + 1); phi[i] = dphi[i] = ITensor(s(i + 1)); for (int j = 1; j <= dim(s(i + 1)); ++j) { phi[i].set(s(i + 1) = j, basis[i](elements[i], j, conv)); @@ -111,7 +76,9 @@ bool conv) { double val = elt(L[d]); + // Gradient combination — independent across k std::vector grad(d); + #pragma omp parallel for num_threads(nt) schedule(static) for (int k = 0; k < d; ++k) { grad[k] = elt(L[k] * (cores[k] * dphi[k]) * R[k + 1]); } @@ -119,105 +86,114 @@ bool conv) { return {val, grad}; } +// ttEval delegates to ttEvalAndGrad and discards the gradient. +// Called infrequently (error measurement, well-tempered height) so the +// cost of computing dphi is negligible relative to avoiding code duplication. +double ttEval(const MPS& tt, + const std::vector& basis, + const std::vector& elements, + bool conv) { + return ttEvalAndGrad(tt, basis, elements, conv).first; + // Compute covariance matrix, marginal means, and partition function of the TT distribution. // Precomputes ITensors for the three moment integrals int0/int1/int2 per dimension, // then evaluates expectations as TT contractions with these integral vectors. -std::tuple, std::vector, double> covMat(const MPS& tt, const std::vector& basis) { - int d = length(tt); - auto s = siteInds(tt); - // integral vectors: int0[i][j] = int phi_j(x_i) dx_i, etc. - std::vector basis_int0(d), basis_int1(d), basis_int2(d); - for(int i = 1; i <= d; ++i) { - basis_int0[i - 1] = basis_int1[i - 1] = basis_int2[i - 1] = ITensor(s(i)); - for(int j = 1; j <= dim(s(i)); ++j) { - basis_int0[i - 1].set(s(i) = j, basis[i - 1].int0(j)); - basis_int1[i - 1].set(s(i) = j, basis[i - 1].int1(j)); - basis_int2[i - 1].set(s(i) = j, basis[i - 1].int2(j)); + std::tuple, std::vector, double> covMat(const MPS& tt, const std::vector& basis) { + int d = length(tt); + auto s = siteInds(tt); + // integral vectors: int0[i][j] = int phi_j(x_i) dx_i, etc. + std::vector basis_int0(d), basis_int1(d), basis_int2(d); + for(int i = 1; i <= d; ++i) { + basis_int0[i - 1] = basis_int1[i - 1] = basis_int2[i - 1] = ITensor(s(i)); + for(int j = 1; j <= dim(s(i)); ++j) { + basis_int0[i - 1].set(s(i) = j, basis[i - 1].int0(j)); + basis_int1[i - 1].set(s(i) = j, basis[i - 1].int1(j)); + basis_int2[i - 1].set(s(i) = j, basis[i - 1].int2(j)); + } } - } - // Z = int tt(x) dx (partition function); normalize to get probability density rho - auto Z = tt(1) * basis_int0[0]; - for(int i = 2; i <= d; ++i) { - Z *= tt(i) * basis_int0[i - 1]; - } - auto rho = tt; - rho /= elt(Z); - // ei[k] = E[x_k], eii[k] = E[x_k^2], eij[k][l] = E[x_k * x_l] for k < l - std::vector ei(d), eii(d); - Matrix eij(d, d); - for(int k = 1; k <= d; ++k) { - // replace dimension k's integral vector with int1 (resp. int2) to get E[x_k] (resp. E[x_k^2]) - auto eival = rho(1) * (k == 1 ? basis_int1[0] : basis_int0[0]); - auto eiival = rho(1) * (k == 1 ? basis_int2[0] : basis_int0[0]); + // Z = int tt(x) dx (partition function); normalize to get probability density rho + auto Z = tt(1) * basis_int0[0]; for(int i = 2; i <= d; ++i) { - eival *= rho(i) * (k == i ? basis_int1[i - 1] : basis_int0[i - 1]); - eiival *= rho(i) * (k == i ? basis_int2[i - 1] : basis_int0[i - 1]); + Z *= tt(i) * basis_int0[i - 1]; } - ei[k - 1] = elt(eival); - eii[k - 1] = elt(eiival); - for(int l = k + 1; l <= d; ++l) { - // replace both dimensions k and l with int1 to get E[x_k * x_l] - auto eijval = rho(1) * (k == 1 ? basis_int1[0] : basis_int0[0]); + auto rho = tt; + rho /= elt(Z); + // ei[k] = E[x_k], eii[k] = E[x_k^2], eij[k][l] = E[x_k * x_l] for k < l + std::vector ei(d), eii(d); + Matrix eij(d, d); + for(int k = 1; k <= d; ++k) { + // replace dimension k's integral vector with int1 (resp. int2) to get E[x_k] (resp. E[x_k^2]) + auto eival = rho(1) * (k == 1 ? basis_int1[0] : basis_int0[0]); + auto eiival = rho(1) * (k == 1 ? basis_int2[0] : basis_int0[0]); for(int i = 2; i <= d; ++i) { - eijval *= rho(i) * (k == i || l == i ? basis_int1[i - 1] : basis_int0[i - 1]); + eival *= rho(i) * (k == i ? basis_int1[i - 1] : basis_int0[i - 1]); + eiival *= rho(i) * (k == i ? basis_int2[i - 1] : basis_int0[i - 1]); + } + ei[k - 1] = elt(eival); + eii[k - 1] = elt(eiival); + for(int l = k + 1; l <= d; ++l) { + // replace both dimensions k and l with int1 to get E[x_k * x_l] + auto eijval = rho(1) * (k == 1 ? basis_int1[0] : basis_int0[0]); + for(int i = 2; i <= d; ++i) { + eijval *= rho(i) * (k == i || l == i ? basis_int1[i - 1] : basis_int0[i - 1]); + } + eij(k - 1, l - 1) = elt(eijval); } - eij(k - 1, l - 1) = elt(eijval); } - } - // sigma(k,l) = Cov(x_k, x_l) = E[x_k*x_l] - E[x_k]*E[x_l] - Matrix sigma(d, d); - for(int k = 1; k <= d; ++k) { - for(int l = k; l <= d; ++l) { - sigma(k - 1, l - 1) = sigma(l - 1, k - 1) = k == l ? eii[k - 1] - pow(ei[k - 1], 2) : eij(k - 1, l - 1) - ei[k - 1] * ei[l - 1]; + // sigma(k,l) = Cov(x_k, x_l) = E[x_k*x_l] - E[x_k]*E[x_l] + Matrix sigma(d, d); + for(int k = 1; k <= d; ++k) { + for(int l = k; l <= d; ++l) { + sigma(k - 1, l - 1) = sigma(l - 1, k - 1) = k == l ? eii[k - 1] - pow(ei[k - 1], 2) : eij(k - 1, l - 1) - ei[k - 1] * ei[l - 1]; + } } + return std::make_tuple(sigma, ei, elt(Z)); } - return std::make_tuple(sigma, ei, elt(Z)); -} // Compute the 2D marginal density of the normalized TT distribution for dimensions // pos1 and pos2 on a (bins x bins) grid. All other dimensions are integrated out // analytically using int0, which collapses those TT cores to scalar factors. -void marginal2d(const MPS& tt, const std::vector& basis, int pos1, int pos2, Matrix& grid, bool conv) { - int bins = grid.nrows(); - int d = length(tt); - auto s = siteInds(tt); - std::vector basis_int0(d); - for(int i = 1; i <= d; ++i) { - basis_int0[i - 1] = ITensor(s(i)); - for(int j = 1; j <= dim(s(i)); ++j) { - basis_int0[i - 1].set(s(i) = j, basis[i - 1].int0(j)); - } - } - // normalize to probability density - auto Z = tt(1) * basis_int0[0]; - for(int i = 2; i <= d; ++i) { - Z *= tt(i) * basis_int0[i - 1]; - } - auto rho = tt; - rho /= elt(Z); - for(int i = 0; i < bins; ++i) { - for(int j = 0; j < bins; ++j) { - double x = basis[pos1 - 1].dom().first + i * (basis[pos1 - 1].dom().second - basis[pos1 - 1].dom().first) / bins; - double y = basis[pos2 - 1].dom().first + j * (basis[pos2 - 1].dom().second - basis[pos2 - 1].dom().first) / bins; - ITensor xevals(s(pos1)), yevals(s(pos2)); - for(int k = 1; k <= dim(s(pos1)); ++k) { - xevals.set(s(pos1) = k, basis[pos1 - 1](x, k, conv)); - yevals.set(s(pos2) = k, basis[pos2 - 1](y, k, conv)); + void marginal2d(const MPS& tt, const std::vector& basis, int pos1, int pos2, Matrix& grid, bool conv) { + int bins = grid.nrows(); + int d = length(tt); + auto s = siteInds(tt); + std::vector basis_int0(d); + for(int i = 1; i <= d; ++i) { + basis_int0[i - 1] = ITensor(s(i)); + for(int j = 1; j <= dim(s(i)); ++j) { + basis_int0[i - 1].set(s(i) = j, basis[i - 1].int0(j)); } - auto val = rho(1) * (pos1 == 1 ? xevals : basis_int0[0]); - for(int k = 2; k <= d; ++k) { - if(pos1 == k) { - val *= rho(k) * xevals; - } else if(pos2 == k) { - val *= rho(k) * yevals; - } else { - val *= rho(k) * basis_int0[k - 1]; + } + // normalize to probability density + auto Z = tt(1) * basis_int0[0]; + for(int i = 2; i <= d; ++i) { + Z *= tt(i) * basis_int0[i - 1]; + } + auto rho = tt; + rho /= elt(Z); + for(int i = 0; i < bins; ++i) { + for(int j = 0; j < bins; ++j) { + double x = basis[pos1 - 1].dom().first + i * (basis[pos1 - 1].dom().second - basis[pos1 - 1].dom().first) / bins; + double y = basis[pos2 - 1].dom().first + j * (basis[pos2 - 1].dom().second - basis[pos2 - 1].dom().first) / bins; + ITensor xevals(s(pos1)), yevals(s(pos2)); + for(int k = 1; k <= dim(s(pos1)); ++k) { + xevals.set(s(pos1) = k, basis[pos1 - 1](x, k, conv)); + yevals.set(s(pos2) = k, basis[pos2 - 1](y, k, conv)); + } + auto val = rho(1) * (pos1 == 1 ? xevals : basis_int0[0]); + for(int k = 2; k <= d; ++k) { + if(pos1 == k) { + val *= rho(k) * xevals; + } else if(pos2 == k) { + val *= rho(k) * yevals; + } else { + val *= rho(k) * basis_int0[k - 1]; + } } + grid(i, j) = elt(val); } - grid(i, j) = elt(val); } } -} } } From 8b1e13da451e0df02323b359fb4293694594e95a Mon Sep 17 00:00:00 2001 From: Carlo Camilloni Date: Wed, 1 Jul 2026 14:37:18 +0200 Subject: [PATCH 5/6] missing parenthesis --- src/ttsketch/TTHelper.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ttsketch/TTHelper.cpp b/src/ttsketch/TTHelper.cpp index 409736d255..6f0b180449 100644 --- a/src/ttsketch/TTHelper.cpp +++ b/src/ttsketch/TTHelper.cpp @@ -94,6 +94,7 @@ double ttEval(const MPS& tt, const std::vector& elements, bool conv) { return ttEvalAndGrad(tt, basis, elements, conv).first; +} // Compute covariance matrix, marginal means, and partition function of the TT distribution. // Precomputes ITensors for the three moment integrals int0/int1/int2 per dimension, From e856374f6eb9a7152169297abccd45babd503508 Mon Sep 17 00:00:00 2001 From: Carlo Camilloni Date: Wed, 1 Jul 2026 16:27:08 +0200 Subject: [PATCH 6/6] astyle --- src/ttsketch/TTHelper.cpp | 158 +++++++++++++++++++------------------- 1 file changed, 79 insertions(+), 79 deletions(-) diff --git a/src/ttsketch/TTHelper.cpp b/src/ttsketch/TTHelper.cpp index 6f0b180449..1730c46016 100644 --- a/src/ttsketch/TTHelper.cpp +++ b/src/ttsketch/TTHelper.cpp @@ -99,102 +99,102 @@ double ttEval(const MPS& tt, // Compute covariance matrix, marginal means, and partition function of the TT distribution. // Precomputes ITensors for the three moment integrals int0/int1/int2 per dimension, // then evaluates expectations as TT contractions with these integral vectors. - std::tuple, std::vector, double> covMat(const MPS& tt, const std::vector& basis) { - int d = length(tt); - auto s = siteInds(tt); - // integral vectors: int0[i][j] = int phi_j(x_i) dx_i, etc. - std::vector basis_int0(d), basis_int1(d), basis_int2(d); - for(int i = 1; i <= d; ++i) { - basis_int0[i - 1] = basis_int1[i - 1] = basis_int2[i - 1] = ITensor(s(i)); - for(int j = 1; j <= dim(s(i)); ++j) { - basis_int0[i - 1].set(s(i) = j, basis[i - 1].int0(j)); - basis_int1[i - 1].set(s(i) = j, basis[i - 1].int1(j)); - basis_int2[i - 1].set(s(i) = j, basis[i - 1].int2(j)); - } +std::tuple, std::vector, double> covMat(const MPS& tt, const std::vector& basis) { + int d = length(tt); + auto s = siteInds(tt); + // integral vectors: int0[i][j] = int phi_j(x_i) dx_i, etc. + std::vector basis_int0(d), basis_int1(d), basis_int2(d); + for(int i = 1; i <= d; ++i) { + basis_int0[i - 1] = basis_int1[i - 1] = basis_int2[i - 1] = ITensor(s(i)); + for(int j = 1; j <= dim(s(i)); ++j) { + basis_int0[i - 1].set(s(i) = j, basis[i - 1].int0(j)); + basis_int1[i - 1].set(s(i) = j, basis[i - 1].int1(j)); + basis_int2[i - 1].set(s(i) = j, basis[i - 1].int2(j)); } - // Z = int tt(x) dx (partition function); normalize to get probability density rho - auto Z = tt(1) * basis_int0[0]; + } + // Z = int tt(x) dx (partition function); normalize to get probability density rho + auto Z = tt(1) * basis_int0[0]; + for(int i = 2; i <= d; ++i) { + Z *= tt(i) * basis_int0[i - 1]; + } + auto rho = tt; + rho /= elt(Z); + // ei[k] = E[x_k], eii[k] = E[x_k^2], eij[k][l] = E[x_k * x_l] for k < l + std::vector ei(d), eii(d); + Matrix eij(d, d); + for(int k = 1; k <= d; ++k) { + // replace dimension k's integral vector with int1 (resp. int2) to get E[x_k] (resp. E[x_k^2]) + auto eival = rho(1) * (k == 1 ? basis_int1[0] : basis_int0[0]); + auto eiival = rho(1) * (k == 1 ? basis_int2[0] : basis_int0[0]); for(int i = 2; i <= d; ++i) { - Z *= tt(i) * basis_int0[i - 1]; + eival *= rho(i) * (k == i ? basis_int1[i - 1] : basis_int0[i - 1]); + eiival *= rho(i) * (k == i ? basis_int2[i - 1] : basis_int0[i - 1]); } - auto rho = tt; - rho /= elt(Z); - // ei[k] = E[x_k], eii[k] = E[x_k^2], eij[k][l] = E[x_k * x_l] for k < l - std::vector ei(d), eii(d); - Matrix eij(d, d); - for(int k = 1; k <= d; ++k) { - // replace dimension k's integral vector with int1 (resp. int2) to get E[x_k] (resp. E[x_k^2]) - auto eival = rho(1) * (k == 1 ? basis_int1[0] : basis_int0[0]); - auto eiival = rho(1) * (k == 1 ? basis_int2[0] : basis_int0[0]); + ei[k - 1] = elt(eival); + eii[k - 1] = elt(eiival); + for(int l = k + 1; l <= d; ++l) { + // replace both dimensions k and l with int1 to get E[x_k * x_l] + auto eijval = rho(1) * (k == 1 ? basis_int1[0] : basis_int0[0]); for(int i = 2; i <= d; ++i) { - eival *= rho(i) * (k == i ? basis_int1[i - 1] : basis_int0[i - 1]); - eiival *= rho(i) * (k == i ? basis_int2[i - 1] : basis_int0[i - 1]); - } - ei[k - 1] = elt(eival); - eii[k - 1] = elt(eiival); - for(int l = k + 1; l <= d; ++l) { - // replace both dimensions k and l with int1 to get E[x_k * x_l] - auto eijval = rho(1) * (k == 1 ? basis_int1[0] : basis_int0[0]); - for(int i = 2; i <= d; ++i) { - eijval *= rho(i) * (k == i || l == i ? basis_int1[i - 1] : basis_int0[i - 1]); - } - eij(k - 1, l - 1) = elt(eijval); + eijval *= rho(i) * (k == i || l == i ? basis_int1[i - 1] : basis_int0[i - 1]); } + eij(k - 1, l - 1) = elt(eijval); } - // sigma(k,l) = Cov(x_k, x_l) = E[x_k*x_l] - E[x_k]*E[x_l] - Matrix sigma(d, d); - for(int k = 1; k <= d; ++k) { - for(int l = k; l <= d; ++l) { - sigma(k - 1, l - 1) = sigma(l - 1, k - 1) = k == l ? eii[k - 1] - pow(ei[k - 1], 2) : eij(k - 1, l - 1) - ei[k - 1] * ei[l - 1]; - } + } + // sigma(k,l) = Cov(x_k, x_l) = E[x_k*x_l] - E[x_k]*E[x_l] + Matrix sigma(d, d); + for(int k = 1; k <= d; ++k) { + for(int l = k; l <= d; ++l) { + sigma(k - 1, l - 1) = sigma(l - 1, k - 1) = k == l ? eii[k - 1] - pow(ei[k - 1], 2) : eij(k - 1, l - 1) - ei[k - 1] * ei[l - 1]; } - return std::make_tuple(sigma, ei, elt(Z)); } + return std::make_tuple(sigma, ei, elt(Z)); +} // Compute the 2D marginal density of the normalized TT distribution for dimensions // pos1 and pos2 on a (bins x bins) grid. All other dimensions are integrated out // analytically using int0, which collapses those TT cores to scalar factors. - void marginal2d(const MPS& tt, const std::vector& basis, int pos1, int pos2, Matrix& grid, bool conv) { - int bins = grid.nrows(); - int d = length(tt); - auto s = siteInds(tt); - std::vector basis_int0(d); - for(int i = 1; i <= d; ++i) { - basis_int0[i - 1] = ITensor(s(i)); - for(int j = 1; j <= dim(s(i)); ++j) { - basis_int0[i - 1].set(s(i) = j, basis[i - 1].int0(j)); - } - } - // normalize to probability density - auto Z = tt(1) * basis_int0[0]; - for(int i = 2; i <= d; ++i) { - Z *= tt(i) * basis_int0[i - 1]; +void marginal2d(const MPS& tt, const std::vector& basis, int pos1, int pos2, Matrix& grid, bool conv) { + int bins = grid.nrows(); + int d = length(tt); + auto s = siteInds(tt); + std::vector basis_int0(d); + for(int i = 1; i <= d; ++i) { + basis_int0[i - 1] = ITensor(s(i)); + for(int j = 1; j <= dim(s(i)); ++j) { + basis_int0[i - 1].set(s(i) = j, basis[i - 1].int0(j)); } - auto rho = tt; - rho /= elt(Z); - for(int i = 0; i < bins; ++i) { - for(int j = 0; j < bins; ++j) { - double x = basis[pos1 - 1].dom().first + i * (basis[pos1 - 1].dom().second - basis[pos1 - 1].dom().first) / bins; - double y = basis[pos2 - 1].dom().first + j * (basis[pos2 - 1].dom().second - basis[pos2 - 1].dom().first) / bins; - ITensor xevals(s(pos1)), yevals(s(pos2)); - for(int k = 1; k <= dim(s(pos1)); ++k) { - xevals.set(s(pos1) = k, basis[pos1 - 1](x, k, conv)); - yevals.set(s(pos2) = k, basis[pos2 - 1](y, k, conv)); - } - auto val = rho(1) * (pos1 == 1 ? xevals : basis_int0[0]); - for(int k = 2; k <= d; ++k) { - if(pos1 == k) { - val *= rho(k) * xevals; - } else if(pos2 == k) { - val *= rho(k) * yevals; - } else { - val *= rho(k) * basis_int0[k - 1]; - } + } + // normalize to probability density + auto Z = tt(1) * basis_int0[0]; + for(int i = 2; i <= d; ++i) { + Z *= tt(i) * basis_int0[i - 1]; + } + auto rho = tt; + rho /= elt(Z); + for(int i = 0; i < bins; ++i) { + for(int j = 0; j < bins; ++j) { + double x = basis[pos1 - 1].dom().first + i * (basis[pos1 - 1].dom().second - basis[pos1 - 1].dom().first) / bins; + double y = basis[pos2 - 1].dom().first + j * (basis[pos2 - 1].dom().second - basis[pos2 - 1].dom().first) / bins; + ITensor xevals(s(pos1)), yevals(s(pos2)); + for(int k = 1; k <= dim(s(pos1)); ++k) { + xevals.set(s(pos1) = k, basis[pos1 - 1](x, k, conv)); + yevals.set(s(pos2) = k, basis[pos2 - 1](y, k, conv)); + } + auto val = rho(1) * (pos1 == 1 ? xevals : basis_int0[0]); + for(int k = 2; k <= d; ++k) { + if(pos1 == k) { + val *= rho(k) * xevals; + } else if(pos2 == k) { + val *= rho(k) * yevals; + } else { + val *= rho(k) * basis_int0[k - 1]; } - grid(i, j) = elt(val); } + grid(i, j) = elt(val); } } +} } }