modernising dynamite

.gitignore

@@ -0,0 +1,141 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+pip-wheel-metadata/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+.python-version
+
+# pipenv
+Pipfile.lock
+
+# PEP 582
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/
+
+# IDEs
+.vscode/
+.idea/
+*.swp
+*.swo
+*~
+
+# OS
+.DS_Store
+Thumbs.db
+
+# Version control
+.git/
+.gitignore
+
+# Dynamite specific
+_version.py

README.md

@@ -11,15 +11,13 @@ It follows the method presented in Pinte et al. 2018, with various improvements,
 ```
 git clone https://github.com/cpinte/dynamite.git
 cd dynamite
-python3 setup.py install
+pip install .
 ```
 
-If you don't have the `sudo` rights, use `python3 setup.py install --user`.
-
 To install in developer mode: (i.e. using symlinks to point directly
 at this directory, so that code changes here are immediately available
 without needing to repeat the above step):
 
 ```
- python3 setup.py develop
+pip install -e .
 ```

pyproject.toml

@@ -0,0 +1,47 @@
+[build-system]
+requires = ["setuptools>=70"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "dynamite"
+version = "0.1.0"
+description = "DYNamical Analysis and MultIscale Tomography of line Emission"
+readme = "README.md"
+requires-python = ">=3.8"
+license = {text = "MIT"}
+authors = [
+    {name = "Christophe Pinte", email = "christophe.pinte@univ-grenoble-alpes.fr"}
+]
+keywords = ["astronomy", "protoplanetary-disks", "kinematics", "tomography", "ALMA"]
+classifiers = [
+    "Programming Language :: Python :: 3",
+    "Programming Language :: Python :: 3.8",
+    "Programming Language :: Python :: 3.9",
+    "Programming Language :: Python :: 3.10",
+    "Programming Language :: Python :: 3.11",
+    "Programming Language :: Python :: 3.12",
+    "Programming Language :: Python :: 3.13",
+    "Programming Language :: Python :: 3.14",
+    "Intended Audience :: Science/Research",
+    "Topic :: Scientific/Engineering :: Astronomy",
+    "Topic :: Scientific/Engineering :: Physics",
+    "License :: OSI Approved :: MIT License",
+    "Operating System :: OS Independent",
+]
+dependencies = [
+    "numpy",
+    "astropy",
+    "scipy",
+    "matplotlib",
+    "tqdm",
+    "casa_cube",
+    "celerite2"
+]
+
+[project.urls]
+Homepage = "https://github.com/cpinte/dynamite"
+Repository = "https://github.com/cpinte/dynamite.git"
+
+[tool.setuptools]
+package-dir = {"" = "src"}
+packages = ["dynamite"]

src/dynamite/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "0.1"
+__version__ = "0.1.0"
 
 from .measure_height import *
 from .toy_model import *

src/dynamite/measure_height.py

@@ -14,22 +14,18 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import scipy.constants as sc
-from alive_progress import alive_bar
+from tqdm import tqdm
 from numpy import ndarray
 from scipy import ndimage
 from scipy.ndimage import rotate
 from scipy.optimize import curve_fit, minimize
-from scipy.optimize import minimize
 from scipy.interpolate import interp1d
 from scipy.signal import find_peaks
 from scipy.stats import binned_statistic
-from astropy.convolution import Gaussian2DKernel, convolve, convolve_fft
-import celerite
-from celerite import terms
+from astropy.convolution import Gaussian2DKernel, convolve_fft
 from scipy import signal
 
-
-import casa_cube
+from casa_cube import Cube
 
 sigma_to_FWHM = 2.0 * np.sqrt(2.0 * np.log(2))
 FWHM_to_sigma = 1.0 / sigma_to_FWHM
@@ -108,7 +104,7 @@ def __init__(self,
 
         if isinstance(cube,str):
             print("Reading cube: "+cube)
-            cube = casa_cube.Cube(cube)
+            cube = Cube(cube)
 
         # Truncating the cube is velocity if needed
         if vmin is not None:
@@ -543,10 +539,8 @@ def _rotate_cube(self):
         self.x_star_rot = center[0] + dx * np.cos(angle) + dy * np.sin(angle)
         self.y_star_rot = center[1] - dx * np.sin(angle) + dy * np.cos(angle)
 
-        with alive_bar(int(self.iv_max-self.iv_min), title="Rotating cube") as bar:
-            for iv in range(self.iv_min,self.iv_max):
-                self.cube.image[iv,:,:] = np.array(rotate(self.cube.image[iv,:,:], self.PA - self.inc_sign * 90.0, reshape=False))
-                bar()
+        for iv in tqdm(range(self.iv_min, self.iv_max + 1), desc="Rotating cube"):
+            self.cube.image[iv,:,:] = np.array(rotate(self.cube.image[iv,:,:], self.PA - self.inc_sign * 90.0, reshape=False))
 
         return
 
@@ -650,12 +644,10 @@ def _extract_isovelocity(self):
         self.I = np.zeros([ns,nv,nx,2])
 
         # Loop over the channels
-        with alive_bar(int(self.iv_max-self.iv_min), title="Extracting isovelocity curves") as bar:
-            for iv in range(self.iv_min,self.iv_max):
-                for iscale in range(self.n_scales):
-                    self._extract_isovelocity_1channel(iv,iscale)
-                    #self._refine_isovelocity_1channel(iv,iscale=iscale)
-                    bar()
+        for iv in tqdm(range(self.iv_min, self.iv_max + 1), desc="Extracting isovelocity curves"):
+            for iscale in range(self.n_scales):
+                self._extract_isovelocity_1channel(iv,iscale)
+                #self._refine_isovelocity_1channel(iv,iscale=iscale)
 
         #--  Additional spectral filtering to clean the data ??
 
@@ -708,11 +700,9 @@ def _make_multiscale_cube(self):
                 self.multiscale_std[iscale] = np.nanstd([im,im1])
 
                 # Make the multiscale cube
-                with alive_bar(int(self.iv_max-self.iv_min), title="Making multi-scale cube: scale #"+str(iscale)) as bar:
-                    for iv in range(self.iv_min,self.iv_max):
-                        im = self.rotated_images[0,iv-self.iv_min,:,:]
-                        self.rotated_images[iscale,iv-self.iv_min,:,:] = convolve_fft(im, beam)
-                        bar()
+                for iv in tqdm(range(self.iv_min, self.iv_max + 1), desc=f"Making multi-scale cube: scale #{iscale}"):
+                    im = self.rotated_images[0,iv-self.iv_min,:,:]
+                    self.rotated_images[iscale,iv-self.iv_min,:,:] = convolve_fft(im, beam)
 
                 self.multiscale_bmaj[iscale] = bmaj
                 self.multiscale_bmin[iscale] = bmin
@@ -1099,7 +1089,7 @@ def find_i(self,num=0):
             metric[i] = self.h_std
 
         plt.plot(inc_array,metric, color="red", markersize=1)
-        plt.xlabel("Inclination ($^\mathrm{o}$)")
+        plt.xlabel(r"Inclination ($^\mathrm{o}$)")
         plt.ylabel("altitude dispersion")
 
         # T dispersion
@@ -1116,7 +1106,7 @@ def find_i(self,num=0):
             metric[i] = self.T_std
 
         plt.plot(inc_array,metric, color="red", markersize=1)
-        plt.xlabel("Inclination ($^\mathrm{o}$)")
+        plt.xlabel(r"Inclination ($^\mathrm{o}$)")
         plt.ylabel("T dispersion")
 
         # Velocity dispersion
@@ -1145,7 +1135,7 @@ def find_i(self,num=0):
             metric[i] = self.v_std
 
         plt.plot(inc_array,metric, color="blue", markersize=1)
-        plt.xlabel("Inclination ($^\mathrm{o}$)")
+        plt.xlabel(r"Inclination ($^\mathrm{o}$)")
         plt.ylabel("Velocity dispersion")
         self.inc = inc_array[np.nanargmin(metric)]
         print("Best fit for inclination =", self.inc, "deg")
@@ -1382,7 +1372,7 @@ def plot_channel(self, iv, iscale=0, radius=3.0, ax=None, clear=True):
         #    im = np.array(rotate(im, self.PA - self.inc_sign * 90.0, reshape=False))
 
         ax.imshow(im, origin="lower", cmap='binary_r')
-        ax.set_title(r'v='+"{:.2f}".format(cube.velocity[iv])+' , $\Delta$v='+"{:.2f}".format(cube.velocity[iv] - self.v_syst)+' , id:'+str(iv), color='k')
+        ax.set_title(r'v='+"{:.2f}".format(cube.velocity[iv])+r' , $\Delta$v='+"{:.2f}".format(cube.velocity[iv] - self.v_syst)+' , id:'+str(iv), color='k')
 
         if n_surf[iscale,iv]:
             ax.plot(x[iscale,iv,:n_surf[iscale,iv]],y[iscale,iv,:n_surf[iscale,iv],0],"o",color="red",markersize=1)
@@ -1768,6 +1758,8 @@ def to_mcfost(self, planet_r=0., planet_PA=0.):
 
 
     def fit_surface_height_gp(self):
+        import celerite2
+        from celerite2 import terms
 
         x = np.array(self.r.ravel().compressed())
         y = self.h.ravel().compressed()
@@ -1778,13 +1770,16 @@ def fit_surface_height_gp(self):
         y=y[order]
         yerr=yerr[order]
 
-        # Set up the GP model
-        k0 = terms.JitterTerm(log_sigma=np.log(np.var(y)))
-        k1 = terms.RealTerm(log_a=np.log(np.var(y)), log_c=np.log(50))
-        k2 = terms.RealTerm(log_a=np.log(np.var(y)), log_c=np.log(5))
-        kernel = k0+k1+k2
-        kernel = k0+k2
-        gp = celerite.GP(kernel, fit_mean=False)
+        # set up the GP model
+        # in celerite2, we use GaussianProcess with a mean function
+        # the JitterTerm and RealTerm are replaced with appropriate terms
+        sigma = np.sqrt(np.var(y))
+        k0 = terms.SHOTerm(sigma=sigma, rho=1.0, Q=1/np.sqrt(2))
+        k1 = terms.SHOTerm(sigma=sigma, rho=50.0, Q=1/np.sqrt(2))
+        k2 = terms.SHOTerm(sigma=sigma, rho=5.0, Q=1/np.sqrt(2))
+        kernel = k0 + k1 + k2
+        kernel = k0 + k2
+        gp = celerite2.GaussianProcess(kernel, mean=0.0)
 
         # Define a cost function
         def neg_log_like(params, y, gp):
@@ -1795,7 +1790,7 @@ def grad_neg_log_like(params, y, gp):
             gp.set_parameter_vector(params)
             return -gp.grad_log_likelihood(y)[1]
 
-        gp.compute(x, yerr)
+        gp.compute(x, yerr=yerr)
 
         initial_params = gp.get_parameter_vector()
         bounds = gp.get_parameter_bounds()
@@ -1807,7 +1802,7 @@ def grad_neg_log_like(params, y, gp):
         # Make the maximum likelihood prediction
         t = np.linspace(np.min(x), np.max(x), 500)
 
-        mu, var = gp.predict(y, t, return_var=True)
+        mu, var = gp.predict(y, t=t, return_var=True)
         std = np.sqrt(var)
 
         return t, mu, std