Skip to content

Analysis

This module defines abstract base classes for commensurability analysis of a galactic potential's phase space. The "dimensionality" associated with the class corresponds with the dimensionality of the relevant orbits.

This module also defines user-facing analysis classes for 2D and 3D commensurability analysis using the tessellation subpackage.

AnalysisBase

Bases: Mapping

Base class for analyzing commensurate orbits within galactic potentials.

This class provides methods for evaluating orbits, constructing images of phase space slices, and saving/loading analysis data.

Attributes:

Name Type Description
ic_function Callable[..., SkyCoord]

Function to generate initial conditions for orbits.
ic_values Mapping[str, Sequence[float]]

Values to pass into ic_function.
axis_names List[str]

Names of the parameters for the initial condition function.
shape Tuple[int]

Shape of the measures array.
size int

Size of the measures array.
potential Any

Potential generated by potential_function.
dt Quantity

Time step for orbit integration.
steps int

Number of integration steps.
pattern_speed Quantity

Pattern speed for orbit integration.
measures ndarray

Array to store orbit measures.
Source code in commensurability/analysis.py 
 41
 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
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
class AnalysisBase(MappingABC):
    """
    Base class for analyzing commensurate orbits within galactic potentials.

    This class provides methods for evaluating orbits, constructing images of
    phase space slices, and saving/loading analysis data.

    Attributes:
        ic_function (Callable[..., c.SkyCoord]): Function to generate initial conditions for orbits.
        ic_values (Mapping[str, Sequence[float]]): Values to pass into ic_function.
        axis_names (List[str]): Names of the parameters for the initial condition function.
        shape (Tuple[int]): Shape of the measures array.
        size (int): Size of the measures array.
        potential (Any): Potential generated by potential_function.
        dt (u.Quantity): Time step for orbit integration.
        steps (int): Number of integration steps.
        pattern_speed (u.Quantity): Pattern speed for orbit integration.
        measures (np.ndarray): Array to store orbit measures.
    """

    @staticmethod
    @abstractmethod
    def evaluate(orbit: c.SkyCoord) -> Evaluation:
        """
        Evaluate an orbit and return an evaluation object.

        Args:
            orbit (c.SkyCoord): Integrated orbit coordinates.

        Returns:
            Evaluation: Evaluation object containing analysis results.
        """
        pass

    def __init__(
        self,
        ic_function: Callable[..., c.SkyCoord],
        values: Mapping[str, Sequence[float]],
        /,
        potential_function: Callable[[], Any],
        dt: Union[float, np.ndarray, u.Quantity],
        steps: Union[int, np.ndarray],
        *,
        pattern_speed: Union[float, u.Quantity] = 0.0,
        backend: Optional[Union[str, Backend]] = None,
        progressbar: bool = True,
        pidgey_chunksize: Optional[int] = None,
        _blank_measures: bool = False,
    ) -> None:
        """
        Initialize AnalysisBase instance.

        Args:
            ic_function (Callable[..., c.SkyCoord]): Function to generate initial conditions for orbits.
            values (Mapping[str, Sequence[float]]): Values to pass into ic_function.
            potential_function (Callable[[], Any]): Function to generate the potential for orbit integration.
            dt (Union[float, np.ndarray, u.Quantity]): Time step for orbit integration.
            steps (Union[int, np.ndarray]): Number of integration steps.
            pattern_speed (Union[float, u.Quantity], optional): Pattern speed for orbit integration (default 0.0).
            backend (Optional[Union[str, Backend]], optional): Backend for orbit computation.
            progressbar (bool, optional): Whether to show progress bar during image construction (default True).
            pidgey_chunksize (int, optional): Chunk size for orbit integration (default 1).
        """
        self.ic_function = ic_function
        self.ic_values = values
        argspec = inspect.getfullargspec(ic_function)
        self.axis_names = argspec.args
        if len(self.ic_values) != len(self.axis_names) or not all(
            name in self.axis_names for name in self.ic_values
        ):
            raise ValueError("values and ic_function signature do not correspond")
        self.shape = tuple(len(self.ic_values[name]) for name in self.axis_names)
        self.size = prod(self.shape)

        self.potential_function = potential_function
        self.potential = potential_function()
        self.dt = make_quantity(dt, u.Gyr)
        if self.dt.value <= 0:
            raise ValueError("dt must be greater than 0")
        self.steps = steps
        if self.steps <= 0:
            raise ValueError("steps must be greater than 0")

        self.pattern_speed = make_quantity(pattern_speed, u.km / u.s / u.kpc)
        if isinstance(backend, str):
            if backend == "agama":
                backend = AgamaBackend()
            elif backend == "gala":
                backend = GalaBackend()
            elif backend == "galpy":
                backend = GalpyBackend()
            else:
                raise ValueError(f"Unrecognized backend: {backend}")
        self.backend = backend or get_backend_from(self.potential)
        if not self.backend:
            raise TypeError(f"Unrecognized potential: {self.potential}")

        if pidgey_chunksize is None:
            pidgey_chunksize = min(int(self.size**0.5), DEFAULT_MAX_CHUNKSIZE)
        if pidgey_chunksize <= 0:
            raise ValueError("pidgey_chunksize must be greater than 0")
        if pidgey_chunksize >= self.size:
            pidgey_chunksize = self.size
            # raise ValueError("chunksize must be less than total number of starting coordinates")

        self.measures = np.zeros(self.shape)
        if not _blank_measures:
            self._construct_image(pidgey_chunksize, progressbar)

    def __len__(self) -> int:
        """
        Returns the size of the measures array.

        Returns:
            int: The size of the measures array.
        """
        return self.size

    def __iter__(self) -> Generator[tuple[int, ...], None, None]:
        """
        Iterate over the measures array's indices.

        Yields:
            Tuple[float]: Tuple of the indices for each grid point.
        """
        for pixel in np.ndindex(self.shape):
            yield pixel

    def __getitem__(self, key: tuple[int, ...]) -> tuple[c.SkyCoord, float]:
        """
        Get the measure at a specific grid point.

        Args:
            key: Tuple of indices for the grid point.

        Returns:
            Tuple[SkyCoord, float]: Tuple of the initial condition and measure at the grid point.
        """
        if not isinstance(key, tuple):
            raise KeyError("key must be a tuple")
        args = [self.ic_values[ax][i] for i, ax in zip(key, self.axis_names)]
        return self.ic_function(*args), self.measures[key]

    def _construct_image(self, pidgey_chunksize: int = 1, progressbar: bool = True):
        """
        Construct an image of a slice of phase space by integrating orbits and evaluating them.

        Args:
            pidgey_chunksize (int, optional): Chunk size for batching orbit integration (default 1).
            progressbar (bool, optional): Whether to show progress bar during construction (default True).
        """
        for pixels in tqdm(
            chunked(np.ndindex(self.shape), pidgey_chunksize),
            desc=f"with {pidgey_chunksize=}",
            total=self.size // pidgey_chunksize,
            disable=not progressbar,
        ):
            coords = []
            for pixel in pixels:
                params = [self.ic_values[ax][i] for i, ax in zip(pixel, self.axis_names)]
                coord = self.ic_function(*params)
                coords.append(coord)
            coords = collapse_coords(coords)

            orbits = self.backend.compute_orbit(
                coords,
                self.potential,
                self.dt,
                self.steps,
                pattern_speed=self.pattern_speed,
            )
            for pixel, orbit in tqdm(
                zip(pixels, orbits),
                desc="commensurability evaluation",
                total=pidgey_chunksize,
                disable=not progressbar,
                leave=False,
            ):
                self.measures[pixel] = self.evaluate(orbit).measure

    def save(self, path: Any):
        """
        Save the analysis data to an HDF5 file.

        Args:
            path: Path to the HDF5 file.
        """
        path = Path(path)
        if not path.parent.exists():
            print("Parent directory does not exist; creating directory.")
            path.parent.mkdir(parents=True)

        # store image mapping function source
        icsource = inspect.getsource(self.ic_function)
        icsource = textwrap.dedent(icsource)
        icsource = icsource.replace(self.ic_function.__name__, "ic_function", 1)

        # store potential function source
        potsource = inspect.getsource(self.potential_function)
        potsource = textwrap.dedent(potsource)
        potsource = potsource.replace(self.potential_function.__name__, "potential_function", 1)

        attrs = dict(
            icfunc=np.void(icsource.encode("utf8")),
            potfunc=np.void(potsource.encode("utf8")),
            dt=self.dt,
            steps=self.steps,
            pattern_speed=self.pattern_speed,
            backend=np.void(self.backend.__class__.__name__.encode("utf8")),
        )
        with h5py.File(path, "w") as f:
            dset = f.create_dataset(self.__class__.__name__, data=self.measures)
            for attr, value in attrs.items():
                dset.attrs[attr] = value
            for attr, value in self.ic_values.items():
                dset.attrs[attr] = value

    @classmethod
    def read_from_hdf5(cls, path: Any, backend_cls: Optional[Backend] = None) -> AnalysisBase:
        """
        Read analysis data from an HDF5 file.

        Args:
            path: Path to the HDF5 file.

        Returns:
            AnalysisBase: Instance of AnalysisBase class with loaded data.
        """
        with h5py.File(path, "r") as f:
            dset = f[cls.__name__]

            if "icfunc" in dset.attrs:
                icsource = dset.attrs["icfunc"].tobytes().decode("utf8")
                namespace: dict[str, Any] = {}
                exec(icsource, {"u": u, "c": c}, namespace)
                ic_function = namespace["ic_function"]
            else:
                warnings.warn("No potential function defined.")

                def ic_function() -> None:
                    pass

            axis_names = inspect.getfullargspec(ic_function).args
            values = {name: dset.attrs[name] for name in axis_names}

            if "potfunc" in dset.attrs:
                potsource = dset.attrs["potfunc"].tobytes().decode("utf8")
                namespace = {}
                print(potsource)
                exec(potsource, {"u": u, "c": c}, namespace)
                potential_function = namespace["potential_function"]
            else:
                warnings.warn("No potential function defined.")

                def potential_function() -> None:
                    pass

            backend_cls = backend_cls or getattr(
                pidgey, dset.attrs["backend"].tobytes().decode("utf8")
            )
            analysis = cls(
                ic_function,
                values,
                potential_function=potential_function,
                dt=dset.attrs["dt"],
                steps=dset.attrs["steps"],
                pattern_speed=dset.attrs["pattern_speed"],
                backend=backend_cls(),
                _blank_measures=True,
            )
            analysis.measures = dset[()]
        return analysis

