Merge pull request #160 from discsim/log_normal

Log-normal fits

HISTORY.rst

@@ -3,6 +3,42 @@
 Changelog
 +++++++++
 
+v.1.2.0
++++++++
+*Introduction of log-normal fits, large amount of code refactoring to support both 'Gaussian' and 'LogNormal' methods*
+
+- default_parameters.json, parameter_descriptions.json:
+    - Adds parameters: 'rescale_flux', 'method'
+- debris_fitters.py
+    - Adds support for fitting optically thin but geometrically thick disks with a known Gaussian scale height.
+- filter.py:
+    - New module that now hosts the routine for optimizing for power spectrum priors, 'CriticalFilter'
+- fit.py:
+    - Adds ability to run either a standard or log-normal fit
+- geometry.py:
+    - Adds routine to rescale the total flux according to the source geometry, 'rescale_total_flux'
+    - Adds 'rescale_factor' property to 'SourceGeometry'
+    - Adds the option to keep all three Fourier components (u,v,w) when deprojecting
+- minimizer.py:
+    - New module that hosts routines for solving non-linear minimization problems: 'BaseLineSearch', 'LineSearch', 'MinimizeNewton'
+- radial_fitters.py:
+    - Code refactoring:
+        * Removes '_HankelRegressor' class
+        * Adds 'FrankGaussianFit' and 'FrankLogNormalFit' classes
+        * Adds 'FrankRadialFit' class
+        * Moves some core functionalities to 'frank.filter' and 'frank.statistical_models'
+    - Adds 'MAP', 'I', 'info' properties, and 'assume_optically_thick' parameter, to 'FrankRadialFit'
+- statistical_models.py:
+   - New module that now hosts 'GaussianModel' class (containing much of the functionality of the now deprecated '_HankelRegressor'), and adds analogous 'LogNormalMAPModel' class
+   - Adds a VisibilityMapping that abstracts the mapping between the brightness profile and visibilities. Handles optically thick (default), optically thin, and debris disk models.
+- tests.py:
+    - Adds test for a log-normal fit
+- Docs:
+    - Updates API
+- Miscellaneous:
+    - Minor bug and typo fixes
+
+
 v.1.1.0
 +++++++
 

docs/py_API.rst

@@ -25,14 +25,27 @@ the deprojected visibilities.
 .. autoclass:: frank.radial_fitters.FrankFitter
   :members: fit, MAP_solution, MAP_spectrum, MAP_spectrum_covariance, r, Rmax, q, Qmax, size, geometry
 
-.. autoclass:: frank.radial_fitters._HankelRegressor
-  :members: mean, covariance, power_spectrum, r, q, Rmax, Qmax, size, geometry, predict, predict_deprojected, log_likelihood, solve_non_negative
+.. autoclass:: frank.radial_fitters.FrankRadialFit
+  :members: I, MAP, r, q, Rmax, Qmax, size, geometry, predict, predict_deprojected
+
+.. autoclass:: frank.radial_fitters.FrankGaussianFit
+  :members: I, mean, MAP, covariance, power_spectrum, r, q, Rmax, Qmax, size, geometry, predict, predict_deprojected, log_likelihood, solve_non_negative
+
+.. autoclass:: frank.radial_fitters.FrankLogNormalFit
+  :members: I, MAP, covariance, power_spectrum, r, q, Rmax, Qmax, size, geometry, predict, predict_deprojected, log_likelihood
+
+.. autoclass:: frank.debris_fitters.FrankDebrisFitter
+  :members: fit, MAP_solution, MAP_spectrum, MAP_spectrum_covariance, r, Rmax, q, Qmax, size, geometry
+
 
 Utility functions and classes
 -----------------------------
 
 These are some useful functions and classes for various aspects of fitting and analysis.
 
+.. autoclass:: frank.hankel.DiscreteHankelTransform
+  :members: r, Rmax, q, Qmax, size, order, transform, coefficients
+
 .. autofunction:: frank.utilities.arcsec_baseline
 
 .. autofunction:: frank.utilities.convolve_profile
@@ -43,6 +56,8 @@ These are some useful functions and classes for various aspects of fitting and a
 
 .. autofunction:: frank.utilities.estimate_weights
 
