import mmf_setup

mmf_setup.nbinit()
from gpe.imports import *

This cell adds /home/docs/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:

  • Choose "Trust Notebook" from the "File" menu.
  • Re-execute this cell.
  • Reload the notebook.

[I 20:01:33 numexpr.utils] NumExpr defaulting to 2 threads.
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[1], line 4
      1 import mmf_setup
      2 
      3 mmf_setup.nbinit()
----> 4 from gpe.imports import *

File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/imports.py:28
     26 from mmfutils.plot import imcontourf  # noqa: E402
     27 from gpe.minimize import MinimizeState  # noqa: E402
---> 28 from gpe.utils import evolve_to, evolve, evolves  # noqa: E402
     29 from gpe.plot_utils import MPLGrid  # noqa: E402
     30 from pytimeode.evolvers import EvolverSplit, EvolverABM  # noqa: E402

File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/utils.py:35
     32 from persist.objects import Archivable
     33 from persist.archive import Archive
---> 35 from pytimeode.evolvers import EvolverABM
     36 from pytimeode.mixins import ArrayStateMixin
     37 from pytimeode.interfaces import implementer, IStateForABMEvolvers

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/__init__.py:2
      1 from . import interfaces
----> 2 from . import mixins
      3 from . import evolvers
      5 __all__ = ["interfaces", "mixins", "evolvers"]

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/mixins.py:39
     30 from .interfaces import (
     31     IState,
     32     IStateApply,
   (...)     36     implementer,
     37 )
     38 from . import interfaces
---> 39 from .utils import expr
     41 __all__ = [
     42     "StateMixin",
     43     "StatesMixin",
   (...)     47     "ArraysStateWithBraketMixin",
     48 ]
     51 @implementer(IState)
     52 class StateMixin:

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/utils/expr.py:18
     15 import numpy as np
     17 try:
---> 18     import sympy
     19 except ImportError:
     20     sympy = None

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/__init__.py:77
     70 from .logic import (to_cnf, to_dnf, to_nnf, And, Or, Not, Xor, Nand, Nor,
     71         Implies, Equivalent, ITE, POSform, SOPform, simplify_logic, bool_map,
     72         true, false, satisfiable)
     74 from .assumptions import (AppliedPredicate, Predicate, AssumptionsContext,
     75         assuming, Q, ask, register_handler, remove_handler, refine)
---> 77 from .polys import (Poly, PurePoly, poly_from_expr, parallel_poly_from_expr,
     78         degree, total_degree, degree_list, LC, LM, LT, pdiv, prem, pquo,
     79         pexquo, div, rem, quo, exquo, half_gcdex, gcdex, invert,
     80         subresultants, resultant, discriminant, cofactors, gcd_list, gcd,
     81         lcm_list, lcm, terms_gcd, trunc, monic, content, primitive, compose,
     82         decompose, sturm, gff_list, gff, sqf_norm, sqf_part, sqf_list, sqf,
     83         factor_list, factor, intervals, refine_root, count_roots, all_roots,
     84         real_roots, nroots, ground_roots, nth_power_roots_poly, cancel,
     85         reduced, groebner, is_zero_dimensional, GroebnerBasis, poly,
     86         symmetrize, horner, interpolate, rational_interpolate, viete, together,
     87         BasePolynomialError, ExactQuotientFailed, PolynomialDivisionFailed,
     88         OperationNotSupported, HeuristicGCDFailed, HomomorphismFailed,
     89         IsomorphismFailed, ExtraneousFactors, EvaluationFailed,
     90         RefinementFailed, CoercionFailed, NotInvertible, NotReversible,
     91         NotAlgebraic, DomainError, PolynomialError, UnificationFailed,
     92         GeneratorsError, GeneratorsNeeded, ComputationFailed,
     93         UnivariatePolynomialError, MultivariatePolynomialError,
     94         PolificationFailed, OptionError, FlagError, minpoly,
     95         minimal_polynomial, primitive_element, field_isomorphism,
     96         to_number_field, isolate, round_two, prime_decomp, prime_valuation,
     97         galois_group, itermonomials, Monomial, lex, grlex,
     98         grevlex, ilex, igrlex, igrevlex, CRootOf, rootof, RootOf,
     99         ComplexRootOf, RootSum, roots, Domain, FiniteField, IntegerRing,
    100         RationalField, RealField, ComplexField, PythonFiniteField,
    101         GMPYFiniteField, PythonIntegerRing, GMPYIntegerRing, PythonRational,
    102         GMPYRationalField, AlgebraicField, PolynomialRing, FractionField,
    103         ExpressionDomain, FF_python, FF_gmpy, ZZ_python, ZZ_gmpy, QQ_python,
    104         QQ_gmpy, GF, FF, ZZ, QQ, ZZ_I, QQ_I, RR, CC, EX, EXRAW,
    105         construct_domain, swinnerton_dyer_poly, cyclotomic_poly,
    106         symmetric_poly, random_poly, interpolating_poly, jacobi_poly,
    107         chebyshevt_poly, chebyshevu_poly, hermite_poly, hermite_prob_poly,
    108         legendre_poly, laguerre_poly, apart, apart_list, assemble_partfrac_list,
    109         Options, ring, xring, vring, sring, field, xfield, vfield, sfield)
    111 from .series import (Order, O, limit, Limit, gruntz, series, approximants,
    112         residue, EmptySequence, SeqPer, SeqFormula, sequence, SeqAdd, SeqMul,
    113         fourier_series, fps, difference_delta, limit_seq)
    115 from .functions import (factorial, factorial2, rf, ff, binomial,
    116         RisingFactorial, FallingFactorial, subfactorial, carmichael,
    117         fibonacci, lucas, motzkin, tribonacci, harmonic, bernoulli, bell, euler,
   (...)    138         Znm, elliptic_k, elliptic_f, elliptic_e, elliptic_pi, beta, mathieus,
    139         mathieuc, mathieusprime, mathieucprime, riemann_xi, betainc, betainc_regularized)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/polys/__init__.py:79
      3 __all__ = [
      4     'Poly', 'PurePoly', 'poly_from_expr', 'parallel_poly_from_expr', 'degree',
      5     'total_degree', 'degree_list', 'LC', 'LM', 'LT', 'pdiv', 'prem', 'pquo',
   (...)     65     'field', 'xfield', 'vfield', 'sfield'
     66 ]
     68 from .polytools import (Poly, PurePoly, poly_from_expr,
     69         parallel_poly_from_expr, degree, total_degree, degree_list, LC, LM,
     70         LT, pdiv, prem, pquo, pexquo, div, rem, quo, exquo, half_gcdex, gcdex,
   (...)     76         nth_power_roots_poly, cancel, reduced, groebner, is_zero_dimensional,
     77         GroebnerBasis, poly)