__getitem__(key)

Get the measure at a specific grid point.

Parameters:

Name Type Description Default
key tuple[int, ...]

Tuple of indices for the grid point.

 required

Returns:

Type Description
tuple[SkyCoord, float]

Tuple[SkyCoord, float]: Tuple of the initial condition and measure at the grid point.
Source code in commensurability/analysis.py 
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def __getitem__(self, key: tuple[int, ...]) -> tuple[c.SkyCoord, float]:
    """
    Get the measure at a specific grid point.

    Args:
        key: Tuple of indices for the grid point.

    Returns:
        Tuple[SkyCoord, float]: Tuple of the initial condition and measure at the grid point.
    """
    if not isinstance(key, tuple):
        raise KeyError("key must be a tuple")
    args = [self.ic_values[ax][i] for i, ax in zip(key, self.axis_names)]
    return self.ic_function(*args), self.measures[key]

__init__(ic_function, values, /, potential_function, dt, steps, *, pattern_speed=0.0, backend=None, progressbar=True, pidgey_chunksize=None, _blank_measures=False)

Initialize AnalysisBase instance.

Parameters:

Name Type Description Default
ic_function Callable[..., SkyCoord]

Function to generate initial conditions for orbits.

 required
values Mapping[str, Sequence[float]]

Values to pass into ic_function.

 required