+.. autofunction:: frank.utilities.make_image
+
 .. autofunction:: frank.utilities.normalize_uv
 
 .. autofunction:: frank.utilities.sweep_profile

docs/tutorials/model_limitations.rst

@@ -19,7 +19,7 @@ is useful to assess and potentially suppress this behavior.
 Allowed regions of negative brightness
 --------------------------------------
 The fitted brightness profile can have negative regions corresponding to spatial scales un- or underconstrained by the visibilities.
-You can perform a fit in which the solution is forced to be positive (given the maximum a posteriori power spectrum) by using the `solve_non_negative <../py_API.rst#frank.radial_fitters._HankelRegressor.solve_non_negative>`_ method provided by the solution returned by `FrankFitter` (if running frank from the terminal, set `hyperparameters : nonnegative` to `true` in your parameter file).
+You can perform a fit in which the solution is forced to be positive (given the maximum a posteriori power spectrum) by using the `solve_non_negative <../py_API.rst#frank.radial_fitters.FrankGaussianFit.solve_non_negative>`_ method provided by the solution returned by `FrankFitter` (if running frank from the terminal, set `hyperparameters : nonnegative` to `true` in your parameter file).
 In tests we've seen the effect on the recovered brightness profile to typically be localized to the regions of negative flux,
 with otherwise minor differences.
 Since enforcing the profile to be non-negative requires some extrapolation beyond the data's longest baseline,

frank/__init__.py

@@ -16,13 +16,14 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>
 #
-__version__ = "1.1.0"
+__version__ = "1.2.0"
 
 from frank import constants
 from frank import geometry
 from frank import hankel
 from frank import io
 from frank import radial_fitters
+from frank import debris_fitters
 from frank import utilities
 
 def enable_logging(log_file=None):

frank/debris_fitters.py

