diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..3e25aac --- /dev/null +++ b/.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 \ No newline at end of file diff --git a/README.md b/README.md index 07c94c4..efd109a 100644 --- a/README.md +++ b/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 . ``` diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..5ec7731 --- /dev/null +++ b/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"] \ No newline at end of file diff --git a/requirements.txt b/requirements.txt deleted file mode 100644 index 7961e23..0000000 --- a/requirements.txt +++ /dev/null @@ -1,5 +0,0 @@ -numpy -astropy -scipy -matplotlib -alive-progress diff --git a/setup.py b/setup.py deleted file mode 100644 index d187e44..0000000 --- a/setup.py +++ /dev/null @@ -1,9 +0,0 @@ -from setuptools import setup - -setup(name='dynamite', - description='DYNamical Analysis and MultIscale Tomography of line Emission', - url='http://github.com/cpinte/dynamite', - author='Christophe Pinte', - license='MIT', - packages=['dynamite'], - zip_safe=False) diff --git a/src/dynamite/__init__.py b/src/dynamite/__init__.py index 12b197a..18c6970 100644 --- a/src/dynamite/__init__.py +++ b/src/dynamite/__init__.py @@ -1,4 +1,4 @@ -__version__ = "0.1" +__version__ = "0.1.0" from .measure_height import * from .toy_model import * diff --git a/src/dynamite/measure_height.py b/src/dynamite/measure_height.py index b55ff99..873bfe4 100644 --- a/src/dynamite/measure_height.py +++ b/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