potential_function Callable[[], Any]

Function to generate the potential for orbit integration.

 required
dt Union[float, ndarray, Quantity]

Time step for orbit integration.

 required
steps Union[int, ndarray]

Number of integration steps.

 required
pattern_speed Union[float, Quantity]

Pattern speed for orbit integration (default 0.0).

 0.0
backend Optional[Union[str, Backend]]

Backend for orbit computation.

 None
progressbar bool

Whether to show progress bar during image construction (default True).

 True
pidgey_chunksize int

Chunk size for orbit integration (default 1).

 None
Source code in commensurability/analysis.py 
 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
def __init__(
    self,
    ic_function: Callable[..., c.SkyCoord],
    values: Mapping[str, Sequence[float]],
    /,
    potential_function: Callable[[], Any],
    dt: Union[float, np.ndarray, u.Quantity],
    steps: Union[int, np.ndarray],
    *,
    pattern_speed: Union[float, u.Quantity] = 0.0,
    backend: Optional[Union[str, Backend]] = None,
    progressbar: bool = True,
    pidgey_chunksize: Optional[int] = None,
    _blank_measures: bool = False,
) -> None:
    """
    Initialize AnalysisBase instance.

    Args:
        ic_function (Callable[..., c.SkyCoord]): Function to generate initial conditions for orbits.
        values (Mapping[str, Sequence[float]]): Values to pass into ic_function.
        potential_function (Callable[[], Any]): Function to generate the potential for orbit integration.
        dt (Union[float, np.ndarray, u.Quantity]): Time step for orbit integration.
        steps (Union[int, np.ndarray]): Number of integration steps.
        pattern_speed (Union[float, u.Quantity], optional): Pattern speed for orbit integration (default 0.0).
        backend (Optional[Union[str, Backend]], optional): Backend for orbit computation.
        progressbar (bool, optional): Whether to show progress bar during image construction (default True).
        pidgey_chunksize (int, optional): Chunk size for orbit integration (default 1).
    """
    self.ic_function = ic_function
    self.ic_values = values
    argspec = inspect.getfullargspec(ic_function)
    self.axis_names = argspec.args
    if len(self.ic_values) != len(self.axis_names) or not all(
        name in self.axis_names for name in self.ic_values
    ):
        raise ValueError("values and ic_function signature do not correspond")
    self.shape = tuple(len(self.ic_values[name]) for name in self.axis_names)
    self.size = prod(self.shape)

    self.potential_function = potential_function
    self.potential = potential_function()
    self.dt = make_quantity(dt, u.Gyr)
    if self.dt.value <= 0:
        raise ValueError("dt must be greater than 0")
    self.steps = steps
    if self.steps <= 0:
        raise ValueError("steps must be greater than 0")

    self.pattern_speed = make_quantity(pattern_speed, u.km / u.s / u.kpc)
    if isinstance(backend, str):
        if backend == "agama":
            backend = AgamaBackend()
        elif backend == "gala":
            backend = GalaBackend()
        elif backend == "galpy":
            backend = GalpyBackend()
        else:
            raise ValueError(f"Unrecognized backend: {backend}")
    self.backend = backend or get_backend_from(self.potential)
    if not self.backend:
        raise TypeError(f"Unrecognized potential: {self.potential}")

    if pidgey_chunksize is None:
        pidgey_chunksize = min(int(self.size**0.5), DEFAULT_MAX_CHUNKSIZE)
    if pidgey_chunksize <= 0:
        raise ValueError("pidgey_chunksize must be greater than 0")
    if pidgey_chunksize >= self.size:
        pidgey_chunksize = self.size
        # raise ValueError("chunksize must be less than total number of starting coordinates")

    self.measures = np.zeros(self.shape)
    if not _blank_measures:
        self._construct_image(pidgey_chunksize, progressbar)

