Coverage for starry/maps.py : 90%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# -*- coding: utf-8 -*-
2from . import config
3from ._constants import *
4from ._core import OpsYlm, OpsLD, OpsReflected, OpsRV, linalg, math
5from ._indices import integers, get_ylm_inds, get_ul_inds, get_ylmw_inds
6from ._plotting import (
7 get_ortho_latitude_lines,
8 get_ortho_longitude_lines,
9 get_projection,
10)
11from ._sht import image2map, healpix2map, array2map
12import numpy as np
13import matplotlib
14import matplotlib.pyplot as plt
15from matplotlib import colors
16from matplotlib.animation import FuncAnimation
17from mpl_toolkits.axes_grid1 import make_axes_locatable
18from IPython.display import HTML
19from astropy import units
20import os
21import logging
23logger = logging.getLogger("starry.maps")
26__all__ = [
27 "Map",
28 "MapBase",
29 "YlmBase",
30 "LimbDarkenedBase",
31 "RVBase",
32 "ReflectedBase",
33]
36class Amplitude(object):
37 def __get__(self, instance, owner):
38 return instance._amp
40 def __set__(self, instance, value):
41 instance._amp = math.cast(np.ones(instance.nw) * value)
44class MapBase(object):
45 """The base class for all `starry` maps."""
47 # The map amplitude (just an attribute)
48 amp = Amplitude()
50 def _no_spectral(self):
51 if self.nw is not None: # pragma: no cover
52 raise NotImplementedError(
53 "Method not yet implemented for spectral maps."
54 )
56 def __init__(self, ydeg, udeg, fdeg, drorder, nw, **kwargs):
57 # Instantiate the Theano ops class
58 self.ops = self._ops_class_(ydeg, udeg, fdeg, drorder, nw, **kwargs)
60 # Dimensions
61 self._ydeg = ydeg
62 self._Ny = (ydeg + 1) ** 2
63 self._udeg = udeg
64 self._Nu = udeg + 1
65 self._fdeg = fdeg
66 self._Nf = (fdeg + 1) ** 2
67 self._deg = ydeg + udeg + fdeg
68 self._N = (ydeg + udeg + fdeg + 1) ** 2
69 self._nw = nw
70 self._drorder = drorder
72 # Basic properties
73 self._inc = math.cast(0.5 * np.pi)
74 self._obl = math.cast(0.0)
75 self._alpha = math.cast(0.0)
77 # Units
78 self.angle_unit = kwargs.pop("angle_unit", units.degree)
80 # Initialize
81 self.reset(**kwargs)
83 @property
84 def angle_unit(self):
85 """An ``astropy.units`` unit defining the angle metric for this map."""
86 return self._angle_unit
88 @angle_unit.setter
89 def angle_unit(self, value):
90 assert value.physical_type == "angle"
91 self._angle_unit = value
92 self._angle_factor = value.in_units(units.radian)
94 @property
95 def ydeg(self):
96 """Spherical harmonic degree of the map. *Read-only*"""
97 return self._ydeg
99 @property
100 def Ny(self):
101 r"""Number of spherical harmonic coefficients. *Read-only*
103 This is equal to :math:`(y_\mathrm{deg} + 1)^2`.
104 """
105 return self._Ny
107 @property
108 def udeg(self):
109 """Limb darkening degree. *Read-only*"""
110 return self._udeg
112 @property
113 def Nu(self):
114 r"""Number of limb darkening coefficients, including :math:`u_0`. *Read-only*
116 This is equal to :math:`u_\mathrm{deg} + 1`.
117 """
118 return self._Nu
120 @property
121 def fdeg(self):
122 """Degree of the multiplicative filter. *Read-only*"""
123 return self._fdeg
125 @property
126 def Nf(self):
127 r"""Number of spherical harmonic coefficients in the filter. *Read-only*
129 This is equal to :math:`(f_\mathrm{deg} + 1)^2`.
130 """
131 return self._Nf
133 @property
134 def deg(self):
135 r"""Total degree of the map. *Read-only*
137 This is equal to :math:`y_\mathrm{deg} + u_\mathrm{deg} + f_\mathrm{deg}`.
138 """
139 return self._deg
141 @property
142 def N(self):
143 r"""Total number of map coefficients. *Read-only*
145 This is equal to :math:`N_\mathrm{y} + N_\mathrm{u} + N_\mathrm{f}`.
146 """
147 return self._N
149 @property
150 def nw(self):
151 """Number of wavelength bins. *Read-only*"""
152 return self._nw
154 @property
155 def drorder(self):
156 """Differential rotation order. *Read-only*"""
157 return self._drorder
159 @property
160 def y(self):
161 """The spherical harmonic coefficient vector. *Read-only*
163 To set this vector, index the map directly using two indices:
164 ``map[l, m] = ...`` where ``l`` is the spherical harmonic degree and
165 ``m`` is the spherical harmonic order. These may be integers or
166 arrays of integers. Slice notation may also be used.
167 """
168 return self._y
170 @property
171 def u(self):
172 """The vector of limb darkening coefficients. *Read-only*
174 To set this vector, index the map directly using one index:
175 ``map[n] = ...`` where ``n`` is the degree of the limb darkening
176 coefficient. This may be an integer or an array of integers.
177 Slice notation may also be used.
178 """
179 return self._u
181 def __getitem__(self, idx):
182 if isinstance(idx, integers) or isinstance(idx, slice):
183 # User is accessing a limb darkening index
184 inds = get_ul_inds(self.udeg, idx)
185 return self._u[inds]
186 elif isinstance(idx, tuple) and len(idx) == 2 and self.nw is None:
187 # User is accessing a Ylm index
188 inds = get_ylm_inds(self.ydeg, idx[0], idx[1])
189 return self._y[inds]
190 elif isinstance(idx, tuple) and len(idx) == 3 and self.nw:
191 # User is accessing a Ylmw index
192 inds = get_ylmw_inds(self.ydeg, self.nw, idx[0], idx[1], idx[2])
193 return self._y[inds]
194 else:
195 raise ValueError("Invalid map index.")
197 def __setitem__(self, idx, val):
198 if isinstance(idx, integers) or isinstance(idx, slice):
199 # User is accessing a limb darkening index
200 inds = get_ul_inds(self.udeg, idx)
201 if 0 in inds:
202 raise ValueError("The u_0 coefficient cannot be set.")
203 if config.lazy:
204 self._u = self.ops.set_map_vector(self._u, inds, val)
205 else:
206 self._u[inds] = val
207 elif isinstance(idx, tuple) and len(idx) == 2 and self.nw is None:
208 # User is accessing a Ylm index
209 inds = get_ylm_inds(self.ydeg, idx[0], idx[1])
210 if 0 in inds:
211 raise ValueError("The Y_{0,0} coefficient cannot be set.")
212 if config.lazy:
213 self._y = self.ops.set_map_vector(self._y, inds, val)
214 else:
215 self._y[inds] = val
216 elif isinstance(idx, tuple) and len(idx) == 3 and self.nw:
217 # User is accessing a Ylmw index
218 inds = get_ylmw_inds(self.ydeg, self.nw, idx[0], idx[1], idx[2])
219 if 0 in inds[0]:
220 raise ValueError("The Y_{0,0} coefficients cannot be set.")
221 if config.lazy:
222 self._y = self.ops.set_map_vector(self._y, inds, val)
223 else:
224 old_shape = self._y[inds].shape
225 new_shape = np.atleast_2d(val).shape
226 if old_shape == new_shape:
227 self._y[inds] = val
228 elif old_shape == new_shape[::-1]:
229 self._y[inds] = np.atleast_2d(val).T
230 else:
231 self._y[inds] = val
232 else:
233 raise ValueError("Invalid map index.")
235 def _check_kwargs(self, method, kwargs):
236 if not config.quiet:
237 for key in kwargs.keys():
238 message = "Invalid keyword `{0}` in call to `{1}()`. Ignoring."
239 message = message.format(key, method)
240 logger.warning(message)
242 def _get_flux_kwargs(self, kwargs):
243 xo = kwargs.pop("xo", 0.0)
244 yo = kwargs.pop("yo", 0.0)
245 zo = kwargs.pop("zo", 1.0)
246 ro = kwargs.pop("ro", 0.0)
247 theta = kwargs.pop("theta", 0.0)
248 theta, xo, yo, zo = math.vectorize(theta, xo, yo, zo)
249 theta, xo, yo, zo, ro = math.cast(theta, xo, yo, zo, ro)
250 theta *= self._angle_factor
251 return theta, xo, yo, zo, ro
253 def reset(self, **kwargs):
254 """Reset all map coefficients and attributes.
256 .. note::
257 Does not reset custom unit settings.
259 """
260 if self.nw is None:
261 y = np.zeros(self.Ny)
262 y[0] = 1.0
263 else:
264 y = np.zeros((self.Ny, self.nw))
265 y[0, :] = 1.0
266 self._y = math.cast(y)
268 u = np.zeros(self.Nu)
269 u[0] = -1.0
270 self._u = math.cast(u)
272 f = np.zeros(self.Nf)
273 f[0] = np.pi
274 self._f = math.cast(f)
276 self._amp = math.cast(kwargs.pop("amp", np.ones(self.nw)))
278 # Reset data and priors
279 self._flux = None
280 self._C = None
281 self._mu = None
282 self._L = None
283 self._solution = None
285 self._check_kwargs("reset", kwargs)
287 def show(self, **kwargs):
288 """
289 Display an image of the map, with optional animation. See the
290 docstring of :py:meth:`render` for more details and additional
291 keywords accepted by this method.
293 Args:
294 ax (optional): A matplotlib axis instance to use. Default is
295 to create a new figure.
296 cmap (string or colormap instance, optional): The matplotlib colormap
297 to use. Defaults to ``plasma``.
298 figsize (tuple, optional): Figure size in inches. Default is
299 (3, 3) for orthographic maps and (7, 3.5) for rectangular
300 maps.
301 projection (string, optional): The map projection. Accepted
302 values are ``ortho``, corresponding to an orthographic
303 projection (as seen on the sky), and ``rect``, corresponding
304 to an equirectangular latitude-longitude projection.
305 Defaults to ``ortho``.
306 colorbar (bool, optional): Display a colorbar? Default is False.
307 grid (bool, optional): Show latitude/longitude grid lines?
308 Defaults to True.
309 interval (int, optional): Interval between frames in milliseconds
310 (animated maps only). Defaults to 75.
311 file (string, optional): The file name (including the extension)
312 to save the figure or animation to. Defaults to None.
313 html5_video (bool, optional): If rendering in a Jupyter notebook,
314 display as an HTML5 video? Default is True. If False, displays
315 the animation using Javascript (file size will be larger.)
316 dpi (int, optional): Image resolution in dots per square inch.
317 Defaults to the value defined in ``matplotlib.rcParams``.
318 bitrate (int, optional): Bitrate in kbps (animations only).
319 Defaults to the value defined in ``matplotlib.rcParams``.
320 norm (optional): The color normalization passed to
321 ``matplotlib.pyplot.imshow``. Default is None.
323 .. note::
324 Pure limb-darkened maps do not accept a ``projection`` keyword.
326 """
327 # Get kwargs
328 cmap = kwargs.pop("cmap", "plasma")
329 grid = kwargs.pop("grid", True)
330 interval = kwargs.pop("interval", 75)
331 file = kwargs.pop("file", None)
332 html5_video = kwargs.pop("html5_video", True)
333 norm = kwargs.pop("norm", None)
334 dpi = kwargs.pop("dpi", None)
335 figsize = kwargs.pop("figsize", None)
336 bitrate = kwargs.pop("bitrate", None)
337 colorbar = kwargs.pop("colorbar", False)
338 ax = kwargs.pop("ax", None)
339 if ax is None:
340 custom_ax = False
341 else:
342 custom_ax = True
344 # Ylm-base maps only
345 if not self.__props__["limbdarkened"]:
347 projection = get_projection(kwargs.get("projection", "ortho"))
349 # Get the map orientation
350 if config.lazy:
351 inc = self._inc.eval()
352 obl = self._obl.eval()
353 else:
354 inc = self._inc
355 obl = self._obl
357 # Get the rotational phase
358 if config.lazy:
359 theta = math.vectorize(
360 math.cast(kwargs.pop("theta", 0.0)) * self._angle_factor
361 ).eval()
362 else:
363 theta = np.atleast_1d(
364 np.array(kwargs.pop("theta", 0.0)) * self._angle_factor
365 )
367 else:
369 inc = np.array(0.5 * np.pi)
370 obl = np.array(0)
371 theta = np.array([0])
373 # Render the map if needed
374 image = kwargs.pop("image", None)
375 if image is None:
377 # We need to evaluate the variables so we can plot the map!
378 if config.lazy:
380 # Get kwargs
381 res = kwargs.pop("res", 300)
383 # Evaluate the variables
384 u = self._u.eval()
386 if not self.__props__["limbdarkened"]:
388 inc = self._inc.eval()
389 obl = self._obl.eval()
390 y = self._y.eval()
391 f = self._f.eval()
392 alpha = self._alpha.eval()
394 # Explicitly call the compiled version of `render`
395 image = self.amp.eval().reshape(
396 -1, 1, 1
397 ) * self.ops.render(
398 res, projection, theta, inc, obl, y, u, f, alpha
399 )
401 else:
403 # Explicitly call the compiled version of `render`
404 image = self.amp.eval().reshape(
405 -1, 1, 1
406 ) * self.ops.render_ld(res, u)
408 else:
410 # Easy!
411 if not self.__props__["limbdarkened"]:
412 image = self.render(
413 theta=theta / self._angle_factor, **kwargs
414 )
415 else:
416 image = self.render(**kwargs)
417 kwargs.pop("res", None)
419 if len(image.shape) == 3:
420 nframes = image.shape[0]
421 else:
422 nframes = 1
423 image = np.reshape(image, (1,) + image.shape)
425 # Animation
426 animated = nframes > 1
427 borders = []
428 latlines = []
429 lonlines = []
431 if (
432 not self.__props__["limbdarkened"]
433 and projection == STARRY_RECTANGULAR_PROJECTION
434 ):
435 # Set up the plot
436 if figsize is None:
437 figsize = (7, 3.75)
438 if ax is None:
439 fig, ax = plt.subplots(1, figsize=figsize)
440 else:
441 fig = ax.figure
442 extent = (-180, 180, -90, 90)
444 # Grid lines
445 if grid:
446 lats = np.linspace(-90, 90, 7)[1:-1]
447 lons = np.linspace(-180, 180, 13)
448 latlines = [None for n in lats]
449 for n, lat in enumerate(lats):
450 latlines[n] = ax.axhline(
451 lat, color="k", lw=0.5, alpha=0.5, zorder=100
452 )
453 lonlines = [None for n in lons]
454 for n, lon in enumerate(lons):
455 lonlines[n] = ax.axvline(
456 lon, color="k", lw=0.5, alpha=0.5, zorder=100
457 )
458 ax.set_xticks(lons)
459 ax.set_yticks(lats)
460 ax.set_xlabel("Longitude [deg]")
461 ax.set_ylabel("Latitude [deg]")
463 else:
464 # Set up the plot
465 if figsize is None:
466 figsize = (3, 3)
467 if ax is None:
468 fig, ax = plt.subplots(1, figsize=figsize)
469 else:
470 fig = ax.figure
471 ax.axis("off")
472 ax.set_xlim(-1.05, 1.05)
473 ax.set_ylim(-1.05, 1.05)
474 extent = (-1, 1, -1, 1)
476 # Grid lines
477 if grid:
478 x = np.linspace(-1, 1, 10000)
479 y = np.sqrt(1 - x ** 2)
480 borders = [None, None]
481 (borders[0],) = ax.plot(x, y, "k-", alpha=1, lw=1)
482 (borders[1],) = ax.plot(x, -y, "k-", alpha=1, lw=1)
483 lats = get_ortho_latitude_lines(inc=inc, obl=obl)
484 latlines = [None for n in lats]
485 for n, l in enumerate(lats):
486 (latlines[n],) = ax.plot(
487 l[0], l[1], "k-", lw=0.5, alpha=0.5, zorder=100
488 )
489 lons = get_ortho_longitude_lines(
490 inc=inc, obl=obl, theta=theta[0]
491 )
492 lonlines = [None for n in lons]
493 for n, l in enumerate(lons):
494 (lonlines[n],) = ax.plot(
495 l[0], l[1], "k-", lw=0.5, alpha=0.5, zorder=100
496 )
498 # Plot the first frame of the image
499 if norm is None or norm == "rv":
500 vmin = np.nanmin(image)
501 vmax = np.nanmax(image)
502 if vmin == vmax:
503 vmin -= 1e-15
504 vmax += 1e-15
505 if norm is None:
506 norm = colors.Normalize(vmin=vmin, vmax=vmax)
507 elif norm == "rv":
508 try:
509 norm = colors.DivergingNorm(
510 vmin=vmin, vcenter=0, vmax=vmax
511 )
512 except AttributeError: # pragma: no cover
513 # DivergingNorm was introduced in matplotlib 3.1
514 norm = colors.Normalize(vmin=vmin, vmax=vmax)
515 img = ax.imshow(
516 image[0],
517 origin="lower",
518 extent=extent,
519 cmap=cmap,
520 norm=norm,
521 interpolation="none",
522 animated=animated,
523 )
525 # Add a colorbar
526 if colorbar:
527 if not custom_ax:
528 fig.subplots_adjust(right=0.85)
529 divider = make_axes_locatable(ax)
530 cax = divider.append_axes("right", size="5%", pad=0.05)
531 fig.colorbar(img, cax=cax, orientation="vertical")
533 # Display or save the image / animation
534 if animated:
536 def updatefig(i):
537 img.set_array(image[i])
538 if (
539 not self.__props__["limbdarkened"]
540 and projection == STARRY_ORTHOGRAPHIC_PROJECTION
541 and grid
542 and len(theta) > 1
543 and self.nw is None
544 ):
545 lons = get_ortho_longitude_lines(
546 inc=inc, obl=obl, theta=theta[i]
547 )
548 for n, l in enumerate(lons):
549 lonlines[n].set_xdata(l[0])
550 lonlines[n].set_ydata(l[1])
551 return (img,) + tuple(lonlines + latlines + borders)
553 ani = FuncAnimation(
554 fig,
555 updatefig,
556 interval=interval,
557 blit=True,
558 frames=image.shape[0],
559 )
561 # Business as usual
562 if (file is not None) and (file != ""):
563 if file.endswith(".mp4"):
564 ani.save(file, writer="ffmpeg", dpi=dpi, bitrate=bitrate)
565 elif file.endswith(".gif"):
566 ani.save(
567 file, writer="imagemagick", dpi=dpi, bitrate=bitrate
568 )
569 else:
570 # Try and see what happens!
571 ani.save(file, dpi=dpi, bitrate=bitrate)
572 if not custom_ax:
573 plt.close()
574 else: # if not custom_ax:
575 try:
576 if "zmqshell" in str(type(get_ipython())):
577 plt.close()
578 with matplotlib.rc_context(
579 {
580 "savefig.dpi": dpi
581 if dpi is not None
582 else "figure",
583 "animation.bitrate": bitrate
584 if bitrate is not None
585 else -1,
586 }
587 ):
588 if html5_video:
589 display(HTML(ani.to_html5_video()))
590 else:
591 display(HTML(ani.to_jshtml()))
592 else:
593 raise NameError("")
594 except NameError:
595 plt.show()
596 plt.close()
598 # Matplotlib generates an annoying empty
599 # file when producing an animation. Delete it.
600 try:
601 os.remove("None0000000.png")
602 except FileNotFoundError:
603 pass
605 else:
606 if (file is not None) and (file != ""):
607 if (
608 not self.__props__["limbdarkened"]
609 and projection == STARRY_ORTHOGRAPHIC_PROJECTION
610 ):
611 fig.subplots_adjust(
612 left=0.01, right=0.99, bottom=0.01, top=0.99
613 )
614 fig.savefig(file)
615 if not custom_ax:
616 plt.close()
617 elif not custom_ax:
618 plt.show()
620 # Check for invalid kwargs
621 if self.__props__["rv"]:
622 kwargs.pop("rv", None)
623 if not self.__props__["limbdarkened"]:
624 kwargs.pop("projection", None)
625 if self.__props__["reflected"]:
626 kwargs.pop("xs", None)
627 kwargs.pop("ys", None)
628 kwargs.pop("zs", None)
629 self._check_kwargs("show", kwargs)
631 def limbdark_is_physical(self):
632 """Check whether the limb darkening profile (if any) is physical.
634 This method uses Sturm's theorem to ensure that the limb darkening
635 intensity is positive everywhere and decreases monotonically toward
636 the limb.
638 Returns:
639 bool: Whether or not the limb darkening profile is physical.
640 """
641 result = self.ops.limbdark_is_physical(self.u)
642 if config.lazy:
643 return result
644 else:
645 return bool(result)
647 def set_data(self, flux, C=None, cho_C=None):
648 """Set the data vector and covariance matrix.
650 This method is required by the :py:meth:`solve` method, which
651 analytically computes the posterior over surface maps given a
652 dataset and a prior, provided both are described as multivariate
653 Gaussians.
655 Args:
656 flux (vector): The observed light curve.
657 C (scalar, vector, or matrix): The data covariance. This may be
658 a scalar, in which case the noise is assumed to be
659 homoscedastic, a vector, in which case the covariance
660 is assumed to be diagonal, or a matrix specifying the full
661 covariance of the dataset. Default is None. Either `C` or
662 `cho_C` must be provided.
663 cho_C (matrix): The lower Cholesky factorization of the data
664 covariance matrix. Defaults to None. Either `C` or
665 `cho_C` must be provided.
666 """
667 self._flux = math.cast(flux)
668 self._C = linalg.Covariance(C, cho_C, N=self._flux.shape[0])
670 def set_prior(self, *, mu=None, L=None, cho_L=None):
671 """Set the prior mean and covariance of the spherical harmonic coefficients.
673 This method is required by the :py:meth:`solve` method, which
674 analytically computes the posterior over surface maps given a
675 dataset and a prior, provided both are described as multivariate
676 Gaussians.
678 Note that the prior is placed on the **amplitude-weighted** coefficients,
679 i.e., the quantity ``x = map.amp * map.y``. Because the first spherical
680 harmonic coefficient is fixed at unity, ``x[0]`` is
681 the amplitude of the map. The actual spherical harmonic coefficients
682 are given by ``x / map.amp``.
684 This convention allows one to linearly fit for an arbitrary map normalization
685 at the same time as the spherical harmonic coefficients, while ensuring
686 the ``starry`` requirement that the coefficient of the :math:`Y_{0,0}`
687 harmonic is always unity.
689 Args:
690 mu (scalar or vector): The prior mean on the amplitude-weighted
691 spherical harmonic coefficients. Default is unity for the
692 first term and zero for the remaining terms. If this is a vector,
693 it must have length equal to :py:attr:`Ny`.
694 L (scalar, vector, or matrix): The prior covariance. This may be
695 a scalar, in which case the covariance is assumed to be
696 homoscedastic, a vector, in which case the covariance
697 is assumed to be diagonal, or a matrix specifying the full
698 prior covariance. Default is None. Either `L` or
699 `cho_L` must be provided.
700 cho_L (matrix): The lower Cholesky factorization of the prior
701 covariance matrix. Defaults to None. Either `L` or
702 `cho_L` must be provided.
704 """
705 if mu is None:
706 mu = np.zeros(self.Ny)
707 mu[0] = 1.0
708 mu = math.cast(mu)
709 self._mu = math.cast(mu) * math.cast(np.ones(self.Ny))
710 self._L = linalg.Covariance(L, cho_L, N=self.Ny)
712 def remove_prior(self):
713 """Remove the prior on the map coefficients."""
714 self._mu = None
715 self._L = None
717 def solve(self, *, design_matrix=None, **kwargs):
718 """Solve the linear least-squares problem for the posterior over maps.
720 This method solves the generalized least squares problem given a
721 light curve and its covariance (set via the :py:meth:`set_data` method)
722 and a Gaussian prior on the spherical harmonic coefficients
723 (set via the :py:meth:`set_prior` method). The map amplitude and
724 coefficients are set to the maximum a posteriori (MAP) solution.
726 Args:
727 design_matrix (matrix, optional): The flux design matrix, the
728 quantity returned by :py:meth:`design_matrix`. Default is
729 None, in which case this is computed based on ``kwargs``.
730 kwargs (optional): Keyword arguments to be passed directly to
731 :py:meth:`design_matrix`, if a design matrix is not provided.
733 Returns:
734 A tuple containing the posterior mean for the amplitude-weighted \
735 spherical harmonic coefficients (a vector) and the Cholesky factorization \
736 of the posterior covariance (a lower triangular matrix).
738 .. note::
739 Users may call :py:meth:`draw` to draw from the
740 posterior after calling this method.
741 """
742 # Not implemented for spectral
743 self._no_spectral()
745 if self._flux is None or self._C is None:
746 raise ValueError("Please provide a dataset with `set_data()`.")
747 elif self._mu is None or self._L is None:
748 raise ValueError("Please provide a prior with `set_prior()`.")
750 # Get the design matrix & remove any amplitude weighting
751 if design_matrix is None:
752 design_matrix = self.design_matrix(**kwargs)
753 X = math.cast(design_matrix)
755 # Compute the MAP solution
756 self._solution = linalg.solve(
757 X, self._flux, self._C.cholesky, self._mu, self._L.inverse
758 )
760 # Set the amplitude and coefficients
761 x, _ = self._solution
762 self.amp = x[0]
763 self[1:, :] = x[1:] / self.amp
765 # Return the mean and covariance
766 return self._solution
768 def lnlike(self, *, design_matrix=None, woodbury=True, **kwargs):
769 """Returns the log marginal likelihood of the data given a design matrix.
771 This method computes the marginal likelihood (marginalized over the
772 spherical harmonic coefficients) given a
773 light curve and its covariance (set via the :py:meth:`set_data` method)
774 and a Gaussian prior on the spherical harmonic coefficients
775 (set via the :py:meth:`set_prior` method).
777 Args:
778 design_matrix (matrix, optional): The flux design matrix, the
779 quantity returned by :py:meth:`design_matrix`. Default is
780 None, in which case this is computed based on ``kwargs``.
781 woodbury (bool, optional): Solve the linear problem using the
782 Woodbury identity? Default is True. The
783 `Woodbury identity <https://en.wikipedia.org/wiki/Woodbury_matrix_identity>`_
784 is used to speed up matrix operations in the case that the
785 number of data points is much larger than the number of
786 spherical harmonic coefficients. In this limit, it can
787 speed up the code by more than an order of magnitude. Keep
788 in mind that the numerical stability of the Woodbury identity
789 is not great, so if you're getting strange results try
790 disabling this. It's also a good idea to disable this in the
791 limit of few data points and large spherical harmonic degree.
792 kwargs (optional): Keyword arguments to be passed directly to
793 :py:meth:`design_matrix`, if a design matrix is not provided.
795 Returns:
796 The log marginal likelihood, a scalar.
797 """
798 # Not implemented for spectral
799 self._no_spectral()
801 if self._flux is None or self._C is None:
802 raise ValueError("Please provide a dataset with `set_data()`.")
803 elif self._mu is None or self._L is None:
804 raise ValueError("Please provide a prior with `set_prior()`.")
806 # Get the design matrix & remove any amplitude weighting
807 if design_matrix is None:
808 design_matrix = self.design_matrix(**kwargs)
809 X = math.cast(design_matrix)
811 # Compute the likelihood
812 if woodbury:
813 return linalg.lnlike_woodbury(
814 X,
815 self._flux,
816 self._C.inverse,
817 self._mu,
818 self._L.inverse,
819 self._C.lndet,
820 self._L.lndet,
821 )
822 else:
823 return linalg.lnlike(
824 X, self._flux, self._C.value, self._mu, self._L.value
825 )
827 @property
828 def solution(self):
829 r"""The posterior probability distribution for the map.
831 This is a tuple containing the mean and lower Cholesky factorization of the
832 covariance of the amplitude-weighted spherical harmonic coefficient vector,
833 obtained by solving the regularized least-squares problem
834 via the :py:meth:`solve` method.
836 Note that to obtain the actual covariance matrix from the lower Cholesky
837 factorization :math:`L`, simply compute :math:`L L^\top`.
839 Note also that this is the posterior for the **amplitude-weighted**
840 map vector. Under this convention, the map amplitude is equal to the
841 first term of the vector and the spherical harmonic coefficients are
842 equal to the vector normalized by the first term.
843 """
844 if self._solution is None:
845 raise ValueError("Please call `solve()` first.")
846 return self._solution
848 def draw(self):
849 """Draw a map from the posterior distribution.
851 This method draws a random map from the posterior distribution and
852 sets the :py:attr:`y` map vector and :py:attr:`amp` map amplitude
853 accordingly. Users should call :py:meth:`solve` to enable this
854 attribute.
855 """
856 if self._solution is None:
857 raise ValueError("Please call `solve()` first.")
859 # Fast multivariate sampling using the Cholesky factorization
860 yhat, cho_ycov = self._solution
861 u = math.cast(np.random.randn(self.Ny))
862 x = yhat + math.dot(cho_ycov, u)
863 self.amp = x[0]
864 self[1:, :] = x[1:] / self.amp
867class YlmBase(object):
868 """The default ``starry`` map class.
870 This class handles light curves and phase curves of objects in
871 emitted light.
873 .. note::
874 Instantiate this class by calling :py:func:`starry.Map` with
875 ``ydeg > 0`` and both ``rv`` and ``reflected`` set to False.
876 """
878 _ops_class_ = OpsYlm
880 def reset(self, **kwargs):
881 if kwargs.get("inc", None) is not None:
882 self.inc = kwargs.pop("inc")
883 else:
884 self._inc = math.cast(0.5 * np.pi)
886 if kwargs.get("obl", None) is not None:
887 self.obl = kwargs.pop("obl")
888 else:
889 self._obl = math.cast(0.0)
891 if kwargs.get("alpha", None) is not None:
892 self.alpha = kwargs.pop("alpha")
893 else:
894 self._alpha = math.cast(0.0)
896 super(YlmBase, self).reset(**kwargs)
898 @property
899 def inc(self):
900 """The inclination of the rotation axis in units of :py:attr:`angle_unit`."""
901 return self._inc / self._angle_factor
903 @inc.setter
904 def inc(self, value):
905 self._inc = math.cast(value) * self._angle_factor
907 @property
908 def obl(self):
909 """The obliquity of the rotation axis in units of :py:attr:`angle_unit`."""
910 return self._obl / self._angle_factor
912 @obl.setter
913 def obl(self, value):
914 self._obl = math.cast(value) * self._angle_factor
916 @property
917 def alpha(self):
918 """The rotational shear coefficient, a number in the range ``[0, 1]``.
920 The parameter :math:`\\alpha` is used to model linear differential
921 rotation. The angular velocity at a given latitude :math:`\\theta`
922 is
924 :math:`\\omega = \\omega_{eq}(1 - \\alpha \\sin^2\\theta)`
926 where :math:`\\omega_{eq}` is the equatorial angular velocity of
927 the object.
929 .. note ::
931 This parameter is only used if :py:attr:`drorder` is greater
932 than zero and/or radial velocity mode is enabled.
933 """
934 return self._alpha
936 @alpha.setter
937 def alpha(self, value):
938 if (self._drorder == 0) and not hasattr(
939 self, "rv"
940 ): # pragma: no cover
941 logger.warning(
942 "Parameter `drorder` is zero, so setting `alpha` has no effect."
943 )
944 else:
945 self._alpha = math.cast(value)
947 def design_matrix(self, **kwargs):
948 r"""Compute and return the light curve design matrix :math:`A`.
950 This matrix is used to compute the flux :math:`f` from a vector of spherical
951 harmonic coefficients :math:`y` and the map amplitude :math:`a`:
952 :math:`f = a A y`.
954 Args:
955 xo (scalar or vector, optional): x coordinate of the occultor
956 relative to this body in units of this body's radius.
957 yo (scalar or vector, optional): y coordinate of the occultor
958 relative to this body in units of this body's radius.
959 zo (scalar or vector, optional): z coordinate of the occultor
960 relative to this body in units of this body's radius.
961 ro (scalar, optional): Radius of the occultor in units of
962 this body's radius.
963 theta (scalar or vector, optional): Angular phase of the body
964 in units of :py:attr:`angle_unit`.
965 """
966 # Orbital kwargs
967 theta, xo, yo, zo, ro = self._get_flux_kwargs(kwargs)
969 # Check for invalid kwargs
970 self._check_kwargs("design_matrix", kwargs)
972 # Compute & return
973 return self.ops.X(
974 theta,
975 xo,
976 yo,
977 zo,
978 ro,
979 self._inc,
980 self._obl,
981 self._u,
982 self._f,
983 self._alpha,
984 )
986 def intensity_design_matrix(self, lat=0, lon=0):
987 """Compute and return the pixelization matrix ``P``.
989 This matrix is used to compute the intensity :math:`I` from a vector of spherical
990 harmonic coefficients :math:`y` and the map amplitude :math:`a`:
991 :math:`I = a P y`.
993 Args:
994 lat (scalar or vector, optional): latitude at which to evaluate
995 the design matrix in units of :py:attr:`angle_unit`.
996 lon (scalar or vector, optional): longitude at which to evaluate
997 the design matrix in units of :py:attr:`angle_unit`.
999 .. note::
1000 This method ignores any filters (such as limb darkening
1001 or velocity weighting) and illumination (for reflected light
1002 maps).
1004 """
1005 # Get the Cartesian points
1006 lat, lon = math.vectorize(*math.cast(lat, lon))
1007 lat *= self._angle_factor
1008 lon *= self._angle_factor
1010 # Compute & return
1011 return self.ops.P(lat, lon)
1013 def flux(self, **kwargs):
1014 """
1015 Compute and return the light curve.
1017 Args:
1018 xo (scalar or vector, optional): x coordinate of the occultor
1019 relative to this body in units of this body's radius.
1020 yo (scalar or vector, optional): y coordinate of the occultor
1021 relative to this body in units of this body's radius.
1022 zo (scalar or vector, optional): z coordinate of the occultor
1023 relative to this body in units of this body's radius.
1024 ro (scalar, optional): Radius of the occultor in units of
1025 this body's radius.
1026 theta (scalar or vector, optional): Angular phase of the body
1027 in units of :py:attr:`angle_unit`.
1028 """
1029 # Orbital kwargs
1030 theta, xo, yo, zo, ro = self._get_flux_kwargs(kwargs)
1032 # Check for invalid kwargs
1033 self._check_kwargs("flux", kwargs)
1035 # Compute & return
1036 return self.amp * self.ops.flux(
1037 theta,
1038 xo,
1039 yo,
1040 zo,
1041 ro,
1042 self._inc,
1043 self._obl,
1044 self._y,
1045 self._u,
1046 self._f,
1047 self._alpha,
1048 )
1050 def intensity(self, lat=0, lon=0, **kwargs):
1051 """
1052 Compute and return the intensity of the map.
1054 Args:
1055 lat (scalar or vector, optional): latitude at which to evaluate
1056 the intensity in units of :py:attr:`angle_unit`.
1057 lon (scalar or vector, optional): longitude at which to evaluate
1058 the intensity in units of :py:attr:`angle_unit``.
1059 theta (scalar, optional): For differentially rotating maps only,
1060 the angular phase at which to evaluate the intensity.
1061 Default 0.
1062 limbdarken (bool, optional): Apply limb darkening
1063 (only if :py:attr:`udeg` > 0)? Default True.
1065 """
1066 # Get the Cartesian points
1067 lat, lon = math.vectorize(*math.cast(lat, lon))
1068 lat *= self._angle_factor
1069 lon *= self._angle_factor
1071 # If differentially rotating, allow a `theta` keyword
1072 if self.drorder > 0:
1073 alpha_theta = math.cast(kwargs.get("theta", 0.0)) * self.alpha
1074 alpha_theta *= self._angle_factor
1075 else:
1076 alpha_theta = math.cast(0.0)
1078 # If limb-darkened, allow a `limbdarken` keyword
1079 if self.udeg > 0 and kwargs.pop("limbdarken", True):
1080 ld = np.array(True)
1081 else:
1082 ld = np.array(False)
1084 # Check for invalid kwargs
1085 self._check_kwargs("intensity", kwargs)
1087 # Compute & return
1088 return self.amp * self.ops.intensity(
1089 lat, lon, self._y, self._u, self._f, alpha_theta, ld
1090 )
1092 def render(self, res=300, projection="ortho", theta=0.0):
1093 """Compute and return the intensity of the map on a grid.
1095 Returns an image of shape ``(res, res)``, unless ``theta`` is a vector,
1096 in which case returns an array of shape ``(nframes, res, res)``, where
1097 ``nframes`` is the number of values of ``theta``. However, if this is
1098 a spectral map, ``nframes`` is the number of wavelength bins and
1099 ``theta`` must be a scalar.
1101 Args:
1102 res (int, optional): The resolution of the map in pixels on a
1103 side. Defaults to 300.
1104 projection (string, optional): The map projection. Accepted
1105 values are ``ortho``, corresponding to an orthographic
1106 projection (as seen on the sky), and ``rect``, corresponding
1107 to an equirectangular latitude-longitude projection.
1108 Defaults to ``ortho``.
1109 theta (scalar or vector, optional): The map rotation phase in
1110 units of :py:attr:`angle_unit`. If this is a vector, an
1111 animation is generated. Defaults to ``0.0``.
1112 """
1113 # Multiple frames?
1114 if self.nw is not None:
1115 animated = True
1116 else:
1117 if config.lazy:
1118 animated = hasattr(theta, "ndim") and theta.ndim > 0
1119 else:
1120 animated = hasattr(theta, "__len__")
1122 # Convert
1123 projection = get_projection(projection)
1124 theta = math.vectorize(math.cast(theta) * self._angle_factor)
1126 # Compute
1127 if self.nw is None or config.lazy:
1128 amp = self.amp
1129 else:
1130 # The intensity has shape `(nw, res, res)`
1131 # so we must reshape `amp` to take the product correctly
1132 amp = self.amp[:, np.newaxis, np.newaxis]
1133 image = amp * self.ops.render(
1134 res,
1135 projection,
1136 theta,
1137 self._inc,
1138 self._obl,
1139 self._y,
1140 self._u,
1141 self._f,
1142 self._alpha,
1143 )
1145 # Squeeze?
1146 if animated:
1147 return image
1148 else:
1149 return math.reshape(image, [res, res])
1151 def load(
1152 self,
1153 image,
1154 healpix=False,
1155 nside=32,
1156 max_iter=3,
1157 sigma=None,
1158 force_psd=False,
1159 **kwargs
1160 ):
1161 """Load an image, array, or ``healpix`` map.
1163 This routine uses various routines in ``healpix`` to compute the
1164 spherical harmonic expansion of the input image and sets the map's
1165 :py:attr:`y` coefficients accordingly.
1167 Args:
1168 image: A path to an image file, a two-dimensional ``numpy``
1169 array, or a ``healpix`` map array (if ``healpix`` is True).
1170 healpix (bool, optional): Treat ``image`` as a ``healpix`` array?
1171 Default is False.
1172 sigma (float, optional): If not None, apply gaussian smoothing
1173 with standard deviation ``sigma`` to smooth over
1174 spurious ringing features. Smoothing is performed with
1175 the ``healpix.sphtfunc.smoothalm`` method.
1176 Default is None.
1177 force_psd (bool, optional): Force the map to be positive
1178 semi-definite? Default is False.
1179 nside (int, optional): The ``NSIDE`` argument to a Healpix
1180 map. This controls the angular resolution of the image; increase
1181 this for maps with higher fidelity, at the expense of extra
1182 compute time. Default is 32.
1183 max_iter (int, optional): Maximum number of iterations in trying
1184 to convert the map to a Healpix array. Each iteration, the
1185 input array is successively doubled in size in order to
1186 increase the angular coverage of the input. Default is 3. If
1187 the number of iterations exceeds this value, an error will
1188 be raised.
1189 kwargs (optional): Any other kwargs passed directly to
1190 :py:meth:`minimize` (only if ``psd`` is True).
1191 """
1192 # Not implemented for spectral
1193 self._no_spectral()
1195 # Is this a file name?
1196 if type(image) is str:
1197 y = image2map(
1198 image,
1199 lmax=self.ydeg,
1200 sigma=sigma,
1201 nside=nside,
1202 max_iter=max_iter,
1203 )
1204 # or is it an array?
1205 elif type(image) is np.ndarray:
1206 if healpix:
1207 y = healpix2map(
1208 image,
1209 lmax=self.ydeg,
1210 sigma=sigma,
1211 nside=nside,
1212 max_iter=max_iter,
1213 )
1214 else:
1215 y = array2map(
1216 image,
1217 lmax=self.ydeg,
1218 sigma=sigma,
1219 nside=nside,
1220 max_iter=max_iter,
1221 )
1222 else:
1223 raise ValueError("Invalid `image` value.")
1225 # Ingest the coefficients w/ appropriate normalization
1226 # This ensures the map intensity will have the same normalization
1227 # as that of the input image
1228 y /= 2 * np.sqrt(np.pi)
1229 self._y = math.cast(y / y[0])
1230 self._amp = math.cast(y[0] * np.pi)
1232 # Ensure positive semi-definite?
1233 if force_psd:
1235 # Find the minimum
1236 _, _, I = self.minimize(**kwargs)
1237 if config.lazy:
1238 I = I.eval()
1240 # Scale the coeffs?
1241 if I < 0:
1242 fac = self._amp / (self._amp - np.pi * I)
1243 if config.lazy:
1244 self._y *= fac
1245 self._y = self.ops.set_map_vector(self._y, 0, 1.0)
1246 else:
1247 self._y[1:] *= fac
1249 def rotate(self, axis, theta):
1250 """Rotate the current map vector an angle ``theta`` about ``axis``.
1252 Args:
1253 axis (vector): The axis about which to rotate the map.
1254 theta (scalar): The angle of (counter-clockwise) rotation.
1255 """
1256 axis = math.cast(axis)
1257 axis /= math.sqrt(math.sum(axis ** 2))
1258 # Note that we rotate by -theta since
1259 # this is the *RHS* rotation operator
1260 y = self.ops.dotR(
1261 math.transpose(
1262 math.reshape(self.y, (-1, 1 if self.nw is None else self.nw))
1263 ),
1264 axis[0],
1265 axis[1],
1266 axis[2],
1267 -math.cast(theta * self._angle_factor),
1268 )
1269 if self.nw is None:
1270 self._y = y[0]
1271 else:
1272 self._y = math.transpose(y)
1274 def add_spot(
1275 self,
1276 amp=None,
1277 intensity=None,
1278 relative=True,
1279 sigma=0.1,
1280 lat=0.0,
1281 lon=0.0,
1282 ):
1283 r"""Add the expansion of a gaussian spot to the map.
1285 This function adds a spot whose functional form is the spherical
1286 harmonic expansion of a gaussian in the quantity
1287 :math:`\cos\Delta\theta`, where :math:`\Delta\theta`
1288 is the angular separation between the center of the spot and another
1289 point on the surface. The spot brightness is controlled by either the
1290 parameter ``amp``, defined as the fractional change in the
1291 total luminosity of the object due to the spot, or the parameter
1292 ``intensity``, defined as the fractional change in the
1293 intensity at the center of the spot.
1295 Args:
1296 amp (scalar or vector, optional): The amplitude of the spot. This
1297 is equal to the fractional change in the luminosity of the map
1298 due to the spot. If the map has more than one wavelength bin,
1299 this must be a vector of length equal to the number of
1300 wavelength bins. Default is None.
1301 Either ``amp`` or ``intensity`` must be given.
1302 intensity (scalar or vector, optional): The intensity of the spot.
1303 This is equal to the fractional change in the intensity of the
1304 map at the *center* of the spot. If the map has more than one
1305 wavelength bin, this must be a vector of length equal to the
1306 number of wavelength bins. Default is None.
1307 Either ``amp`` or ``intensity`` must be given.
1308 relative (bool, optional): If True, computes the spot expansion
1309 assuming the fractional `amp` or `intensity` change is relative
1310 to the **current** map amplitude/intensity. If False, computes
1311 the spot expansion assuming the fractional change is relative
1312 to the **original** map amplitude/intensity (i.e., that of
1313 a featureless map). Defaults to True. Note that if True,
1314 adding two spots with the same values of `amp` or `intensity`
1315 will generally result in *different* intensities at their
1316 centers, since the first spot will have changed the map
1317 intensity everywhere! Defaults to True.
1318 sigma (scalar, optional): The standard deviation of the gaussian.
1319 Defaults to 0.1.
1320 lat (scalar, optional): The latitude of the spot in units of
1321 :py:attr:`angle_unit`. Defaults to 0.0.
1322 lon (scalar, optional): The longitude of the spot in units of
1323 :py:attr:`angle_unit`. Defaults to 0.0.
1325 """
1326 # Parse the amplitude
1327 if (amp is None and intensity is None) or (
1328 amp is not None and intensity is not None
1329 ):
1330 raise ValueError("Please provide either `amp` or `intensity`.")
1331 elif amp is not None:
1332 amp, _ = math.vectorize(math.cast(amp), np.ones(self.nw))
1333 # Normalize?
1334 if not relative:
1335 amp /= self.amp
1336 else:
1337 # Vectorize if needed
1338 intensity, _ = math.vectorize(
1339 math.cast(intensity), np.ones(self.nw)
1340 )
1341 # Normalize?
1342 if not relative:
1343 baseline = 1.0 / np.pi
1344 else:
1345 baseline = self.intensity(lat=lat, lon=lon, limbdarken=False)
1346 DeltaI = baseline * intensity
1347 # The integral of the gaussian in cos(Delta theta) over the
1348 # surface of the sphere is sigma * sqrt(2 * pi^3). Combining
1349 # this with the normalization convention of starry (a factor of 4),
1350 # the corresponding spot amplitude is...
1351 amp = sigma * np.sqrt(2 * np.pi ** 3) * DeltaI / 4
1352 if not relative:
1353 amp /= self.amp
1355 # Parse remaining kwargs
1356 sigma, lat, lon = math.cast(sigma, lat, lon)
1358 # Get the Ylm expansion of the spot. Note that yspot[0] is not
1359 # unity, so we'll need to normalize it before setting self._y
1360 yspot = self.ops.expand_spot(
1361 amp, sigma, lat * self._angle_factor, lon * self._angle_factor
1362 )
1363 y_new = self._y + yspot
1364 amp_new = self._amp * y_new[0]
1365 y_new /= y_new[0]
1367 # Update the map and the normalizing amplitude
1368 self._y = y_new
1369 self._amp = amp_new
1371 def minimize(self, oversample=1, ntries=1, return_info=False):
1372 """Find the global minimum of the map intensity.
1374 Args:
1375 oversample (int): Factor by which to oversample the initial
1376 grid on which the brute force search is performed. Default 1.
1377 ntries (int): Number of times the nonlinear minimizer is called.
1378 Default 1.
1379 return_info (bool): Return the info from the minimization call?
1380 Default is False.
1382 Returns:
1383 A tuple of the latitude, longitude, and the value of the intensity \
1384 at the minimum. If ``return_info`` is True, also returns the detailed \
1385 solver information.
1386 """
1387 # Not implemented for spectral
1388 self._no_spectral()
1390 self.ops._minimize.setup(oversample=oversample, ntries=ntries)
1391 lat, lon, I = self.ops.get_minimum(self.y)
1392 if return_info: # pragma: no cover
1393 return (
1394 lat / self._angle_factor,
1395 lon / self._angle_factor,
1396 self._amp * I,
1397 self.ops._minimize.result,
1398 )
1399 else:
1400 return (
1401 lat / self._angle_factor,
1402 lon / self._angle_factor,
1403 self._amp * I,
1404 )
1407class LimbDarkenedBase(object):
1408 """The ``starry`` map class for purely limb-darkened maps.
1410 This class handles light curves of purely limb-darkened objects in
1411 emitted light.
1413 .. note::
1414 Instantiate this class by calling :py:func:`starry.Map` with
1415 ``ydeg`` set to zero and both ``rv`` and ``reflected`` set to False.
1416 """
1418 _ops_class_ = OpsLD
1420 def flux(self, **kwargs):
1421 """
1422 Compute and return the light curve.
1424 Args:
1425 xo (scalar or vector, optional): x coordinate of the occultor
1426 relative to this body in units of this body's radius.
1427 yo (scalar or vector, optional): y coordinate of the occultor
1428 relative to this body in units of this body's radius.
1429 zo (scalar or vector, optional): z coordinate of the occultor
1430 relative to this body in units of this body's radius.
1431 ro (scalar, optional): Radius of the occultor in units of
1432 this body's radius.
1433 """
1434 # Orbital kwargs
1435 theta = kwargs.pop("theta", None)
1436 _, xo, yo, zo, ro = self._get_flux_kwargs(kwargs)
1438 # Check for invalid kwargs
1439 if theta is not None:
1440 # If the user passed in `theta`, make sure a warning is raised
1441 kwargs["theta"] = theta
1442 self._check_kwargs("flux", kwargs)
1444 # Compute & return
1445 return self.amp * self.ops.flux(xo, yo, zo, ro, self._u)
1447 def intensity(self, mu=None, x=None, y=None):
1448 r"""
1449 Compute and return the intensity of the map.
1451 Args:
1452 mu (scalar or vector, optional): the radial parameter :math:`\mu`,
1453 equal to the cosine of the angle between the line of sight and
1454 the normal to the surface. Default is None.
1455 x (scalar or vector, optional): the Cartesian x position on the
1456 surface in units of the body's radius. Default is None.
1457 y (scalar or vector, optional): the Cartesian y position on the
1458 surface in units of the body's radius. Default is None.
1460 .. note::
1461 Users must provide either `mu` **or** `x` and `y`.
1462 """
1463 # Get the Cartesian points
1464 if mu is not None:
1465 mu = math.vectorize(math.cast(mu))
1466 assert (
1467 x is None and y is None
1468 ), "Please provide either `mu` or `x` and `y`, but not both."
1469 else:
1470 assert (
1471 x is not None and y is not None
1472 ), "Please provide either `mu` or `x` and `y`."
1473 x, y = math.vectorize(*math.cast(x, y))
1474 mu = (1 - x ** 2 - y ** 2) ** 0.5
1476 # Compute & return
1477 return self.amp * self.ops.intensity(mu, self._u)
1479 def render(self, res=300):
1480 """Compute and return the intensity of the map on a grid.
1482 Returns an image of shape ``(res, res)``.
1484 Args:
1485 res (int, optional): The resolution of the map in pixels on a
1486 side. Defaults to 300.
1487 """
1488 # Multiple frames?
1489 if self.nw is not None:
1490 animated = True
1491 else:
1492 animated = False
1494 # Compute
1495 image = self.amp * self.ops.render_ld(res, self._u)
1497 # Squeeze?
1498 if animated:
1499 return image
1500 else:
1501 return math.reshape(image, [res, res])
1504class RVBase(object):
1505 """The radial velocity ``starry`` map class.
1507 This class handles velocity-weighted intensities for use in
1508 Rossiter-McLaughlin effect investigations. It has all the same
1509 attributes and methods as :py:class:`starry.maps.YlmBase`, with the
1510 additions and modifications listed below.
1512 .. note::
1513 Instantiate this class by calling :py:func:`starry.Map` with
1514 ``ydeg > 0`` and ``rv`` set to True.
1515 """
1517 _ops_class_ = OpsRV
1519 def reset(self, **kwargs):
1520 self.velocity_unit = kwargs.pop("velocity_unit", units.m / units.s)
1521 self.veq = kwargs.pop("veq", 0.0)
1522 super(RVBase, self).reset(**kwargs)
1524 @property
1525 def velocity_unit(self):
1526 """An ``astropy.units`` unit defining the velocity metric for this map."""
1527 return self._velocity_unit
1529 @velocity_unit.setter
1530 def velocity_unit(self, value):
1531 assert value.physical_type == "speed"
1532 self._velocity_unit = value
1533 self._velocity_factor = value.in_units(units.m / units.s)
1535 @property
1536 def veq(self):
1537 """The equatorial velocity of the body in units of :py:attr:`velocity_unit`.
1539 .. warning::
1540 If this map is associated with a :py:class:`starry.Body`
1541 instance in a Keplerian system, changing the body's
1542 radius and rotation period does not currently affect this
1543 value. The user must explicitly change this value to affect
1544 the map's radial velocity.
1546 """
1547 return self._veq / self._velocity_factor
1549 @veq.setter
1550 def veq(self, value):
1551 self._veq = math.cast(value) * self._velocity_factor
1553 def _unset_RV_filter(self):
1554 f = np.zeros(self.Nf)
1555 f[0] = np.pi
1556 self._f = math.cast(f)
1558 def _set_RV_filter(self):
1559 self._f = self.ops.compute_rv_filter(
1560 self._inc, self._obl, self._veq, self._alpha
1561 )
1563 def rv(self, **kwargs):
1564 """Compute the net radial velocity one would measure from the object.
1566 The radial velocity is computed as the ratio
1568 :math:`\\Delta RV = \\frac{\\int Iv \\mathrm{d}A}{\\int I \\mathrm{d}A}`
1570 where both integrals are taken over the visible portion of the
1571 projected disk. :math:`I` is the intensity field (described by the
1572 spherical harmonic and limb darkening coefficients) and :math:`v`
1573 is the radial velocity field (computed based on the equatorial velocity
1574 of the star, its orientation, etc.)
1576 Args:
1577 xo (scalar or vector, optional): x coordinate of the occultor
1578 relative to this body in units of this body's radius.
1579 yo (scalar or vector, optional): y coordinate of the occultor
1580 relative to this body in units of this body's radius.
1581 zo (scalar or vector, optional): z coordinate of the occultor
1582 relative to this body in units of this body's radius.
1583 ro (scalar, optional): Radius of the occultor in units of
1584 this body's radius.
1585 theta (scalar or vector, optional): Angular phase of the body
1586 in units of :py:attr:`angle_unit`.
1587 """
1588 # Orbital kwargs
1589 theta, xo, yo, zo, ro = self._get_flux_kwargs(kwargs)
1591 # Check for invalid kwargs
1592 self._check_kwargs("rv", kwargs)
1594 # Compute
1595 return self.ops.rv(
1596 theta,
1597 xo,
1598 yo,
1599 zo,
1600 ro,
1601 self._inc,
1602 self._obl,
1603 self._y,
1604 self._u,
1605 self._veq,
1606 self._alpha,
1607 )
1609 def intensity(self, **kwargs):
1610 """
1611 Compute and return the intensity of the map.
1613 Args:
1614 lat (scalar or vector, optional): latitude at which to evaluate
1615 the intensity in units of :py:attr:`angle_unit`.
1616 lon (scalar or vector, optional): longitude at which to evaluate
1617 the intensity in units of :py:attr:`angle_unit`.
1618 rv (bool, optional): If True, computes the velocity-weighted
1619 intensity instead. Defaults to True.
1620 theta (scalar, optional): For differentially rotating maps only,
1621 the angular phase at which to evaluate the intensity.
1622 Default 0.
1623 limbdarken (bool, optional): Apply limb darkening
1624 (only if :py:attr:`udeg` > 0)? Default True.
1626 """
1627 # Compute the velocity-weighted intensity if `rv==True`
1628 rv = kwargs.pop("rv", True)
1629 if rv:
1630 self._set_RV_filter()
1631 res = super(RVBase, self).intensity(**kwargs)
1632 if rv:
1633 self._unset_RV_filter()
1634 return res
1636 def render(self, **kwargs):
1637 """
1638 Compute and return the intensity of the map on a grid.
1640 Returns an image of shape ``(res, res)``, unless ``theta`` is a vector,
1641 in which case returns an array of shape ``(nframes, res, res)``, where
1642 ``nframes`` is the number of values of ``theta``. However, if this is
1643 a spectral map, ``nframes`` is the number of wavelength bins and
1644 ``theta`` must be a scalar.
1646 Args:
1647 res (int, optional): The resolution of the map in pixels on a
1648 side. Defaults to 300.
1649 projection (string, optional): The map projection. Accepted
1650 values are ``ortho``, corresponding to an orthographic
1651 projection (as seen on the sky), and ``rect``, corresponding
1652 to an equirectangular latitude-longitude projection.
1653 Defaults to ``ortho``.
1654 theta (scalar or vector, optional): The map rotation phase in
1655 units of :py:attr:`angle_unit`. If this is a vector, an
1656 animation is generated. Defaults to ``0.0``.
1657 rv (bool, optional): If True, computes the velocity-weighted
1658 intensity instead. Defaults to True.
1659 """
1660 # Render the velocity map if `rv==True`
1661 # Override the `projection` kwarg if we're
1662 # plotting the radial velocity.
1663 rv = kwargs.pop("rv", True)
1664 if rv:
1665 kwargs.pop("projection", None)
1666 self._set_RV_filter()
1667 res = super(RVBase, self).render(**kwargs)
1668 if rv:
1669 self._unset_RV_filter()
1670 return res
1672 def show(self, rv=True, **kwargs):
1673 # Show the velocity map if `rv==True`
1674 # Override some kwargs if we're
1675 # plotting the radial velocity.
1676 if rv:
1677 kwargs.pop("projection", None)
1678 self._set_RV_filter()
1679 kwargs["cmap"] = kwargs.pop("cmap", "RdBu_r")
1680 kwargs["norm"] = kwargs.pop("norm", "rv")
1681 res = super(RVBase, self).show(rv=rv, **kwargs)
1682 if rv:
1683 self._unset_RV_filter()
1684 return res
1687class ReflectedBase(object):
1688 """The reflected light ``starry`` map class.
1690 This class handles light curves and phase curves of objects viewed
1691 in reflected light. It has all the same attributes and methods as
1692 :py:class:`starry.maps.YlmBase`, with the
1693 additions and modifications listed below.
1695 The spherical harmonic coefficients of a map in reflected light are
1696 an expansion of the object's *albedo* (instead of its emissivity, in
1697 the default case).
1699 The illumination source is currently assumed to be a point source for
1700 the purposes of computing the illumination profile on the surface of the
1701 body. However, if the illumination source occults the body, the flux
1702 *is* computed correctly (i.e., the occulting body has a finite radius).
1703 This approximation holds if the distance between the occultor and the source
1704 is large compared to the size of the source. It fails, for example, in the
1705 case of an extremely short-period planet, in which case signficantly more
1706 than half the planet surface is illuminated by the star at any given time.
1707 We plan to account for this effect in the future, so stay tuned.
1709 The ``xo``, ``yo``, and ``zo`` parameters in several of the methods below
1710 specify the position of the illumination source in units of this body's
1711 radius. The flux returned by the :py:meth:`flux` method is normalized such
1712 that when the distance between the occultor and the illumination source is
1713 unity, a uniform unit-amplitude map will emit a flux of unity when viewed
1714 at noon.
1716 This class does not currently support occultations. If an occultation
1717 does occur, a ``ValueError`` will be raised. Support for occultations in
1718 reflected light will be added in an upcoming version, so stay tuned.
1720 .. note::
1721 Instantiate this class by calling
1722 :py:func:`starry.Map` with ``ydeg > 0`` and ``reflected`` set to True.
1723 """
1725 _ops_class_ = OpsReflected
1727 def _get_flux_kwargs(self, kwargs):
1728 xo = kwargs.pop("xo", 0.0)
1729 yo = kwargs.pop("yo", 0.0)
1730 zo = kwargs.pop("zo", 1.0)
1731 ro = kwargs.pop("ro", 0.0)
1732 xs = kwargs.pop("xs", 0.0)
1733 ys = kwargs.pop("ys", 0.0)
1734 zs = kwargs.pop("zs", 1.0)
1735 theta = kwargs.pop("theta", 0.0)
1736 theta, xs, ys, zs, xo, yo, zo = math.vectorize(
1737 theta, xs, ys, zs, xo, yo, zo
1738 )
1739 theta, xs, ys, zs, xo, yo, zo, ro = math.cast(
1740 theta, xs, ys, zs, xo, yo, zo, ro
1741 )
1742 theta *= self._angle_factor
1743 return theta, xs, ys, zs, xo, yo, zo, ro
1745 def design_matrix(self, **kwargs):
1746 r"""
1747 Compute and return the light curve design matrix, :math:`A`.
1749 This matrix is used to compute the flux :math:`f` from a vector of spherical
1750 harmonic coefficients :math:`y` and the map amplitude :math:`a`:
1751 :math:`f = a A y`.
1753 Args:
1754 xs (scalar or vector, optional): x coordinate of the illumination
1755 source relative to this body in units of this body's radius.
1756 ys (scalar or vector, optional): y coordinate of the illumination
1757 source relative to this body in units of this body's radius.
1758 zs (scalar or vector, optional): z coordinate of the illumination
1759 source relative to this body in units of this body's radius.
1760 xo (scalar or vector, optional): x coordinate of the occultor
1761 relative to this body in units of this body's radius.
1762 yo (scalar or vector, optional): y coordinate of the occultor
1763 relative to this body in units of this body's radius.
1764 zo (scalar or vector, optional): z coordinate of the occultor
1765 relative to this body in units of this body's radius.
1766 ro (scalar, optional): Radius of the occultor in units of
1767 this body's radius.
1768 theta (scalar or vector, optional): Angular phase of the body
1769 in units of :py:attr:`angle_unit`.
1771 .. note::
1772 ``starry`` does not yet support occultations in reflected light.
1774 """
1775 # Orbital kwargs
1776 theta, xs, ys, zs, xo, yo, zo, ro = self._get_flux_kwargs(kwargs)
1778 # Check for invalid kwargs
1779 self._check_kwargs("X", kwargs)
1781 # Compute & return
1782 return self.ops.X(
1783 theta,
1784 xs,
1785 ys,
1786 zs,
1787 xo,
1788 yo,
1789 zo,
1790 ro,
1791 self._inc,
1792 self._obl,
1793 self._u,
1794 self._f,
1795 self._alpha,
1796 )
1798 def flux(self, **kwargs):
1799 """
1800 Compute and return the reflected flux from the map.
1802 Args:
1803 xs (scalar or vector, optional): x coordinate of the illumination
1804 source relative to this body in units of this body's radius.
1805 ys (scalar or vector, optional): y coordinate of the illumination
1806 source relative to this body in units of this body's radius.
1807 zs (scalar or vector, optional): z coordinate of the illumination
1808 source relative to this body in units of this body's radius.
1809 xo (scalar or vector, optional): x coordinate of the occultor
1810 relative to this body in units of this body's radius.
1811 yo (scalar or vector, optional): y coordinate of the occultor
1812 relative to this body in units of this body's radius.
1813 zo (scalar or vector, optional): z coordinate of the occultor
1814 relative to this body in units of this body's radius.
1815 ro (scalar, optional): Radius of the occultor in units of
1816 this body's radius.
1817 theta (scalar or vector, optional): Angular phase of the body
1818 in units of :py:attr:`angle_unit`.
1820 .. note::
1821 ``starry`` does not yet support occultations in reflected light.
1823 """
1824 # Orbital kwargs
1825 theta, xs, ys, zs, xo, yo, zo, ro = self._get_flux_kwargs(kwargs)
1827 # Check for invalid kwargs
1828 self._check_kwargs("flux", kwargs)
1830 # Compute & return
1831 return self.amp * self.ops.flux(
1832 theta,
1833 xs,
1834 ys,
1835 zs,
1836 xo,
1837 yo,
1838 zo,
1839 ro,
1840 self._inc,
1841 self._obl,
1842 self._y,
1843 self._u,
1844 self._f,
1845 self._alpha,
1846 )
1848 def intensity(self, lat=0, lon=0, xs=0, ys=0, zs=1, **kwargs):
1849 """
1850 Compute and return the intensity of the map.
1852 Args:
1853 lat (scalar or vector, optional): latitude at which to evaluate
1854 the intensity in units of :py:attr:`angle_unit`.
1855 lon (scalar or vector, optional): longitude at which to evaluate
1856 the intensity in units of :py:attr:`angle_unit`.
1857 xs (scalar or vector, optional): x coordinate of the illumination
1858 source relative to this body in units of this body's radius.
1859 ys (scalar or vector, optional): y coordinate of the illumination
1860 source relative to this body in units of this body's radius.
1861 zs (scalar or vector, optional): z coordinate of the illumination
1862 source relative to this body in units of this body's radius.
1863 theta (scalar, optional): For differentially rotating maps only,
1864 the angular phase at which to evaluate the intensity.
1865 Default 0.
1866 limbdarken (bool, optional): Apply limb darkening
1867 (only if :py:attr:`udeg` > 0)? Default True.
1868 """
1869 # Get the Cartesian points
1870 lat, lon = math.vectorize(*math.cast(lat, lon))
1871 lat *= self._angle_factor
1872 lon *= self._angle_factor
1874 # Get the source position
1875 xs, ys, zs = math.vectorize(*math.cast(xs, ys, zs))
1877 # Get the amplitude
1878 if self.nw is None or config.lazy:
1879 amp = self.amp
1880 else:
1881 # The intensity has shape `(nsurf_pts, nw, nsource_pts)`
1882 # so we must reshape `amp` to take the product correctly
1883 amp = self.amp[np.newaxis, :, np.newaxis]
1885 # If differentially rotating, allow a `theta` keyword
1886 if self.drorder > 0:
1887 alpha_theta = math.cast(kwargs.get("theta", 0.0)) * self.alpha
1888 alpha_theta *= self._angle_factor
1889 else:
1890 alpha_theta = math.cast(0.0)
1891 self._check_kwargs("intensity", kwargs)
1893 # If limb-darkened, allow a `limbdarken` keyword
1894 if self.udeg > 0:
1895 ld = np.array(True)
1896 else:
1897 ld = np.array(False)
1899 # Compute & return
1900 return amp * self.ops.intensity(
1901 lat, lon, self._y, self._u, self._f, xs, ys, zs, alpha_theta, ld
1902 )
1904 def render(
1905 self,
1906 res=300,
1907 projection="ortho",
1908 illuminate=True,
1909 theta=0.0,
1910 xs=0,
1911 ys=0,
1912 zs=1,
1913 ):
1914 """
1915 Compute and return the intensity of the map on a grid.
1917 Returns an image of shape ``(res, res)``, unless ``theta`` is a vector,
1918 in which case returns an array of shape ``(nframes, res, res)``, where
1919 ``nframes`` is the number of values of ``theta``. However, if this is
1920 a spectral map, ``nframes`` is the number of wavelength bins and
1921 ``theta`` must be a scalar.
1923 Args:
1924 res (int, optional): The resolution of the map in pixels on a
1925 side. Defaults to 300.
1926 projection (string, optional): The map projection. Accepted
1927 values are ``ortho``, corresponding to an orthographic
1928 projection (as seen on the sky), and ``rect``, corresponding
1929 to an equirectangular latitude-longitude projection.
1930 Defaults to ``ortho``.
1931 illuminate (bool, optional): Illuminate the map? Default is True.
1932 theta (scalar or vector, optional): The map rotation phase in
1933 units of :py:attr:`angle_unit`. If this is a vector, an
1934 animation is generated. Defaults to ``0.0``.
1935 xs (scalar or vector, optional): x coordinate of the illumination
1936 source relative to this body in units of this body's radius.
1937 ys (scalar or vector, optional): y coordinate of the illumination
1938 source relative to this body in units of this body's radius.
1939 zs (scalar or vector, optional): z coordinate of the illumination
1940 source relative to this body in units of this body's radius.
1942 """
1943 # Multiple frames?
1944 if self.nw is not None:
1945 animated = True
1946 else:
1947 if config.lazy:
1948 animated = hasattr(theta, "ndim") and theta.ndim > 0
1949 else:
1950 animated = hasattr(theta, "__len__")
1952 # Convert stuff as needed
1953 projection = get_projection(projection)
1954 theta = math.cast(theta) * self._angle_factor
1955 xs = math.cast(xs)
1956 ys = math.cast(ys)
1957 zs = math.cast(zs)
1958 theta, xs, ys, zs = math.vectorize(theta, xs, ys, zs)
1959 illuminate = int(illuminate)
1961 # Compute
1962 if self.nw is None or config.lazy:
1963 amp = self.amp
1964 else:
1965 # The intensity has shape `(nw, res, res)`
1966 # so we must reshape `amp` to take the product correctly
1967 amp = self.amp[:, np.newaxis, np.newaxis]
1969 image = amp * self.ops.render(
1970 res,
1971 projection,
1972 illuminate,
1973 theta,
1974 self._inc,
1975 self._obl,
1976 self._y,
1977 self._u,
1978 self._f,
1979 self._alpha,
1980 xs,
1981 ys,
1982 zs,
1983 )
1985 # Squeeze?
1986 if animated:
1987 return image
1988 else:
1989 return math.reshape(image, [res, res])
1991 def show(self, **kwargs):
1992 # We need to evaluate the variables so we can plot the map!
1993 if config.lazy and kwargs.get("image", None) is None:
1995 # Get kwargs
1996 res = kwargs.pop("res", 300)
1997 projection = get_projection(kwargs.get("projection", "ortho"))
1998 theta = math.cast(kwargs.pop("theta", 0.0)) * self._angle_factor
1999 xs = math.cast(kwargs.pop("xs", 0))
2000 ys = math.cast(kwargs.pop("ys", 0))
2001 zs = math.cast(kwargs.pop("zs", 1))
2002 theta, xs, ys, zs = math.vectorize(theta, xs, ys, zs)
2003 illuminate = int(kwargs.pop("illuminate", True))
2005 # Evaluate the variables
2006 theta = theta.eval()
2007 xs = xs.eval()
2008 ys = ys.eval()
2009 zs = zs.eval()
2010 inc = self._inc.eval()
2011 obl = self._obl.eval()
2012 y = self._y.eval()
2013 u = self._u.eval()
2014 f = self._f.eval()
2015 alpha = self._alpha.eval()
2017 # Explicitly call the compiled version of `render`
2018 kwargs["image"] = self.ops.render(
2019 res,
2020 projection,
2021 illuminate,
2022 theta,
2023 inc,
2024 obl,
2025 y,
2026 u,
2027 f,
2028 alpha,
2029 xs,
2030 ys,
2031 zs,
2032 )
2033 kwargs["theta"] = theta / self._angle_factor
2035 return super(ReflectedBase, self).show(**kwargs)
2038def Map(
2039 ydeg=0, udeg=0, drorder=0, nw=None, rv=False, reflected=False, **kwargs
2040):
2041 """A generic ``starry`` surface map.
2043 This function is a class factory that returns either
2044 a :doc:`spherical harmonic map <SphericalHarmonicMap>`,
2045 a :doc:`limb darkened map <LimbDarkenedMap>`,
2046 a :doc:`radial velocity map <RadialVelocityMap>`, or
2047 a :doc:`reflected light map <ReflectedLightMap>`,
2048 depending on the arguments provided by the user. The default is
2049 a :doc:`spherical harmonic map <SphericalHarmonicMap>`. If ``rv`` is True,
2050 instantiates a :doc:`radial velocity map <RadialVelocityMap>` map, and
2051 if ``reflected`` is True, instantiates a :doc:`reflected light map
2052 <ReflectedLightMap>`. Otherwise, if ``ydeg`` is zero, instantiates a
2053 :doc:`limb darkened map <LimbDarkenedMap>`.
2055 Args:
2056 ydeg (int, optional): Degree of the spherical harmonic map.
2057 Defaults to 0.
2058 udeg (int, optional): Degree of the limb darkening filter.
2059 Defaults to 0.
2060 drorder (int, optional): Order of the differential rotation
2061 approximation. Defaults to 0.
2062 nw (int, optional): Number of wavelength bins. Defaults to None
2063 (for monochromatic light curves).
2064 rv (bool, optional): If True, enable computation of radial velocities
2065 for modeling the Rossiter-McLaughlin effect. Defaults to False.
2066 reflected (bool, optional): If True, models light curves in reflected
2067 light. Defaults to False.
2068 """
2069 # Check args
2070 ydeg = int(ydeg)
2071 assert ydeg >= 0, "Keyword `ydeg` must be positive."
2072 udeg = int(udeg)
2073 assert udeg >= 0, "Keyword `udeg` must be positive."
2074 if nw is not None:
2075 nw = int(nw)
2076 assert nw > 0, "Number of wavelength bins must be positive."
2077 drorder = int(drorder)
2078 assert (drorder >= 0) and (
2079 drorder <= 3
2080 ), "Differential rotation orders above 3 are not supported."
2081 if drorder > 0:
2082 assert ydeg > 0, "Differential rotation requires `ydeg` >= 1."
2084 # TODO: phase this next warning out
2085 logger.warning(
2086 "Differential rotation is still an experimental feature. "
2087 "Use it with care."
2088 )
2090 Ddeg = (4 * drorder + 1) * ydeg
2091 if Ddeg >= 50:
2092 logger.warning(
2093 "The degree of the differential rotation operator "
2094 "is currently {0}, ".format(Ddeg)
2095 + "which will likely cause the code to run very slowly. "
2096 "Consider decreasing the degree of the map or the order "
2097 "of differential rotation."
2098 )
2100 # Limb-darkened?
2101 if (ydeg == 0) and (rv is False) and (reflected is False):
2103 # TODO: Add support for wavelength-dependent limb darkening
2104 if nw is not None:
2105 raise NotImplementedError(
2106 "Multi-wavelength limb-darkened maps are not yet supported."
2107 )
2109 Bases = (LimbDarkenedBase, MapBase)
2110 else:
2111 Bases = (YlmBase, MapBase)
2113 # Radial velocity / reflected light?
2114 if rv:
2115 Bases = (RVBase,) + Bases
2116 fdeg = 3
2117 elif reflected:
2118 Bases = (ReflectedBase,) + Bases
2119 fdeg = 1
2120 else:
2121 fdeg = 0
2123 # Ensure we're not doing both
2124 if rv and reflected:
2125 raise NotImplementedError(
2126 "Radial velocity maps not implemented in reflected light."
2127 )
2129 # Construct the class
2130 class Map(*Bases):
2132 # Tags
2133 __props__ = dict(
2134 limbdarkened=LimbDarkenedBase in Bases,
2135 reflected=ReflectedBase in Bases,
2136 rv=RVBase in Bases,
2137 spectral=nw is not None,
2138 differential_rotation=drorder > 0,
2139 )
2141 def __init__(self, *args, **kwargs):
2142 # Once a map has been instantiated, no changes
2143 # to the config are allowed.
2144 config.freeze()
2145 super(Map, self).__init__(*args, **kwargs)
2147 return Map(ydeg, udeg, fdeg, drorder, nw, **kwargs)