@@ -0,0 +1,153 @@
+# Frankenstein: 1D disc brightness profile reconstruction from Fourier data
+# using non-parametric Gaussian Processes
+#
+# Copyright (C) 2019-2020  R. Booth, J. Jennings, M. Tazzari
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>
+#
+"""This module contains methods for fitting a radial brightness profile to a set
+  of deprojected visibities. Routines in this file assume that the emission is
+  optically thin with a Gaussian vertical structure.
+"""
+import abc
+from collections import defaultdict
+import logging
+import numpy as np
+
+from frank.radial_fitters import FourierBesselFitter, FrankFitter
+
+class FourierBesselDebrisFitter(FourierBesselFitter):
+    """
+    Fourier-Bessel series optically-thin model for fitting visibilities.
+
+    The brightness model is :math:`I(R, z) = I(R) exp(-z^2/2H(R)^2)`, where
+    :math:`H(R)` is the (known) scale-height.
+
+    Parameters
+    ----------
+    Rmax : float, unit = arcsec
+        Radius of support for the functions to transform, i.e.,
+            f(r) = 0 for R >= Rmax
+    N : int
+        Number of collocation points
+    geometry : SourceGeometry object
+        Geometry used to deproject the visibilities before fitting
+    scale_height : function R --> H
+        Specifies the thickness of disc as a function of radius. Both
+        units should be in arcsec.
+    nu : int, default = 0
+        Order of the discrete Hankel transform (DHT)
+    block_data : bool, default = True
+        Large temporary matrices are needed to set up the data. If block_data
+        is True, we avoid this, limiting the memory requirement to block_size
+        elements.
+    block_size : int, default = 10**5
+        Size of the matrices if blocking is used
+    verbose : bool, default = False
+        Whether to print notification messages
+    """
+
+    def __init__(self, Rmax, N, geometry, scale_height, nu=0, block_data=True,
+                 block_size=10 ** 5, verbose=True):
+
+        # All functionality is provided by the base class. 
+        # FourierBesselDebrisFitter is just a sub-set of FourierBesselFitter
+        super(FourierBesselDebrisFitter, self).__init__(
+            Rmax, N, geometry, nu=nu, block_data=block_data,
+            assume_optically_thick=False, scale_height=scale_height,
+            block_size=block_size, verbose=verbose
+        )
+
+class FrankDebrisFitter(FrankFitter):
+    """
+    Fit a Gaussian process model using the Discrete Hankel Transform of
+    Baddour & Chouinard (2015).
+
+    The brightness model is :math:`I(R, z) = I(R) exp(-z^2/2H(R)^2)`, where
+    :math:`H(R)` is the (known) scale-height.
+
+    The GP model is based upon Oppermann et al. (2013), which use a maximum
+    a posteriori estimate for the power spectrum as the GP prior for the
+    real-space coefficients.
+
+    Parameters
+    ----------
+    Rmax : float, unit = arcsec
+        Radius of support for the functions to transform, i.e., f(r) = 0 for
+        R >= Rmax.
+    N : int
+        Number of collaction points
+    geometry : SourceGeometry object
+        Geometry used to deproject the visibilities before fitting
+    scale_height : function R --> H
+        Specifies the thickness of disc as a function of radius. Both
+        units should be in arcsec.
+    nu : int, default = 0
+        Order of the discrete Hankel transform, given by J_nu(r)
+    block_data : bool, default = True
+        Large temporary matrices are needed to set up the data. If block_data
+        is True, we avoid this, limiting the memory requirement to block_size
+        elements
+    block_size : int, default = 10**5
+        Size of the matrices if blocking is used
+    alpha : float >= 1, default = 1.05
+        Order parameter of the inverse gamma prior for the power spectrum
+        coefficients
+    p_0 : float >= 0, default = None, unit=Jy^2
+        Scale parameter of the inverse gamma prior for the power spectrum
+        coefficients. If not provided p_0 = 1e-15 (method="Normal") or 
+        1e-35 (method="LogNormal") will be used.
+    weights_smooth : float >= 0, default = 1e-4
+        Spectral smoothness prior parameter. Zero is no smoothness prior
+    tol : float > 0, default = 1e-3
+        Tolerence for convergence of the power spectrum iteration
+    method : string, default="Normal"
+        Model used for the brightness reconstrution. This must be one of
+        "Normal" of "LogNormal".
+    I_scale : float, default = 1e5, unit= Jy/Sr
+        Brightness scale. Only used in the LogNormal model. Notet the 
+        LogNormal model produces I(Rmax) =  I_scale.
+    max_iter: int, default = 2000
+        Maximum number of fit iterations
+    check_qbounds: bool, default = True
+        Whether to check if the first (last) collocation point is smaller
+        (larger) than the shortest (longest) deprojected baseline in the dataset
+    store_iteration_diagnostics: bool, default = False
+        Whether to store the power spectrum parameters and brightness profile
+        for each fit iteration
+    verbose:
+        Whether to print notification messages
+
+    References
+    ----------
+        Baddour & Chouinard (2015)
+            DOI: https://doi.org/10.1364/JOSAA.32.000611
+        Oppermann et al. (2013)
+            DOI:  https://doi.org/10.1103/PhysRevE.87.032136
+    """
+
+    def __init__(self, Rmax, N, geometry, scale_height, nu=0, block_data=True,
+                 block_size=10 ** 5, alpha=1.05, p_0=None, weights_smooth=1e-4,
+                 tol=1e-3, method='Normal', I_scale=1e5, max_iter=2000, check_qbounds=True,
+                 store_iteration_diagnostics=False, verbose=True):
+
+        # All functionality is provided by the base class. FrankDebrisFitter is just a 
+        # sub-set of FrankFitter
+        super(FrankDebrisFitter, self).__init__(
+            Rmax, N, geometry, nu=nu, block_data=block_data,
+            block_size=block_size, alpha=alpha, p_0=p_0, weights_smooth=weights_smooth,
+            tol=tol, method=method, I_scale=I_scale, max_iter=max_iter, 
+            check_qbounds=check_qbounds, store_iteration_diagnostics=store_iteration_diagnostics, 
+            assume_optically_thick=False, scale_height=scale_height, verbose=verbose)

frank/default_parameters.json

@@ -25,7 +25,8 @@
     "inc"              : 0.0,
     "pa"               : 0.0,
     "dra"              : 0.0,
-    "ddec"             : 0.0
+    "ddec"             : 0.0,
+    "rescale_flux"     : true
   },
 
   "hyperparameters" : {
@@ -36,7 +37,8 @@
     "wsmooth"          : 1e-4,
     "iter_tol"         : 1e-3,
     "max_iter"         : 2000,
-    "nonnegative"      : false
+    "nonnegative"      : false,
+    "method"           : "Normal"
   },
 
   "plotting" : {