__iter__()

Iterate over the measures array's indices.

Yields:

Type Description
tuple[int, ...]

Tuple[float]: Tuple of the indices for each grid point.
Source code in commensurability/analysis.py 
159
160
161
162
163
164
165
166
167
def __iter__(self) -> Generator[tuple[int, ...], None, None]:
    """
    Iterate over the measures array's indices.

    Yields:
        Tuple[float]: Tuple of the indices for each grid point.
    """
    for pixel in np.ndindex(self.shape):
        yield pixel

__len__()

Returns the size of the measures array.

Returns:

Name Type Description
int int

The size of the measures array.
Source code in commensurability/analysis.py 
150
151
152
153
154
155
156
157
def __len__(self) -> int:
    """
    Returns the size of the measures array.

    Returns:
        int: The size of the measures array.
    """
    return self.size

evaluate(orbit) abstractmethod staticmethod

Evaluate an orbit and return an evaluation object.

Parameters:

Name Type Description Default
orbit SkyCoord

Integrated orbit coordinates.

 required

Returns:

Name Type Description
Evaluation Evaluation

Evaluation object containing analysis results.
Source code in commensurability/analysis.py 
61
62
63
64
65
66
67
68
69
70
71
72
73
@staticmethod
@abstractmethod
def evaluate(orbit: c.SkyCoord) -> Evaluation:
    """
    Evaluate an orbit and return an evaluation object.

    Args:
        orbit (c.SkyCoord): Integrated orbit coordinates.

    Returns:
        Evaluation: Evaluation object containing analysis results.
    """
    pass

read_from_hdf5(path, backend_cls=None) classmethod

Read analysis data from an HDF5 file.

Parameters:

Name Type Description Default
path Any

Path to the HDF5 file.

 required

Returns:

Name Type Description
AnalysisBase AnalysisBase

Instance of AnalysisBase class with loaded data.
Source code in commensurability/analysis.py 
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
@classmethod
def read_from_hdf5(cls, path: Any, backend_cls: Optional[Backend] = None) -> AnalysisBase:
    """
    Read analysis data from an HDF5 file.

    Args:
        path: Path to the HDF5 file.

    Returns:
        AnalysisBase: Instance of AnalysisBase class with loaded data.
    """
    with h5py.File(path, "r") as f:
        dset = f[cls.__name__]

        if "icfunc" in dset.attrs:
            icsource = dset.attrs["icfunc"].tobytes().decode("utf8")
            namespace: dict[str, Any] = {}
            exec(icsource, {"u": u, "c": c}, namespace)
            ic_function = namespace["ic_function"]
        else:
            warnings.warn("No potential function defined.")

            def ic_function() -> None:
                pass

        axis_names = inspect.getfullargspec(ic_function).args
        values = {name: dset.attrs[name] for name in axis_names}

        if "potfunc" in dset.attrs:
            potsource = dset.attrs["potfunc"].tobytes().decode("utf8")
            namespace = {}
            print(potsource)
            exec(potsource, {"u": u, "c": c}, namespace)
            potential_function = namespace["potential_function"]
        else:
            warnings.warn("No potential function defined.")

            def potential_function() -> None:
                pass

        backend_cls = backend_cls or getattr(
            pidgey, dset.attrs["backend"].tobytes().decode("utf8")
        )
        analysis = cls(
            ic_function,
            values,
            potential_function=potential_function,
            dt=dset.attrs["dt"],
            steps=dset.attrs["steps"],
            pattern_speed=dset.attrs["pattern_speed"],
            backend=backend_cls(),
            _blank_measures=True,
        )
        analysis.measures = dset[()]
    return analysis