---> 79 from .polyfuncs import (symmetrize, horner, interpolate,
     80         rational_interpolate, viete)
     82 from .rationaltools import together
     84 from .polyerrors import (BasePolynomialError, ExactQuotientFailed,
     85         PolynomialDivisionFailed, OperationNotSupported, HeuristicGCDFailed,
     86         HomomorphismFailed, IsomorphismFailed, ExtraneousFactors,
   (...)     91         MultivariatePolynomialError, PolificationFailed, OptionError,
     92         FlagError)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/polys/polyfuncs.py:10
      8 from sympy.polys.polyoptions import allowed_flags, build_options
      9 from sympy.polys.polytools import poly_from_expr, Poly
---> 10 from sympy.polys.specialpolys import (
     11     symmetric_poly, interpolating_poly)
     12 from sympy.polys.rings import sring
     13 from sympy.utilities import numbered_symbols, take, public

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/polys/specialpolys.py:298
    294     return dmp_mul(f, h, n, K), dmp_mul(g, h, n, K), h
    296 # A few useful polynomials from Wang's paper ('78).
--> 298 from sympy.polys.rings import ring
    300 def _f_0():
    301     R, x, y, z = ring("x,y,z", ZZ)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/polys/rings.py:31
     27 from sympy.polys.polyoptions import (Domain as DomainOpt,
     28                                      Order as OrderOpt, build_options)
     29 from sympy.polys.polyutils import (expr_from_dict, _dict_reorder,
     30                                    _parallel_dict_from_expr)
---> 31 from sympy.printing.defaults import DefaultPrinting
     32 from sympy.utilities import public, subsets
     33 from sympy.utilities.iterables import is_sequence

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/printing/__init__.py:11
      7 from .mathml import mathml, print_mathml
      9 from .python import python, print_python
---> 11 from .pycode import pycode
     13 from .codeprinter import print_ccode, print_fcode
     15 from .codeprinter import ccode, fcode, cxxcode, rust_code # noqa:F811

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/printing/pycode.py:11
      9 from sympy.core.mod import Mod
     10 from .precedence import precedence
---> 11 from .codeprinter import CodePrinter
     13 _kw = {
     14     'and', 'as', 'assert', 'break', 'class', 'continue', 'def', 'del', 'elif',
     15     'else', 'except', 'finally', 'for', 'from', 'global', 'if', 'import', 'in',
     16     'is', 'lambda', 'not', 'or', 'pass', 'raise', 'return', 'try', 'while',
     17     'with', 'yield', 'None', 'False', 'nonlocal', 'True'
     18 }
     20 _known_functions = {
     21     'Abs': 'abs',
     22     'Min': 'min',
     23     'Max': 'max',
     24 }

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/printing/codeprinter.py:13
     11 from sympy.core.sorting import default_sort_key
     12 from sympy.core.symbol import Symbol
---> 13 from sympy.functions.elementary.complexes import re
     14 from sympy.printing.str import StrPrinter
     15 from sympy.printing.precedence import precedence, PRECEDENCE

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/functions/__init__.py:28
     25 from sympy.functions.elementary.integers import floor, ceiling, frac
     26 from sympy.functions.elementary.piecewise import (Piecewise, piecewise_fold,
     27                                                   piecewise_exclusive)
