diff --git a/docs/source/api/methods.rst b/docs/source/api/methods.rst index 08b4185d..df65752b 100644 --- a/docs/source/api/methods.rst +++ b/docs/source/api/methods.rst @@ -60,6 +60,15 @@ Motion profiles :no-special-members: +Steady State Detection +^^^^^^^^^^^^^^^^^^^^^^ + +.. autoapimodule:: steady_state_detector + :members: + :no-private-members: + :no-special-members: + + Utilities ------ diff --git a/notebooks/demo-data/bottleneck/rho_v_Voronoi_traj_AO_b240.txt_id_1.dat b/notebooks/demo-data/bottleneck/rho_v_Voronoi_traj_AO_b240.txt_id_1.dat new file mode 100644 index 00000000..4856600e --- /dev/null +++ b/notebooks/demo-data/bottleneck/rho_v_Voronoi_traj_AO_b240.txt_id_1.dat @@ -0,0 +1,1074 @@ +#Frame Voronoi density(m^(-2)) Voronoi velocity(m/s) +29 1.021 0.034 +30 1.054 0.043 +31 1.058 0.056 +32 1.064 0.067 +33 1.071 0.079 +34 1.076 0.065 +35 1.081 0.083 +36 1.090 0.108 +37 1.093 0.135 +38 1.101 0.174 +39 1.111 0.224 +40 1.129 0.287 +41 1.150 0.363 +42 1.180 0.448 +43 1.219 0.556 +44 1.266 0.683 +45 1.325 0.831 +46 1.376 0.960 +47 1.413 1.077 +48 1.439 1.177 +49 1.480 1.247 +50 1.542 1.299 +51 1.599 1.337 +52 1.663 1.361 +53 1.732 1.367 +54 1.808 1.366 +55 1.894 1.363 +56 1.975 1.359 +57 2.050 1.356 +58 2.111 1.350 +59 2.151 1.343 +60 2.175 1.330 +61 2.186 1.312 +62 2.196 1.301 +63 2.216 1.287 +64 2.212 1.264 +65 2.214 1.249 +66 2.210 1.235 +67 2.193 1.222 +68 2.167 1.223 +69 2.156 1.221 +70 2.168 1.212 +71 2.197 1.201 +72 2.248 1.185 +73 2.307 1.173 +74 2.359 1.165 +75 2.398 1.159 +76 2.424 1.158 +77 2.442 1.157 +78 2.457 1.155 +79 2.469 1.154 +80 2.482 1.150 +81 2.506 1.140 +82 2.538 1.127 +83 2.569 1.114 +84 2.592 1.101 +85 2.595 1.101 +86 2.618 1.100 +87 2.657 1.094 +88 2.605 1.084 +89 2.645 1.072 +90 2.694 1.059 +91 2.751 1.046 +92 2.819 1.031 +93 2.905 1.017 +94 2.993 1.003 +95 3.069 0.990 +96 3.130 0.979 +97 3.132 0.974 +98 3.148 0.974 +99 3.164 0.978 +100 3.187 0.983 +101 3.206 0.988 +102 3.213 0.992 +103 3.208 0.994 +104 3.209 0.992 +105 3.205 0.985 +106 3.191 0.975 +107 3.182 0.961 +108 3.174 0.941 +109 3.159 0.917 +110 3.138 0.899 +111 3.095 0.894 +112 3.057 0.895 +113 3.045 0.899 +114 3.053 0.907 +115 3.078 0.917 +116 3.098 0.922 +117 3.102 0.925 +118 3.108 0.924 +119 3.130 0.918 +120 3.158 0.909 +121 3.184 0.900 +122 3.206 0.890 +123 3.163 0.879 +124 3.121 0.876 +125 3.084 0.871 +126 3.045 0.870 +127 3.006 0.871 +128 2.983 0.878 +129 2.966 0.883 +130 2.954 0.889 +131 2.944 0.895 +132 2.931 0.901 +133 2.923 0.903 +134 2.906 0.903 +135 2.898 0.895 +136 2.887 0.881 +137 2.880 0.869 +138 2.878 0.859 +139 2.898 0.859 +140 2.936 0.858 +141 2.973 0.861 +142 3.002 0.869 +143 2.971 0.882 +144 2.933 0.897 +145 2.919 0.906 +146 2.922 0.908 +147 2.932 0.903 +148 2.955 0.896 +149 2.985 0.888 +150 2.989 0.881 +151 2.983 0.875 +152 2.991 0.872 +153 2.999 0.870 +154 3.011 0.870 +155 3.030 0.871 +156 3.060 0.874 +157 3.085 0.881 +158 3.114 0.889 +159 3.138 0.893 +160 3.171 0.896 +161 3.215 0.896 +162 3.274 0.899 +163 3.347 0.902 +164 3.430 0.901 +165 3.512 0.897 +166 3.580 0.890 +167 3.655 0.883 +168 3.705 0.877 +169 3.733 0.875 +170 3.742 0.873 +171 3.725 0.872 +172 3.729 0.869 +173 3.739 0.867 +174 3.775 0.865 +175 3.830 0.867 +176 3.904 0.865 +177 3.982 0.861 +178 4.047 0.856 +179 4.090 0.846 +180 4.136 0.836 +181 4.169 0.827 +182 4.173 0.818 +183 4.158 0.810 +184 4.132 0.805 +185 4.099 0.797 +186 4.056 0.787 +187 3.996 0.779 +188 3.927 0.774 +189 3.864 0.775 +190 3.811 0.775 +191 3.771 0.779 +192 3.762 0.787 +193 3.752 0.802 +194 3.707 0.817 +195 3.652 0.829 +196 3.611 0.838 +197 3.583 0.835 +198 3.560 0.826 +199 3.536 0.818 +200 3.466 0.810 +201 3.416 0.800 +202 3.388 0.795 +203 3.380 0.792 +204 3.369 0.796 +205 3.355 0.804 +206 3.350 0.810 +207 3.331 0.821 +208 3.300 0.825 +209 3.277 0.825 +210 3.259 0.823 +211 3.245 0.817 +212 3.248 0.806 +213 3.265 0.794 +214 3.289 0.785 +215 3.324 0.778 +216 3.359 0.778 +217 3.391 0.782 +218 3.412 0.791 +219 3.422 0.797 +220 3.416 0.801 +221 3.404 0.803 +222 3.405 0.806 +223 3.416 0.810 +224 3.443 0.809 +225 3.483 0.807 +226 3.482 0.803 +227 3.498 0.798 +228 3.530 0.791 +229 3.569 0.779 +230 3.595 0.763 +231 3.607 0.745 +232 3.597 0.730 +233 3.545 0.723 +234 3.503 0.719 +235 3.473 0.720 +236 3.451 0.722 +237 3.441 0.723 +238 3.435 0.725 +239 3.439 0.728 +240 3.444 0.730 +241 3.447 0.729 +242 3.449 0.728 +243 3.451 0.726 +244 3.435 0.724 +245 3.419 0.722 +246 3.409 0.720 +247 3.403 0.718 +248 3.412 0.709 +249 3.431 0.698 +250 3.457 0.689 +251 3.406 0.684 +252 3.384 0.681 +253 3.369 0.680 +254 3.352 0.683 +255 3.341 0.682 +256 3.318 0.687 +257 3.301 0.702 +258 3.289 0.720 +259 3.280 0.738 +260 3.268 0.749 +261 3.275 0.756 +262 3.293 0.752 +263 3.321 0.744 +264 3.358 0.736 +265 3.394 0.733 +266 3.428 0.729 +267 3.464 0.724 +268 3.495 0.721 +269 3.512 0.720 +270 3.509 0.721 +271 3.466 0.719 +272 3.424 0.719 +273 3.390 0.716 +274 3.379 0.713 +275 3.384 0.711 +276 3.391 0.712 +277 3.400 0.717 +278 3.422 0.724 +279 3.455 0.732 +280 3.422 0.738 +281 3.416 0.743 +282 3.390 0.744 +283 3.371 0.739 +284 3.355 0.731 +285 3.333 0.727 +286 3.312 0.724 +287 3.298 0.721 +288 3.293 0.717 +289 3.299 0.712 +290 3.309 0.712 +291 3.326 0.715 +292 3.339 0.721 +293 3.363 0.726 +294 3.392 0.726 +295 3.420 0.720 +296 3.453 0.711 +297 3.497 0.702 +298 3.536 0.692 +299 3.548 0.684 +300 3.534 0.681 +301 3.536 0.678 +302 3.553 0.674 +303 3.587 0.675 +304 3.625 0.677 +305 3.657 0.684 +306 3.695 0.692 +307 3.749 0.701 +308 3.806 0.710 +309 3.864 0.718 +310 3.922 0.724 +311 3.976 0.726 +312 4.006 0.726 +313 4.022 0.722 +314 4.029 0.724 +315 4.034 0.724 +316 4.035 0.722 +317 4.043 0.727 +318 4.039 0.739 +319 4.033 0.750 +320 4.020 0.754 +321 4.012 0.755 +322 3.982 0.749 +323 3.946 0.743 +324 3.888 0.735 +325 3.828 0.730 +326 3.785 0.723 +327 3.758 0.715 +328 3.737 0.707 +329 3.719 0.702 +330 3.703 0.702 +331 3.695 0.701 +332 3.702 0.702 +333 3.706 0.706 +334 3.718 0.709 +335 3.727 0.714 +336 3.727 0.716 +337 3.688 0.719 +338 3.653 0.720 +339 3.627 0.719 +340 3.596 0.714 +341 3.575 0.707 +342 3.562 0.700 +343 3.559 0.698 +344 3.562 0.698 +345 3.569 0.700 +346 3.545 0.700 +347 3.529 0.699 +348 3.518 0.696 +349 3.511 0.695 +350 3.486 0.695 +351 3.478 0.701 +352 3.460 0.707 +353 3.434 0.713 +354 3.400 0.720 +355 3.360 0.725 +356 3.329 0.730 +357 3.312 0.730 +358 3.307 0.732 +359 3.321 0.733 +360 3.354 0.733 +361 3.393 0.733 +362 3.439 0.731 +363 3.474 0.730 +364 3.477 0.729 +365 3.478 0.732 +366 3.468 0.738 +367 3.446 0.743 +368 3.392 0.750 +369 3.356 0.755 +370 3.333 0.759 +371 3.328 0.758 +372 3.342 0.754 +373 3.362 0.751 +374 3.399 0.750 +375 3.449 0.748 +376 3.515 0.748 +377 3.569 0.751 +378 3.605 0.750 +379 3.621 0.747 +380 3.628 0.736 +381 3.618 0.727 +382 3.595 0.721 +383 3.571 0.715 +384 3.551 0.712 +385 3.548 0.709 +386 3.545 0.711 +387 3.540 0.709 +388 3.530 0.709 +389 3.521 0.709 +390 3.518 0.711 +391 3.507 0.715 +392 3.498 0.720 +393 3.486 0.725 +394 3.470 0.732 +395 3.453 0.733 +396 3.434 0.732 +397 3.395 0.729 +398 3.352 0.722 +399 3.331 0.716 +400 3.326 0.714 +401 3.331 0.713 +402 3.342 0.714 +403 3.361 0.713 +404 3.384 0.715 +405 3.420 0.720 +406 3.465 0.722 +407 3.505 0.722 +408 3.532 0.724 +409 3.551 0.728 +410 3.506 0.726 +411 3.501 0.725 +412 3.500 0.729 +413 3.504 0.737 +414 3.523 0.741 +415 3.559 0.744 +416 3.605 0.743 +417 3.657 0.745 +418 3.719 0.747 +419 3.784 0.744 +420 3.823 0.740 +421 3.797 0.737 +422 3.784 0.736 +423 3.779 0.733 +424 3.788 0.730 +425 3.812 0.724 +426 3.841 0.722 +427 3.869 0.719 +428 3.898 0.716 +429 3.922 0.715 +430 3.930 0.715 +431 3.916 0.717 +432 3.894 0.717 +433 3.873 0.721 +434 3.858 0.723 +435 3.838 0.721 +436 3.770 0.718 +437 3.719 0.711 +438 3.680 0.706 +439 3.666 0.710 +440 3.670 0.717 +441 3.686 0.724 +442 3.716 0.731 +443 3.774 0.736 +444 3.842 0.739 +445 3.890 0.738 +446 3.931 0.735 +447 3.958 0.729 +448 3.946 0.726 +449 3.901 0.717 +450 3.843 0.710 +451 3.792 0.700 +452 3.724 0.692 +453 3.660 0.689 +454 3.617 0.689 +455 3.614 0.693 +456 3.634 0.698 +457 3.661 0.706 +458 3.696 0.711 +459 3.733 0.713 +460 3.770 0.711 +461 3.814 0.709 +462 3.845 0.707 +463 3.850 0.705 +464 3.834 0.701 +465 3.815 0.696 +466 3.793 0.690 +467 3.770 0.681 +468 3.757 0.672 +469 3.749 0.663 +470 3.740 0.660 +471 3.721 0.659 +472 3.696 0.666 +473 3.671 0.675 +474 3.652 0.687 +475 3.647 0.697 +476 3.650 0.706 +477 3.652 0.712 +478 3.649 0.718 +479 3.643 0.722 +480 3.616 0.722 +481 3.593 0.719 +482 3.573 0.710 +483 3.530 0.704 +484 3.483 0.700 +485 3.424 0.700 +486 3.375 0.703 +487 3.366 0.709 +488 3.398 0.716 +489 3.430 0.735 +490 3.461 0.753 +491 3.508 0.768 +492 3.553 0.778 +493 3.600 0.782 +494 3.647 0.781 +495 3.682 0.776 +496 3.704 0.769 +497 3.726 0.761 +498 3.741 0.752 +499 3.747 0.746 +500 3.745 0.740 +501 3.733 0.736 +502 3.718 0.730 +503 3.671 0.721 +504 3.605 0.716 +505 3.535 0.709 +506 3.478 0.701 +507 3.440 0.693 +508 3.421 0.687 +509 3.409 0.681 +510 3.398 0.676 +511 3.396 0.671 +512 3.392 0.665 +513 3.383 0.660 +514 3.345 0.660 +515 3.322 0.662 +516 3.308 0.663 +517 3.309 0.663 +518 3.322 0.667 +519 3.314 0.673 +520 3.292 0.679 +521 3.271 0.688 +522 3.265 0.691 +523 3.260 0.690 +524 3.253 0.693 +525 3.252 0.692 +526 3.263 0.693 +527 3.291 0.690 +528 3.323 0.685 +529 3.374 0.679 +530 3.439 0.675 +531 3.495 0.674 +532 3.555 0.678 +533 3.616 0.685 +534 3.687 0.688 +535 3.763 0.687 +536 3.833 0.688 +537 3.900 0.689 +538 3.966 0.689 +539 4.032 0.689 +540 4.095 0.686 +541 4.136 0.685 +542 4.145 0.687 +543 4.131 0.688 +544 4.112 0.691 +545 4.092 0.692 +546 4.071 0.689 +547 4.043 0.686 +548 4.019 0.685 +549 4.006 0.686 +550 4.021 0.692 +551 4.053 0.697 +552 4.072 0.706 +553 3.963 0.719 +554 3.935 0.729 +555 3.890 0.738 +556 3.838 0.740 +557 3.778 0.741 +558 3.729 0.736 +559 3.682 0.730 +560 3.654 0.727 +561 3.639 0.726 +562 3.639 0.723 +563 3.648 0.715 +564 3.666 0.710 +565 3.683 0.702 +566 3.693 0.696 +567 3.708 0.690 +568 3.721 0.689 +569 3.724 0.688 +570 3.724 0.685 +571 3.715 0.685 +572 3.714 0.685 +573 3.720 0.687 +574 3.741 0.693 +575 3.779 0.702 +576 3.824 0.711 +577 3.811 0.719 +578 3.837 0.728 +579 3.887 0.733 +580 3.950 0.734 +581 4.013 0.736 +582 4.067 0.736 +583 4.113 0.738 +584 4.135 0.739 +585 4.146 0.738 +586 4.144 0.736 +587 4.140 0.732 +588 4.134 0.726 +589 4.116 0.721 +590 4.099 0.717 +591 4.076 0.712 +592 4.047 0.710 +593 3.992 0.711 +594 3.930 0.713 +595 3.876 0.714 +596 3.833 0.717 +597 3.809 0.716 +598 3.798 0.717 +599 3.796 0.720 +600 3.794 0.725 +601 3.789 0.726 +602 3.781 0.726 +603 3.772 0.726 +604 3.762 0.727 +605 3.742 0.727 +606 3.713 0.726 +607 3.657 0.721 +608 3.577 0.717 +609 3.502 0.711 +610 3.431 0.709 +611 3.368 0.707 +612 3.316 0.706 +613 3.278 0.708 +614 3.256 0.707 +615 3.252 0.710 +616 3.257 0.714 +617 3.274 0.716 +618 3.292 0.715 +619 3.309 0.714 +620 3.316 0.707 +621 3.317 0.701 +622 3.301 0.695 +623 3.289 0.695 +624 3.278 0.697 +625 3.275 0.697 +626 3.286 0.695 +627 3.304 0.693 +628 3.332 0.690 +629 3.368 0.686 +630 3.408 0.685 +631 3.436 0.681 +632 3.451 0.679 +633 3.461 0.674 +634 3.473 0.672 +635 3.460 0.671 +636 3.448 0.669 +637 3.433 0.672 +638 3.411 0.679 +639 3.376 0.687 +640 3.337 0.693 +641 3.296 0.698 +642 3.264 0.703 +643 3.241 0.707 +644 3.235 0.711 +645 3.234 0.712 +646 3.239 0.713 +647 3.242 0.715 +648 3.246 0.713 +649 3.235 0.707 +650 3.195 0.704 +651 3.171 0.701 +652 3.147 0.698 +653 3.146 0.698 +654 3.112 0.700 +655 3.083 0.706 +656 3.068 0.709 +657 3.070 0.712 +658 3.088 0.717 +659 3.119 0.724 +660 3.155 0.727 +661 3.190 0.730 +662 3.224 0.736 +663 3.260 0.743 +664 3.295 0.748 +665 3.329 0.749 +666 3.354 0.751 +667 3.365 0.747 +668 3.371 0.742 +669 3.382 0.739 +670 3.403 0.740 +671 3.426 0.741 +672 3.447 0.740 +673 3.454 0.738 +674 3.454 0.737 +675 3.455 0.739 +676 3.445 0.743 +677 3.444 0.746 +678 3.438 0.747 +679 3.405 0.751 +680 3.362 0.753 +681 3.330 0.756 +682 3.305 0.758 +683 3.282 0.759 +684 3.258 0.760 +685 3.239 0.760 +686 3.239 0.762 +687 3.255 0.763 +688 3.279 0.763 +689 3.301 0.764 +690 3.319 0.767 +691 3.336 0.769 +692 3.359 0.772 +693 3.363 0.778 +694 3.362 0.787 +695 3.373 0.797 +696 3.363 0.803 +697 3.367 0.806 +698 3.365 0.811 +699 3.355 0.814 +700 3.353 0.818 +701 3.358 0.823 +702 3.379 0.824 +703 3.420 0.819 +704 3.470 0.812 +705 3.524 0.807 +706 3.576 0.801 +707 3.632 0.794 +708 3.623 0.786 +709 3.615 0.778 +710 3.607 0.773 +711 3.590 0.771 +712 3.562 0.769 +713 3.533 0.765 +714 3.511 0.755 +715 3.485 0.740 +716 3.457 0.719 +717 3.437 0.704 +718 3.428 0.696 +719 3.419 0.693 +720 3.406 0.691 +721 3.382 0.692 +722 3.353 0.696 +723 3.272 0.701 +724 3.204 0.708 +725 3.158 0.716 +726 3.125 0.725 +727 3.104 0.728 +728 3.089 0.724 +729 3.072 0.719 +730 3.054 0.713 +731 3.046 0.707 +732 3.053 0.700 +733 3.070 0.695 +734 3.020 0.690 +735 3.009 0.689 +736 3.001 0.688 +737 2.993 0.690 +738 2.977 0.695 +739 2.952 0.699 +740 2.932 0.701 +741 2.924 0.704 +742 2.928 0.703 +743 2.932 0.702 +744 2.939 0.700 +745 2.951 0.699 +746 2.960 0.697 +747 2.960 0.686 +748 2.955 0.683 +749 2.958 0.680 +750 2.964 0.682 +751 2.975 0.680 +752 2.644 0.660 +753 2.685 0.660 +754 2.725 0.667 +755 2.759 0.673 +756 2.791 0.679 +757 2.807 0.688 +758 2.819 0.694 +759 2.829 0.701 +760 2.835 0.708 +761 2.846 0.714 +762 2.855 0.721 +763 2.869 0.728 +764 2.885 0.732 +765 2.910 0.735 +766 2.952 0.734 +767 3.007 0.728 +768 3.068 0.726 +769 3.132 0.727 +770 3.195 0.728 +771 3.253 0.727 +772 3.288 0.730 +773 3.286 0.733 +774 3.277 0.743 +775 3.270 0.757 +776 3.270 0.774 +777 3.283 0.790 +778 3.308 0.799 +779 3.348 0.801 +780 3.398 0.797 +781 3.443 0.788 +782 3.475 0.774 +783 3.490 0.762 +784 3.496 0.746 +785 3.495 0.730 +786 3.486 0.715 +787 3.470 0.703 +788 3.444 0.699 +789 3.404 0.701 +790 3.364 0.707 +791 3.330 0.716 +792 3.304 0.723 +793 3.292 0.730 +794 3.178 0.738 +795 3.127 0.745 +796 3.088 0.750 +797 3.054 0.751 +798 3.021 0.749 +799 3.003 0.748 +800 2.999 0.748 +801 3.013 0.746 +802 3.035 0.750 +803 3.060 0.754 +804 3.085 0.753 +805 3.102 0.750 +806 3.109 0.740 +807 3.112 0.728 +808 3.127 0.715 +809 3.148 0.702 +810 3.197 0.693 +811 3.190 0.692 +812 3.202 0.698 +813 3.232 0.705 +814 3.150 0.710 +815 3.176 0.713 +816 3.162 0.715 +817 3.127 0.716 +818 3.088 0.718 +819 3.054 0.723 +820 3.023 0.726 +821 3.013 0.728 +822 3.022 0.726 +823 3.036 0.721 +824 3.049 0.716 +825 3.068 0.713 +826 3.080 0.711 +827 3.097 0.710 +828 3.122 0.712 +829 3.147 0.720 +830 3.113 0.725 +831 3.088 0.728 +832 3.071 0.733 +833 3.082 0.741 +834 3.106 0.750 +835 3.146 0.758 +836 3.193 0.760 +837 3.249 0.758 +838 3.303 0.752 +839 3.350 0.739 +840 3.391 0.728 +841 3.427 0.719 +842 3.463 0.713 +843 3.489 0.707 +844 3.097 0.699 +845 3.085 0.697 +846 3.077 0.695 +847 3.078 0.690 +848 3.087 0.682 +849 3.111 0.674 +850 3.138 0.667 +851 3.164 0.662 +852 3.182 0.658 +853 3.191 0.659 +854 3.204 0.662 +855 3.219 0.665 +856 3.239 0.672 +857 3.269 0.682 +858 3.300 0.694 +859 3.321 0.706 +860 3.322 0.721 +861 3.318 0.740 +862 3.302 0.759 +863 3.258 0.775 +864 3.205 0.785 +865 3.148 0.789 +866 3.085 0.787 +867 3.028 0.780 +868 2.978 0.771 +869 2.936 0.755 +870 2.778 0.743 +871 2.731 0.734 +872 2.697 0.730 +873 2.674 0.727 +874 2.662 0.726 +875 2.677 0.731 +876 2.657 0.741 +877 2.632 0.750 +878 2.604 0.754 +879 2.584 0.759 +880 2.565 0.760 +881 2.533 0.759 +882 2.491 0.757 +883 2.446 0.755 +884 2.411 0.752 +885 2.379 0.751 +886 2.348 0.750 +887 2.311 0.755 +888 2.273 0.763 +889 2.232 0.768 +890 2.185 0.774 +891 2.144 0.781 +892 2.110 0.785 +893 2.086 0.785 +894 2.066 0.780 +895 2.053 0.772 +896 2.033 0.762 +897 2.012 0.751 +898 1.981 0.741 +899 1.950 0.739 +900 1.927 0.737 +901 1.918 0.738 +902 1.924 0.741 +903 1.943 0.750 +904 1.973 0.761 +905 2.015 0.771 +906 2.059 0.779 +907 2.099 0.788 +908 2.135 0.797 +909 2.161 0.803 +910 2.184 0.808 +911 2.203 0.810 +912 2.221 0.810 +913 2.238 0.804 +914 2.258 0.796 +915 2.235 0.786 +916 2.242 0.779 +917 2.250 0.778 +918 2.252 0.781 +919 2.251 0.786 +920 2.247 0.784 +921 2.196 0.779 +922 2.152 0.772 +923 2.117 0.762 +924 2.097 0.750 +925 2.084 0.738 +926 2.067 0.726 +927 2.056 0.715 +928 2.046 0.708 +929 2.031 0.707 +930 2.022 0.709 +931 2.022 0.711 +932 2.029 0.712 +933 2.036 0.707 +934 2.048 0.699 +935 2.056 0.691 +936 2.064 0.684 +937 2.066 0.684 +938 2.073 0.685 +939 2.082 0.680 +940 2.105 0.679 +941 2.112 0.684 +942 2.179 0.695 +943 2.183 0.711 +944 2.194 0.729 +945 2.206 0.745 +946 2.209 0.756 +947 2.205 0.764 +948 2.189 0.770 +949 2.165 0.773 +950 2.135 0.775 +951 2.099 0.774 +952 2.068 0.770 +953 2.042 0.767 +954 2.022 0.762 +955 2.014 0.756 +956 2.017 0.753 +957 2.021 0.756 +958 2.025 0.764 +959 2.030 0.775 +960 2.040 0.784 +961 2.073 0.790 +962 2.076 0.793 +963 2.078 0.794 +964 2.071 0.790 +965 2.055 0.788 +966 2.035 0.788 +967 2.012 0.790 +968 1.982 0.791 +969 1.935 0.790 +970 1.888 0.788 +971 1.845 0.784 +972 1.805 0.781 +973 1.768 0.780 +974 1.740 0.779 +975 1.712 0.777 +976 1.687 0.772 +977 1.668 0.766 +978 1.652 0.762 +979 1.633 0.759 +980 1.607 0.756 +981 1.577 0.755 +982 1.545 0.753 +983 1.517 0.749 +984 1.493 0.746 +985 1.475 0.745 +986 1.465 0.751 +987 1.460 0.758 +988 1.458 0.764 +989 1.456 0.774 +990 1.459 0.779 +991 1.452 0.776 +992 1.430 0.768 +993 1.399 0.758 +994 1.354 0.750 +995 1.307 0.747 +996 1.258 0.743 +997 1.208 0.745 +998 1.160 0.750 +999 1.114 0.754 +1000 1.069 0.760 +1001 1.025 0.769 +1002 0.984 0.778 +1003 0.950 0.785 +1004 0.916 0.784 +1005 0.885 0.782 +1006 0.855 0.779 +1007 0.826 0.775 +1008 0.799 0.772 +1009 0.771 0.777 +1010 0.746 0.781 +1011 0.723 0.780 +1012 0.704 0.777 +1013 0.688 0.773 +1014 0.675 0.765 +1015 0.660 0.750 +1016 0.640 0.735 +1017 0.615 0.718 +1018 0.590 0.701 +1019 0.567 0.686 +1020 0.544 0.673 +1021 0.516 0.661 +1022 0.489 0.642 +1023 0.460 0.622 +1024 0.431 0.598 +1025 0.396 0.565 +1026 0.361 0.528 +1027 0.329 0.494 +1028 0.295 0.455 +1029 0.260 0.416 +1030 0.228 0.377 +1031 0.197 0.341 +1032 0.168 0.302 +1033 0.140 0.260 +1034 0.113 0.215 +1035 0.087 0.169 +1036 0.063 0.124 +1037 0.041 0.084 +1038 0.028 0.056 +1039 0.016 0.033 +1040 0.007 0.013 +1041 0.000 0.000 +1042 0.000 0.000 +1043 0.000 0.000 +1044 0.000 0.000 +1045 0.000 0.000 +1046 0.000 0.000 +1047 0.000 0.000 +1048 0.000 0.000 +1049 0.000 0.000 +1050 0.000 0.000 +1051 0.000 0.000 +1052 0.000 0.000 +1053 0.000 0.000 +1054 0.000 0.000 +1055 0.000 0.000 +1056 0.000 0.000 +1057 0.000 0.000 +1058 0.000 0.000 +1059 0.000 0.000 +1060 0.000 0.000 +1061 0.000 0.000 +1062 0.000 0.000 +1063 0.000 0.000 +1064 0.000 0.000 +1065 0.000 0.000 +1066 0.000 0.000 +1067 0.000 0.000 +1068 0.000 0.000 +1069 0.000 0.000 +1070 0.000 0.000 +1071 0.000 0.000 +1072 0.000 0.000 +1073 0.000 0.000 +1074 0.000 0.000 +1075 0.000 0.000 +1076 0.000 0.000 +1077 0.000 0.000 +1078 0.000 0.000 +1079 0.000 0.000 +1080 0.000 0.000 +1081 0.000 0.000 +1082 0.000 0.000 +1083 0.000 0.000 +1084 0.000 0.000 +1085 0.000 0.000 +1086 0.000 0.000 +1087 0.000 0.000 +1088 0.000 0.000 +1089 0.000 0.000 +1090 0.000 0.000 +1091 0.000 0.000 +1092 0.000 0.000 +1093 0.000 0.000 +1094 0.000 0.000 +1095 0.000 0.000 +1096 0.000 0.000 +1097 0.000 0.000 +1098 0.000 0.000 +1099 0.000 0.000 +1100 0.000 0.000 +1101 0.000 0.000 diff --git a/notebooks/user_guide.ipynb b/notebooks/user_guide.ipynb index 242e4eef..89718900 100644 --- a/notebooks/user_guide.ipynb +++ b/notebooks/user_guide.ipynb @@ -4529,6 +4529,187 @@ "plt.show()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Steady State Detection\n", + "\n", + "When analysing pedestrian experiments, it is important to distinguish between transient and steady states.\n", + "The transient phase (at the beginning and end of an experiment) can strongly influence the results.\n", + "*PedPy* provides a steady state detection method based on a modified CUSUM (Cumulative Sum Control Chart)\n", + "algorithm from [Liao et al. (2016)](https://doi.org/10.1016/j.physa.2016.05.051).\n", + "\n", + "The algorithm works on time series of density and speed:\n", + "\n", + "1. Select a **reference process** — a sub-range that visually looks like steady state\n", + "2. The algorithm computes CUSUM statistics measuring deviations from this reference\n", + "3. Intervals where the statistics drop below a calibrated threshold are detected as steady state\n", + "4. For flow analysis, steady states of density and speed are **combined** (overlapping intervals)\n", + "\n", + "#### Load density/speed time series\n", + "\n", + "We use Voronoi density and speed data from a bottleneck experiment (experiment AO, b=240cm).\n", + "The data file has three columns: frame, density, and speed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "data = np.loadtxt(\n", + " \"demo-data/bottleneck/rho_v_Voronoi_traj_AO_b240.txt_id_1.dat\",\n", + " usecols=[0, 1, 2],\n", + ")\n", + "# Remove rows where density is zero (no measurement)\n", + "data = data[data[:, 1] != 0]\n", + "\n", + "frames = data[:, 0]\n", + "density = data[:, 1]\n", + "speed = data[:, 2]\n", + "\n", + "print(f\"Frames: {frames[0]:.0f} to {frames[-1]:.0f} ({len(frames)} observations)\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Detect steady state for density and speed\n", + "\n", + "We select a reference process from frames 240 to 640 (a visually steady region)\n", + "and run the detection independently for density and speed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pedpy import combine_steady_states, detect_steady_state\n", + "\n", + "result_density = detect_steady_state(\n", + " frames=frames,\n", + " values=density,\n", + " ref_start=240,\n", + " ref_end=640,\n", + ")\n", + "\n", + "result_speed = detect_steady_state(\n", + " frames=frames,\n", + " values=speed,\n", + " ref_start=240,\n", + " ref_end=640,\n", + ")\n", + "\n", + "print(\n", + " f\"Density - theta: {result_density.theta}, \"\n", + " f\"acf: {result_density.acf:.3f}, \"\n", + " f\"mean: {result_density.mean:.3f}, std: {result_density.std:.3f}\"\n", + ")\n", + "for s, e in zip(result_density.frame_start, result_density.frame_end):\n", + " print(f\" Steady interval: frame {s:.0f} to {e:.0f}\")\n", + "\n", + "print(\n", + " f\"\\nSpeed - theta: {result_speed.theta}, \"\n", + " f\"acf: {result_speed.acf:.3f}, \"\n", + " f\"mean: {result_speed.mean:.3f}, std: {result_speed.std:.3f}\"\n", + ")\n", + "for s, e in zip(result_speed.frame_start, result_speed.frame_end):\n", + " print(f\" Steady interval: frame {s:.0f} to {e:.0f}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Combine steady states\n", + "\n", + "The final steady state of the flow is the **overlap** of the density and speed steady states." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "combined = combine_steady_states([result_density, result_speed])\n", + "\n", + "for start, end in combined:\n", + " ratio = (end - start) / len(frames) * 100\n", + " print(f\"Combined steady state: frame {start:.0f} to {end:.0f} (ratio: {ratio:.1f}% of total observations)\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Visualise the CUSUM statistics and detected steady states" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.patches as mpatches\n", + "import matplotlib.pyplot as plt\n", + "\n", + "fps = 16 # frames per second\n", + "ref_start, ref_end = 240, 640\n", + "\n", + "fig, axes = plt.subplots(2, 2, figsize=(14, 10))\n", + "\n", + "for col, (result, label, unit) in enumerate(\n", + " [\n", + " (result_density, \"Density\", r\"$\\rho\\;/\\mathrm{m^{-2}}$\"),\n", + " (result_speed, \"Speed\", r\"$v\\;/\\mathrm{m/s}$\"),\n", + " ]\n", + "):\n", + " values_col = density if col == 0 else speed\n", + "\n", + " # Top row: CUSUM statistics\n", + " ax = axes[0, col]\n", + " ax.plot(frames / fps, result.cusum, \"b--\", lw=1.5, label=r\"$S_k$\")\n", + " ax.axhline(y=result.theta, color=\"r\", lw=2, label=rf\"$\\theta = {result.theta:.0f}$\")\n", + " ax.set_xlabel(\"t / s\")\n", + " ax.set_ylabel(f\"Statistics ({label})\")\n", + " ax.set_ylim(-5, 110)\n", + " ax.legend(fontsize=12)\n", + " ax.set_title(f\"{label}: CUSUM Statistics\")\n", + "\n", + " # Bottom row: time series with steady state highlighted\n", + " ax = axes[1, col]\n", + " ax.plot(frames / fps, values_col, \"b-\", lw=1, alpha=0.8)\n", + " ax.axvline(x=ref_start / fps, color=\"g\", ls=\"--\", lw=1.5, label=\"Reference\")\n", + " ax.axvline(x=ref_end / fps, color=\"g\", ls=\"--\", lw=1.5)\n", + "\n", + " # shade individual steady intervals\n", + " for i, (s, e) in enumerate(zip(result.frame_start, result.frame_end)):\n", + " lbl = f\"Steady ({label})\" if i == 0 else None\n", + " ax.axvspan(s / fps, e / fps, alpha=0.15, color=\"red\", label=lbl)\n", + "\n", + " # shade combined steady interval\n", + " for i, (s, e) in enumerate(combined):\n", + " lbl = \"Combined steady\" if i == 0 else None\n", + " ax.axvspan(s / fps, e / fps, alpha=0.2, color=\"gold\", label=lbl)\n", + "\n", + " ax.set_xlabel(\"t / s\")\n", + " ax.set_ylabel(unit)\n", + " ax.legend(fontsize=10, loc=\"upper right\")\n", + " ax.set_title(f\"{label}: Time Series with Steady State\")\n", + "\n", + "fig.tight_layout()\n", + "plt.show()" + ] + }, { "cell_type": "markdown", "metadata": { diff --git a/pedpy/__init__.py b/pedpy/__init__.py index 50a37d79..b9663ad3 100644 --- a/pedpy/__init__.py +++ b/pedpy/__init__.py @@ -126,6 +126,11 @@ compute_species, compute_voronoi_speed, ) +from .methods.steady_state_detector import ( + SteadyStateResult, + combine_steady_states, + detect_steady_state, +) from .plotting.plotting import ( PEDPY_BLUE, PEDPY_GREEN, @@ -196,6 +201,9 @@ "is_species_valid", "is_trajectory_valid", "compute_pair_distribution_function", + "SteadyStateResult", + "combine_steady_states", + "detect_steady_state", "DensityMethod", "RsetMethod", "SpeedMethod", diff --git a/pedpy/methods/steady_state_detector.py b/pedpy/methods/steady_state_detector.py new file mode 100644 index 00000000..4524c61a --- /dev/null +++ b/pedpy/methods/steady_state_detector.py @@ -0,0 +1,419 @@ +"""Steady state detection using the modified CUSUM algorithm. + +Implements the method from: + Liao et al., "Measuring the steady state of pedestrian flow in + bottleneck experiments", Physica A 461 (2016) 248-261. + +The algorithm detects intervals in density/speed time series where +the process is stationary ("steady state"), using a modified +Cumulative Sum Control Chart (CUSUM) with a detection threshold +calibrated via an autoregressive model. +""" + +import math +from dataclasses import dataclass +from typing import List, Tuple + +import numpy as np +from scipy.stats import norm, pearsonr + + +@dataclass(frozen=True) +class SteadyStateResult: + """Result of steady state detection on a single time series. + + Attributes: + frame_start: start frames of detected steady intervals + frame_end: end frames of detected steady intervals + cusum: the CUSUM statistics array (one value per observation) + theta: the calibrated detection threshold + mean: mean of the reference process + std: standard deviation of the reference process + acf: lag-1 autocorrelation of the reference process + """ + + frame_start: np.ndarray + frame_end: np.ndarray + cusum: np.ndarray + theta: float + mean: float + std: float + acf: float + + +def detect_steady_state( + *, + frames: np.ndarray, + values: np.ndarray, + ref_start: int, + ref_end: int, + alpha: float = 0.99, + gamma: float = 0.99, + s_max: int = 100, + grid_size: int = 500, + grid_half_width: float = 3.2, +) -> SteadyStateResult: + """Detect steady state intervals in a time series. + + Uses the modified CUSUM algorithm from Liao et al. (2016) to find + intervals where the time series is stationary relative to a + user-selected reference process. + + Args: + frames: array of frame numbers (monotonically increasing) + values: array of observed values (density or speed), same length + as frames + ref_start: start frame of the reference process (inclusive) + ref_end: end frame of the reference process (inclusive) + alpha: critical probability for the step function F + (default 0.99) + gamma: confidence level for the detection threshold theta + (default 0.99) + s_max: upper boundary of the CUSUM statistics (default 100) + grid_size: number of grid points K for the numerical theta + computation (default 500) + grid_half_width: half-width of the discretisation grid + (default 3.2) + + Returns: + SteadyStateResult with detected steady intervals and diagnostics + + Raises: + ValueError: if reference range is invalid or data is too short. + """ + frames = np.asarray(frames, dtype=float) + values = np.asarray(values, dtype=float) + + if len(frames) != len(values): + raise ValueError("frames and values must have the same length.") + + if ref_start > ref_end: + raise ValueError(f"ref_start ({ref_start}) must be <= ref_end ({ref_end}).") + + # --- reference statistics --- + ref_mask = (frames >= ref_start) & (frames <= ref_end) + ref_series = values[ref_mask] + if len(ref_series) < 3: + raise ValueError("Reference range must contain at least 3 observations.") + + ref_mean = float(np.mean(ref_series)) + ref_std = float(np.std(ref_series)) + ref_acf = float(pearsonr(ref_series[:-1], ref_series[1:])[0]) + + if ref_std == 0: + raise ValueError("Reference series has zero standard deviation.") + + # --- calibrate detection threshold theta --- + q = norm.ppf(alpha) + theta = _compute_theta( + acf=ref_acf, + q=q, + gamma=gamma, + s_max=s_max, + grid_size=grid_size, + grid_half_width=grid_half_width, + ) + + # --- compute CUSUM statistics --- + cusum = _compute_cusum(values, ref_mean, ref_std, q, s_max) + + # --- extract steady intervals --- + starts, ends = _find_steady_intervals( + frames=frames, + cusum=cusum, + theta=theta, + s_max=s_max, + ) + + return SteadyStateResult( + frame_start=starts, + frame_end=ends, + cusum=cusum, + theta=theta, + mean=ref_mean, + std=ref_std, + acf=ref_acf, + ) + + +def combine_steady_states( + results: List[SteadyStateResult], +) -> List[Tuple[float, float]]: + """Find the overlapping steady intervals across multiple series. + + For pedestrian flow analysis, the final steady state is the + intersection of the steady states detected independently for + density and speed. + + Args: + results: list of SteadyStateResult objects (e.g. one for + density, one for speed) + + Returns: + list of (start_frame, end_frame) tuples for the combined + steady intervals + """ + if not results: + return [] + + # collect all intervals from all results + all_intervals = [] + for r in results: + for s, e in zip(r.frame_start, r.frame_end, strict=True): + all_intervals.append((s, e)) + + if not all_intervals: + return [] + + if len(results) == 1: + return all_intervals + + # find pairwise overlaps across all results + # start with intervals from the first result, intersect with each + # subsequent result + current = list(zip(results[0].frame_start, results[0].frame_end, strict=True)) + for r in results[1:]: + next_intervals = list(zip(r.frame_start, r.frame_end, strict=True)) + current = _intersect_interval_lists(current, next_intervals) + + return current + + +# --------------- internal helpers --------------- + + +def _compute_cusum( + values: np.ndarray, + ref_mean: float, + ref_std: float, + q: float, + s_max: int, +) -> np.ndarray: + """Compute the modified CUSUM statistics (Eq. 7 in Liao 2016). + + s_i = min(max(0, s_{i-1} + F(x_tilde_i)), s_max) + with s_0 = s_max and F(x) = 1 if |x| > q else -1. + """ + x_tilde = (values - ref_mean) / ref_std + f_vals = np.where(np.abs(x_tilde) > q, 1.0, -1.0) + s = np.empty(len(values), dtype=float) + s[0] = s_max + for i in range(1, len(s)): + s[i] = min(max(0.0, s[i - 1] + f_vals[i]), s_max) + return s + + +def _find_steady_intervals( + *, + frames: np.ndarray, + cusum: np.ndarray, + theta: float, + s_max: int, +) -> Tuple[np.ndarray, np.ndarray]: + """Extract steady state frame ranges from CUSUM statistics. + + A point is in steady state when cusum < theta. The raw intervals + are then adjusted for the reaction times: + t_reaching = (s_max - theta) / fps (entering steady) + t_leaving = theta / fps (leaving steady) + Since we work in frame units (fps=1 frame), the corrections are + applied directly in frame counts. + """ + frame_step = 1.0 + if len(frames) > 1: + frame_step = float(np.median(np.diff(frames))) + + steady_mask = cusum < theta + if not np.any(steady_mask): + return np.array([]), np.array([]) + + # The legacy implementation updates the statistic after processing the + # current observation and records it on the next frame boundary. + steady_frames = frames[steady_mask] + frame_step + + # split into contiguous runs + starts = [] + ends = [] + run_start = steady_frames[0] + for i in range(1, len(steady_frames)): + # detect gap (non-consecutive frames) + if not np.isclose(steady_frames[i] - steady_frames[i - 1], frame_step): + run_end = steady_frames[i - 1] + starts.append(run_start) + ends.append(run_end) + run_start = steady_frames[i] + # last run + starts.append(run_start) + ends.append(steady_frames[-1]) + + # apply reaction time corrections (Eq. 8 in Liao 2016) + # t_reaching = (s_max - theta): entering steady state + # t_leaving = theta: leaving steady state + reaching = s_max - theta + leaving = theta + + adjusted_starts = [] + adjusted_ends = [] + for s, e in zip(starts, ends, strict=True): + adj_s = s - reaching + adj_e = e - leaving + if adj_s < adj_e: + adjusted_starts.append(adj_s) + adjusted_ends.append(adj_e) + + return np.array(adjusted_starts), np.array(adjusted_ends) + + +def _compute_theta( + *, + acf: float, + q: float, + gamma: float, + s_max: int, + grid_size: int, + grid_half_width: float, +) -> float: + """Calibrate the detection threshold theta via the autoregressive model. + + Uses the block-tridiagonal Thomas Algorithm to solve for the + stationary distribution of the Markov chain (y_i, s_i) as + described in Appendix A of Liao et al. (2016). + + The autoregressive process is: + y_i = c * y_{i-1} + sqrt(1 - c^2) * epsilon_i + + theta is the upper gamma-percentile of the marginal distribution + of s from the stationary distribution. + """ + c = acf + n_blocks = s_max # number of s-blocks (0..s_max) + + d = 2 * grid_half_width / grid_size + xi = np.arange(-grid_half_width, grid_half_width + 0.0001, d) + kp1 = grid_size + 1 # grid points per block + + # find grid indices ia, ib corresponding to -q and +q + ia, ib, dnorm = _init_grid_params(xi, q) + + # build block matrices b1 (outside [-q,q]) and b2 (inside [-q,q]) + b1, b2 = _build_block_matrices(xi, ia, ib, kp1, c, d) + eye = np.eye(kp1) + + # solve block-tridiagonal system via Thomas algorithm + td = _solve_thomas(b1, b2, eye, kp1, n_blocks, dnorm) + + # assemble full solution vector and normalise + tms = np.vstack(td) + tms = tms / (d * np.sum(tms)) + + # marginal distribution of s: integrate over x for each s + tm = np.zeros(n_blocks + 1) + for j in range(n_blocks + 1): + begin = j * kp1 + end = begin + kp1 + tm[j] = d * np.sum(tms[begin:end]) + + # find theta as upper gamma-percentile + cumulative = tm[0] + theta = 1 + while theta + 1 < len(tm) and cumulative + tm[theta] < gamma: + cumulative += tm[theta] + theta += 1 + + return float(theta) + + +def _init_grid_params(xi: np.ndarray, q: float) -> Tuple[int, int, float]: + """Find grid indices for -q and +q, and compute normalisation constant.""" + ia = ib = 0 + dnorm = 0.0 + for i in range(len(xi)): + if ia == 0 and i + 1 < len(xi) and xi[i + 1] > -q: + ia = i + if ib == 0 and xi[i] > q: + ib = i + dnorm += (1.0 / math.sqrt(2 * math.pi)) * math.exp(-xi[i] * xi[i] / 2.0) + return ia, ib, dnorm + + +def _build_block_matrices( + xi: np.ndarray, + ia: int, + ib: int, + kp1: int, + c: float, + d: float, +) -> Tuple[np.ndarray, np.ndarray]: + """Build block transition matrices b1 and b2.""" + + def func_b(i: int, k: int) -> float: + if abs(c) >= 1.0: + return d / math.sqrt(2 * math.pi) * math.exp(-((xi[i] - xi[k]) ** 2) / 2.0) + return ( + d * math.exp(-((xi[i] - xi[k] * c) ** 2) / (2.0 * (1.0 - c * c))) / math.sqrt(2 * math.pi * (1.0 - c * c)) + ) + + b1 = np.zeros((kp1, kp1)) + b2 = np.zeros((kp1, kp1)) + for i in range(kp1): + for j in range(kp1): + if i < ia or i >= ib: + b1[i, j] = func_b(i, j) + if ia <= i < ib: + b2[i, j] = func_b(i, j) + return b1, b2 + + +def _solve_thomas( + b1: np.ndarray, + b2: np.ndarray, + eye: np.ndarray, + kp1: int, + n_blocks: int, + dnorm: float, +) -> List[np.ndarray]: + """Solve the block-tridiagonal system using the Thomas algorithm.""" + td: List[np.ndarray] = [np.zeros((kp1, 1)) for _ in range(n_blocks + 1)] + td[n_blocks][-1, 0] = dnorm + + tc_arr: List[np.ndarray] = [np.zeros((kp1, kp1)) for _ in range(n_blocks + 1)] + + # block 0: Tb = b2 - I, Tc = b2 + tb0 = b2 - eye + tc_arr[0] = np.linalg.solve(tb0, b2.copy()) + td[0] = np.linalg.solve(tb0, td[0]) + + # blocks 1..n_blocks + for i in range(1, n_blocks + 1): + if i < n_blocks: + tb_i = -eye + tc_i = b2.copy() + else: + tb_i = b1 - eye + tc_i = np.zeros((kp1, kp1)) + + aa = tb_i - b1 @ tc_arr[i - 1] + if i < n_blocks: + tc_arr[i] = np.linalg.solve(aa, tc_i) + td[i] = np.linalg.solve(aa, td[i] - b1 @ td[i - 1]) + + # backward sweep + for i in range(n_blocks - 1, -1, -1): + td[i] = td[i] - tc_arr[i] @ td[i + 1] + + return td + + +def _intersect_interval_lists( + a: List[Tuple[float, float]], + b: List[Tuple[float, float]], +) -> List[Tuple[float, float]]: + """Compute pairwise intersections of two lists of intervals.""" + result = [] + for a_s, a_e in a: + for b_s, b_e in b: + start = max(a_s, b_s) + end = min(a_e, b_e) + if start <= end: + result.append((start, end)) + return result diff --git a/tests/unit_tests/methods/test_steady_state_detector.py b/tests/unit_tests/methods/test_steady_state_detector.py new file mode 100644 index 00000000..d2c940ce --- /dev/null +++ b/tests/unit_tests/methods/test_steady_state_detector.py @@ -0,0 +1,346 @@ +"""Tests for steady state detection (modified CUSUM algorithm).""" + +import numpy as np +import pytest + +from pedpy.methods.steady_state_detector import ( + SteadyStateResult, + _compute_cusum, + _compute_theta, + _find_steady_intervals, + _intersect_interval_lists, + combine_steady_states, + detect_steady_state, +) + + +class TestComputeCusum: + """Tests for the CUSUM statistics computation.""" + + def test_constant_series_at_reference_mean_decreases(self): + """CUSUM should decrease when observations match reference.""" + values = np.full(200, 5.0) + q = 2.326 # norm.ppf(0.99) + cusum = _compute_cusum(values, ref_mean=5.0, ref_std=1.0, q=q, s_max=100) + assert cusum[0] == 100 # starts at s_max + # should decrease since all x_tilde = 0, F = -1 + assert cusum[-1] == 0 + + def test_outlier_series_stays_high(self): + """CUSUM should stay high when observations are far from reference.""" + values = np.full(50, 100.0) + q = 2.326 + cusum = _compute_cusum(values, ref_mean=5.0, ref_std=1.0, q=q, s_max=100) + # all values are outliers, F = +1 always, stays at s_max + assert cusum[-1] == 100 + + def test_cusum_bounded(self): + """CUSUM should always be in [0, s_max].""" + rng = np.random.default_rng(42) + values = rng.normal(0, 3, size=500) + q = 2.326 + cusum = _compute_cusum(values, ref_mean=0.0, ref_std=1.0, q=q, s_max=100) + assert np.all(cusum >= 0) + assert np.all(cusum <= 100) + + +class TestComputeTheta: + """Tests for theta calibration via autoregressive model.""" + + def test_theta_matches_original_implementation(self): + """Theta values should match the original SteadyState.py code.""" + # These values were verified against the original implementation + expected = {0.95: 84, 0.97: 78, 0.99: 53} + for acf, expected_theta in expected.items(): + theta = _compute_theta( + acf=acf, + q=2.326, + gamma=0.99, + s_max=100, + grid_size=500, + grid_half_width=3.2, + ) + assert theta == expected_theta, f"acf={acf}: got theta={theta}, expected {expected_theta}" + + def test_theta_in_valid_range(self): + """Theta should be between 1 and s_max.""" + theta = _compute_theta( + acf=0.9, + q=2.326, + gamma=0.99, + s_max=100, + grid_size=500, + grid_half_width=3.2, + ) + assert 1 <= theta <= 100 + + +class TestFindSteadyIntervals: + """Tests for extracting steady intervals from CUSUM.""" + + def test_no_steady_state(self): + """When CUSUM never drops below theta, no intervals found.""" + frames = np.arange(100, dtype=float) + cusum = np.full(100, 90.0) + starts, ends = _find_steady_intervals( + frames=frames, + cusum=cusum, + theta=50.0, + s_max=100, + ) + assert len(starts) == 0 + assert len(ends) == 0 + + def test_single_steady_interval(self): + """Detect a single steady interval.""" + frames = np.arange(200, dtype=float) + cusum = np.full(200, 80.0) + # create a dip below theta in the middle + cusum[50:150] = 10.0 + starts, ends = _find_steady_intervals( + frames=frames, + cusum=cusum, + theta=50.0, + s_max=100, + ) + assert len(starts) == 1 + # the interval should be adjusted by reaction times + + +class TestDetectSteadyState: + """Integration tests for the full detection pipeline.""" + + def test_stationary_signal(self): + """A stationary signal should have a large steady interval.""" + rng = np.random.default_rng(123) + n = 1000 + frames = np.arange(n, dtype=float) + values = rng.normal(5.0, 0.5, size=n) + + result = detect_steady_state( + frames=frames, + values=values, + ref_start=200, + ref_end=400, + ) + # should detect at least one steady interval + assert len(result.frame_start) >= 1 + assert result.theta > 0 + assert result.mean != 0 + + def test_step_function_detects_transition(self): + """A step function should yield a steady interval only in flat parts.""" + n = 600 + frames = np.arange(n, dtype=float) + values = np.concatenate( + [ + np.full(200, 3.0) + np.random.default_rng(1).normal(0, 0.1, 200), + np.full(200, 6.0) + np.random.default_rng(2).normal(0, 0.1, 200), + np.full(200, 3.0) + np.random.default_rng(3).normal(0, 0.1, 200), + ] + ) + # reference from first plateau + result = detect_steady_state( + frames=frames, + values=values, + ref_start=50, + ref_end=150, + ) + # steady intervals should NOT include the middle plateau (around 200-400) + for s, e in zip(result.frame_start, result.frame_end): + # no steady interval should be fully contained in the step region + assert not (s >= 210 and e <= 390), f"Unexpected steady interval {s}-{e} in step region" + + def test_invalid_reference_range(self): + """Should raise ValueError for invalid reference range.""" + frames = np.arange(100, dtype=float) + values = np.ones(100) + with pytest.raises(ValueError): + detect_steady_state( + frames=frames, + values=values, + ref_start=50, + ref_end=10, + ) + + def test_mismatched_lengths(self): + """Should raise ValueError for mismatched input lengths.""" + with pytest.raises(ValueError): + detect_steady_state( + frames=np.arange(10), + values=np.ones(5), + ref_start=0, + ref_end=5, + ) + + def test_zero_std_reference(self): + """Should raise ValueError when reference has zero std.""" + frames = np.arange(100, dtype=float) + values = np.ones(100) # all identical + with pytest.raises(ValueError, match="zero standard deviation"): + detect_steady_state( + frames=frames, + values=values, + ref_start=10, + ref_end=50, + ) + + +class TestCombineSteadyStates: + """Tests for combining steady states across multiple series.""" + + def test_overlapping_intervals(self): + """Overlapping intervals from two series should be intersected.""" + r1 = SteadyStateResult( + frame_start=np.array([100.0]), + frame_end=np.array([500.0]), + cusum=np.array([]), + theta=50.0, + mean=3.0, + std=0.5, + acf=0.9, + ) + r2 = SteadyStateResult( + frame_start=np.array([200.0]), + frame_end=np.array([600.0]), + cusum=np.array([]), + theta=70.0, + mean=0.7, + std=0.1, + acf=0.95, + ) + combined = combine_steady_states([r1, r2]) + assert len(combined) == 1 + assert combined[0][0] == 200.0 + assert combined[0][1] == 500.0 + + def test_no_overlap(self): + """Non-overlapping intervals should yield empty result.""" + r1 = SteadyStateResult( + frame_start=np.array([100.0]), + frame_end=np.array([200.0]), + cusum=np.array([]), + theta=50.0, + mean=3.0, + std=0.5, + acf=0.9, + ) + r2 = SteadyStateResult( + frame_start=np.array([300.0]), + frame_end=np.array([400.0]), + cusum=np.array([]), + theta=70.0, + mean=0.7, + std=0.1, + acf=0.95, + ) + combined = combine_steady_states([r1, r2]) + assert len(combined) == 0 + + def test_single_result(self): + """Single result should return its intervals directly.""" + r1 = SteadyStateResult( + frame_start=np.array([100.0, 300.0]), + frame_end=np.array([200.0, 400.0]), + cusum=np.array([]), + theta=50.0, + mean=3.0, + std=0.5, + acf=0.9, + ) + combined = combine_steady_states([r1]) + assert len(combined) == 2 + + def test_empty_results(self): + """Empty result list should return empty.""" + assert combine_steady_states([]) == [] + + +class TestIntersectIntervalLists: + """Tests for interval list intersection.""" + + def test_basic_overlap(self): + """Two overlapping intervals should produce their intersection.""" + result = _intersect_interval_lists( + [(1, 10)], + [(5, 15)], + ) + assert result == [(5, 10)] + + def test_multiple_overlaps(self): + """Multiple intervals should produce all pairwise intersections.""" + result = _intersect_interval_lists( + [(0, 10), (20, 30)], + [(5, 25)], + ) + assert len(result) == 2 + assert (5, 10) in result + assert (20, 25) in result + + +class TestReferenceData: + """Validation against the original SteadyState.py reference results.""" + + @pytest.fixture + def reference_data(self): + """Load the reference dataset.""" + import os + + data_path = os.path.join( + os.path.dirname(__file__), + "..", + "..", + "..", + "materials", + "SteadyState", + "rho_v_Voronoi_traj_AO_b240.txt_id_1.dat", + ) + if not os.path.exists(data_path): + pytest.skip("Reference data not available") + data = np.loadtxt(data_path, usecols=[0, 1, 2]) + data = data[data[:, 1] != 0] + return data + + def test_density_steady_state(self, reference_data): + """Density steady state should match reference within tolerance.""" + data = reference_data + result = detect_steady_state( + frames=data[:, 0], + values=data[:, 1], + ref_start=240, + ref_end=640, + ) + # Reference from the legacy SteadyState.py implementation. + np.testing.assert_array_equal(result.frame_start, np.array([161.0])) + np.testing.assert_array_equal(result.frame_end, np.array([866.0])) + + def test_speed_steady_state(self, reference_data): + """Speed steady state should match reference within tolerance.""" + data = reference_data + result = detect_steady_state( + frames=data[:, 0], + values=data[:, 2], + ref_start=240, + ref_end=640, + ) + # Reference from the legacy SteadyState.py implementation. + np.testing.assert_array_equal(result.frame_start, np.array([231.0])) + np.testing.assert_array_equal(result.frame_end, np.array([965.0])) + + def test_combined_steady_state(self, reference_data): + """Combined steady state should match reference within tolerance.""" + data = reference_data + result_rho = detect_steady_state( + frames=data[:, 0], + values=data[:, 1], + ref_start=240, + ref_end=640, + ) + result_v = detect_steady_state( + frames=data[:, 0], + values=data[:, 2], + ref_start=240, + ref_end=640, + ) + combined = combine_steady_states([result_rho, result_v]) + assert combined == [(231.0, 866.0)]