save(path)

Save the analysis data to an HDF5 file.

Parameters:

Name Type Description Default
path Any

Path to the HDF5 file.

 required
Source code in commensurability/analysis.py 
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
def save(self, path: Any):
    """
    Save the analysis data to an HDF5 file.

    Args:
        path: Path to the HDF5 file.
    """
    path = Path(path)
    if not path.parent.exists():
        print("Parent directory does not exist; creating directory.")
        path.parent.mkdir(parents=True)

    # store image mapping function source
    icsource = inspect.getsource(self.ic_function)
    icsource = textwrap.dedent(icsource)
    icsource = icsource.replace(self.ic_function.__name__, "ic_function", 1)

    # store potential function source
    potsource = inspect.getsource(self.potential_function)
    potsource = textwrap.dedent(potsource)
    potsource = potsource.replace(self.potential_function.__name__, "potential_function", 1)

    attrs = dict(
        icfunc=np.void(icsource.encode("utf8")),
        potfunc=np.void(potsource.encode("utf8")),
        dt=self.dt,
        steps=self.steps,
        pattern_speed=self.pattern_speed,
        backend=np.void(self.backend.__class__.__name__.encode("utf8")),
    )
    with h5py.File(path, "w") as f:
        dset = f.create_dataset(self.__class__.__name__, data=self.measures)
        for attr, value in attrs.items():
            dset.attrs[attr] = value
        for attr, value in self.ic_values.items():
            dset.attrs[attr] = value

AnalysisBase2D

Bases: MPAnalysisBase

Base class for commensurability analysis on 2D orbits.

This class extends MPAnalysisBase and provides additional methods for launching interactive plots.

Attributes:

Name Type Description
ic_function Callable[..., SkyCoord]

Function to generate initial conditions for orbits.
ic_values Mapping[str, Sequence[float]]

Values to pass into ic_function.
axis_names List[str]

Names of the parameters for the initial condition function.
shape Tuple[int]

Shape of the measures array.
size int

Size of the measures array.
potential Any

Potential generated by potential_function.
dt Quantity

Time step for orbit integration.
steps int

Number of integration steps.
pattern_speed Quantity

Pattern speed for orbit integration.
measures ndarray

Array to store orbit measures.
Source code in commensurability/analysis.py 
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
class AnalysisBase2D(MPAnalysisBase):
    """
    Base class for commensurability analysis on 2D orbits.

    This class extends MPAnalysisBase and provides additional methods for launching interactive plots.

    Attributes:
        ic_function (Callable[..., c.SkyCoord]): Function to generate initial conditions for orbits.
        ic_values (Mapping[str, Sequence[float]]): Values to pass into ic_function.
        axis_names (List[str]): Names of the parameters for the initial condition function.
        shape (Tuple[int]): Shape of the measures array.
        size (int): Size of the measures array.
        potential (Any): Potential generated by potential_function.
        dt (u.Quantity): Time step for orbit integration.
        steps (int): Number of integration steps.
        pattern_speed (u.Quantity): Pattern speed for orbit integration.
        measures (np.ndarray): Array to store orbit measures.
    """

    def launch_interactive_plot(
        self,
        x_axis: Optional[str] = None,
        y_axis: Optional[str] = None,
        var_axis: Optional[str] = None,
    ):
        """
        Launch an interactive plot for 2D orbits.

        Args:
            x_axis (str): Name of the x-axis parameter.
            y_axis (str): Name of the y-axis parameter.
            var_axis (Optional[str]): Name of the axis varied by scrolling (optional).
        """
        iplot: Viewer = AnalysisViewer2D(self, x_axis, y_axis, var_axis)
        iplot.show()

launch_interactive_plot(x_axis=None, y_axis=None, var_axis=None)

Launch an interactive plot for 2D orbits.

Parameters:

Name Type Description Default
x_axis str

Name of the x-axis parameter.

 None
y_axis str

Name of the y-axis parameter.

 None
var_axis Optional[str]

Name of the axis varied by scrolling (optional).

 None
Source code in commensurability/analysis.py 
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
def launch_interactive_plot(
    self,
    x_axis: Optional[str] = None,
    y_axis: Optional[str] = None,
    var_axis: Optional[str] = None,
):
    """
    Launch an interactive plot for 2D orbits.

    Args:
        x_axis (str): Name of the x-axis parameter.
        y_axis (str): Name of the y-axis parameter.
        var_axis (Optional[str]): Name of the axis varied by scrolling (optional).
    """
    iplot: Viewer = AnalysisViewer2D(self, x_axis, y_axis, var_axis)
    iplot.show()

