API reference

pyEDW.model

The Exo-Daisy World model and its stochastic integrator.

This module ports updateExoDaisyWorld.m (the single time step) and eDW_BHsim.m / exoDaisyWorld.m (the ensemble drivers) to Python. The dynamics are the dimensionless system of Appendix B of Sowinski, Ghoshal & Frank, Planet. Sci. J. 6, 176 (2025):

• Eq. (B8) d f_B / dt = w(T_B - 1)(f - f_B - f_W) f_B - (gamma_D/gamma_G) f_B

• Eq. (B9) d f_W / dt = w(T_W - 1)(f - f_B - f_W) f_W - (gamma_D/gamma_G) f_W

• Eq. (B10) d T   / dt = (1/gamma_G tau_E)[(1 - (dA_B f_B + dA_W f_W)/(1-A_G)) L - T^4]

• Eq. (B11) d L   / dt = (1/gamma_G tau_S)(1 + lambda - L) + sqrt(2/gamma_G tau_S) delta (1+lambda) eta

with growth window w(x) = exp(-alpha x^4) (Eq. B7), scaled daisy temperatures T_alpha^4 = T^4 + Q dA_B f_B + Q dA_W f_W - Q dA_alpha (Eq. B12), and dA_alpha = A_alpha - A_G. Temperatures are in units of the optimal temperature T_opt and luminosity in units of L_opt.

The noise is a single scalar Wiener process entering only the luminosity, so the SDE is scalar-driven and the strong-order-1 stochastic Runge-Kutta scheme of A. Roberts (2012) applies. Parameters packs the physical constants into the dimensionless theta vector consumed by the integrator, whose layout mirrors the MATLAB comment block exactly:

idx	theta component
0	f
1	1 - A_G
2	A_B - A_G
3	A_W - A_G
4	Q
5	gamma_D / gamma_G
6	1 / (gamma_G tau_E)
7	1 / (gamma_G tau_S)
8	8 (T_opt / Delta T)^4
9	delta
10	lambda

class pyEDW.model.Parameters(f: float = 0.88, A_G: float = 0.3, A_B: float = 0.1, A_W: float = 0.6, T_opt: float = 300.0, Q: float = 0.1, gamma_G: float = 1.0, gamma_D: float = 0.2, tau_S: float = 3.0, tau_E: float = 5.0, dT: float = 30.0, delta: float = 0.05, lam: float = 0.0)

Bases: object

Physical parameters of Exo-Daisy World.

All temperatures are expressed in units of the optimal temperature T_opt and luminosities in units of L_opt once dimensionalized; the fields below are the dimensional constants of Table 1, from which to_theta() builds the dimensionless theta vector.

Parameters:

• f (float) – Habitable land fraction.

• A_G (float) – Ground, black-daisy, and white-daisy albedos (A_B < A_G < A_W).

• A_B (float) – Ground, black-daisy, and white-daisy albedos (A_B < A_G < A_W).

• A_W (float) – Ground, black-daisy, and white-daisy albedos (A_B < A_G < A_W).

• T_opt (float) – Optimal growth temperature (sets the temperature scale).

• Q (float) – Temperature-shift heat constant coupling daisy albedo to local temperature.

• gamma_G (float) – Maximal growth rate and daisy decay rate.

• gamma_D (float) – Maximal growth rate and daisy decay rate.

• tau_S (float) – Stellar-fluctuation and environmental (thermal) timescales, in units of 1 / gamma_G.

• tau_E (float) – Stellar-fluctuation and environmental (thermal) timescales, in units of 1 / gamma_G.

• dT (float) – Growth-rate bandwidth Delta T (range of bearable temperatures).

• delta (float) – Relative amplitude of the stellar-luminosity fluctuations.

• lam (float) – Luminosity offset lambda; the mean luminosity is <L> = (1 + lambda) L_opt.

f: float = 0.88

A_G: float = 0.3

A_B: float = 0.1

A_W: float = 0.6

T_opt: float = 300.0

Q: float = 0.1

gamma_G: float = 1.0

gamma_D: float = 0.2

tau_S: float = 3.0

tau_E: float = 5.0

dT: float = 30.0

delta: float = 0.05

lam: float = 0.0

to_theta()

Return the 11-element dimensionless theta vector.

Returns:

theta in the layout documented in the module docstring, matching eDW_BHsim.m.

Return type:

numpy.ndarray

with_bandwidth(dT)

Return a copy with a new growth-rate bandwidth Delta T.

with_lambda(lam)

Return a copy with a new luminosity offset lambda.

class pyEDW.model.ExoDaisyWorld(params, x0=None, dt=0.1, rng=None)

Bases: object

The coupled agent+environment stochastic Daisy World.

Wraps the compiled integrator behind a small stateful interface that keeps the MATLAB usage: build with Parameters, advance with evolve() (one step) or run() (many), and reach for the static ensemble() to sweep a luminosity axis in parallel.

Parameters:

• params (Parameters) – Physical parameters.

• x0 (array_like, optional) – Initial [f_B, f_W, T, L]. Defaults to a small equal daisy seed at the equilibrium temperature T = L**0.25 for L = 1 + lambda.

