Reference
Analysis¶
The quantifiers that consume any System. Prose-first
treatments live in the Analysis section.
Orbit diagrams¶
orbit_diagram
¶
orbit_diagram(
system: Any,
param: str,
values: Any,
*,
points_per_value: int = 200,
transient: int = 500,
carry_state: bool = True,
components: int | str | tuple[Any, ...] = 0,
section: Any | None = None,
ic: Any | None = None,
seed: int | None = 0,
dt: float = _FLOW_DT,
max_time: float = _FLOW_MAX_TIME,
) -> OrbitDiagram
Sweep a parameter and record the asymptotic orbit at each value.
Pass any system — this is the one-liner::
ts.analysis.orbit_diagram(ts.systems.Lorenz(), "rho", np.linspace(0.0, 50.0, 200))
A :class:~tsdynamics.families.DiscreteMap is swept directly. A flow is
reduced to a discrete view automatically, because a bifurcation diagram of a
flow is the common case and a section always has to be chosen:
- by default, successive maxima of the recorded component — the next-amplitude (peak) map, the view Lorenz (1963) used and the one every textbook bifurcation diagram of a flow is drawn from. It needs no parameter-dependent constant, so it keeps producing points as the attractor moves across the sweep, and a value whose flow has settled on an equilibrium records that single state (the fixed-point branch);
- or the section you name with
section=(any :class:~tsdynamics.derived.PoincareMapplanespelling, e.g.section=("z", 27.0, "up")).
The choice is never silent: it is recorded in meta["section"] and
printed under the figure's title. A flow already wrapped in a
:class:~tsdynamics.derived.PoincareMap /
:class:~tsdynamics.derived.StroboscopicMap is used as given.
orbit_diagram is the one spelling: bifurcation_diagram was the
same object under a second name and is gone (a shared implementation can name
only one of its spellings in a message, so half of all callers were answered
about a function they had never typed). The picture is still called a
bifurcation diagram — PlotKind.BIFURCATION_DIAGRAM is untouched.
ODE parameter changes reuse the cached lowered tape (control parameters), so flow sweeps stay cheap; DDE sweeps re-lower per value (their structure depends on all parameters).
Notes
How each value is run. A genuine :class:~tsdynamics.families.DiscreteMap
sweeps the whole parameter array in one Rust kernel call. A
:class:~tsdynamics.derived.PoincareMap collects each value's section in one
engine call through :meth:~tsdynamics.derived.PoincareMap.trajectory (the
wired Rust event march) — ~70x faster than the per-dt step() loop it
replaces. Everything else (StroboscopicMap, a map whose _step will not
lower, a wheel-free environment) drives the per-step protocol loop.
The Poincaré fast path marches the fixed-step rk4 kernel at the map's
detection dt (see :mod:tsdynamics.derived._crossings), which is the same
discretisation :func:~tsdynamics.analysis.poincare_section and
:func:~tsdynamics.analysis.return_map already use — so an orbit diagram and
a section of the same flow are now consistent with each other. It is a
different discretisation from the flow's adaptive default that the old
step() loop used. Measured over a 40-value Rössler c-cascade at the
default dt=0.01: in a periodic window the branch values shift by ~5e-8
(up to ~3e-3 for a value sitting right at a window edge, where a tiny shift
in the orbit is amplified), and on a chaotic band the two are the same
attractor rather than the same points (support bounds agree to ~6e-2). The
detected branch structure — OrbitDiagram.periods() and
bifurcation_points() — is unchanged, at both a short and a long
transient.
rk4 is fixed-step, so its stability region is bounded (|λh| ≲ 2.785).
A sweep that runs into a regime where the flow is stiff relative to the
detection dt therefore records an empty point set with a
RuntimeWarning for those values, where the old adaptive step() loop
would have shrunk its step and carried on (on the Rössler c-sweep at
dt=0.01 this starts at c ≈ 200). Pass a smaller dt to
:class:~tsdynamics.derived.PoincareMap when sweeping into a fast-timescale
regime.
The flow path. With no section=, each parameter value is integrated
(output step dt, horizon grown from 100 time units up to max_time)
until transient + n maxima of the recorded component have been collected;
the last n are the column. transient and n therefore count
peaks, exactly as they count iterates for a map and crossings for a
section. Each peak is sharpened by the same parabolic interpolation
:func:~tsdynamics.analysis.return_map uses, so a coarse dt still gives a
crisp period-1 branch.
| PARAMETER | DESCRIPTION |
|---|---|
system
|
The system to sweep — a map, a flow, or a discrete view of a flow. Never
mutated: each value gets a fresh
TYPE:
|
param
|
Parameter name to sweep.
TYPE:
|
values
|
Parameter values, in sweep order.
TYPE:
|
points_per_value
|
Points recorded per parameter value (map iterates / section crossings / peaks, according to the view) — the height of one vertical stripe of the diagram. .. versionchanged:: 6.0
Named
TYPE:
|
transient
|
Points discarded before recording, at every value — a count, in the
unit the sweep advances in: map iterations for a map, crossings for a
Poincaré/stroboscopic view, recorded peaks for the flow peak map. (The
subject of an orbit diagram is always a discrete view, so this word never
means time here, unlike
TYPE:
|
carry_state
|
Start each value from the previous value's final state (follows the
attractor branch; the classic way to draw clean diagrams). When
False, every value starts from
TYPE:
|
components
|
Which state component(s) to record (names allowed when the system
declares
TYPE:
|
section
|
Slice a flow with this Poincaré section instead of the default peak map.
Any :class:
TYPE:
|
ic
|
Initial state for the first value (and every value when
TYPE:
|
seed
|
Seed for the random initial condition when the system has none, so the
diagram is reproducible. Pass .. versionchanged:: 6.0
Was
TYPE:
|
dt
|
Output sampling step of the flow path's integration, in time units. Ignored for a map / an already-discrete view, which iterate.
TYPE:
|
max_time
|
Ceiling on the flow path's integration horizon per parameter value.
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
OrbitDiagram
|
The swept |
| RAISES | DESCRIPTION |
|---|---|
InvalidInputError
|
If |
InvalidParameterError
|
If |
ValueError
|
If a named |
| WARNS | DESCRIPTION |
|---|---|
RuntimeWarning
|
When a parameter value diverges (that value records an empty set and the
sweep continues), or when the flow path exhausts |
References
May, R. M. (1976). Simple mathematical models with very complicated dynamics. Nature, 261, 459--467.
Lorenz, E. N. (1963). Deterministic nonperiodic flow. Journal of the Atmospheric Sciences, 20, 130--141. (Successive maxima as the discrete view of a flow.)
Examples:
>>> # a flow, straight up — the section is chosen for you and reported:
>>> od = orbit_diagram(Lorenz(), "rho", np.linspace(0.0, 50.0, 200))
>>> od.meta["section"]
'successive maxima of x'
>>> # ... or name the section yourself:
>>> od = orbit_diagram(Rossler(), "c", np.linspace(2, 6, 80),
... section=("y", 0.0, "up"))
>>> # a map:
>>> od = orbit_diagram(Logistic(), "r", np.linspace(2.5, 4.0, 600), points_per_value=120)
>>> x, y = od.flat()
Source code in src/tsdynamics/analysis/orbits/orbit_diagram.py
884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 | |
OrbitDiagram
dataclass
¶
OrbitDiagram(
param: str = "",
values: ndarray = (lambda: empty(0))(),
points: list[ndarray] = list(),
components: tuple[int, ...] = (),
equilibria: ndarray = (lambda: empty(0, dtype=bool))(),
*,
meta: Mapping[str, Any] = dict(),
)
Bases: AnalysisResult
Result of :func:orbit_diagram.
An :class:~tsdynamics.analysis._result.AnalysisResult, so it carries
.meta / the readout repr / .to_dict() / the .plot seam. Iterate to
get (value, points) pairs, or use :meth:flat for the scatter-ready
arrays.
| ATTRIBUTE | DESCRIPTION |
|---|---|
param |
Name of the swept control parameter.
TYPE:
|
values |
The swept parameter axis — see the note below.
TYPE:
|
points |
The asymptotic orbit recorded at each value,
TYPE:
|
.. note |
TYPE:
|
flat
¶
Flatten to scatter-plot arrays (x, y).
x repeats each parameter value once per recorded point; y is
the chosen recorded component.
Source code in src/tsdynamics/analysis/orbits/orbit_diagram.py
periods
¶
Return the detected period at each swept parameter value.
Counts the distinct asymptotic branches in the recorded orbit — the
period of a periodic window — by clustering the points of one component
with a scale-free gap test: a new branch starts where the sorted-value
gap exceeds rtol times the orbit's range. A finite period p >= 2
is only reported when the recorded iterate sequence actually revisits
its values cyclically (v[i] ≈ v[i + p] to rtol); a chaotic band
whose finite-sample points merely cluster into p bins fails this
repeat test and is reported as aperiodic (0).
| PARAMETER | DESCRIPTION |
|---|---|
component
|
Which recorded component to count branches in.
TYPE:
|
max_period
|
Periods above this are reported as
TYPE:
|
rtol
|
Relative gap (fraction of the per-value range) separating branches.
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
numpy.ndarray of int
|
One entry per parameter value: the period |
Source code in src/tsdynamics/analysis/orbits/orbit_diagram.py
bifurcation_points
¶
bifurcation_points(
*,
component: int = 0,
max_period: int = 16,
rtol: float = 0.01,
labelled: bool = False,
) -> ndarray
Parameter values where the detected period changes.
Locates the boundaries of the period-doubling cascade (and other
bifurcations) as the midpoints between consecutive swept values across
which :meth:periods differs. Transitions touching a diverged value
(-1) are skipped.
| PARAMETER | DESCRIPTION |
|---|---|
component
|
Which recorded component to count branches in.
TYPE:
|
max_period
|
Periods above this are treated as aperiodic when detecting changes.
TYPE:
|
rtol
|
Relative gap separating branches in :meth:
TYPE:
|
labelled
|
Return a structured array carrying the period on each side of every transition, so a genuine period-doubling can be told from a band split inside the chaotic regime. Measured on the logistic map, the plain array holds 20 values of
which the first ( .. versionadded:: 6.0
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
ndarray
|
|
Examples:
>>> import numpy as np, tsdynamics as ts
>>> od = ts.analysis.orbit_diagram(
... ts.systems.Logistic(), "r", np.linspace(2.8, 3.6, 120)
... )
>>> rows = od.bifurcation_points(labelled=True)
>>> float(rows["value"][rows["kind"] == "period-doubling"][0])
3.0033613445378155
References
Feigenbaum, M. J. (1978). Quantitative universality for a class of nonlinear transformations. Journal of Statistical Physics, 19, 25--52.
Source code in src/tsdynamics/analysis/orbits/orbit_diagram.py
Poincaré sections¶
poincare_section
¶
poincare_section(
system: Any,
plane: tuple[Any, ...] | None = None,
*,
direction: int | str = +1,
crossings: int = 1000,
skip_crossings: int = 0,
ic: Any | None = None,
dt: float = 0.01,
max_time: float = 10000.0,
seed: int | None = 0,
) -> PoincareSection
Poincaré surface of section.
Two input modes:
- System → wraps it in a :class:
~tsdynamics.derived.PoincareMapand collectscrossingsroot-refined crossings on the fast Rust event engine (stream WS-CROSSKERNEL). - Trajectory → finds the plane crossings between consecutive samples by linear interpolation (pure data path; accuracy limited by the trajectory's sampling interval).
| PARAMETER | DESCRIPTION |
|---|---|
system
|
A flow to section, or measured trajectory data (the
TYPE:
|
plane
|
The section, in any of three spellings:
For example Omit it and one is chosen, from a short probe orbit: the
widest-spread component (largest interquartile range) among those the
orbit repeatedly returns to, sectioned at its median — an offset the
probe orbit provably straddles (so the auto section cannot miss it; on
a system with no
TYPE:
|
direction
|
Crossing direction filter (
TYPE:
|
crossings
|
Number of crossings to collect (system mode) — the number of points the returned section holds. .. versionchanged:: 6.0
Named
TYPE:
|
skip_crossings
|
Number of leading crossings to discard before recording. (A section
transient is a count of crossings, deliberately distinct from the
time/step
TYPE:
|
ic
|
Where to start the section run — .. versionadded:: 6.0
Every sibling that runs a system took
TYPE:
|
dt
|
Crossing-detection step, in time units (system mode). The march
is fixed-step
TYPE:
|
max_time
|
Ceiling on the integration horizon, in time units (system mode):
how long to march before giving up on collecting
TYPE:
|
seed
|
Seed for the random initial condition when the system has none
(system mode), so the section is reproducible. Pass .. versionchanged:: 6.0
Was
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
PoincareSection
|
A :class: |
Examples:
>>> section = poincare_section(Rossler(), plane=("y", 0.0, "up"), crossings=500)
>>> section = poincare_section(Rossler(), crossings=500) # section chosen + recorded
>>> section = poincare_section(traj, plane=("z", 25.0)) # from data
Source code in src/tsdynamics/analysis/orbits/poincare.py
42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 | |
Return maps¶
return_map
¶
return_map(
system: Any,
components: int | str = 0,
*,
kind: str = "max",
plane: tuple[Any, ...] | None = None,
direction: int = +1,
n: int = 2000,
final_time: float = 200.0,
dt: float = 0.01,
transient: float = 0.0,
skip_crossings: int = 0,
ic: Any | None = None,
seed: int | None = 0,
**integrate_kwargs: Any,
) -> ReturnMap
First-return map of a recurring observable.
Records a sequence of scalar values :math:v_0, v_1, \dots from the motion
and pairs each with its successor, giving the one-dimensional map
:math:v_{n+1} = F(v_n) that organises the dynamics.
| PARAMETER | DESCRIPTION |
|---|---|
system
|
What to read the observable from. A continuous
:class:
TYPE:
|
components
|
Which state component to record (names allowed when the system /
trajectory declares |
kind
|
TYPE:
|
plane
|
TYPE:
|
direction
|
Crossing-direction filter (
TYPE:
|
n
|
Number of section crossings to collect when integrating a system in
TYPE:
|
final_time
|
Integration horizon, in time units, when
TYPE:
|
dt
|
Step, in time units, when
TYPE:
|
transient
|
Dynamics discarded before recording extrema, in time units — the unit
TYPE:
|
skip_crossings
|
Number of leading crossings discarded before recording
(
TYPE:
|
ic
|
Initial state when
TYPE:
|
seed
|
Seed for the random initial condition when the system has none, so the
map is reproducible. Pass .. versionchanged:: 6.0
Was
TYPE:
|
**integrate_kwargs
|
Forwarded to
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
ReturnMap
|
The recorded |
Examples:
>>> rm = return_map(Lorenz(), "z", kind="max", final_time=400.0, transient=40.0)
>>> x, y = rm.flat() # the cusp map z_n -> z_{n+1}
>>> rm = return_map(Rossler(), 0, kind="poincare", plane=(0, 0.0), n=400)
Source code in src/tsdynamics/analysis/orbits/return_map.py
261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 | |
ReturnMap
dataclass
¶
ReturnMap(
current: ndarray = (lambda: empty(0))(),
successor: ndarray = (lambda: empty(0))(),
values: ndarray = (lambda: empty(0))(),
times: ndarray = (lambda: empty(0))(),
observable: int = 0,
kind: str = "max",
*,
meta: Mapping[str, Any] = dict(),
)
Bases: AnalysisResult
Result of :func:return_map.
An :class:~tsdynamics.analysis._result.AnalysisResult, so it carries
.meta / the readout repr / .to_dict() / the .plot seam. Iterate for
(current, successor) pairs, or use :meth:flat for the scatter-ready
arrays.
| ATTRIBUTE | DESCRIPTION |
|---|---|
current, successor |
The return map itself — :math:
TYPE:
|
values |
The raw observable series the map was built from: every recorded
:math:
TYPE:
|
times |
When each entry of :attr:
TYPE:
|
kind |
TYPE:
|
flat
¶
cobweb
¶
Describe the return map's cobweb (staircase) as a :class:PlotSpec.
The cobweb diagram traces the iteration :math:v_{n+1} = F(v_n) as a
staircase: from a point on the diagonal it steps vertically to the return
curve, then horizontally back to the diagonal, and repeats. This emits a
COBWEB spec carrying the scatter of the return points
:math:(v_n, v_{n+1}), the diagonal :math:v_{n+1} = v_n, and the
staircase LINE itself built from the recorded sequence. The
:mod:tsdynamics.viz.spec import is lazy, so building a spec never pulls
a plotting library.
| PARAMETER | DESCRIPTION |
|---|---|
kind
|
Override the semantic kind.
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
PlotSpec
|
|
Source code in src/tsdynamics/analysis/orbits/return_map.py
Lyapunov quantifiers¶
lyapunov_spectrum
¶
lyapunov_spectrum(
system: Any,
*,
k: int | None = None,
final_time: float | None = None,
n: int | None = None,
transient: float | None = None,
dt: float | None = None,
ic: Any | None = None,
solver: str | None = None,
rtol: float | None = None,
atol: float | None = None,
backend: str | None = None,
reortho_interval: int | None = None,
d0: float | None = None,
steps_per: int | None = None,
**aliases: Any,
) -> LyapunovSpectrum
Lyapunov spectrum — lambda_1 > 0 is the chaos test.
The one entry point for every family, and since v6 the one entry point
for the maximal exponent too: ask for k=1.
Dispatches to the family implementation (QR tangent dynamics for maps, the
extended variational system on the engine for ODEs, the engine
function-space estimator for DDEs), translating this one signature to each
family's native keywords. The exponents are obtained by Benettin
renormalisation of an evolving orthonormal frame (Benettin et al. 1980).
A system with no right-hand side to differentiate — a
:class:~tsdynamics.families.WrappedSystem around an opaque external
stepper — is answered instead by Benettin two-trajectory rescaling, which
forms no Jacobian; it resolves one exponent, so k > 1 is refused there.
.. versionchanged:: 6.0
Absorbed max_lyapunov, which was a second public door onto the same
question and answered it with a different number: on Hénon at one
nominal horizon it returned 0.4233 where this function returned 0.4160.
The two things it had that this one lacked came with it — the map
burn-in (see transient, now on by default) and the
two-trajectory machine for a system with no Jacobian.
| PARAMETER | DESCRIPTION |
|---|---|
system
|
A flow (ODE/DDE) or a discrete map.
TYPE:
|
k
|
Number of exponents to compute. Defaults to
TYPE:
|
final_time
|
Averaging-window length for a flow (after the transient). Mutually
exclusive with
TYPE:
|
n
|
Number of iterations for a map. Mutually exclusive with
TYPE:
|
**aliases
|
TYPE:
|
rtol
|
Solver tolerances for the flow / DDE variational integration (they have no meaning for a map, whose tangent iteration is exact arithmetic).
TYPE:
|
atol
|
Solver tolerances for the flow / DDE variational integration (they have no meaning for a map, whose tangent iteration is exact arithmetic).
TYPE:
|
transient
|
Dynamics discarded before the exponents are accumulated — time units
for a flow (the same unit as A Lyapunov exponent is a property of the attractor, so the orbit has
to be on it first; a map's default is
:data: .. versionchanged:: 6.0
A map used to refuse this word — "its QR iteration
reorthonormalises from the initial condition, so there is nothing to
discard", which confuses the tangent frame with the base orbit. The
burn-in is folded in from the retired
TYPE:
|
dt
|
Integration step for the variational system, in time units (flows only; a map iterates). Defaults to the family's own step.
TYPE:
|
ic
|
Initial condition. Falls back to
TYPE:
|
solver
|
Numerical kernel for the variational integration (continuous flows
only), e.g. .. versionchanged:: 6.0
Was
TYPE:
|
backend
|
Which engine runs the tangent dynamics. Named here because this free function is the ONLY door in v6 (ruling A2): a keyword the removed method accepted and this one refuses is a signature bug, not a curation (C1).
TYPE:
|
reortho_interval
|
Iterations between Gram--Schmidt reorthonormalisations (maps only).
TYPE:
|
d0
|
Separation the two trajectories are reset to after every rescaling —
only for a system with no right-hand side to differentiate, which is
the one case answered that way. Default It is a real knob, not a formality: on an exactly-linear expanding flow
the growth rate is independent of the perturbation size, and a .. versionadded:: 6.0
Carried over from the retired
TYPE:
|
steps_per
|
Protocol steps between rescalings of that same two-trajectory loop
(default 10). Refused, by name, for a system that has a tangent frame —
a frame is reorthonormalised on its own cadence
( .. versionadded:: 6.0
Carried over from the retired
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
LyapunovSpectrum
|
The exponents (largest first), a drop-in for the bare |
| RAISES | DESCRIPTION |
|---|---|
InvalidInputError
|
If |
ValueError
|
If |
Examples:
>>> lyapunov_spectrum(Lorenz(), final_time=300.0) # [0.91, ~0, -14.57]
>>> lyapunov_spectrum(Henon(), k=2, n=5000) # [0.42, -1.62]
References
G. Benettin, L. Galgani, A. Giorgilli & J.-M. Strelcyn, "Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them", Meccanica 15 (1980) 9--20 (Part 1) and 21--30 (Part 2).
Source code in src/tsdynamics/analysis/lyapunov/__init__.py
445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 | |
lyapunov_from_data
¶
lyapunov_from_data(
data: ndarray,
*,
dt: float | None = None,
dimension: int = 5,
delay: int | None = None,
theiler: int | None = None,
k_max: int | None = None,
eps: float | None = None,
n_neighbors: int = 1,
method: str = "kantz",
fit: tuple[int, int] | None = None,
) -> LyapunovFromData
Estimate the maximal Lyapunov exponent from a time series.
Reconstructs an m-dimensional delay embedding of series, measures how
fast nearby trajectories diverge as a function of the look-ahead k, and
reads the exponent off the slope of the resulting log-divergence curve.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
1-D scalar series, or 2-D
TYPE:
|
dt
|
Sampling interval (time between consecutive samples); the exponent is
reported per unit of
TYPE:
|
dimension
|
Embedding dimension. Should be large enough to unfold the attractor
(Takens' sufficient condition is
TYPE:
|
delay
|
Embedding delay, in samples.
TYPE:
|
theiler
|
Theiler window (Theiler 1986): neighbours with
TYPE:
|
k_max
|
Number of forward samples over which divergence is tracked; the curve
spans
TYPE:
|
eps
|
Neighbour-ball radius for
TYPE:
|
n_neighbors
|
Minimum neighbours a reference point needs to contribute (Kantz only).
TYPE:
|
method
|
Divergence estimator (see module docstring).
TYPE:
|
fit
|
Inclusive sample range |
| RETURNS | DESCRIPTION |
|---|---|
LyapunovFromData
|
The estimated exponent, the full divergence curve, and the parameters
used. Casts to |
| WARNS | DESCRIPTION |
|---|---|
ScalingRegionWarning
|
When the requested The other untrustworthy outcome, a curve with no stable-slope plateau
(all transient and/or all saturation, so the slope is a fallback fit over
the whole curve), is reported through |
| RAISES | DESCRIPTION |
|---|---|
InvalidParameterError
|
If a reconstruction parameter is invalid ( |
ConvergenceError
|
If no usable neighbour is found (no reference point has a neighbour
within |
Notes
The estimate is only as good as the embedding and the chosen scaling region.
Always look at result.times vs result.divergence: a trustworthy
estimate comes from a clear straight segment before the curve saturates.
Measured at the defaults against independent truth (logistic r = 4, where
the exponent is exactly :math:\ln 2; Hénon, Lorenz, Chen and Rössler
against variational spectra), the automatic region lands within 1 % on a
map and within −22 % to +9 % on a flow, depending on the observable and
the sampling rate — a Kantz/Rosenstein reading of a scalar flow observable is
a ~20 %-accurate quantity, not a precision one. Treat it as "is there a
positive exponent, and roughly how large", and use
:func:~tsdynamics.analysis.lyapunov.lyapunov_spectrum when the system (not
just a measured series) is in hand.
Examples:
>>> import tsdynamics as ts
>>> traj = ts.systems.Henon().run(6000, transient=500, ic=[0.1, 0.1])
>>> res = ts.analysis.lyapunov_from_data(traj.y[:, 0], dimension=4, k_max=12, fit=(0, 6))
>>> 0.30 < float(res) < 0.55 # ≈ 0.42
True
References
H. Kantz, "A robust method to estimate the maximal Lyapunov exponent of a time series", Physics Letters A 185 (1994) 77--87.
M. T. Rosenstein, J. J. Collins & C. J. De Luca, "A practical method for calculating largest Lyapunov exponents from small data sets", Physica D 65 (1993) 117--134.
Source code in src/tsdynamics/analysis/lyapunov/from_data.py
758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 | |
LyapunovFromData
dataclass
¶
LyapunovFromData(
estimate: float = 0.0,
stderr: float = 0.0,
abscissa: ndarray = (lambda: empty(0))(),
ordinate: ndarray = (lambda: empty(0))(),
fit_region: tuple[int, int] = (0, 0),
intercept: float = 0.0,
embedding_dim: int = 0,
delay: int = 0,
theiler: int = 0,
n_reference: int = 0,
method: str = "kantz",
decorrelation: int = 0,
trusted: bool = True,
*,
meta: Mapping[str, Any] = dict(),
)
Bases: ScalingResult
Outcome of :func:lyapunov_from_data: the divergence curve and its slope.
A :class:~tsdynamics.analysis._result.ScalingResult — the maximal Lyapunov
exponent is read off the slope of the stretching curve, the same shape every
fractal dimension and embedding diagnostic share — so it inherits the canonical
estimate / abscissa / ordinate / fit_region schema, the result
surface (.meta / the readout repr / .to_dict() / the .plot seam) and
float(result) (the exponent). The domain names for the curve are
:attr:lyapunov, :attr:times and :attr:divergence.
| ATTRIBUTE | DESCRIPTION |
|---|---|
lyapunov |
Estimated maximal Lyapunov exponent (per unit time), the slope of
:attr:
TYPE:
|
times |
Relative times
TYPE:
|
divergence |
The stretching curve
TYPE:
|
fit_region |
Inclusive index range into the curve used for the slope. |
embedding_dim, delay, theiler |
Reconstruction parameters actually used.
TYPE:
|
method |
TYPE:
|
trusted |
TYPE:
|
Notes
Four further quantities are carried and readable but kept off dir()
(contract §11.3): abscissa / ordinate (R2 — :attr:times and
:attr:divergence are the names this estimator's curve is known by, and the
ones the docs use), and n_reference / decorrelation, the two
estimator-internal counts whose readable combination is the public
:attr:independent_windows. All four remain in to_dict().
independent_windows
property
¶
independent_windows: float
How many independent embedding windows the record holds.
n_reference / max((m - 1) * τ, τ_dec) — the reference points divided
by the span one embedding vector covers, floored at the series' own
decorrelation time. This is the record-length half of :attr:trusted,
and it is the half that was missing: measured on 2000 samples of
Lorenz x (truth λ = 0.9056) the estimator returned 1.2414 —
37% high — with R² = 0.99872 and trusted = True, because
:math:R^2 measures how straight the fitted line is and says nothing at
all about whether there was enough data behind it.
Why the floor. Without it the denominator is entirely the caller's
to choose, so the one thing the library hands a reader whose call was
refused — "reduce dimension/delay" — also silently restored the
trust flag: on 500 Lorenz samples, dimension=2, delay=1 takes the
span from 48 to 1 and the count from 2.2 to 358.0, a 160x jump in
apparent record length bought by shrinking the reconstruction rather
than by measuring anything. :attr:decorrelation is read off the data,
so it cannot be bought. Where the embedding is auto-chosen the floor
never binds (τ ≈ τ_dec and m = 5, so (m-1)τ is already
4 τ_dec) and no shipped number moves.
fit_efoldings
property
¶
fit_efoldings: float
How many e-folds of separation the fitted window covers.
:math:\lambda \cdot \Delta t_{fit} — the exponent times the width of
the window it was fitted over. A Lyapunov exponent is a growth rate,
so a window the separation barely grows across carries no rate to read;
this is the diagnostic neither :math:R^2 nor the record length can be.
See :data:_MIN_FIT_EFOLDINGS for the measurement that sets the floor.
lyapunov
property
¶
lyapunov: float
The estimated maximal Lyapunov exponent (alias of :attr:estimate).
chaotic
property
¶
chaotic: bool | None
Is this chaotic? None when the fit cannot support a verdict.
The one adjective-named verdict (contract §4.2 rule 10), three-valued on
purpose: None exactly when :attr:trusted is False, so an
unreadable scaling region can never be read as a confident "no".
This estimator is billed as the data-side twin of the "one number that
says chaotic", and it used to be the only one of the three that refused
to say it — verdict was None on every result, while
zero_one_test, which carries no trust flag at all, answered
chaotic (K ~ 1). For a user who came because they could not make
the call themselves, that is backwards.
kaplan_yorke_dimension
¶
kaplan_yorke_dimension(spectrum: Any) -> ScalarResult
Kaplan--Yorke (Lyapunov) dimension from a Lyapunov spectrum.
.. math::
D_{KY} = j + \frac{\lambda_1 + \cdots + \lambda_j}{|\lambda_{j+1}|},
where :math:j is the largest index whose cumulative exponent sum is
non-negative (Kaplan & Yorke 1979). Interpolating between the :math:j-th
and :math:(j+1)-th exponents gives a fractional estimate of the attractor's
information dimension.
When you already hold a :class:LyapunovSpectrum, spectrum.kaplan_yorke
is the same number and the shorter spelling — it calls this. Reach for the
function when the exponents came from somewhere else: a literature spectrum,
a hand-written array, another library's output. That is what it takes, and
why both exist::
ts.analysis.kaplan_yorke_dimension([0.906, 0.0, -14.572]) # Lorenz, from the paper
| PARAMETER | DESCRIPTION |
|---|---|
spectrum
|
Lyapunov exponents (any order; sorted descending internally). Must be
non-empty and finite — an empty spectrum has no dimension to report and a
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
ScalarResult
|
The dimension, a drop-in for its
:math: |
| RAISES | DESCRIPTION |
|---|---|
InvalidParameterError
|
If |
Examples:
>>> float(kaplan_yorke_dimension([0.906, 0.0, -14.57])) # Lorenz
2.062...
>>> float(kaplan_yorke_dimension([-0.5, -1.0])) # stable fixed point
0.0
>>> float(kaplan_yorke_dimension([1.0, -1.0])) # conservative: saturates
2.0
References
J. L. Kaplan & J. A. Yorke, "Chaotic behavior of multidimensional difference equations", in Functional Differential Equations and Approximation of Fixed Points, Lecture Notes in Mathematics 730, Springer (1979) 204--227.
Source code in src/tsdynamics/analysis/lyapunov/__init__.py
327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 | |
Fixed points & periodic orbits¶
fixed_points
¶
fixed_points(
system: Any,
region: Any = None,
*,
n_seeds: int = 200,
tol: float = 1e-12,
max_iter: int = 60,
dedup_tol: float = 1e-06,
method: str = "newton",
lam: float = 0.05,
beta: float = 1.0,
max_c: int | None = None,
seed: int | None = 0,
) -> FixedPointSet
Find the fixed points of a map, or the equilibria of a flow.
Seeds are drawn uniformly from region plus points sampled from a short
orbit; each runs the chosen root finder, and converged roots are deduplicated
and classified by the Jacobian spectrum (maps: |lambda| < 1; flows:
Re lambda < 0).
| PARAMETER | DESCRIPTION |
|---|---|
system
|
A discrete map (fixed points) or a continuous flow (equilibria). Delay and stochastic systems are not supported.
TYPE:
|
region
|
Search region. An explicit region is a hard search domain: seeds are
drawn in it and converged roots outside it are discarded.
TYPE:
|
n_seeds
|
Random seeds in the hull box; the padded box takes half as many again,
and 20 orbit points are added on top. This is a
completeness knob, and the cheap one: multi-start Newton is a heuristic
that can always miss a root, so raise it when the count matters. For a
guaranteed complete set use
TYPE:
|
tol
|
Residual tolerance (
TYPE:
|
max_iter
|
Root-finding iterations per seed.
TYPE:
|
dedup_tol
|
Distance below which two roots are merged.
TYPE:
|
method
|
TYPE:
|
lam
|
Step size of the Schmelcher--Diakonos iteration (
TYPE:
|
beta
|
Regularisation strength of the Davidchack--Lai iteration
(
TYPE:
|
max_c
|
Cap on the number of stabilising matrices tried (
TYPE:
|
seed
|
RNG seed for the multi-start sampling. All randomness (the box seeds and
the burn-in orbit's starting state) is drawn from a local
:class: .. versionchanged:: 6.0
Was
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
FixedPointSet
|
A list-like |
| RAISES | DESCRIPTION |
|---|---|
NotImplementedError
|
If |
ValueError
|
If |
InvalidInputError
|
If |
Examples:
>>> fixed_points(Henon()) # two saddles of the Hénon map
>>> fixed_points(Lorenz()) # the origin and the two C± equilibria
>>> fixed_points(Henon(), region=[(-3, 3), (-3, 3)], method="interval") # rigorous
References
Schmelcher & Diakonos (1997), Phys. Rev. Lett. 78, 4733. Davidchack & Lai (1999), Phys. Rev. E 60, 6172. Krawczyk (1969), Computing 4, 187. Neumaier (1990), Interval Methods for Systems of Equations, CUP.
Source code in src/tsdynamics/analysis/fixedpoints/fixed.py
470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 | |
FixedPoint
dataclass
¶
FixedPoint(
x: ndarray = (lambda: empty(0))(),
eigenvalues: ndarray = (lambda: empty(0))(),
stable: bool = False,
continuous: bool = False,
*,
meta: Mapping[str, Any] = dict(),
)
Bases: _ArrayBacked, AnalysisResult
A fixed point (map) or equilibrium (flow) with its linear stability data.
It IS its point (v6): np.asarray(fp) is the (dim,) state,
fp[0] is a float, fp - other and np.linalg.norm(fp) work, and
len(fp) is the state-space dimension — so fixed_points(sys)[0] hands
back numbers and the class stays invisible, showing itself only in the repr,
which states the answer. The stability data is one dot away and never in the
way: :attr:eigenvalues, :attr:stable, :attr:continuous.
It is also an :class:~tsdynamics.analysis._result.AnalysisResult, so it
carries .meta / the readout repr / .to_dict() / .to_frame() /
the .plot seam.
| ATTRIBUTE | DESCRIPTION |
|---|---|
x |
The point, shape
TYPE:
|
eigenvalues |
Eigenvalues of the Jacobian at
TYPE:
|
stable |
For a map,
TYPE:
|
continuous |
TYPE:
|
eigenvalue_plane
¶
Describe the Jacobian spectrum as an :class:EIGENVALUE_PLANE spec.
Plots the eigenvalues / multipliers of :attr:eigenvalues in the complex
plane. The stability boundary is drawn as a reference geometry: the unit
circle for a map (|λ| = 1) or the imaginary axis for a flow
(Re λ = 0), per :attr:continuous. The :mod:tsdynamics.viz.spec
import is lazy, so building a spec never pulls a plotting library.
| PARAMETER | DESCRIPTION |
|---|---|
kind
|
Override the semantic kind.
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
PlotSpec
|
|
Source code in src/tsdynamics/analysis/fixedpoints/fixed.py
periodic_orbits
¶
periodic_orbits(
system: Any,
period: int | float | None = None,
region: Any = None,
*,
n_seeds: int = 300,
method: str = "newton",
lam: float = 0.05,
beta: float = 1.0,
tol: float | None = None,
max_iter: int | None = None,
dedup_tol: float = 1e-06,
prime: bool = True,
max_c: int | None = None,
seed: int | None = 0,
ic: Any = None,
transient: float = 0.0,
steps_per_period: int = 2000,
n_points: int = 400,
min_amplitude: float = 1e-06,
) -> OrbitSet
Periodic orbits of a map, or a flow's limit cycle.
Solves :math:f^{p}(x) = x by multi-start root finding, recovers each orbit
by forward iteration, filters orbits whose minimal period properly divides
period (prime=True), and merges the cyclic shifts of one orbit.
The default is plain Newton on :math:f^{p}. It was Davidchack--Lai
until v6, on the reasoning that the stabilising transformations reach
unstable orbits Newton misses; measured, that is not what happens at the
periods users actually explore. On the logistic map at r = 4, where the
prime-cycle counts are known exactly (2, 1, 2, 3, 6, 9, 18 for
:math:p = 1 \ldots 7), Newton at the default n_seeds recovers all of
them, and so do "sd" and "dl" — at up to 19x the cost. On
Henon(1.4, 0.3) the two agree on the count at every :math:p \le 7 while
Newton is 112-323x faster (p=7: 0.09 s against 29.0 s), because
"dl" runs every seed against each of the :math:2^d d! stabilising
matrices and each residual is a p-fold monodromy sweep. The cheap axis
for completeness is n_seeds, not the transformation: raise it first, and
reach for "dl" when a high period still comes up short.
| PARAMETER | DESCRIPTION |
|---|---|
system
|
The map, or the flow whose limit cycle to shoot for.
TYPE:
|
period
|
A map: the period |
region
|
Search region, one
TYPE:
|
seed
|
RNG seed for the multi-start / burn-in sampling, so the orbits found are
reproducible. Pass .. versionchanged:: 6.0
Was
TYPE:
|
n_seeds
|
Seeding controls (see :func:
TYPE:
|
dedup_tol
|
Seeding controls (see :func:
TYPE:
|
method
|
Root finder.
TYPE:
|
lam
|
Stabilising-transformation controls (see
TYPE:
|
beta
|
Stabilising-transformation controls (see
TYPE:
|
max_c
|
Stabilising-transformation controls (see
TYPE:
|
tol
|
Residual tolerance
TYPE:
|
max_iter
|
Iterations per seed/matrix.
TYPE:
|
prime
|
Keep only orbits of minimal period
TYPE:
|
| PARAMETER | DESCRIPTION |
|---|---|
ic |
Starting point. A flow: the shooting guess for
TYPE:
|
transient |
Dynamics discarded before the shooting burn-in lands, in time units
(flows only — the unit
TYPE:
|
steps_per_period |
Integration steps per period guess in the shooting/monodromy march (flows only). Higher is more accurate and slower.
TYPE:
|
n_points |
How many points of the converged cycle to return in
:attr:
TYPE:
|
min_amplitude |
Below this peak-to-peak amplitude a converged "cycle" is an equilibrium the shooting collapsed onto, and is rejected (flows only). .. versionchanged:: 6.0 These five were on the signature and in no docstring at all.
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
OrbitSet
|
A list-like |
| RAISES | DESCRIPTION |
|---|---|
InvalidInputError
|
If |
ValueError
|
If |
Examples:
>>> periodic_orbits(Logistic(params={"r": 3.2}), 2) # the stable 2-cycle
>>> periodic_orbits(Logistic(params={"r": 3.83}), 3) # stable node + saddle
Source code in src/tsdynamics/analysis/fixedpoints/periodic.py
455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 | |
PeriodicOrbit
dataclass
¶
PeriodicOrbit(
points: ndarray = (lambda: empty(0))(),
period: int | float = 0,
multipliers: ndarray = (lambda: empty(0))(),
stable: bool = False,
continuous: bool = False,
residual: float = 0.0,
*,
meta: Mapping[str, Any] = dict(),
)
Bases: _ArrayBacked, AnalysisResult
A periodic orbit with its Floquet/multiplier stability data.
It IS its orbit (v6): np.asarray(orbit) is the (n_points, dim)
point array, orbit[0] is the first point, len(orbit) is how many
points it is stored as, and iteration walks them — so
periodic_orbits(sys)[0] hands back numbers and the class stays invisible,
showing itself only in the repr. The stability data is one dot away:
:attr:period, :attr:multipliers, :attr:stable, :attr:residual.
| ATTRIBUTE | DESCRIPTION |
|---|---|
points |
The orbit, shape
TYPE:
|
period |
The (minimal) period — an integer iteration count for a map, the time
|
multipliers |
Stability multipliers: eigenvalues of :math:
TYPE:
|
stable |
TYPE:
|
continuous |
TYPE:
|
residual |
Closure residual
TYPE:
|
eigenvalue_plane
¶
Describe the multiplier spectrum as an :class:EIGENVALUE_PLANE spec.
Plots the stability multipliers in the complex plane against the unit
circle (a map's eigenvalues of :math:Df^{p}, or a flow's Floquet
multipliers, both judged by |μ| < 1). For a flow the trivial
multiplier ≈ 1 along the flow direction is split into its own
distinctly-marked layer — located here for the plot by argmin|μ − 1|
(a presentation heuristic; the stability flag itself uses the more robust
eigenvector-alignment test). The :mod:tsdynamics.viz.spec import is
lazy.
| PARAMETER | DESCRIPTION |
|---|---|
kind
|
Override the semantic kind.
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
PlotSpec
|
|
Source code in src/tsdynamics/analysis/fixedpoints/periodic.py
estimate_period
¶
estimate_period(
data: Any,
*,
dt: float | None = None,
components: int | str | None = None,
method: str = "autocorrelation",
max_delay: int | None = None,
detrend: bool = True,
) -> ScalarResult
Estimate the dominant period of a sampled signal.
Accepts a :class:~tsdynamics.data.Trajectory (the sampling step is read from
its time grid), a 1-D array, or a 2-D array (one row per sample); for
multi-component input, components selects the channel (default: the
highest-variance one).
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The signal.
TYPE:
|
dt
|
Sampling step. For a bare array it sets the time unit (default
TYPE:
|
components
|
Which column of a multi-component input to read, by index or by
name (default: the highest-variance one). It selects a channel, never
a sample: |
method
|
TYPE:
|
max_delay
|
Largest lag considered (
TYPE:
|
detrend
|
Subtract the mean before estimating (default
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
ScalarResult
|
The estimated period in time units ( |
| RAISES | DESCRIPTION |
|---|---|
ValueError
|
If fewer than 8 samples are given, the signal is constant, the
autocorrelation is degenerate or has no zero-crossing, the spectrum has
no dominant frequency, or |
Examples:
References
Box, G. E. P. & Jenkins, G. M. (1970). Time Series Analysis: Forecasting and Control. Holden-Day (autocorrelation method).
Source code in src/tsdynamics/analysis/fixedpoints/periodic.py
1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 | |
Chaos indicators¶
Three literature-validated answers to "is this orbit chaotic?": the Generalized Alignment Index (GALI), the 0--1 test, and Hunt--Ott expansion entropy.
gali
¶
gali(
system: Any,
k: int = 2,
*,
n: int | None = None,
final_time: float | None = None,
dt: float | None = None,
ic: Any | None = None,
transient: float | None = None,
seed: int | None = 0,
n_internal: int = 10,
) -> GALIResult
Compute the GALI\ :sub:k time series of a map or flow.
| PARAMETER | DESCRIPTION |
|---|---|
system
|
The system whose tangent dynamics to track. Delay/stochastic systems are not supported (their tangent space is not finite-dimensional here).
TYPE:
|
k
|
Number of deviation vectors,
TYPE:
|
n
|
Number of iterations (maps). Default 1000.
TYPE:
|
final_time
|
Integration time (flows). Default 100.0.
TYPE:
|
dt
|
Sampling/recording step for flows. Default 0.1. Not valid for maps.
TYPE:
|
ic
|
Initial condition (defaults to the system's resolved IC).
TYPE:
|
transient
|
Dynamics discarded before tracking, in the unit the family advances in:
iterations for a map, time units for a flow — the same rule
TYPE:
|
seed
|
Seed for the random orthonormal initial deviation frame.
TYPE:
|
n_internal
|
RK4 sub-steps per
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
GALIResult
|
|
| RAISES | DESCRIPTION |
|---|---|
InvalidInputError
|
If If an explicit |
InvalidParameterError
|
If |
ConvergenceError
|
Only when |
Notes
For a chaotic orbit GALI\ :sub:k decays as
:math:\exp\!\big[-\sum_{i=2}^{k}(\lambda_1-\lambda_i)\,t\big], so its
log-slope measures the leading Lyapunov-exponent gaps; a regular orbit keeps
it bounded. See the module docstring for the full statement.
Examples:
>>> gali(Henon(), k=2, n=60).is_chaotic() # exponential collapse
True
>>> float(gali(Lorenz(), k=2, final_time=25.0)) # ~0: Lorenz is chaotic
0.0...
References
Skokos, Bountis & Antonopoulos, "Geometrical properties of local dynamics in Hamiltonian systems: The Generalized Alignment Index (GALI) method", Physica D 231 (2007) 30--54.
Source code in src/tsdynamics/analysis/chaos/gali.py
308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 | |
GALIResult
dataclass
¶
GALIResult(
k: int = 0,
times: ndarray = (lambda: empty(0))(),
values: ndarray = (lambda: empty(0))(),
is_discrete: bool = False,
*,
meta: Mapping[str, Any] = dict(),
)
Bases: AnalysisResult
A GALI\ :sub:k time series with the tools to read order vs chaos off it.
An :class:~tsdynamics.analysis._result.AnalysisResult, so it carries
.meta / the readout repr / .to_dict() / the .plot seam.
float(result) is the final value — :math:\approx 1 for a regular orbit,
:math:\to 0 for a chaotic one — so the result drops straight into a
threshold test.
| ATTRIBUTE | DESCRIPTION |
|---|---|
k |
Number of deviation vectors.
TYPE:
|
times |
Time (flows) or iteration index (maps) at each sample.
TYPE:
|
values |
GALI\ :sub:
TYPE:
|
is_discrete |
Whether the underlying system is a map.
TYPE:
|
final
property
¶
final: float
The last GALI value.
| RAISES | DESCRIPTION |
|---|---|
InvalidParameterError
|
If the series is empty (no samples were recorded). |
applicable
property
¶
applicable: bool
Whether enough of the GALI\ :sub:k curve was recorded to read a verdict off.
A verdict is a reading of where the curve went, so it needs a curve:
a single sample reported GALI_2 = 1 at the end regular (GALI
bounded) from one point, which is the initial condition, not a
measurement.
| RETURNS | DESCRIPTION |
|---|---|
bool
|
|
chaotic
property
¶
chaotic: bool | None
Whether the orbit is chaotic, at the shipped collapse threshold.
The one adjective-named spelling of the verdict (contract §4.2 rule 10);
:meth:is_chaotic is the same test with a tunable threshold. None —
never False — when the curve is not :attr:applicable or ends in the
undecided band between the two thresholds.
| RETURNS | DESCRIPTION |
|---|---|
bool or None
|
|
decay_rate
¶
Exponential decay rate of GALI\ :sub:k (a positive number for chaos).
Fits a line to :math:\ln \mathrm{GALI}_k over the samples that lie
above floor (i.e. before the index reaches the floating-point noise
floor) and returns -slope. For a chaotic orbit this estimates the
Lyapunov gap sum :math:(\lambda_1-\lambda_2)+\cdots+(\lambda_1-\lambda_k)
(Skokos et al. 2008); for a regular orbit it is ~0.
| PARAMETER | DESCRIPTION |
|---|---|
floor
|
Ignore samples at or below this value (numerical underflow region).
TYPE:
|
t_min
|
Also ignore samples before this time/iteration (skip the initial alignment transient).
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
float
|
Estimated decay rate; |
Source code in src/tsdynamics/analysis/chaos/gali.py
is_chaotic
¶
Whether GALI\ :sub:k collapsed below your threshold by the end.
This is the one thing :attr:chaotic cannot do, and the reason a
second spelling of one verdict is allowed to exist here. Everywhere
else in the library a classification has exactly one name (contract §4.2
rule 10) and that name answers with the library's own judgement. But
GALI\ :sub:k is a decay test, and how far the index must fall before
you call an orbit chaotic depends on how long you integrated and how
cleanly the system separates — so the threshold is genuinely the
caller's to set, and a closed-vocabulary property has nowhere to put it.
Prefer :attr:chaotic (which knows about :attr:applicable and hedges
with None in the undecided band) and reach for this only when you
are sweeping the threshold or reproducing a published cut::
gali = ts.analysis.gali(system, k=2, final_time=500)
gali.chaotic # the library's verdict, or None
gali.is_chaotic(threshold=1e-8) # yours, on the same curve
Unlike :attr:chaotic this always returns a plain bool: you asked a
sharp question, so it gives a sharp answer and does not consult
:attr:applicable. On a curve that never resolved, that answer is
False because nothing collapsed — check :attr:applicable first if
that distinction matters.
| PARAMETER | DESCRIPTION |
|---|---|
threshold
|
Collapse level for :attr:
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
bool
|
|
See Also
chaotic : the verdict at the shipped thresholds, hedging when undecided. decay_rate : the fitted exponential rate the collapse happened at.
Source code in src/tsdynamics/analysis/chaos/gali.py
zero_one_test
¶
zero_one_test(
system: Any,
*,
components: int | str = 0,
final_time: float | None = None,
n: int | None = None,
dt: float | None = None,
transient: float | None = None,
ic: Any | None = None,
n_c: int = 100,
c_range: tuple[float, float] = (
pi / 5.0,
4.0 * pi / 5.0,
),
n_cut: int | None = None,
seed: int | None = 0,
oversampling: str = "resample",
) -> ZeroOneResult | tuple[ZeroOneResult, ndarray]
Run the 0--1 test for chaos on a system or a measured observable.
| PARAMETER | DESCRIPTION |
|---|---|
system
|
A dynamical system (integrated / iterated internally to produce the
observable, like :func:
TYPE:
|
components
|
Which state component is the scalar observable, by index or by name.
The test reads ONE observable; the default is the first component, so
|
final_time
|
Integration horizon for a flow (system input). Default 1000.0.
TYPE:
|
n
|
Number of iterations for a map / discrete view (system input). Default 5000.
TYPE:
|
dt
|
Sampling / integration step for a flow (system input). Default 0.1.
TYPE:
|
transient
|
Dynamics discarded before recording, in the unit the family advances in:
time units for a flow, iterations for a map (system input) — the
same rule
TYPE:
|
ic
|
Initial condition (system input).
TYPE:
|
n_c
|
Number of random frequencies :math:
TYPE:
|
c_range
|
Interval the frequencies are drawn from. The default avoids the
resonances near :math: |
n_cut
|
Largest displacement lag used in the mean-square displacement. Default
TYPE:
|
seed
|
Seed for the frequency draw (makes :math:
TYPE:
|
oversampling
|
What to do when the observable is oversampled — the classic misuse of
this test, and one it fails silently: an observable sampled many times
per oscillation makes :math:
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
ZeroOneResult
|
The median correlation growth indicator :math: |
| WARNS | DESCRIPTION |
|---|---|
OversamplingWarning
|
When the observable is oversampled and the guard could not (or was asked
not to) fix it — see |
| RAISES | DESCRIPTION |
|---|---|
InvalidParameterError
|
If the observable is shorter than 200 points (too short for the test to
be meaningful); if |
Examples:
>>> zero_one_test(Logistic(params={"r": 4.0}), n=5000) > 0.9 # chaotic
True
>>> x = Logistic(params={"r": 4.0}).run(steps=5000)["x"]
>>> zero_one_test(x) > 0.9 # the data overload
True
>>> lorenz = Lorenz(ic=[1.0, 1.0, 1.0]) # a flow, sampled fine
>>> zero_one_test(lorenz, components=0, dt=0.02, final_time=600.0) > 0.9
True
References
Gottwald & Melbourne, "A new test for chaos in deterministic systems", Proc. R. Soc. Lond. A 460 (2004) 603--611.
Gottwald & Melbourne, "On the implementation of the 0--1 test for chaos", SIAM J. Appl. Dyn. Syst. 8 (2009) 129--145.
Source code in src/tsdynamics/analysis/chaos/zero_one.py
332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 | |
expansion_entropy
¶
expansion_entropy(
system: Any,
region: Any = None,
*,
n_samples: int = 1000,
n: int | None = None,
final_time: float | None = None,
dt: float | None = None,
fit_range: tuple[int, int] | None = None,
seed: int | None = 0,
n_internal: int = 5,
) -> ExpansionEntropyResult
Estimate the expansion entropy :math:H of a map or flow on a region.
| PARAMETER | DESCRIPTION |
|---|---|
system
|
The system whose expansion to measure.
TYPE:
|
region
|
The restricting region :math:
TYPE:
|
n_samples
|
Number of initial conditions sampled uniformly in the region — the
same quantity
TYPE:
|
n
|
Number of iterations (maps). Default 15. (Kept modest: the raw tangent product is not renormalised, so very long horizons overflow.)
TYPE:
|
final_time
|
Integration time (flows). Default 5.0.
TYPE:
|
dt
|
Recording step for flows. Default 0.1. Not valid for maps.
TYPE:
|
fit_range
|
Inclusive index range in the :math: |
seed
|
Seed for the initial-condition sampling.
TYPE:
|
n_internal
|
RK4 sub-steps per
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
ExpansionEntropyResult
|
|
| RAISES | DESCRIPTION |
|---|---|
InvalidInputError
|
If |
InvalidParameterError
|
If |
Examples:
References
Hunt & Ott, "Defining chaos", Chaos 25 (2015) 097618.
Source code in src/tsdynamics/analysis/chaos/expansion.py
176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 | |
ExpansionEntropyResult
dataclass
¶
ExpansionEntropyResult(
estimate: float = 0.0,
stderr: float = 0.0,
abscissa: ndarray = (lambda: empty(0))(),
ordinate: ndarray = (lambda: empty(0))(),
fit_region: tuple[int, int] = (0, 0),
intercept: float = 0.0,
n_samples: int = 0,
n_survivors: int = 0,
*,
meta: Mapping[str, Any] = dict(),
)
Bases: ScalingResult
An expansion-entropy estimate with the growth curve it was read from.
A :class:~tsdynamics.analysis._result.ScalingResult — the entropy is the
slope of :math:\ln E(t) against :math:t — so it inherits the canonical
estimate / abscissa / ordinate / fit_region schema, the result
surface (.meta / the readout repr / .to_dict() / the .plot seam) and
float(result) (the entropy :math:H). The domain names for the curve
are :attr:entropy, :attr:times and ordinate.
| ATTRIBUTE | DESCRIPTION |
|---|---|
entropy |
The estimated expansion entropy :math:
TYPE:
|
times |
The :math:
TYPE:
|
ordinate |
:math:
TYPE:
|
fit_region |
Inclusive |
n_samples |
Number of initial conditions sampled in the region.
TYPE:
|
n_survivors |
How many of them stayed in the region for the whole run.
TYPE:
|
.. versionchanged |
TYPE:
|
applicable
property
¶
applicable: bool
Whether enough trajectories survived the region for the fit to mean anything.
False on a degenerate fit — fewer than two survivors, or a
stderr of exactly zero, which is what an empty survivor set produces.
Before v6 the "indistinguishable from zero" guard was written
se > 0.0 and abs(h) < 2*se, so 0 ± 0 fell straight through to the
most confident branch and a run where nothing survived printed
H0 = 0 ± 0 non-chaotic (H0 ≤ 0).
| RETURNS | DESCRIPTION |
|---|---|
bool
|
|
Fractal dimensions¶
correlation_dimension
¶
correlation_dimension(
data: Any,
*,
theiler: int | str = "auto",
metric: str | float = "euclidean",
radii: ndarray | None = None,
n_radii: int = 24,
c_lo: float | None = None,
c_hi: float = _DEFAULT_C_HI,
min_window: int = 5,
flatness: float = 1.5,
) -> DimensionResult
Grassberger--Procaccia correlation dimension :math:D_2.
Computes the correlation sum, then reads :math:D_2 off the slope of
:math:\log C(r) vs :math:\log r in the automatically selected scaling
region (:func:~tsdynamics.analysis.dimensions._scaling.fit_scaling_region).
Not the same estimator as generalized_dimension(data, 2.0), despite
both being called :math:D_2: this one counts pairwise distances and the
Rényi one counts box occupancies, and they disagree by more than round-off
(measured on a Lorenz orbit: 2.09 vs 1.86). Prefer this one — it is the
standard :math:D_2 estimator, it is far less sensitive to the box
alignment, and it works down to smaller scales on a modest point set. Reach
for :func:generalized_dimension when you want the whole :math:D_q
spectrum on one consistent box partition.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The point set.
TYPE:
|
theiler
|
Theiler window — exclude pairs with :math:
TYPE:
|
metric
|
Distance metric. |
radii
|
Explicit radii. The default grid spans the informative window of
:math:
TYPE:
|
n_radii
|
Number of radii when
TYPE:
|
c_lo
|
Target :math:
TYPE:
|
c_hi
|
Target :math:
TYPE:
|
min_window
|
Minimum number of radii in the fitted scaling region.
TYPE:
|
flatness
|
How flat the fitted scaling region has to be: a window is admitted when
its straight-line residual is within this factor of the flattest window
found. Larger admits wider but less straight regions. (It is not a
solver tolerance — that is what the v6 rename off
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
DimensionResult
|
|
Notes
Two defaults changed in v6, both because the old ones needed overriding to
be right. The radius grid used to span the 1st to the 50th percentile of the
pair-distance distribution — i.e. up to :math:C = 0.5, deep inside the
saturation bend — and the Theiler window used to default to 0, which
counts temporally adjacent samples of a flow as genuine neighbours. On the
Hénon map the two together returned :math:D_2 = 1.169 against the published
:math:1.220 \pm 0.005; with the current defaults the same call returns
1.206, and — unlike before — stays there as the radius grid is refined.
References
P. Grassberger and I. Procaccia, "Characterization of strange attractors", Phys. Rev. Lett. 50, 346 (1983).
J. Theiler, "Spurious dimension from correlation algorithms applied to limited time-series data", Phys. Rev. A 34, 2427 (1986).
Examples:
| RAISES | DESCRIPTION |
|---|---|
InvalidParameterError
|
If fewer than
:data: |
Source code in src/tsdynamics/analysis/dimensions/correlation.py
212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 | |
correlation_sum
¶
correlation_sum(
data: Any,
radii: ndarray | None = None,
*,
theiler: int | str = "auto",
metric: str | float = "euclidean",
n_radii: int = 24,
c_lo: float | None = None,
c_hi: float = _DEFAULT_C_HI,
) -> CorrelationSum
Correlation sum :math:C(r) over a grid of radii — the raw curve.
This is the curve :func:correlation_dimension fits a slope to, exposed
on its own for anyone who wants to choose the scaling region by eye rather
than let the fitter choose it. correlation_dimension(data) carries the
same data in log form (.abscissa / .ordinate), so reach for this one
when you want :math:r and :math:C(r) themselves — a plot, a custom fit, a
comparison against an analytic :math:C(r) — and for the dimension reach
for :func:correlation_dimension.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The point set (a :class:
TYPE:
|
radii
|
Radii at which to evaluate :math:
TYPE:
|
theiler
|
Exclude pairs with :math:
TYPE:
|
metric
|
Distance metric ( |
n_radii
|
Number of radii when
TYPE:
|
c_lo
|
Lower target value of :math:
TYPE:
|
c_hi
|
Upper target value of :math:
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
CorrelationSum
|
A named 2-tuple |
| RAISES | DESCRIPTION |
|---|---|
ValueError
|
If the Theiler window leaves no valid pairs, or |
Source code in src/tsdynamics/analysis/dimensions/correlation.py
generalized_dimension
¶
generalized_dimension(
data: Any,
q: float = 2.0,
*,
scales: ndarray | None = None,
n_scales: int = 20,
min_occupancy: float = _DEFAULT_MIN_OCCUPANCY,
min_resolution: float = _DEFAULT_MIN_RESOLUTION,
min_window: int = 5,
flatness: float = 1.5,
offsets: tuple[float, ...] = _DEFAULT_OFFSETS,
) -> DimensionResult
Generalized (Rényi) dimension :math:D_q by box counting.
The one Rényi estimator: :func:box_counting_dimension is this at
:math:q = 0 and :func:information_dimension is this at :math:q = 1,
both by delegation, so those names cost nothing to learn and can never give a
different answer. :math:D_2 from this function is the Rényi
:math:q = 2, which is a different estimator from
:func:correlation_dimension (box occupancies vs pairwise distances);
prefer that one for :math:D_2 alone.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The point set.
TYPE:
|
q
|
Rényi order.
TYPE:
|
scales
|
Box sizes :math:
TYPE:
|
n_scales
|
Number of box sizes when
TYPE:
|
min_occupancy
|
Minimum mean number of points per occupied box for a scale to be used.
This is the saturation guard: at fewer than a few points per box the box
count is capped by the sample size rather than set by the geometry, and
the curve flattens. (It replaces the pre-v6
TYPE:
|
min_resolution
|
Minimum linear resolution
TYPE:
|
min_window
|
Minimum number of box sizes in the fitted scaling region.
TYPE:
|
flatness
|
How flat the fitted scaling region has to be: a window is admitted when its straight-line residual is within this factor of the flattest window found. (It is not a solver tolerance — hence the v6 rename.)
TYPE:
|
offsets
|
Grid-origin offsets (as fractions of each box side) swept per scale; the
offset with the fewest occupied boxes (the minimal cover) is kept, which
removes the alignment bias of a single fixed grid origin. Pass
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
DimensionResult
|
|
Notes
This is the raw single-order estimator: it cannot check itself, because the
only cheap internal consistency test on a box-counting dimension is that
:math:D_q is non-increasing in :math:q, which needs more than one order.
:func:box_counting_dimension and :func:dimension_spectrum do run that
check. Small :math:q — and :math:q = 0 above all — converge slowly,
because they weight the rarely-visited boxes a finite sample resolves worst.
References
H. G. E. Hentschel and I. Procaccia, "The infinite number of generalized dimensions of fractals and strange attractors", Physica D 8, 435 (1983).
| RAISES | DESCRIPTION |
|---|---|
InvalidParameterError
|
If |
Source code in src/tsdynamics/analysis/dimensions/generalized.py
524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 | |
box_counting_dimension
¶
box_counting_dimension(
data: Any, *, strict: bool = False, **kwargs: Any
) -> DimensionResult
Box-counting (capacity) dimension :math:D_0, with a self-consistency check.
Numerically identical to generalized_dimension(data, 0.0) — it calls
it — so there is one estimator here, under the name the literature uses for
:math:q = 0. Choose this one when :math:D_0 is the quantity you want;
choose :func:generalized_dimension when :math:q is a variable you are
sweeping, and :func:dimension_spectrum when you want the whole curve.
:math:D_0 is the count of occupied boxes,
:math:N(\epsilon) \sim \epsilon^{-D_0} — the :math:q = 0 member of
:func:generalized_dimension, and the hardest order to estimate: weighting
every occupied box equally makes it dominated by the sparsely-visited parts
of the support, which are exactly what a finite sample resolves last. An
under-resolved box count is biased down, and silently so.
So this function does not just return :math:D_0. It also computes
:math:D_1 and :math:D_2 from the same box partition (free — the
occupancies are already there) and checks the exact inequality
:math:D_0 \ge D_1 \ge D_2. A violation means the estimate has not
resolved: the number is still returned, but marked trusted=False and
accompanied by a :class:NonMonotoneSpectrumWarning naming the violation. On an 8000-point Lorenz trajectory, for
instance, this fires: :math:D_0 = 1.75 against :math:D_2 = 2.00, and the
:math:D_0 end is the one that is wrong -- an independent box count of the
same points puts the :math:\log N(\epsilon) slope at ~1.75 and the
:math:q = 2 slope at ~2.00, so it is box counting's notoriously slow
convergence at :math:q = 0, not a property of the attractor.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The point set.
TYPE:
|
strict
|
TYPE:
|
**kwargs
|
Forwarded to :func:
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
DimensionResult
|
|
| RAISES | DESCRIPTION |
|---|---|
ConvergenceError
|
Only when |
Source code in src/tsdynamics/analysis/dimensions/generalized.py
information_dimension
¶
information_dimension(
data: Any, *, strict: bool = False, **kwargs: Any
) -> DimensionResult
Information dimension :math:D_1, with a self-consistency check.
Numerically identical to generalized_dimension(data, 1.0) — it calls
it — under the name the literature uses for :math:q = 1. See
:func:box_counting_dimension for when to reach for which spelling.
The :math:q \to 1 limit of the Rényi family: the slope of the Shannon
information :math:\sum_i p_i \log p_i against :math:\log \epsilon.
(The Rényi formula is a 0/0 at :math:q = 1; the entropy form is its exact
limit and is what is evaluated — see :func:_partition_ordinate.)
Like :func:box_counting_dimension this also computes the neighbouring
orders off the same partition and checks :math:D_0 \ge D_1 \ge D_2.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The point set.
TYPE:
|
strict
|
Raise on a non-monotone spectrum;
TYPE:
|
**kwargs
|
Forwarded to :func:
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
DimensionResult
|
|
| RAISES | DESCRIPTION |
|---|---|
ConvergenceError
|
If the shared-partition spectrum increases with |
Source code in src/tsdynamics/analysis/dimensions/generalized.py
dimension_spectrum
¶
dimension_spectrum(
data: Any,
qs: Any = None,
*,
scales: ndarray | None = None,
n_scales: int = 20,
min_occupancy: float = _DEFAULT_MIN_OCCUPANCY,
min_resolution: float = _DEFAULT_MIN_RESOLUTION,
min_window: int = 5,
flatness: float = 1.5,
offsets: tuple[float, ...] = _DEFAULT_OFFSETS,
strict: bool = False,
) -> dict[float, DimensionResult]
Compute the :math:D_q spectrum over several Rényi orders.
Computes box occupancies once per scale and reuses them across every q,
so the whole spectrum costs barely more than a single :math:D_q. For a
monofractal the spectrum is flat; a decreasing :math:D_q signals
multifractality.
An increasing :math:D_q signals nothing about the attractor at all —
:math:D_q is non-increasing in :math:q for every measure — so it means the
estimate has not resolved.
strict defaults to False here and to True on the single-value
estimators (:func:box_counting_dimension, :func:information_dimension), and
the asymmetry is deliberate. Asking for one number leaves the caller no way to
see that it is unresolved, so those refuse rather than hand back a figure that
looks authoritative. Asking for the whole spectrum returns the very evidence
of the failure — the rising :math:D_q values are right there in the result —
so it warns (:class:NonMonotoneSpectrumWarning) and returns them, which is
what makes the diagnosis possible. Box counting is slow to converge at small
:math:q, so this fires on ordinary inputs (a Lorenz trajectory, a Gaussian
sample); :func:correlation_dimension is the right tool when it does.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The point set.
TYPE:
|
qs
|
Rényi orders. Default:
TYPE:
|
scales
|
Box sizes; default as in :func:
TYPE:
|
n_scales
|
As in :func:
TYPE:
|
min_occupancy
|
As in :func:
TYPE:
|
min_resolution
|
As in :func:
TYPE:
|
min_window
|
As in :func:
TYPE:
|
flatness
|
As in :func:
TYPE:
|
offsets
|
As in :func:
TYPE:
|
strict
|
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
dict[float, DimensionResult]
|
|
| RAISES | DESCRIPTION |
|---|---|
ConvergenceError
|
If the computed spectrum increases with |
Source code in src/tsdynamics/analysis/dimensions/generalized.py
771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 | |
fixed_mass_dimension
¶
fixed_mass_dimension(
data: Any,
*,
ks: ndarray | None = None,
theiler: int | str = "auto",
metric: str | float = "euclidean",
n_ref: int | None = 1500,
n_ks: int = 16,
min_window: int = 5,
flatness: float = 1.5,
seed: int = 0,
) -> DimensionResult
Fixed-mass (nearest-neighbour) dimension.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The point set.
TYPE:
|
ks
|
Neighbour counts (masses) to probe. Default: a log-spaced integer grid
from 1 to
TYPE:
|
theiler
|
Exclude neighbours with :math:
TYPE:
|
metric
|
Distance metric. |
n_ref
|
Number of reference points averaged over (randomly sub-sampled, seeded).
TYPE:
|
n_ks
|
Number of masses when
TYPE:
|
min_window
|
Minimum number of masses in the fitted scaling region.
TYPE:
|
flatness
|
How flat the fitted scaling region has to be: a window is admitted when its straight-line residual is within this factor of the flattest window found. (It is not a solver tolerance — hence the v6 rename.)
TYPE:
|
seed
|
Seed for the reference sub-sample (keeps the estimate reproducible).
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
DimensionResult
|
|
Notes
The ordinate is the digamma :math:\psi(k) rather than :math:\log k. For
a fixed mass the enclosing radius :math:r_k is a random order statistic, and
:math:\langle\log r_k\rangle = D^{-1}\psi(k) + \text{const} holds without the
:math:O(1/k) bias that :math:\log k carries at small :math:k (Grassberger
1985; the digamma correction of the Kozachenko–Leonenko nearest-neighbour
estimators).
References
R. Badii and A. Politi, "Statistical description of chaotic attractors: The dimension function", J. Stat. Phys. 40, 725 (1985).
P. Grassberger, "Generalizations of the Hausdorff dimension of fractal measures", Phys. Lett. A 107, 101 (1985).
Source code in src/tsdynamics/analysis/dimensions/fixedmass.py
73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 | |
DimensionResult
dataclass
¶
DimensionResult(
estimate: float = 0.0,
stderr: float = 0.0,
abscissa: ndarray = (lambda: empty(0))(),
ordinate: ndarray = (lambda: empty(0))(),
fit_region: tuple[int, int] = (0, 0),
intercept: float = 0.0,
kind: str = "",
q: float | None = None,
trusted: bool = True,
*,
meta: Mapping[str, Any] = dict(),
)
Bases: ScalingResult
A fractal-dimension estimate with the log--log curve it was read from.
Returned by every estimator in this subpackage. A
:class:~tsdynamics.analysis._result.ScalingResult — the dimension is the
fitted slope of a log--log curve — so it inherits the canonical estimate /
abscissa / ordinate / fit_region schema, the result surface
(.meta / the readout repr / .to_dict() / the .plot seam) and behaves
as the dimension number (float(result) and comparisons). The
domain-named :attr:dimension is the spelling for the number.
| ATTRIBUTE | DESCRIPTION |
|---|---|
dimension |
The estimated dimension (the fitted slope).
TYPE:
|
stderr |
Standard error of the slope over the selected scaling region.
TYPE:
|
kind |
Which estimator produced it (
TYPE:
|
abscissa, ordinate |
The log--log curve the slope was fitted to (log-radius vs log-C for the correlation sum; log-scale vs partition ordinate for the generalized dimensions; mean-log-radius vs log-mass for fixed mass).
TYPE:
|
fit_region |
Inclusive |
intercept |
Intercept of the fitted line.
TYPE:
|
q |
Rényi order, for the generalized dimensions (
TYPE:
|
.. versionchanged |
TYPE:
|
unembedded
property
¶
unembedded: bool
Whether this dimension was measured on a single coordinate.
See :class:UnembeddedSeriesWarning: the estimate is then ~1 by
construction, so it is reported untrusted and the repr says why.
x
property
¶
x: ndarray
The abscissa of the scaling curve (alias of :attr:abscissa; see :attr:kind).
y
property
¶
y: ndarray
The ordinate of the scaling curve (alias of :attr:ordinate; see :attr:kind).
local_slopes
property
¶
local_slopes: ndarray
Pointwise local slope of the log--log curve (the diagnostic plateau).
scaling_window
property
¶
The (x_lo, x_hi) abscissa span of the selected scaling region.
Choosing the window yourself
Reading this and wanting to set it is the usual next step, and the
keyword is named here rather than guessed at, because it differs from the
one :func:~tsdynamics.analysis.lyapunov_from_data uses for the same job:
- here —
c_lo=/c_hi=bound the correlation-sum fraction the window is drawn from, andflatness=says how flat the fitted log–log window has to be; - there —
fit=(lo, hi)takes index bounds straight into the stretching curveS(k).
They are spelled differently because they select on different things: a dimension window is chosen by where the correlation sum is well-sampled, a Lyapunov window by which stretching-curve samples are in the linear regime.
Delay embeddings¶
State-space reconstruction from a scalar (or multivariate) measurement (Takens, 1981): the time-delay map, plus the delay- and dimension-selection heuristics that parameterise it.
embed
¶
embed(
data: Any,
dimension: int | Sequence[int] | None = None,
delay: int | Sequence[int] | None = None,
*,
components: int | str | None = None,
**renamed: Any,
) -> Embedding
Time-delay embedding of a scalar series (or a multivariate bundle).
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The source signal. A 1-D series (or a single selected
TYPE:
|
dimension
|
Embedding dimension :math:
TYPE:
|
delay
|
Delay :math:
TYPE:
|
components
|
Select a single channel from a multi-component |
**renamed
|
Not a real parameter: it exists so that a near-miss spelling of one of
the two above —
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
Embedding
|
The delay-coordinate matrix (behaves as an |
| RAISES | DESCRIPTION |
|---|---|
ValueError
|
If |
InvalidParameterError
|
If a keyword is a banned spelling of |
Notes
The matrix is consumable directly by the point-set analyses — e.g.
correlation_dimension(embed(x, m, tau)) estimates :math:D_2 of the
reconstructed attractor. Rows are index-ordered with the original sampling,
so an index-based Theiler window still removes temporally-correlated pairs.
Estimated parameters are recorded, never silent. meta["delay_auto"]
/ meta["dimension_auto"] say which of the two the caller supplied, so a
reconstruction always carries how it was parameterised. Estimating both
costs a mutual-information scan and a Cao sweep; pass either one to skip its
estimate. Auto-selection is for the univariate case (one series has one
:math:\tau and one :math:m) — a multivariate bundle needs its per-channel
values named.
References
F. Takens, "Detecting strange attractors in turbulence", in Dynamical Systems and Turbulence, Lecture Notes in Mathematics 898, 366 (1981).
Examples:
>>> import numpy as np
>>> from tsdynamics.analysis.embedding import embed
>>> x = np.arange(10.0)
>>> y = embed(x, dimension=3, delay=2)
>>> y.shape
(6, 3)
>>> y[0]
array([0., 2., 4.])
Source code in src/tsdynamics/analysis/embedding/embed.py
174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 | |
Delay selection¶
optimal_delay
¶
optimal_delay(
data: Any,
*,
method: str = "mi",
max_delay: int = 50,
bins: int | None = None,
components: int | str | None = None,
) -> CountResult
Recommend an embedding delay :math:\tau (in samples).
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The scalar series (or a selected
TYPE:
|
method
|
TYPE:
|
max_delay
|
Largest lag considered.
TYPE:
|
bins
|
Histogram bins for the mutual-information estimate (
TYPE:
|
components
|
Component selector for a multi-component input. |
| RETURNS | DESCRIPTION |
|---|---|
CountResult
|
The recommended delay (behaves as an |
| RAISES | DESCRIPTION |
|---|---|
ValueError
|
If |
Notes
If no first minimum / crossing is found within max_delay (e.g. a
slowly-decaying curve), the criterion's global fallback is used — the
location of the smallest mutual information, or max_delay for the
autocorrelation rules — so a usable delay is always returned.
The mutual-information criterion is the more widely recommended nonlinear measure (Fraser & Swinney 1986); the autocorrelation rules are the classic linear alternatives.
The first local minimum is taken at the onset of the dip: on a
flat-bottomed (quantised) mutual-information valley the recommended
:math:\tau is the first lag of the plateau, not its far edge — so the
delay is not over-estimated on coarsely sampled curves.
References
A. M. Fraser and H. L. Swinney, "Independent coordinates for strange attractors from mutual information", Phys. Rev. A 33, 1134 (1986).
Examples:
>>> import numpy as np
>>> from tsdynamics.analysis.embedding import optimal_delay
>>> t = np.linspace(0.0, 100.0, 4000)
>>> x = np.sin(2.0 * np.pi * 0.05 * t)
>>> int(optimal_delay(x, method="mi", max_delay=60)) >= 1
True
Source code in src/tsdynamics/analysis/embedding/delay.py
394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 | |
mutual_information
¶
mutual_information(
data: Any,
*,
max_delay: int = 50,
bins: int | None = None,
base: float = e,
components: int | str | None = None,
) -> MutualInformation
Time-delayed mutual information :math:I(\tau) up to max_delay.
The histogram estimator of
.. math::
I(\tau) = \sum_{a,b} p_{ab}(\tau)\,
\log\frac{p_{ab}(\tau)}{p_a\, p_b},
where :math:p_{ab} is the joint distribution of :math:(x_i, x_{i+\tau})
over an equal-width 2-D histogram and :math:p_a, p_b its marginals.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The scalar series (or a selected
TYPE:
|
max_delay
|
Largest lag returned. Clamped to
TYPE:
|
bins
|
Number of histogram bins per axis. Default: a sample-size-dependent
rule,
TYPE:
|
base
|
Logarithm base —
TYPE:
|
components
|
Component selector for a multi-component input. |
| RETURNS | DESCRIPTION |
|---|---|
MutualInformation
|
Behaves as an |
| RAISES | DESCRIPTION |
|---|---|
ValueError
|
If |
References
A. M. Fraser and H. L. Swinney, "Independent coordinates for strange attractors from mutual information", Phys. Rev. A 33, 1134 (1986).
Examples:
>>> import numpy as np
>>> from tsdynamics.analysis.embedding import mutual_information
>>> t = np.linspace(0.0, 100.0, 2000)
>>> x = np.sin(2.0 * np.pi * 0.1 * t)
>>> mi = mutual_information(x, max_delay=40)
>>> int(mi.optimal_lag) >= 1
True
Source code in src/tsdynamics/analysis/embedding/delay.py
186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 | |
autocorrelation
¶
autocorrelation(
data: Any,
*,
max_delay: int = 50,
components: int | str | None = None,
) -> ndarray
Normalised autocorrelation function up to max_delay.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The scalar series (or a selected
TYPE:
|
max_delay
|
Largest lag returned. Clamped to
TYPE:
|
components
|
Component selector for a multi-component input. |
| RETURNS | DESCRIPTION |
|---|---|
(ndarray, shape(max_delay + 1))
|
|
| RAISES | DESCRIPTION |
|---|---|
ValueError
|
If |
Notes
Computed via FFT (Wiener--Khinchin) on the mean-subtracted series. This is the standard biased autocorrelation estimator — each lag is normalised by the full zero-lag variance (not by the shrinking overlap count at that lag) — which is positive-definite and the conventional choice for delay selection.
Source code in src/tsdynamics/analysis/embedding/delay.py
Dimension selection¶
embedding_dimension
¶
embedding_dimension(
data: Any,
*,
method: str = "cao",
delay: int | None = None,
max_dim: int = 10,
components: int | str | None = None,
**kwargs: Any,
) -> EmbeddingDimension
Estimate the minimum embedding dimension by the chosen method.
The one door for the question "how many coordinates does this attractor
need?", with method= naming the estimator — the shape every other
estimator-choosing analysis in the library has (fixed_points(method=),
optimal_delay(method=), lyapunov_from_data(method=)). The two
estimators keep their own names and their own citations, so
embedding_dimension(x, method="cao") and cao_dimension(x) return the
identical object; type whichever reads better at your call site. Reach for
:func:cao_dimension first if you have no reason to prefer otherwise — Cao's
:math:E_1 needs no threshold to be chosen for it, which is the usual
difficulty with :func:false_nearest_neighbors.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The scalar series (or a selected
TYPE:
|
method
|
TYPE:
|
delay
|
Embedding delay in samples. .. versionchanged:: 6.0
Was a hard
TYPE:
|
max_dim
|
Largest dimension evaluated.
TYPE:
|
components
|
Component selector for a multi-component input. |
**kwargs
|
Forwarded to the selected estimator (
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
EmbeddingDimension
|
|
| RAISES | DESCRIPTION |
|---|---|
ValueError
|
If |
See Also
cao_dimension : Cao's averaged false-neighbour estimator. false_nearest_neighbors : Kennel's false-nearest-neighbour estimator.
Source code in src/tsdynamics/analysis/embedding/dimension.py
cao_dimension
¶
cao_dimension(
data: Any,
*,
delay: int | None = None,
max_dim: int = 10,
threshold: float = 0.9,
theiler: int = 0,
components: int | str | None = None,
) -> EmbeddingDimension
Cao's averaged-false-neighbour minimum embedding dimension.
The recommended estimator, and what embedding_dimension(x, method="cao")
(the generic door) calls: Cao's :math:E_1 saturates on its own, with no
rejection threshold to pick, where :func:false_nearest_neighbors needs
rtol / atol chosen for the data.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The scalar series (or a selected
TYPE:
|
delay
|
Embedding delay :math:
TYPE:
|
max_dim
|
Largest dimension :math:
TYPE:
|
threshold
|
Saturation threshold: the estimate is the smallest :math:
TYPE:
|
theiler
|
Exclude temporally-close neighbours with :math:
TYPE:
|
components
|
Component selector for a multi-component input. |
| RETURNS | DESCRIPTION |
|---|---|
EmbeddingDimension
|
|
| RAISES | DESCRIPTION |
|---|---|
ValueError
|
If |
Notes
The d = 1 point of the returned :math:E_1 / :math:E_2 curves
(afn_e1[0] / afn_e2[0]) is diagnostically unreliable. At
:math:d = 1 the base distance :math:R_d is the 1-D nearest-neighbour gap
:math:|x_i - x_n|; on a densely sampled series, coincident-but-distinct
samples make :math:R_d arbitrarily small, so the per-point ratio
:math:a(i, 1) = R_{d+1}/R_d can blow up and inflate the mean by orders of
magnitude. This corrupts only the first point of the diagnostic curve and
does not affect the recommended dimension :math:m, which is read off the
saturation plateau past :math:d = 1 (use a Theiler window of a few
autocorrelation times to mitigate the artefact).
References
L. Cao, "Practical method for determining the minimum embedding dimension of a scalar time series", Physica D 110, 43 (1997).
Examples:
>>> import numpy as np
>>> from tsdynamics.analysis.embedding import cao_dimension
>>> t = np.linspace(0.0, 200.0, 4000)
>>> x = np.sin(t) + 0.5 * np.sin(2.0 * t)
>>> int(cao_dimension(x, delay=10, max_dim=8)) >= 2
True
Source code in src/tsdynamics/analysis/embedding/dimension.py
323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 | |
false_nearest_neighbors
¶
false_nearest_neighbors(
data: Any,
*,
delay: int | None = None,
max_dim: int = 10,
rtol: float = 15.0,
atol: float = 2.0,
threshold: float = 0.01,
theiler: int = 0,
components: int | str | None = None,
) -> EmbeddingDimension
Kennel's false-nearest-neighbour minimum embedding dimension.
The classic estimator, and what embedding_dimension(x, method="fnn")
calls. Choose it over :func:cao_dimension when you want an explicit,
interpretable rejection criterion (the fraction of neighbours that are false)
rather than a saturation curve — at the cost of choosing rtol / atol
for your data.
A neighbour found in dimension :math:d is false when extending to
:math:d+1 either stretches the pair by more than rtol relative to their
:math:d-dimensional distance, or pushes them apart by more than atol
relative to the attractor size :math:R_A (the series' standard deviation) —
the second test catching neighbours whose :math:d-dimensional distance is
essentially zero.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The scalar series (or a selected
TYPE:
|
delay
|
Embedding delay :math:
TYPE:
|
max_dim
|
Largest dimension evaluated.
TYPE:
|
rtol
|
First-criterion tolerance :math:
TYPE:
|
atol
|
Second-criterion tolerance :math:
TYPE:
|
threshold
|
The estimate is the smallest :math:
TYPE:
|
theiler
|
Exclude temporally-close neighbours with :math:
TYPE:
|
components
|
Component selector for a multi-component input. |
| RETURNS | DESCRIPTION |
|---|---|
EmbeddingDimension
|
|
| RAISES | DESCRIPTION |
|---|---|
ValueError
|
If |
References
M. B. Kennel, R. Brown and H. D. I. Abarbanel, "Determining embedding dimension for phase-space reconstruction using a geometrical construction", Phys. Rev. A 45, 3403 (1992).
Examples:
>>> import numpy as np
>>> from tsdynamics.analysis.embedding import false_nearest_neighbors
>>> t = np.linspace(0.0, 200.0, 4000)
>>> x = np.sin(t) + 0.5 * np.sin(2.0 * t)
>>> int(false_nearest_neighbors(x, delay=10, max_dim=8)) >= 2
True
Source code in src/tsdynamics/analysis/embedding/dimension.py
442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 | |
EmbeddingDimension
dataclass
¶
EmbeddingDimension(
dimension: int = 0,
dims: ndarray = (lambda: empty(0))(),
method: str = "",
delay: int = 0,
afn_e1: ndarray | None = None,
afn_e2: ndarray | None = None,
fnn_fraction: ndarray | None = None,
*,
meta: Mapping[str, Any] = dict(),
)
Bases: AnalysisResult
A minimum-embedding-dimension estimate with the curve it was read from.
An :class:~tsdynamics.analysis._result.AnalysisResult, so it carries
.meta / the readout repr / .to_dict() / the .plot seam. It also
behaves as the dimension integer (int(result) drops it straight into
:func:~tsdynamics.analysis.embedding.embed.embed), while carrying the
per-dimension diagnostic so the saturation/decay can be inspected.
| ATTRIBUTE | DESCRIPTION |
|---|---|
dimension |
The recommended minimum embedding dimension :math:
TYPE:
|
dims |
The dimensions :math:
TYPE:
|
method |
TYPE:
|
delay |
The delay (in samples) used to build the reconstructions.
TYPE:
|
afn_e1, afn_e2 |
Cao's :math:
TYPE:
|
fnn_fraction |
Kennel's false-nearest-neighbour fraction per dimension (decays to 0).
TYPE:
|
saturated
property
¶
saturated: bool
Whether the criterion was met inside the scanned range.
False when the selected dimension is the largest one tried, which
means the search hit its ceiling rather than finding saturation — and
the repr then says so instead of presenting the ceiling as the answer
(measured: m = 10 reported with :math:E_1 = 0.6, nowhere near
saturation, with no flag).
| RETURNS | DESCRIPTION |
|---|---|
bool
|
|
Recurrence & RQA¶
Recurrence plots (Eckmann, Kamphorst & Ruelle, 1987) record when a trajectory revisits its own past; recurrence quantification analysis (Marwan et al., 2007) reduces that structure to scalar measures of determinism and laminarity, run globally or in a sliding window.
recurrence_matrix
¶
recurrence_matrix(
data: Any,
*,
threshold: float | None = None,
recurrence_rate: float | None = None,
metric: str | float = "euclidean",
theiler: int = 0,
) -> RecurrenceMatrix
Build a recurrence matrix from a trajectory or point set.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The state points (a :class:
TYPE:
|
threshold
|
Fixed distance threshold :math:
TYPE:
|
recurrence_rate
|
Target matrix density in
TYPE:
|
metric
|
Distance metric ( |
theiler
|
Exclude the near-diagonal band :math:
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
RecurrenceMatrix
|
|
| RAISES | DESCRIPTION |
|---|---|
InvalidInputError
|
If |
InvalidParameterError
|
If neither or both of |
ValueError
|
If the Theiler window leaves no valid pairs. |
References
J.-P. Eckmann, S. O. Kamphorst and D. Ruelle, "Recurrence plots of dynamical systems", Europhys. Lett. 4, 973 (1987).
Source code in src/tsdynamics/analysis/recurrence/matrix.py
260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 | |
RecurrenceMatrix
dataclass
¶
RecurrenceMatrix(
matrix: Any = None,
epsilon: float = 0.0,
metric: str | float = "euclidean",
theiler_window: int = 0,
*,
meta: Mapping[str, Any] = dict(),
)
Bases: AnalysisResult
A binary recurrence matrix with the parameters it was built from.
The matrix is symmetric, sparse and excludes the line of identity (plus the
Theiler band when theiler_window > 0). Pass it straight to
:func:~tsdynamics.analysis.rqa for quantification, or read
:attr:recurrence_rate / densify with :meth:toarray for inspection.
| ATTRIBUTE | DESCRIPTION |
|---|---|
matrix |
The :math:
TYPE:
|
epsilon |
The distance threshold actually used.
TYPE:
|
metric |
The metric the threshold is measured in. |
theiler_window |
Excluded near-diagonal band :math:
TYPE:
|
size
property
¶
size: int
Number of states N (the matrix is N x N); 0 when unpopulated.
The field defaults to None (it carries no matrix until an estimator
fills it in), so this reports 0 rather than raising — the repr reads
it, and a repr must never be the thing that raises in a console.
recurrence_rate
property
¶
recurrence_rate: float
Matrix density :math:RR = \#\{R_{ij}=1\}/N^2 (0 when unpopulated).
to_frame
¶
to_frame() -> Any
Return the recurrent pairs as a tidy (nnz, 2) i / j table.
The base :meth:~tsdynamics.analysis._result.AnalysisResult.to_frame
would drop the matrix and tabulate ['epsilon', 'metric',
'theiler_window'] — the settings, with the answer missing (contract
§4.2 rule 8). The matrix is never densified: the stored COO coordinates
are exactly the recurrences, so the frame is nnz rows, not
:math:N^2. The settings ride on frame.attrs.
| RETURNS | DESCRIPTION |
|---|---|
DataFrame
|
Columns |
| RAISES | DESCRIPTION |
|---|---|
ImportError
|
If :mod: |
Source code in src/tsdynamics/analysis/recurrence/matrix.py
rqa
¶
rqa(
data: Any,
*,
threshold: float | None = None,
recurrence_rate: float | None = None,
metric: str | float = "euclidean",
theiler: int = 0,
min_diagonal: int = 2,
min_vertical: int = 2,
) -> RQAResult
Recurrence quantification of a trajectory, series, or recurrence matrix.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
A prebuilt :class:
TYPE:
|
threshold
|
Threshold or target recurrence rate when
TYPE:
|
recurrence_rate
|
Threshold or target recurrence rate when
TYPE:
|
metric
|
Distance metric (ignored when |
theiler
|
Excluded near-diagonal band (ignored when
TYPE:
|
min_diagonal
|
Shortest diagonal line counted toward
TYPE:
|
min_vertical
|
Shortest vertical line counted toward
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
RQAResult
|
The recurrence-quantification measures (RR, DET, LAM, L, L_max, DIV, ENTR, TT, V_max) plus the raw diagonal / vertical line-length histograms. |
| RAISES | DESCRIPTION |
|---|---|
ValueError
|
If |
| WARNS | DESCRIPTION |
|---|---|
UserWarning
|
When |
Notes
L_max and V_max are the longest diagonal / vertical lines over the
full line-length histograms (no min_* filter), matching Marwan et al.
(2007); DIV = 1 / L_max. The min_diagonal / min_vertical cut
applies only to the ratio and mean measures (DET / LAM / L / TT / ENTR).
Line of identity. The diagonal measures are read off the matrix without
:math:R_{ii} (one trivial infinite line). At theiler == 0 the vertical
measures are read off the matrix with it, as Marwan et al. (2007) §3.5
define them — a laminar sojourn is one vertical line through the diagonal,
and dropping the diagonal point bisects it (measured on 60 identical states
followed by 140 separated ones: TT 30.5 and V_max 59 without the
restoration, against the closed-form 60 and 60).
Theiler window and the vertical measures. theiler > 0 removes the
band :math:|i - j| \le w from the matrix, which is what the diagonal
statistics need but punches a hole through every vertical line: LAM /
TT / V_max are then computed on a truncated matrix and are not
comparable to published values. Quantify laminarity at theiler=0 (the
default) and use a nonzero window only for the diagonal measures. The
excluded band contains the line of identity, so :math:R_{ii} is not
restored in that case — fabricating those N points inside the removed
band would make LAM / TT / V_max describe a matrix that is
neither the textbook one nor the one the caller asked for.
References
N. Marwan, M. C. Romano, M. Thiel and J. Kurths, "Recurrence plots for the analysis of complex systems", Phys. Rep. 438, 237 (2007).
Examples:
>>> import numpy as np
>>> import tsdynamics as ts
>>> t = np.linspace(0.0, 100.0, 1000)
>>> emb = ts.analysis.embed(np.sin(t), dimension=2, delay=5)
>>> res = ts.analysis.rqa(emb, recurrence_rate=0.05)
>>> 0.0 <= res.determinism <= 1.0
True
Source code in src/tsdynamics/analysis/recurrence/rqa.py
443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 | |
RQAResult
dataclass
¶
RQAResult(
recurrence_rate: float,
determinism: float,
laminarity: float,
avg_diagonal_length: float,
max_diagonal_length: int,
divergence: float,
diagonal_entropy: float,
trapping_time: float,
max_vertical_length: int,
size: int,
epsilon: float,
theiler_window: int,
min_diagonal: int,
min_vertical: int,
diagonal_lengths: ndarray,
vertical_lengths: ndarray,
*,
meta: Mapping[str, Any] = dict(),
)
Bases: AnalysisResult
Recurrence-quantification measures of one recurrence matrix.
| ATTRIBUTE | DESCRIPTION |
|---|---|
recurrence_rate |
TYPE:
|
determinism |
TYPE:
|
laminarity |
TYPE:
|
avg_diagonal_length |
TYPE:
|
max_diagonal_length |
TYPE:
|
divergence |
TYPE:
|
diagonal_entropy |
TYPE:
|
trapping_time |
TYPE:
|
max_vertical_length |
TYPE:
|
size |
Number of states
TYPE:
|
epsilon |
Threshold the matrix was built with.
TYPE:
|
theiler_window |
Excluded near-diagonal band.
TYPE:
|
min_diagonal, min_vertical |
Minimum line lengths counted as diagonal / vertical lines.
TYPE:
|
diagonal_lengths, vertical_lengths |
Raw line-length histograms (every run, before the
TYPE:
|
applicable
property
¶
applicable: bool
Whether there were any recurrence points to quantify at all.
False when the matrix holds no recurrence point (RR = 0) or no
diagonal line was counted: DET/LAM/ENTR are then 0 by
vacuity, not by measurement. Before v6 that vacuum printed the full
verdict — rqa(np.random.normal(size=300), threshold=1e-12) returned
DET = 0.000 … stochastic (few diagonal lines), a confident reading of
a matrix with zero points — the same defect
:attr:~tsdynamics.analysis.results.WadaResult.applicable was introduced
for, one subpackage away.
| RETURNS | DESCRIPTION |
|---|---|
bool
|
|
windowed_rqa
¶
windowed_rqa(
data: Any,
*,
window: int | None = None,
step: int | None = None,
threshold: float | None = None,
recurrence_rate: float | None = None,
metric: str | float = "euclidean",
theiler: int = 0,
min_diagonal: int = 2,
min_vertical: int = 2,
) -> WindowedRQA
RQA in a sliding window: nonstationarity in time.
Runs :func:~tsdynamics.analysis.rqa over overlapping windows.
| PARAMETER | DESCRIPTION |
|---|---|
data
|
The state points (or a 1-D series).
TYPE:
|
window
|
Window length in samples (
TYPE:
|
step
|
Stride between windows (default:
TYPE:
|
threshold
|
Exactly one; passed to each window's matrix (see
:func:
TYPE:
|
recurrence_rate
|
Exactly one; passed to each window's matrix (see
:func:
TYPE:
|
metric
|
Distance metric. |
theiler
|
Excluded near-diagonal band, applied within each window. Same caveats as
:func:
TYPE:
|
min_diagonal
|
Minimum line lengths (see :func:
TYPE:
|
min_vertical
|
Minimum line lengths (see :func:
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
WindowedRQA
|
|
| RAISES | DESCRIPTION |
|---|---|
InvalidInputError
|
If |
InvalidParameterError
|
If |
Source code in src/tsdynamics/analysis/recurrence/windowed.py
225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 | |
WindowedRQA
dataclass
¶
WindowedRQA(
centers: ndarray = (lambda: empty(0))(),
results: tuple[RQAResult, ...] = (),
window: int = 0,
step: int = 0,
*,
meta: Mapping[str, Any] = dict(),
)
Bases: AnalysisResult
RQA measures over a sliding window.
Each scalar measure of :class:~tsdynamics.analysis.RQAResult is available as
a per-window array of the same length as :attr:centers (the window-centre
sample indices), e.g. windowed.determinism or
windowed.measure("laminarity").
| ATTRIBUTE | DESCRIPTION |
|---|---|
centers |
Window-centre positions in sample units (
TYPE:
|
results |
The per-window results, in order. |
window |
Window length in samples.
TYPE:
|
step |
Stride between consecutive windows in samples.
TYPE:
|
details
property
¶
The per-window readouts, in order — the objects [] no longer hands back.
w.details[3] is window 3's full
:class:~tsdynamics.analysis.results.RQAResult, repr, verdict and all.
measures
property
¶
Names of the columns of :meth:table / np.asarray(self) / self[i].
table
¶
table() -> ndarray
Return every measure over every window as one (n_windows, 9) array.
Columns are :attr:measures, in that order.
| RETURNS | DESCRIPTION |
|---|---|
ndarray
|
|
Source code in src/tsdynamics/analysis/recurrence/windowed.py
measure
¶
Return one RQA measure as an array over windows.
| PARAMETER | DESCRIPTION |
|---|---|
name
|
Any scalar attribute of :class:
TYPE:
|
Source code in src/tsdynamics/analysis/recurrence/windowed.py
Attractors & basins¶
The global stability picture of a multistable system. Attractors are located by following trajectories through a cell tessellation until they recurrently revisit cells (Datseris & Wagemakers, 2022); the basin of each is the set of initial conditions reaching it. Basin stability (Menck et al., 2013) is an attractor's share of a sampled region; basin entropy (Daza et al., 2016) and the uncertainty exponent (Grebogi et al., 1983) quantify how fractal the boundaries are; continuation tracks attractors and their basins across a parameter.
attractors
¶
attractors(
system: Any,
region: Grid
| Box
| Ball
| Sequence[tuple[float, ...]]
| None = None,
*,
resolution: int | tuple[int, ...] = 100,
n_seeds: int = 1000,
seed: int | None = 0,
dt: float = 1.0,
max_steps: int = 10000,
merge_tol: float | None = None,
**fsm: Any,
) -> AttractorSet
Find the attractors a system has within a region, via recurrences.
Tessellate region into cells, draw n_seeds random initial conditions
from it, and follow each until it settles into a recurrent cell set (a new
attractor) or inherits one already found (Datseris & Wagemakers, 2022).
| PARAMETER | DESCRIPTION |
|---|---|
system
|
A discrete map or continuous flow. Delay and stochastic systems are not supported (their state is not a finite-dimensional point).
TYPE:
|
region
|
Where to sample initial conditions and the box the recurrence cells cover.
TYPE:
|
resolution
|
Recurrence cells per axis when
TYPE:
|
n_seeds
|
Number of random initial conditions to classify.
TYPE:
|
seed
|
Seed for the initial-condition sampler (reproducible).
TYPE:
|
dt
|
Integration step between cell checks for a flow (ignored for a map).
TYPE:
|
max_steps
|
Per-seed step cap before declaring divergence.
TYPE:
|
merge_tol
|
Merge attractors whose centroids lie within this distance (a split-set
cleanup).
TYPE:
|
**fsm
|
Finite-state-machine thresholds forwarded to :class:
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
AttractorSet
|
The located attractors plus how many seeds diverged. |
| RAISES | DESCRIPTION |
|---|---|
TypeError
|
If |
| WARNS | DESCRIPTION |
|---|---|
UserWarning
|
When a located set fails the invariance audit and is discarded — the
recurrence predicate was fooled by slow motion ( |
Notes
Every located attractor is verified before it is returned
(:func:audit_attractors). The bare recurrence predicate — mx_fnd
consecutive steps into already-visited cells — is also satisfied by an orbit
that merely moves less than one cell per dt, so a bottleneck, a slow
manifold or the crawl past an unstable spiral would otherwise be reported as
an attractor (measured: 2–4 phantom attractors on :math:\dot x = \mu + x^2,
which provably has none; 2–4 fragments of the one Lorenz attractor at
dt <= 0.1). The audit re-marches from each located set, fuses mutually
reachable fragments and drops sets the orbit leaves and never returns to.
The seed march is sequential by design, not a parallelism oversight: each
seed is followed cell-by-cell, and the persistent cell labels (_att_cells /
_bas_cells) accumulated by earlier seeds let later seeds settle cheaply by
reaching an already-labelled cell. That shared, order-dependent labelling state
is what makes the sweep amortise — and is exactly why the march cannot be
parallelised without changing the result or the determinism.
On a supported engine run (an ODE flow or a map whose _step lowers, on the
interp / jit backend) the whole per-IC march now runs in one
sequential Rust kernel call (stream perf/basin-march) — stepping,
cell-binning and the shared-label early-out all in Rust, with no per-dt
Python→FFI round-trip — so it is fast without parallelism. The per-cell-check
numerics are byte-for-byte the released system.step(), so the result is
bit-identical to the pure-Python loop, which stays the fallback (and the oracle)
for reference, a non-lowering _step, and DDE/SDE systems.
Why thread-parallelism is a net loss here, and why the kernel is therefore
serial (measured, so the next reader does not re-derive it). The shared
early-out is the dominant work-saver: on the two-well Duffing 60×60 grid a
serial march takes ~42k engine steps, whereas marching every seed independently
to max_steps (the only way to lift the serial label dependency) takes
~1.4M — a ~34× work inflation that 16 cores cannot recover (the independent
full-march alone clocked ~5× slower than the whole serial run). A
"speculative march in parallel, fold in seed order" scheme would reproduce the
labels bit-for-bit (each seed's state stream is a pure function of its IC) but
pays exactly that 34× over-march, so it is abandoned. Dense-block FFI batching
is not an option either: the FSM checks the cell after every per-dt
step() restart, and an adaptive dense-output block over several dt does
not reproduce those per-dt restart states bit-for-bit (it drifts at ~5e-7),
so it would change the cell sequence and the labels. The march stays serial
because that is the algorithm; the Rust kernel makes that serial march cheap.
References
G. Datseris and A. Wagemakers, "Effortless estimation of basins of attraction", Chaos 32, 023104 (2022).
Source code in src/tsdynamics/analysis/basins/attractors.py
1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1466 1467 1468 1469 1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 | |
AttractorSet
dataclass
¶
AttractorSet(
attractors: dict[int, Attractor] = dict(),
diverged: int = 0,
seeds: int = 0,
*,
meta: Mapping[str, Any] = dict(),
)
Bases: AnalysisResult
The attractors found in a region, indexing to numbers.
aset[0] is the first attractor's (dim,) centre as a plain
:class:numpy.ndarray, equal to np.asarray(aset)[0]; iteration yields
those centres. The :class:Attractor records are :attr:details (by
position) and :meth:by_id (by label), and :attr:centers / :attr:cells /
:attr:ids read the same data column-wise, so nothing here needs a loop::
aset = ts.analysis.attractors(system, region)
aset[0] # array([-1., 0.]) the centre
aset.centers # (n, dim)
aset.cells # (n,) int
aset.details[0] # Attractor #1 centre [-1, 0] · 3 cells
aset.by_id(2) # the attractor LABELLED 2, as a record
| ATTRIBUTE | DESCRIPTION |
|---|---|
attractors |
Located attractors, id → :class: |
diverged |
How many seeds left the region / never settled.
TYPE:
|
seeds |
How many seeds were classified in total.
TYPE:
|
details
property
¶
The :class:Attractor records, in the order [] indexes.
aset[0] is the centre; aset.details[0] is the attractor it came
from, with its repr, its points cloud and its cells count. One
name across every collection in the library, so it is learned once.
centers
property
¶
centers: ndarray
Stack of attractor representatives, shape (n_attractors, dim).
The same array np.asarray(aset) returns: a SET arrays as one
representative per member, while a MEMBER arrays as its own point cloud.
by_id
¶
Return the record of the attractor labelled key.
[] indexes by position and hands back numbers; an id is something you
read off a basin image and then want to know about, so this door hands
back the :class:Attractor itself — the same object :attr:details
holds.
| RAISES | DESCRIPTION |
|---|---|
KeyError
|
If no attractor carries that id. |
Source code in src/tsdynamics/analysis/basins/attractors.py
match
¶
Return the id of the attractor closest to point (or None if empty).
Uses :func:tsdynamics.data.set_distance; point may be a single state
or a point cloud.
Source code in src/tsdynamics/analysis/basins/attractors.py
Attractor
dataclass
¶
Attractor(
id: int = 0,
points: ndarray = (lambda: empty((0, 0)))(),
cells: int = 0,
*,
meta: Mapping[str, Any] = dict(),
)
Bases: _ArrayBacked, AnalysisResult
One located attractor: a point cloud of states sampled on it.
It IS its point cloud (v6): np.asarray(attractor) is the
(m, dim) array of sampled states, attractor[0] is the first of them,
and len(attractor) is how many were sampled — so attractors(sys, …)[0]
hands back numbers and the class stays invisible, showing itself only in the
repr. The label and the coarse size are one dot away: :attr:id,
:attr:center, :attr:cells.
| ATTRIBUTE | DESCRIPTION |
|---|---|
id |
Integer label (
TYPE:
|
points |
States sampled while the trajectory was on the attractor. A fixed point
collapses to one repeated point; a cycle/chaotic set spreads out.
TYPE:
|
cells |
Number of distinct grid cells the attractor occupies (a coarse size).
TYPE:
|
basins
¶
basins(
system: Any,
region: Grid
| Box
| Ball
| Sequence[tuple[float, ...]]
| None = None,
*,
recurrence: Box
| Grid
| Sequence[tuple[float, ...]]
| None = None,
recurrence_resolution: int | tuple[int, ...] = 100,
seed: int | None = 0,
dt: float = 1.0,
max_steps: int = 10000,
merge_tol: float | None = None,
**fsm: Any,
) -> BasinsResult
Classify every point of a grid region by the attractor it converges to.
Each lattice point is followed until it settles into a recurrent cell set (Datseris & Wagemakers, 2022); by default the region doubles as the recurrence tessellation, so labels accumulate and most points settle cheaply by reaching an already-labelled cell.
For a higher-dimensional flow whose basins are viewed on a slice (e.g. a
position grid of a system that also carries velocities), pass a full-dimension
recurrence box covering the whole trajectory range and let region be the
thin slice of initial conditions (free axes pinned with counts == 1).
| PARAMETER | DESCRIPTION |
|---|---|
system
|
A discrete map or continuous flow.
TYPE:
|
region
|
The lattice of initial conditions (full state dimension; a slice pins free
axes with
TYPE:
|
recurrence
|
Full-dimension region whose tessellation recurrences are detected on.
Defaults to |
recurrence_resolution
|
Recurrence cells per axis when
TYPE:
|
seed
|
Accepted for signature uniformity with :func:
TYPE:
|
dt
|
Integration step between cell checks for a flow (ignored for a map).
TYPE:
|
max_steps
|
Per-point step cap before declaring divergence.
TYPE:
|
merge_tol
|
Merge attractors whose centroids lie within this distance.
TYPE:
|
**fsm
|
Finite-state-machine thresholds forwarded to
:class:
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
BasinsResult
|
The labelled basin image and the attractors it refers to. |
| RAISES | DESCRIPTION |
|---|---|
TypeError
|
If |
| WARNS | DESCRIPTION |
|---|---|
UserWarning
|
When a located set fails the invariance audit and is discarded — the
recurrence predicate mistook slow motion for convergence ( |
References
G. Datseris and A. Wagemakers, "Effortless estimation of basins of attraction", Chaos 32, 023104 (2022).
Source code in src/tsdynamics/analysis/basins/basins.py
544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 | |
BasinsResult
dataclass
¶
BasinsResult(
labels: ndarray = (lambda: empty(0, dtype=int))(),
grid: Grid | None = None,
attractors: AttractorSet = AttractorSet(),
*,
meta: Mapping[str, Any] = dict(),
)
Bases: AnalysisResult
A basin diagram: every grid cell labelled by the attractor it reaches.
| ATTRIBUTE | DESCRIPTION |
|---|---|
labels |
Attractor id (
TYPE:
|
grid |
The lattice the labels are laid out on.
TYPE:
|
attractors |
The attractors the labels refer to.
TYPE:
|
fractions
property
¶
Fraction of grid CELLS in each attractor's basin.
A census of this image: every cell of the lattice, counted once, with
the diverged/lost share reported separately by :attr:diverged_fraction.
This is not the same statistic as :func:basin_fractions, despite the
shared word, and the two measurably disagree -- on Van der Pol,
{1: 0.9922} here against {1: 1.0} there. That function is Monte
Carlo basin stability (Menck et al. 2013): an estimate, with a
standard error, of the share of a measure over the region, from
randomly drawn initial conditions. This one is an exact count of a
particular grid, so it depends on the grid and carries no error bar --
and it includes cells that never settled, which the sampler drops.
Use this to read off the picture you just drew; use
:func:basin_fractions to estimate how likely a randomly perturbed state
is to end up on each attractor.
basin_fractions
¶
basin_fractions(
system: Any,
region: Grid
| Box
| Ball
| Sequence[tuple[float, ...]]
| None = None,
*,
n_seeds: int = 10000,
resolution: int | tuple[int, ...] = 100,
seed: int | None = 0,
dt: float = 1.0,
max_steps: int = 10000,
merge_tol: float | None = None,
**fsm: Any,
) -> BasinFractions
Estimate basin stability: each attractor's share of a sampled region.
Draw n random initial conditions from region and classify each; the
fraction converging to an attractor estimates its basin stability (Menck et
al., 2013). The estimate is dimension-free — its standard error
:math:\sqrt{p(1-p)/n} depends only on the fraction and n.
Not the same statistic as basins(system, region).fractions, despite
the shared word: that is an exact census of one grid image (every cell
counted once, diverged cells included in the total), this is a Monte Carlo
estimate over a measure with a standard error. They measurably disagree —
on Van der Pol, {1: 1.0} here against {1: 0.9922} there. Ask this
one "how likely is a random perturbation to land here?"; ask that one "what
does this picture show?".
| PARAMETER | DESCRIPTION |
|---|---|
system
|
A discrete map or continuous flow.
TYPE:
|
region
|
The measure to sample initial conditions from (uniform over a Box/Ball, or a Grid's bounding box).
TYPE:
|
n_seeds
|
Number of random initial conditions drawn from
TYPE:
|
resolution
|
Recurrence cells per axis (a Grid uses its own
TYPE:
|
seed
|
Seed for the sampler (reproducible).
TYPE:
|
dt
|
Integration step between cell checks for a flow (ignored for a map).
TYPE:
|
max_steps
|
Per-sample step cap before declaring divergence.
TYPE:
|
merge_tol
|
Merge attractors whose centroids lie within this distance.
TYPE:
|
**fsm
|
Finite-state-machine thresholds forwarded to
:class:
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
BasinFractions
|
Fractions per attractor, the diverged share, and the attractors. |
| RAISES | DESCRIPTION |
|---|---|
TypeError
|
If |
| WARNS | DESCRIPTION |
|---|---|
UserWarning
|
When a located set fails the invariance audit and is discarded — the
recurrence predicate mistook slow motion for convergence ( |
References
P. J. Menck, J. Heitzig, N. Marwan and J. Kurths, "How basin stability complements the linear-stability paradigm", Nature Physics 9, 89 (2013).
Source code in src/tsdynamics/analysis/basins/basins.py
693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 | |
BasinFractions
dataclass
¶
BasinFractions(
fractions: dict[int, float] = dict(),
diverged: float = 0.0,
n: int = 0,
attractors: AttractorSet = AttractorSet(),
*,
meta: Mapping[str, Any] = dict(),
)
Bases: AnalysisResult
Monte-Carlo basin stability: each attractor's share of a sampled region.
| ATTRIBUTE | DESCRIPTION |
|---|---|
fractions |
Attractor id → fraction of sampled initial conditions converging to it. |
diverged |
Fraction of samples that diverged / never settled.
TYPE:
|
n |
Number of initial conditions sampled.
TYPE:
|
attractors |
The attractors the ids refer to.
TYPE:
|
standard_error
property
¶
Binomial standard error :math:\sqrt{p(1-p)/n} per fraction.
dominant
property
¶
dominant: int | None
Id of the attractor with the largest basin (None if all diverged).
to_frame
¶
to_frame() -> Any
Return a tidy frame: one row per attractor, with its share and error.
The base :meth:~tsdynamics.analysis.results.AnalysisResult.to_frame
builds ONE row of the declared fields, and the answer here lives in a
dict[int, float], which a table cell cannot hold — so it was dropped:
measured, the frame came back (1, 3) of diverged/n/
dominant with no fractions.
| RETURNS | DESCRIPTION |
|---|---|
DataFrame
|
Columns |
Source code in src/tsdynamics/analysis/basins/basins.py
by_id
¶
Return the basin fraction of the attractor labelled key.
[] is positional (the sequence convention every collection here
follows); this is the explicit id lookup, the same pairing
:class:~tsdynamics.analysis.results.AttractorSet uses. The whole
{id: share} mapping is still :attr:fractions.
| RAISES | DESCRIPTION |
|---|---|
KeyError
|
If no attractor carries that id. |
Source code in src/tsdynamics/analysis/basins/basins.py
Boundary structure¶
basin_entropy
¶
basin_entropy(
basins: Any,
*,
box_size: int = 5,
base: float = e,
include_diverged: bool = False,
) -> BasinEntropy
Basin entropy :math:S_b and boundary basin entropy :math:S_{bb}.
Partition the basin image into non-overlapping boxes of box_size cells per
axis. Within box :math:i, with colour fractions :math:p_{ij}, the Gibbs
entropy is :math:S_i = -\sum_j p_{ij}\log p_{ij}. Then
:math:S_b = \langle S_i\rangle over all boxes and
:math:S_{bb} = \langle S_i\rangle over boundary boxes (more than one
colour). :math:S_{bb} > \log 2 is sufficient for a fractal boundary, since
a box straddling a smooth boundary holds at most two colours.
| PARAMETER | DESCRIPTION |
|---|---|
basins
|
The basin image (attractor ids
TYPE:
|
box_size
|
Box side length in cells.
TYPE:
|
base
|
Logarithm base. With the default natural log the fractal threshold is
:math:
TYPE:
|
include_diverged
|
If
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
BasinEntropy
|
|
| RAISES | DESCRIPTION |
|---|---|
ValueError
|
If |
References
A. Daza, A. Wagemakers, B. Georgeot, D. Guéry-Odelin and M. A. F. Sanjuán, "Basin entropy: a new tool to analyze uncertainty in dynamical systems", Scientific Reports 6, 31416 (2016).
Source code in src/tsdynamics/analysis/basins/metrics.py
629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 | |
BasinEntropy
dataclass
¶
BasinEntropy(
sb: float = 0.0,
sbb: float = 0.0,
n_boxes: int = 0,
n_boundary_boxes: int = 0,
box_size: int = 0,
log_base: float = 0.0,
fractal_boundary: bool = False,
*,
meta: Mapping[str, Any] = dict(),
)
Bases: AnalysisResult
Basin entropy and boundary basin entropy of a basin diagram.
| ATTRIBUTE | DESCRIPTION |
|---|---|
sb |
Basin entropy :math:
TYPE:
|
sbb |
Boundary basin entropy :math:
TYPE:
|
n_boxes |
Number of boxes the grid was partitioned into.
TYPE:
|
n_boundary_boxes |
Number of boxes containing more than one basin.
TYPE:
|
fractal_boundary |
TYPE:
|
Notes
box_size (the box side length in cells) and log_base (the base of the
logarithm, e by default) are the caller's own settings echoed back, so
they are carried and exported but kept off dir() — R1, contract §11.3.
Both still resolve and both appear in :meth:to_dict.
uncertainty_exponent
¶
uncertainty_exponent(
basins: Any,
*,
radii: tuple[int, ...] = (1, 2, 3, 4, 5),
cell_size: float | ndarray = 1.0,
include_diverged: bool = False,
) -> UncertaintyExponent
Estimate the uncertainty exponent of a basin boundary.
A cell is :math:\varepsilon-uncertain when at least one axis-neighbour at
distance :math:\varepsilon carries a different basin label. The fraction of
such cells scales as :math:f(\varepsilon)\sim\varepsilon^{\alpha}
(Grebogi et al., 1983); :math:\alpha is the slope of a log-log fit and the
boundary box-counting dimension is :math:D_0 = D - \alpha.
| PARAMETER | DESCRIPTION |
|---|---|
basins
|
The basin image.
TYPE:
|
radii
|
Perturbation radii in cells.
TYPE:
|
cell_size
|
Physical cell spacing (a scalar or per-axis); rescales |
include_diverged
|
If
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
UncertaintyExponent
|
|
| RAISES | DESCRIPTION |
|---|---|
ValueError
|
If fewer than two |
References
C. Grebogi, S. W. McDonald, E. Ott and J. A. Yorke, "Final state sensitivity: an obstruction to predictability", Physics Letters A 99, 415 (1983).
Source code in src/tsdynamics/analysis/basins/metrics.py
752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 | |
UncertaintyExponent
dataclass
¶
UncertaintyExponent(
alpha: float = 0.0,
boundary_dimension: float = 0.0,
state_dimension: int = 0,
epsilons: ndarray = (lambda: empty(0))(),
uncertain_fractions: ndarray = (lambda: empty(0))(),
r_squared: float = 0.0,
*,
meta: Mapping[str, Any] = dict(),
)
Bases: AnalysisResult
The uncertainty exponent of a basin boundary.
| ATTRIBUTE | DESCRIPTION |
|---|---|
alpha |
Uncertainty exponent :math:
TYPE:
|
boundary_dimension |
Box-counting dimension of the boundary, :math:
TYPE:
|
state_dimension |
State-space (grid) dimension :math:
TYPE:
|
epsilons |
Perturbation radii used (in state-space units).
TYPE:
|
uncertain_fractions |
Fraction of :math:
TYPE:
|
r_squared |
Coefficient of determination of the log-log fit.
TYPE:
|
.. versionchanged |
The curve was called
TYPE:
|
D0
property
¶
D0: float
Boundary dimension in Grebogi's casing — the same number as :attr:boundary_dimension.
applicable
property
¶
applicable: bool
Whether enough radii were measured for a power law to be read off.
False below :data:_MIN_UNCERTAINTY_RADII radii: :math:R^2 clears
any acceptance level trivially on two or three points, so the verdict was
gated on a diagnostic that could not fail — measured, a two-radius run
printed final-state sensitive (fractal boundary).
| RETURNS | DESCRIPTION |
|---|---|
bool
|
|
slope_drift
property
¶
slope_drift: float
Ratio of the last local slope to the first, across the window.
A power law has one slope everywhere, so this is 1 for a converged
measurement and falls towards 0 as the uncertain fraction saturates
against its ceiling of 1. It is the diagnostic :math:R^2 cannot be:
:math:R^2 says how straight the fitted line is, not whether the
relation is a power law at all (measured, R² = 0.9926 on a fit whose
local slope had already halved).
| RETURNS | DESCRIPTION |
|---|---|
float
|
|
resolved
property
¶
resolved: bool
Whether the power law is converged enough for :math:\alpha to be read.
False when the local slope decays by more than
:data:_SLOPE_STABLE_RATIO across the probed radii — the uncertain
fraction has saturated, so the fitted :math:\alpha is a lower bound
set by the grid rather than a measurement of the boundary.
final_state_sensitive
property
¶
final_state_sensitive: bool | None
Whether :math:\alpha \ll 1 — a fractal basin boundary.
The one adjective-named spelling of the verdict (contract §4.2 rule 10).
None — never False — when the test does not apply or the power
law itself did not fit.
| RETURNS | DESCRIPTION |
|---|---|
bool or None
|
|
contradicted_by_basin_entropy
property
¶
contradicted_by_basin_entropy: bool
Whether :func:basin_entropy proves fractal where this fit reads smooth.
Only one direction is a contradiction. :math:S_{bb} > \ln 2 is
sufficient for a fractal boundary and not necessary (Daza et al.
2016), so failing it says "not established", never "smooth" — but
passing it while :math:\alpha \approx 1 is two answers to one
question about one image.
wada_property
¶
wada_property(
basins: Any,
*,
radii: tuple[int, ...] = (1, 2, 3, 4, 5),
threshold: float = 0.9,
include_diverged: bool = False,
) -> WadaResult
Test a basin diagram for the Wada property on a grid.
With three or more basins, a boundary cell is "Wada-complete" at radius
:math:r when its Chebyshev-:math:r neighbourhood contains every basin
colour. The fraction :math:W(r) of such boundary cells tends to one for a
genuine Wada boundary (Daza et al., 2015). This is a sufficient grid test,
not a topological proof.
| PARAMETER | DESCRIPTION |
|---|---|
basins
|
The basin image. Diverged cells (
TYPE:
|
radii
|
Chebyshev radii (in cells) at which to grow the neighbourhood.
TYPE:
|
threshold
|
Minimum :math:
TYPE:
|
include_diverged
|
If
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
WadaResult
|
|
| RAISES | DESCRIPTION |
|---|---|
ValueError
|
If |
References
A. Daza, A. Wagemakers, M. A. F. Sanjuán and J. A. Yorke, "Testing for basins of Wada", Scientific Reports 5, 16579 (2015).
Source code in src/tsdynamics/analysis/basins/metrics.py
919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 | |
WadaResult
dataclass
¶
WadaResult(
is_wada: bool = False,
n_basins: int = 0,
radii: ndarray = (lambda: empty(0))(),
fractions: ndarray = (lambda: empty(0))(),
n_boundary_cells: int = 0,
threshold: float = 0.0,
*,
meta: Mapping[str, Any] = dict(),
)
Bases: AnalysisResult
A grid test for the Wada property of a basin diagram.
| ATTRIBUTE | DESCRIPTION |
|---|---|
wada |
The verdict:
TYPE:
|
W |
The measured Wada fraction at the largest radius.
TYPE:
|
applicable |
Whether the test could be run at all (it needs ≥ 3 basins).
TYPE:
|
n_basins |
Number of attractor basins (colours) considered.
TYPE:
|
radii |
Chebyshev radii tested.
TYPE:
|
fractions |
Fraction of boundary cells whose neighbourhood contains every basin, per
radius (:math:
TYPE:
|
n_boundary_cells |
Number of boundary cells.
TYPE:
|
.. versionchanged |
TYPE:
|
applicable
property
¶
applicable: bool
Whether the Wada test could be run on this image at all.
The Wada property is defined as every boundary point being on the
boundary of three or more basins, so the test needs at least three
basins and at least one boundary cell. When either is missing,
:func:wada_property early-returns zeros — and a W = 0 reads as a
measured negative ("we looked, the boundaries are not Wada") when the
truth is that nothing was measured. This flag separates the two, and
:meth:_answer prints the reason instead of the number.
Derived from the fields the result already carries, so no estimator
change: it is exactly the condition wada_property tests before it
early-returns.
wada
property
¶
wada: bool | None
Whether the boundary is Wada — the one adjective-named verdict spelling.
None — never False — when the test does not :attr:apply
<applicable>. This is the spelling: it is the one every other
classifying result uses (contract §4.2 rule 10), and the only one that
distinguishes a measured no from nothing was measured. The raw
stored flag is_wada is still readable (and still exported) but is
off dir(), because it answers False in both cases.
| RETURNS | DESCRIPTION |
|---|---|
bool or None
|
|
resilience
¶
resilience(
result: BasinsResult, attractor_id: int | None = None
) -> ScalarResult
Minimal-fatal-shock resilience of an attractor: its distance to the boundary.
The smallest perturbation that pushes the attractor out of its own basin, estimated as the state-space distance from the attractor's representative to the nearest cell of another basin (Halekotte & Feudel, 2020). Larger means more resilient.
| PARAMETER | DESCRIPTION |
|---|---|
result
|
A basin image carrying its grid and attractors.
TYPE:
|
attractor_id
|
Which attractor (basin label) to measure. A basin image with exactly one attractor needs no choice, so it may be omitted there; with two or more it is required, and the error lists the ids actually present.
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
ScalarResult
|
Distance from the attractor to its basin boundary (behaves as a
|
| RAISES | DESCRIPTION |
|---|---|
TypeError
|
If |
ValueError
|
If the labels are not laid out on the grid (a pre-squeezed slice), or
|
Notes
A sliced basin image (a label cube with one or more degenerate counts == 1
axes, e.g. a position-plane slice of a higher-dimensional flow) is collapsed to
its free axes before the distance transform, so a pinned axis does not inject a
spurious one-cell-away boundary that would cap the reported distance.
References
L. Halekotte and U. Feudel, "Minimal fatal shocks in multistable complex networks", Scientific Reports 10, 11374 (2020). Builds on the integral stability of C. Mitra, J. Kurths and R. V. Donner, Scientific Reports 5, 16196 (2015).
Source code in src/tsdynamics/analysis/basins/metrics.py
1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 | |
Continuation & tipping¶
continuation
¶
continuation(
system: Any,
param: str,
values: Any,
region: Box
| Ball
| Grid
| Sequence[tuple[float, ...]]
| None = None,
*,
n_seeds: int = 2000,
resolution: int | tuple[int, ...] = 100,
seed: int | None = 0,
dt: float = 1.0,
max_steps: int = 10000,
min_fraction: float = 0.0,
match_method: str = "centroid",
match_threshold: float | None = None,
**fsm: Any,
) -> ContinuationResult
Track attractors and basin fractions as param sweeps values.
At each parameter value the attractors are found and their basin fractions
estimated (:func:~tsdynamics.analysis.basins.basin_fractions); consecutive
value's attractors are matched greedily by nearest state-space distance so a
persisting attractor keeps its id and a vanishing one drops out (Datseris,
Rossi & Wagemakers, 2023).
| PARAMETER | DESCRIPTION |
|---|---|
system
|
A discrete map or continuous flow exposing
TYPE:
|
param
|
Name of the parameter to sweep.
TYPE:
|
values
|
Parameter values, in the order to walk them.
TYPE:
|
region
|
The region whose basin fractions are measured at each value.
TYPE:
|
n_seeds
|
Number of random initial conditions sampled at each parameter value
— the word
TYPE:
|
resolution
|
Recurrence cells per axis (a Grid uses its own
TYPE:
|
seed
|
Sampler seed (shared across values for a fair comparison).
TYPE:
|
dt
|
Integration step between cell checks for a flow.
TYPE:
|
max_steps
|
Per-sample step cap.
TYPE:
|
min_fraction
|
Drop attractors whose basin fraction is below this at a given value before
matching — a filter for the tiny spurious sets the recurrence finder can
report near unstable equilibria.
TYPE:
|
match_method
|
Set distance used to match attractors between values.
TYPE:
|
match_threshold
|
Reject a match farther apart than this (so an attractor that jumps is not
spuriously tied to a different one).
TYPE:
|
**fsm
|
Finite-state-machine thresholds forwarded to the recurrence finder.
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
ContinuationResult
|
|
| RAISES | DESCRIPTION |
|---|---|
TypeError
|
If |
References
G. Datseris, K. L. Rossi and A. Wagemakers, "Framework for global stability analysis of dynamical systems", Chaos 33, 073151 (2023).
Source code in src/tsdynamics/analysis/basins/continuation.py
308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 | |
ContinuationResult
dataclass
¶
ContinuationResult(
param: str = "",
values: ndarray = (lambda: empty(0))(),
fractions: dict[int, ndarray] = dict(),
attractors: list[dict[int, Attractor]] = list(),
diverged: ndarray = (lambda: empty(0))(),
*,
meta: Mapping[str, Any] = dict(),
)
Bases: AnalysisResult
Attractors and basin fractions tracked across a parameter sweep.
| ATTRIBUTE | DESCRIPTION |
|---|---|
param |
The swept parameter name.
TYPE:
|
values |
The swept parameter axis, not a measurement — the values of
:attr:
TYPE:
|
fractions |
Global attractor id → basin fraction at each value ( |
attractors |
Per value, the located attractors keyed by their global (matched) id —
the records. For the numbers, use :attr: |
diverged |
Diverged-or-untracked fraction at each value: the diverged share plus
the basin mass of any attractor dropped by
TYPE:
|
centers
property
¶
centers: ndarray
Where each tracked attractor sits at each value: (n_values, n_ids, dim).
The numbers behind :attr:attractors, which is a list of dicts of
records — three containers deep before a coordinate. Row k is the
parameter value :attr:values [k], column j the attractor
:attr:ids [j], and a value where that id is absent is nan, so
the array is rectangular and aligns with np.asarray(self) index for
index::
cont.centers[:, 0, :] # attractor ids[0]'s path through state space
An empty sweep gives shape (0, 0, 0).
to_frame
¶
to_frame() -> Any
Return a tidy frame: one row per (value, attractor), with the share.
The base :meth:~tsdynamics.analysis.results.AnalysisResult.to_frame
builds ONE row of the declared fields, and a continuation's answer — the
basin fractions — lives in a dict[int, ndarray], which a table cell
cannot hold, so it was dropped: measured, the frame came back
(1, 13) with the parameter values spread across columns and no
fractions at all. A frame that silently omits the payload is worse
than no frame.
| RETURNS | DESCRIPTION |
|---|---|
DataFrame
|
Columns |
Source code in src/tsdynamics/analysis/basins/continuation.py
tipping_points
¶
tipping_points(
result: ContinuationResult, *, threshold: float = 0.0
) -> CollectionResult
Read tipping events off a continuation.
A tipping event is a basin fraction crossing threshold between consecutive
parameter values: a disappear event (fraction falls to/through it — an
attractor and its basin annihilating) or an appear event (a new attractor
gaining a basin). With the default threshold=0 only true
appearance/annihilation is reported; raise it to flag basins shrinking past a
safety margin.
| PARAMETER | DESCRIPTION |
|---|---|
result
|
A continuation from :func:
TYPE:
|
threshold
|
Basin-fraction level whose crossing counts as a tipping event.
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
CollectionResult
|
Behaves as a |
Source code in src/tsdynamics/analysis/basins/continuation.py
527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 | |
Sampling¶
Sagitta-based tools for time-ordered samples: choose an output dt (the
largest stride whose mid-point bow off the chord — the sagitta — stays under a
geometric tolerance), or read the per-point sagitta as a local "how sharply it
bends" field.
estimate_dt_from_sagitta
¶
estimate_dt_from_sagitta(
y: ndarray,
dt0: float,
*,
epsilon: float,
percentile: float = 95.0,
coarsen_only: bool = True,
use_relative: bool | None = None,
min_points_per_segment: int = 3,
search_growth: float = 1.5,
) -> SagittaDt
Output dt that keeps the drawn bow under epsilon.
The sagitta-based output-step :math:\Delta t^\ast selector.
A derivative-free, scale-invariant, idempotent heuristic for the output
sampling step. For a range of candidate strides it measures the sagitta —
the perpendicular bow of the midpoint of each (i-span, i, i+span) triple
off its chord — at a robust percentile, then picks the largest stride whose
sagitta still satisfies the geometric tolerance epsilon.
A one-dimensional input is first delay-embedded so the geometric criterion
lives in a reconstructed state space: the delay is the first minimum of the
time-delayed mutual information (Fraser & Swinney 1986) and the dimension is
Kennel's false-nearest-neighbour estimate (Kennel, Brown & Abarbanel 1992),
both via the public :mod:tsdynamics.analysis.embedding estimators
(:func:~tsdynamics.analysis.embedding.optimal_delay and
:func:~tsdynamics.analysis.embedding.embedding_dimension with
method="fnn"). Multivariate input is used directly, per-feature
σ-normalised (ddof=1).
| PARAMETER | DESCRIPTION |
|---|---|
y
|
Time-ordered samples. A 1-D (or
TYPE:
|
dt0
|
Base sampling step :math:
TYPE:
|
epsilon
|
Geometric tolerance :math:
TYPE:
|
percentile
|
Robust percentile :math:
TYPE:
|
coarsen_only
|
If
TYPE:
|
use_relative
|
Whether to use the relative sagitta criterion
(
TYPE:
|
min_points_per_segment
|
Minimum number of triples
TYPE:
|
search_growth
|
Multiplicative growth factor (
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
SagittaDt
|
The chosen :math: |
| RAISES | DESCRIPTION |
|---|---|
InvalidParameterError
|
If |
References
A. M. Fraser and H. L. Swinney, "Independent coordinates for strange attractors from mutual information", Phys. Rev. A 33, 1134 (1986).
M. B. Kennel, R. Brown and H. D. I. Abarbanel, "Determining embedding dimension for phase-space reconstruction using a geometrical construction", Phys. Rev. A 45, 3403 (1992).
Source code in src/tsdynamics/analysis/sampling/sagitta.py
246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 | |
sagitta_profile
¶
sagitta_profile(
samples: ndarray,
*,
span: int = 1,
relative: bool = True,
normalize: bool = True,
) -> ndarray
Per-point sagitta along a sampled curve — its local bow off the chord.
For each interior point i the sagitta is the perpendicular distance of
samples[i] from the chord through its neighbours (i-span, i+span) (the
same geometry the sagitta-based Δt selector uses, here evaluated at every
point rather than a strided percentile). Large where the trajectory bends
sharply, ~0 on straight runs; the first/last span points are 0.
| PARAMETER | DESCRIPTION |
|---|---|
samples
|
Time-ordered points (a 1-D series is treated as a column).
TYPE:
|
span
|
Neighbour offset of the chord
TYPE:
|
relative
|
Divide each bow by its chord length, giving a dimensionless bend ratio
that is (largely) independent of the local speed/step. Default
TYPE:
|
normalize
|
σ-normalize each feature first (per-feature std,
TYPE:
|
| RETURNS | DESCRIPTION |
|---|---|
(ndarray, shape(n))
|
The per-point sagitta, aligned to |