AnalysisBase3D

Bases: MPAnalysisBase

Base class for commensurability analysis on 3D orbits.

This class extends MPAnalysisBase and provides additional methods for launching interactive plots.

Attributes:

Name Type Description
ic_function Callable[..., SkyCoord]

Function to generate initial conditions for orbits.
ic_values Mapping[str, Sequence[float]]

Values to pass into ic_function.
axis_names List[str]

Names of the parameters for the initial condition function.
shape Tuple[int]

Shape of the measures array.
size int

Size of the measures array.
potential Any

Potential generated by potential_function.
dt Quantity

Time step for orbit integration.
steps int

Number of integration steps.
pattern_speed Quantity

Pattern speed for orbit integration.
measures ndarray

Array to store orbit measures.
Source code in commensurability/analysis.py 
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
class AnalysisBase3D(MPAnalysisBase):
    """
    Base class for commensurability analysis on 3D orbits.

    This class extends MPAnalysisBase and provides additional methods for launching interactive plots.

    Attributes:
        ic_function (Callable[..., c.SkyCoord]): Function to generate initial conditions for orbits.
        ic_values (Mapping[str, Sequence[float]]): Values to pass into ic_function.
        axis_names (List[str]): Names of the parameters for the initial condition function.
        shape (Tuple[int]): Shape of the measures array.
        size (int): Size of the measures array.
        potential (Any): Potential generated by potential_function.
        dt (u.Quantity): Time step for orbit integration.
        steps (int): Number of integration steps.
        pattern_speed (u.Quantity): Pattern speed for orbit integration.
        measures (np.ndarray): Array to store orbit measures.
    """

    def launch_interactive_plot(
        self,
        x_axis: Optional[str] = None,
        y_axis: Optional[str] = None,
        var_axis: Optional[str] = None,
    ):
        """
        Launch an interactive plot for 3D orbits.

        Args:
            x_axis (str): Name of the x-axis parameter.
            y_axis (str): Name of the y-axis parameter.
            var_axis (Optional[str]): Name of the variable axis (optional).
        """
        iplot: Viewer = AnalysisViewer3D(self, x_axis, y_axis, var_axis)
        iplot.show()

launch_interactive_plot(x_axis=None, y_axis=None, var_axis=None)

Launch an interactive plot for 3D orbits.

Parameters:

Name Type Description Default
x_axis str

Name of the x-axis parameter.

 None
y_axis str

Name of the y-axis parameter.

 None
var_axis Optional[str]

Name of the variable axis (optional).

 None
Source code in commensurability/analysis.py 
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
def launch_interactive_plot(
    self,
    x_axis: Optional[str] = None,
    y_axis: Optional[str] = None,
    var_axis: Optional[str] = None,
):
    """
    Launch an interactive plot for 3D orbits.

    Args:
        x_axis (str): Name of the x-axis parameter.
        y_axis (str): Name of the y-axis parameter.
        var_axis (Optional[str]): Name of the variable axis (optional).
    """
    iplot: Viewer = AnalysisViewer3D(self, x_axis, y_axis, var_axis)
    iplot.show()

MPAnalysisBase

Bases: AnalysisBase

Base class for analyzing commensurate orbits within galactic potentials. This makes use of the multiprocessing module for parallelization of orbit evaluation.

This class provides methods for evaluating orbits, constructing images of phase space slices, and saving/loading analysis data.

Attributes:

Name Type Description
ic_function Callable[..., SkyCoord]

Function to generate initial conditions for orbits.
ic_values Mapping[str, Sequence[float]]

Values to pass into ic_function.
axis_names List[str]

Names of the parameters for the initial condition function.
shape Tuple[int]

Shape of the measures array.
size int

Size of the measures array.
potential Any

Potential generated by potential_function.
dt Quantity

Time step for orbit integration.
steps int

Number of integration steps.
pattern_speed Quantity

Pattern speed for orbit integration.
measures ndarray