---> 28 from sympy.functions.special.error_functions import (erf, erfc, erfi, erf2,
     29         erfinv, erfcinv, erf2inv, Ei, expint, E1, li, Li, Si, Ci, Shi, Chi,
     30         fresnels, fresnelc)
     31 from sympy.functions.special.gamma_functions import (gamma, lowergamma,
     32         uppergamma, polygamma, loggamma, digamma, trigamma, multigamma)
     33 from sympy.functions.special.zeta_functions import (dirichlet_eta, zeta,
     34         lerchphi, polylog, stieltjes, riemann_xi)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/functions/special/error_functions.py:2710
   2702         return super()._eval_aseries(n, args0, x, logx)
   2705 ###############################################################################
   2706 #################### HELPER FUNCTIONS #########################################
   2707 ###############################################################################
-> 2710 class _erfs(DefinedFunction):
   2711     """
   2712     Helper function to make the $\\mathrm{erf}(z)$ function
   2713     tractable for the Gruntz algorithm.
   2714 
   2715     """
   2716     @classmethod
   2717     def eval(cls, arg):

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/core/basic.py:218, in Basic.__init_subclass__(cls)
    212 def __init_subclass__(cls):
    213     # Initialize the default_assumptions FactKB and also any assumptions
    214     # property methods. This method will only be called for subclasses of
    215     # Basic but not for Basic itself so we call
    216     # _prepare_class_assumptions(Basic) below the class definition.
    217     super().__init_subclass__()
--> 218     _prepare_class_assumptions(cls)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/core/assumptions.py:642, in _prepare_class_assumptions(cls)
    640 cls._prop_handler = {}
    641 for k in _assume_defined:
--> 642     eval_is_meth = getattr(cls, '_eval_is_%s' % k, None)
    643     if eval_is_meth is not None:
    644         cls._prop_handler[k] = eval_is_meth

KeyboardInterrupt: 

Plotting#

Here we demonstrate some tools for plotting, animating results, and making movies.

MPLGrid#

Here as a simple example, we start with a harmonically trapped gas, phase imprint a dark soliton, and then watch the “evolution” as we cool into the soliton state.

from gpe.bec import State, u

s = State(Nxyz=(128,), Lxyz=(20.0,))
s[:] *= np.sign(s.xyz[0])
s.cooling_phase = 1j
e = EvolverABM(s, dt=0.1 * s.t_scale)


def plot(state):
    plt.plot(state.xyz[0] / u.micron, state.get_density())


states = [s]
with NoInterrupt(ignore=True) as interrupted:
    for n in range(6):
        e.evolve(200)
        states.append(e.get_y())
        plt.clf()
        plot(e.y)
        display(plt.gcf())
        clear_output(wait=True)

Our first demonstration is to stack these frames using MPLGrid. We can stack going down or going right.

direction='down'#

import gpe.plot_utils

reload(gpe.plot_utils)
from gpe.plot_utils import MPLGrid

plt.figure(figsize=(10, 5))
grid = MPLGrid(direction="down", share=True)
for s in states:
    ax = grid.next()
    ax.plot(s.xyz[0] / u.micron, s.get_density())
plt.tight_layout()

direction='right'#

The plots can also be stacked to the right:

plt.figure(figsize=(15, 2))
grid = MPLGrid(direction="right", share=True)
for s in states:
    ax = grid.next()
    ax.plot(s.xyz[0] / u.micron, s.get_density())
plt.tight_layout()

right=True#

plt.figure(figsize=(15, 2))
grid = MPLGrid(direction="right", share=True, right=True)
for s in states:
    ax = grid.next()
    ax.plot(s.xyz[0] / u.micron, s.get_density())
plt.tight_layout()

grid()#

More complex layouts can be created using subgrids:

plt.figure(figsize=(10, 3))
grid = MPLGrid(direction="down", share=True)
for s in states[::3]:
    ax = grid.next(size=1.0 / 3)  # Make these 1/3 of the size
    ax.plot(s.xyz[0] / u.micron, s.get_density())
subgrid = grid.grid(size=1, direction="right", share=True)
for s in states:
    ax = subgrid.next(2)
    ax.plot(s.xyz[0] / u.micron, s.get_density())

plt.tight_layout()

Note in this case that we obscure the x-axis of the upper plots. This is a limitation of the current MPLGrid class - all axes must either be shared or not shared. To add a space, the first three panes should be in their own subgrid:

reload(gpe.plot_utils)
from gpe.plot_utils import MPLGrid
from matplotlib.gridspec import GridSpecFromSubplotSpec

plt.figure(figsize=(10, 3))
grid = MPLGrid(direction="down", share=False, space=0.5)
subgrid = grid.grid(direction="down", share=True)
for s in states[::3]:
    ax = subgrid.next(1.0 / 2)  # Make these 1/2 of the size
    ax.plot(s.xyz[0] / u.micron, s.get_density())
ax.set_xlabel("x")

subgrid = grid.grid(
    direction="right", share=True, space=0.1
)  # Some space might help lables not clash
for s in states:
    ax = subgrid.next(2)
    ax.plot(s.xyz[0] / u.micron, s.get_density())

plt.tight_layout()

Note that tight_layout() does not work properly here. Not sure why. Probably the solution is to use the new constrained layout feature of matplotlib once it is settled.