• dt (float, optional) – Integration time step (default 0.1, as in the MATLAB).

• rng (numpy.random.Generator, optional) – Used only to seed NumPy’s global state for the compiled loops (Numba’s random stream). Pass an integer-seeded generator for repeatability.

evolve()

Advance the state by one SRK1 step and return it.

run(nsteps, record=False)

Integrate nsteps steps from the current state.

Parameters:

• nsteps (int) – Number of SRK1 steps.

• record (bool, optional) – If True, return the full (nsteps + 1, 4) trajectory including the current state; otherwise advance in place and return the (4,) endpoint.

Returns:

Endpoint or trajectory, per record.

Return type:

numpy.ndarray

static ensemble(params, Ls=None, nsteps=2000, N=500, dt=0.1, agent_free=False, seed=None)

Parallel endpoint ensemble over a luminosity axis.

Ports eDW_BHsim.m (biotic) and the agent-free driver of exoDaisyWorld.m. For each luminosity, N independent instances are integrated for nsteps steps from freshly sampled initial conditions and their endpoints recorded.

Parameters:

• params (Parameters) – Physical parameters; params.lam is overridden per luminosity so that <L> = 1 + lambda tracks the axis.

• Ls (array_like, optional) – Mean luminosities. Defaults to 1 + linspace(-0.7, 1.4, 400).

• nsteps (int, optional) – Steps per instance (default 2000; the agent-free run used 200).

• N (int, optional) – Instances per luminosity (default 500).

• dt (float, optional) – Time step (default 0.1).

• agent_free (bool, optional) – If True, seed both daisy fractions at zero (environment only).

• seed (int, optional) – Seeds NumPy’s global RNG before the parallel loop.

Returns:

Endpoints of shape (len(Ls), N, 4).

Return type:

numpy.ndarray

pyEDW.metrics

Viability and the information architecture of Exo-Daisy World.

This module ports the entropy pipeline of Daisy_World_Paper_Figures.m (lines 307-392) and the measure assembly of makeFIG4.m / makeFig5.m. Given an endpoint ensemble data of shape (n_lum, N, 4) – the output of pyEDW.model.ExoDaisyWorld.ensemble() – it estimates, at each luminosity, the joint distribution p_AE(f_B, f_W, T, L) by histogramming and reads off the Shannon information measures of Section 4.

Binning follows the paper’s square-root rule: with n living instances the number of bins is N_bins = 1 + ceil(sqrt(n)) and each variable is split by N_bins - 1 equally spaced interior edges spanning its observed range, so a value lands in bin sum(v > edges) (0-indexed) – exactly the MATLAB 1 + sum(v > x) up to the index base.

Entropy table

entropy_table() returns, per luminosity, a 16-column array H whose columns are the entropies of every marginal of p_AE needed downstream (agent dof a_1 = f_B, a_2 = f_W; environment dof e_1 = T, e_2 = L):

col	entropy	col	entropy
0	H[a1]	8	H[a1 e1]
1	H[a2]	9	H[a1 a2]
2	H[e1]	10	H[a2 e1 e2]
3	H[e2]	11	H[a1 e1 e2]
4	H[e1 e2]	12	H[a1 a2 e2]
5	H[a2 e2]	13	H[a1 a2 e1]
6	H[a2 e1]	14	H[a1 a2 e1 e2]
7	H[a1 e2]	15	N_bins

The information measures are then

• I(A:E)  = H[e1 e2] + H[a1 a2] - H[a1 a2 e1 e2] (mutual_information())

• dI      = I(e1:e2) - I0(e1:e2) (delta_I())

• C(a1:a2||E) = H[a1] + H[a2] + H[e1 e2] + H[a1 a2 e1 e2] `` - H[a1 a2] - H[a1 e1 e2] - H[a2 e1 e2]`` (cooperation())

Reference: Sowinski, Ghoshal & Frank, Planet. Sci. J. 6, 176 (2025), Section 4 and Appendix B.3.

pyEDW.metrics.entropy_table(data)

Per-luminosity entropy table for a biotic endpoint ensemble.

Parameters:

data (numpy.ndarray) – Ensemble endpoints of shape (n_lum, N, 4) with the last axis [f_B, f_W, T, L] (from pyEDW.model.ExoDaisyWorld.ensemble()).

Returns:

Array H of shape (n_lum, 16); see the module docstring for the column layout. Luminosities with fewer than two living instances (f_B f_W > 1e-7) are left as NaN except the bin count.

Return type:

numpy.ndarray

pyEDW.metrics.entropy_table_env(data_env)

Environment-only entropy table for the agent-free ensemble.

Ports the data_entropyE block: with the biome removed, only (T, L) vary. The bin count is fixed by the full instance count N (not a living subset), matching Nbins = 1 + ceil(sqrt(N)) in the MATLAB.

Parameters:

data_env (numpy.ndarray) – Agent-free endpoints of shape (n_lum, N, 4).

Returns:

Array of shape (n_lum, 4) with columns [H[e1], H[e2], H[e1 e2], N_bins].

Return type:

numpy.ndarray

pyEDW.metrics.viability(data, f)