Array to store orbit measures.
Source code in commensurability/analysis.py 
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
class MPAnalysisBase(AnalysisBase):
    """
    Base class for analyzing commensurate orbits within galactic potentials.
    This makes use of the multiprocessing module for parallelization of orbit evaluation.

    This class provides methods for evaluating orbits, constructing images of
    phase space slices, and saving/loading analysis data.

    Attributes:
        ic_function (Callable[..., c.SkyCoord]): Function to generate initial conditions for orbits.
        ic_values (Mapping[str, Sequence[float]]): Values to pass into ic_function.
        axis_names (List[str]): Names of the parameters for the initial condition function.
        shape (Tuple[int]): Shape of the measures array.
        size (int): Size of the measures array.
        potential (Any): Potential generated by potential_function.
        dt (u.Quantity): Time step for orbit integration.
        steps (int): Number of integration steps.
        pattern_speed (u.Quantity): Pattern speed for orbit integration.
        measures (np.ndarray): Array to store orbit measures.
    """

    @staticmethod
    @abstractmethod
    def __eval__(orbit: c.SkyCoord) -> float:
        return 0.0

    def __init__(
        self,
        ic_function: Callable[..., c.SkyCoord],
        values: Mapping[str, Sequence[float]],
        /,
        potential_function: Callable[[], Any],
        dt: Union[float, np.ndarray, u.Quantity],
        steps: Union[int, np.ndarray],
        *,
        pattern_speed: Union[float, u.Quantity] = 0.0,
        backend: Optional[Union[str, Backend]] = None,
        progressbar: bool = True,
        pidgey_chunksize: Optional[int] = None,
        mp_chunksize: Optional[int] = None,
        _blank_measures: bool = False,
    ) -> None:
        super().__init__(
            ic_function,
            values,
            potential_function,
            dt,
            steps,
            pattern_speed=pattern_speed,
            backend=backend,
            progressbar=progressbar,
            pidgey_chunksize=pidgey_chunksize,
            _blank_measures=True,
        )
        if pidgey_chunksize is None:
            # set default chunk size for pidgey if not provided
            pidgey_chunksize = min(int(self.size**0.5), DEFAULT_MAX_CHUNKSIZE)
        if pidgey_chunksize <= 0:
            raise ValueError("pidgey_chunksize must be greater than 0")
        if pidgey_chunksize >= self.size:
            pidgey_chunksize = self.size
            # raise ValueError("chunksize must be less than total number of starting coordinates")

        if mp_chunksize is None:
            # set default chunk size for multiprocessing if not provided
            mp_chunksize = int(pidgey_chunksize**0.5)
        if mp_chunksize <= 0:
            raise ValueError("mp_chunksize must be greater than 0")
        if mp_chunksize > pidgey_chunksize:
            raise ValueError("mp_chunksize must not be greater than pidgey_chunksize")

        if not _blank_measures:
            self._construct_image_with_mp(pidgey_chunksize, mp_chunksize, progressbar)

    def _construct_image_with_mp(
        self, pidgey_chunksize: int = 1, mp_chunksize: int = 1, progressbar: bool = True
    ):
        for pixels in tqdm(
            chunked(np.ndindex(self.shape), pidgey_chunksize),
            desc=f"with {pidgey_chunksize=}",
            total=self.size // pidgey_chunksize,
            disable=not progressbar,
        ):
            coords = []
            for pixel in pixels:
                params = [self.ic_values[ax][i] for i, ax in zip(pixel, self.axis_names)]
                coord = self.ic_function(*params)
                coords.append(coord)
            coords = collapse_coords(coords)

            # NOTE: causes weird behavior when velocities are all set to 0
            # params = np.array(
            #     [
            #         [self.ic_values[ax][i] for i, ax in zip(pixel, self.axis_names)]
            #         for pixel in pixels
            #     ]
            # )
            # coords = self.ic_function(*params.T)

            orbits = self.backend.compute_orbit(
                coords,
                self.potential,
                self.dt,
                self.steps,
                pattern_speed=self.pattern_speed,
            )
            with Pool() as p:
                values = tuple(
                    tqdm(
                        p.imap(self.__eval__, orbits, chunksize=mp_chunksize),
                        desc=f"with {mp_chunksize=}",
                        total=pidgey_chunksize,
                        leave=False,
                    )
                )
            for pixel, value in zip(pixels, values):
                self.measures[pixel] = value

TessellationAnalysis

Bases: AnalysisBase3D

Analysis class for tessellation analysis on 3D orbits.

This class extends AnalysisBase3D and implements the evaluate method for tessellation analysis.

Attributes:

Name Type Description
ic_function Callable[..., SkyCoord]

Function to generate initial conditions for orbits.
ic_values Mapping[str, Sequence[float]]

Values to pass into ic_function.
axis_names List[str]

Names of the parameters for the initial condition function.
shape Tuple[int]

Shape of the measures array.
size int

Size of the measures array.
potential Any

Potential generated by potential_function.
dt Quantity

Time step for orbit integration.
steps int

Number of integration steps.
pattern_speed Quantity

Pattern speed for orbit integration.
measures ndarray

Array to store orbit measures.
Source code in commensurability/analysis.py 
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
class TessellationAnalysis(AnalysisBase3D):
    """
    Analysis class for tessellation analysis on 3D orbits.

    This class extends AnalysisBase3D and implements the evaluate method for tessellation analysis.

    Attributes:
        ic_function (Callable[..., c.SkyCoord]): Function to generate initial conditions for orbits.
        ic_values (Mapping[str, Sequence[float]]): Values to pass into ic_function.
        axis_names (List[str]): Names of the parameters for the initial condition function.
        shape (Tuple[int]): Shape of the measures array.
        size (int): Size of the measures array.
        potential (Any): Potential generated by potential_function.
        dt (u.Quantity): Time step for orbit integration.
        steps (int): Number of integration steps.
        pattern_speed (u.Quantity): Pattern speed for orbit integration.
        measures (np.ndarray): Array to store orbit measures.
    """

    @staticmethod
    def evaluate(orbit: c.SkyCoord) -> TessellationBase:
        """
        Evaluate an orbit using the tessellation and trimming algorithm.

        Args:
            orbit (c.SkyCoord): Integrated orbit coordinates.

        Returns:
            TessellationBase: Tessellation object containing tessellation results.
        """
        return Tessellation(orbit, incremental=False)

    @staticmethod
    def __eval__(orbit: c.SkyCoord) -> float:
        return Tessellation(orbit, incremental=False).measure

evaluate(orbit) staticmethod

Evaluate an orbit using the tessellation and trimming algorithm.

Parameters:

Name Type Description Default
orbit SkyCoord

Integrated orbit coordinates.

 required

Returns:

Name Type Description
TessellationBase TessellationBase

Tessellation object containing tessellation results.
Source code in commensurability/analysis.py 
532
533
534
535
536
537
538
539
540
541
542
543
@staticmethod
def evaluate(orbit: c.SkyCoord) -> TessellationBase:
    """
    Evaluate an orbit using the tessellation and trimming algorithm.

    Args:
        orbit (c.SkyCoord): Integrated orbit coordinates.

    Returns:
        TessellationBase: Tessellation object containing tessellation results.
    """
    return Tessellation(orbit, incremental=False)

TessellationAnalysis2D

Bases: AnalysisBase2D

Analysis class for tessellation analysis on 2D orbits.

This class extends AnalysisBase2D and implements the evaluate method for tessellation analysis.

Attributes:

Name Type Description
ic_function Callable[..., SkyCoord]

Function to generate initial conditions for orbits.
ic_values Mapping[str, Sequence[float]]

Values to pass into ic_function.
axis_names List[str]

Names of the parameters for the initial condition function.
shape Tuple[int]

Shape of the measures array.
size int

Size of the measures array.
potential Any

Potential generated by potential_function.
dt Quantity

Time step for orbit integration.
steps int

Number of integration steps.
pattern_speed Quantity

Pattern speed for orbit integration.
measures ndarray

Array to store orbit measures.
Source code in commensurability/analysis.py 
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
class TessellationAnalysis2D(AnalysisBase2D):
    """
    Analysis class for tessellation analysis on 2D orbits.

    This class extends AnalysisBase2D and implements the evaluate method for tessellation analysis.

    Attributes:
        ic_function (Callable[..., c.SkyCoord]): Function to generate initial conditions for orbits.
        ic_values (Mapping[str, Sequence[float]]): Values to pass into ic_function.
        axis_names (List[str]): Names of the parameters for the initial condition function.
        shape (Tuple[int]): Shape of the measures array.
        size (int): Size of the measures array.
        potential (Any): Potential generated by potential_function.
        dt (u.Quantity): Time step for orbit integration.
        steps (int): Number of integration steps.
        pattern_speed (u.Quantity): Pattern speed for orbit integration.
        measures (np.ndarray): Array to store orbit measures.
    """

    @staticmethod
    def evaluate(orbit) -> TessellationBase:
        """
        Evaluate an orbit using the tessellation and trimming algorithm.

        Args:
            orbit (c.SkyCoord): Integrated orbit coordinates.

        Returns:
            TessellationBase: Tessellation object containing tessellation results.
        """
        return Tessellation(orbit.xyz[:2].T, incremental=False)

    @staticmethod
    def __eval__(orbit: c.SkyCoord) -> float:
        return Tessellation(orbit.xyz[:2].T, incremental=False).measure

evaluate(orbit) staticmethod

Evaluate an orbit using the tessellation and trimming algorithm.

Parameters:

Name Type Description Default
orbit SkyCoord

Integrated orbit coordinates.

 required

Returns:

Name Type Description
TessellationBase TessellationBase

Tessellation object containing tessellation results.
Source code in commensurability/analysis.py 
569
570
571
572
573
574
575
576
577
578
579
580
@staticmethod
def evaluate(orbit) -> TessellationBase:
    """
    Evaluate an orbit using the tessellation and trimming algorithm.

    Args:
        orbit (c.SkyCoord): Integrated orbit coordinates.

    Returns:
        TessellationBase: Tessellation object containing tessellation results.
    """
    return Tessellation(orbit.xyz[:2].T, incremental=False)