Viability V = E^A[(f_B + f_W)/f] per luminosity (Eq. 9).

Parameters:

• data (numpy.ndarray) – Ensemble endpoints (n_lum, N, 4).

• f (float) – Habitable land fraction.

Returns:

V of shape (n_lum,), the ensemble-mean occupied fraction of the habitable area (0 <= V <= 1).

Return type:

numpy.ndarray

pyEDW.metrics.efficacy(data)

Thermoregulatory efficacy E = <T> - 1 per luminosity.

The mean planetary temperature’s deviation from the optimum (in units of T_opt); zero means the biome holds the surface exactly at T_opt.

pyEDW.metrics.mutual_information(H)

Agent-environment mutual information I(A:E) (bits).

I(A:E) = H[e1 e2] + H[a1 a2] - H[a1 a2 e1 e2].

pyEDW.metrics.intra_environment_information(H)

Intra-environment correlation I(e1:e2) = H[e1] + H[e2] - H[e1 e2].

pyEDW.metrics.delta_I(H, H_env)

Agent-induced change in intra-environment correlation, dI.

dI = I(e1:e2) - I0(e1:e2) (Eq. 12): the intra-environment mutual information with the biome present minus its agent-free value. Positive means the biome strengthens the temperature-luminosity correlation.

Parameters:

• H (numpy.ndarray) – Biotic entropy table, shape (..., 16) from entropy_table().

• H_env (numpy.ndarray) – Agent-free table, shape (n_lum, 4) from entropy_table_env(), broadcast against H along the luminosity axis.

Returns:

dI in bits.

Return type:

numpy.ndarray

pyEDW.metrics.cooperation(H)

Intra-agent cooperation given the environment, C(a1:a2||E) (Eq. 14).

The interaction information

C = H[a1] + H[a2] + H[e1 e2] + H[a1 a2 e1 e2] `` - H[a1 a2] - H[a1 e1 e2] - H[a2 e1 e2]``,

equal to I(a1:a2) - I(a1:a2 | E). Negative values indicate a synergistic effect where the environment enhances the correlation between the two daisy species.

pyEDW.figures

Plotting helpers that regenerate the figures of the eDW paper.

These port the visual content of makeFIG1.m – makeFig5.m. Each function takes an endpoint ensemble (and, where needed, a precomputed entropy table) and returns a Matplotlib (fig, axes) pair, so the caller controls saving and styling. Colors follow the paper: white daisies in pink, black daisies in blue, temperature in green.

• ensemble_scatter() – Fig. 1: daisy fractions, albedo, and temperature scattered over luminosity for one growth-rate bandwidth.

• corner() – Fig. 2: pairwise joint distribution of the four dof at a single luminosity.

• cooperation_vs_deltaI() – Fig. 4: cooperation against the agent-induced correlation change, split by viability.

• information_vs_viability() – Fig. 5: I(A:E) against viability.

pyEDW.figures.albedo(data, params)

Planetary albedo per instance, A = A_G + dA_B f_B + dA_W f_W (Eq. 2).

pyEDW.figures.ensemble_scatter(data, params, Ls=None, s=3, alpha=0.05)

Fig. 1: one bandwidth column – fractions, albedo, temperature vs L.

Parameters:

• data (numpy.ndarray) – Endpoint ensemble (n_lum, N, 4).

• params (pyEDW.model.Parameters) – Parameters used to generate data (for albedo and T_opt lines).

• Ls (array_like, optional) – Luminosity axis; defaults to 1 + linspace(-0.7, 1.4, n_lum).

• s (float) – Scatter marker size and opacity.

• alpha (float) – Scatter marker size and opacity.

Returns:

Figure and its three stacked axes (fractions, albedo, temperature).

Return type:

(matplotlib.figure.Figure, numpy.ndarray)

pyEDW.figures.corner(data, ll, labels=('$f_B$', '$f_W$', '$T$', '$L$'), bins=40)

Fig. 2: corner plot of the four dof at luminosity index ll.

Diagonal panels show the marginal histograms; lower-triangle panels show the pairwise joint densities, revealing the intra-agent (f_B vs f_W) and intra-environment (T vs L) correlations.

Returns:

Figure and its 4 x 4 axis grid.

Return type:

(matplotlib.figure.Figure, numpy.ndarray)

pyEDW.figures.cooperation_vs_deltaI(H, viability, v_hi=0.75, v_lo=0.3, bins=60)

Fig. 4: cooperation C(a1:a2||E) vs correlation change dI.

Points are pooled over luminosity (and, if 2-D, bandwidth) and split into high- and low-viability clouds.

Parameters:

• H (numpy.ndarray) – Entropy table of shape (..., 16) from pyEDW.metrics().

• viability (numpy.ndarray) – Viability broadcastable to H[..., 0].

• v_hi (float) – Viability thresholds for the two clouds.

• v_lo (float) – Viability thresholds for the two clouds.

Return type:

(matplotlib.figure.Figure, matplotlib.axes.Axes)

pyEDW.figures.information_vs_viability(H, viability, bins=120)

Fig. 5: agent-environment information I(A:E) vs viability V.