sciapy.regress.models_cel.ProxyModel.log_prior() - Code Metrics - Inspection of "deps: Use `regressproxy` package" - st-bender/sciapy - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Push — master ( af9f8c...efc6a8 )

by Stefan

created 2022-10-20 08:21 UTC

sciapy.regress.models_cel.ProxyModel.log_prior() A

↳ Parent: sciapy.regress.models_cel

Complexity

Conditions

Size

Total Lines	6
Code Lines	6

Duplication

Lines	0
Ratio	0 %

Code Coverage

Tests	1
CRAP Score	7.608

Importance

Changes

Metric	Value
cc	3
eloc	6
nop	1
dl	0
loc	6
ccs	1
cts	5
cp	0.2
crap	7.608
rs	10
c	0
b	0
f	0

# -*- coding: utf-8 -*-
# vim:fileencoding=utf-8
#
# Copyright (c) 2017-2018 Stefan Bender
#
# This module is part of sciapy.
# sciapy is free software: you can redistribute it or modify
# it under the terms of the GNU General Public License as published
# by the Free Software Foundation, version 2.
# See accompanying LICENSE file or http://www.gnu.org/licenses/gpl-2.0.html.
"""SCIAMACHY regression models (celerite version)

Model classes for SCIAMACHY data regression fits using the
:mod:`celerite` [#]_ modelling protocol.

.. [#] https://celerite.readthedocs.io
"""
from __future__ import absolute_import, division, print_function

import numpy as np
from scipy.interpolate import interp1d

from celerite.modeling import Model, ModelSet, ConstantModel

__all__ = ["ConstantModel",
		"HarmonicModelCosineSine", "HarmonicModelAmpPhase",
		"ProxyModel", "TraceGasModelSet",
		"setup_proxy_model_with_bounds", "trace_gas_model"]


class HarmonicModelCosineSine(Model):
	"""Model for harmonic terms

	Models harmonic terms using a cosine and sine part.
	The total amplitude and phase can be inferred from that.

	Parameters
	----------
	freq : float
		The frequency in years^-1
	cos : float
		The amplitude of the cosine part
	sin : float
		The amplitude of the sine part
	"""
	parameter_names = ("freq", "cos", "sin")

	def get_value(self, t):
		t = np.atleast_1d(t)
		return (self.cos * np.cos(self.freq * 2 * np.pi * t) +
				self.sin * np.sin(self.freq * 2 * np.pi * t))

	def get_amplitude(self):
		return np.sqrt(self.cos**2 + self.sin**2)

	def get_phase(self):
		return np.arctan2(self.sin, self.cos)

	def compute_gradient(self, t):
		t = np.atleast_1d(t)
		dcos = np.cos(self.freq * 2 * np.pi * t)
		dsin = np.sin(self.freq * 2 * np.pi * t)
		df = 2 * np.pi * t * (self.sin * dcos - self.cos * dsin)
		return np.array([df, dcos, dsin])


class HarmonicModelAmpPhase(Model):
	"""Model for harmonic terms

	Models harmonic terms using a cosine and sine part.
	The total amplitude and phase can be inferred from that.

	Parameters
	----------
	freq : float
		The frequency in years^-1
	amp : float
		The amplitude of the harmonic term
	phase : float
		The phase of the harmonic part
	"""
	parameter_names = ("freq", "amp", "phase")

	def get_value(self, t):
		t = np.atleast_1d(t)
		return self.amp * np.cos(self.freq * 2 * np.pi * t + self.phase)

	def get_amplitude(self):
		return self.amp

	def get_phase(self):
		return self.phase

	def compute_gradient(self, t):
		t = np.atleast_1d(t)
		damp = np.cos(self.freq * 2 * np.pi * t + self.phase)
		dphi = -self.amp * np.sin(self.freq * 2 * np.pi * t + self.phase)
		df = 2 * np.pi * t * dphi
		return np.array([df, damp, dphi])


class ProxyModel(Model):
	"""Model for proxy terms

	Models proxy terms with a finite and (semi-)annually varying life time.

	Parameters
	----------
	proxy_times : (N,) array_like
		The data times of the proxy values
	proxy_vals : (N,) array_like
		The proxy values at `proxy_times`
	amp : float
		The amplitude of the proxy term
	lag : float
		The lag of the proxy value in years.
	tau0 : float
		The base life time of the proxy
	taucos1 : float
		The amplitude of the cosine part of the annual life time variation.
	tausin1 : float
		The amplitude of the sine part of the annual life time variation.
	taucos2 : float
		The amplitude of the cosine part of the semi-annual life time variation.
	tausin2 : float
		The amplitude of the sine part of the semi-annual life time variation.
	ltscan : float
		The number of days to sum the previous proxy values. If it is
		negative, the value will be set to three times the maximal lifetime.
		No lifetime adjustemets are calculated when set to zero.
	center : bool, optional
		Centers the proxy values by subtracting the overall mean. The mean is
		calculated from the whole `proxy_vals` array and is stored in the
		`mean` attribute.
		Default: False
	sza_intp : scipy.interpolate.interp1d() instance, optional
		When not `None`, cos(sza) and sin(sza) are used instead
		of the time to model the annual variation of the lifetime.
		Semi-annual variations are not used in that case.
		Default: None
	fit_phase : bool, optional
		Fit the phase shift directly instead of using sine and cosine
		terms for the (semi-)annual lifetime variations. If True, the fitted
		cosine parameter is the amplitude and the sine parameter the phase.
		Default: False (= fit sine and cosine terms)
	lifetime_prior : str, optional
		The prior probability density for each coefficient of the lifetime.
		Possible types are "flat" or `None` for a flat prior, "exp" for an
		exponential density ~ :math:`\\text{exp}(-|\\tau| / \\text{metric})`,
		and "normal" for a normal distribution
		~ :math:`\\text{exp}(-\\tau^2 / (2 * \\text{metric}^2))`.
		Default: None (= flat prior).
	lifetime_metric : float, optional
		The metric (scale) of the lifetime priors in days, see `prior`.
		Default 1.
	days_per_time_unit : float, optional
		The number of days per time unit, used to normalize the lifetime
		units. Use 365.25 if the times are in fractional years, or 1 if
		they are in days.
		Default: 365.25
	"""
	parameter_names = ("amp", "lag", "tau0",
			"taucos1", "tausin1", "taucos2", "tausin2",
			"ltscan")

	def __init__(self, proxy_times, proxy_vals,
			center=False,
			sza_intp=None, fit_phase=False,
			lifetime_prior=None, lifetime_metric=1.,
			days_per_time_unit=365.25,
			*args, **kwargs):
		self.mean = 0.
		if center:
			self.mean = np.nanmean(proxy_vals)
		self.times = proxy_times
		self.dt = 1.
		self.values = proxy_vals - self.mean
		self.sza_intp = sza_intp
		self.fit_phase = fit_phase
		self.days_per_time_unit = days_per_time_unit
		self.omega = 2 * np.pi * days_per_time_unit / 365.25
		self.lifetime_prior = lifetime_prior
		self.lifetime_metric = lifetime_metric
		# Makes "(m)jd" and "jyear" compatible for the lifetime
		# seasonal variation. The julian epoch (the default)
		# is slightly offset with respect to (modified) julian days.
		self.t_adj = 0.
		if self.days_per_time_unit == 1:
			# discriminate between julian days and modified julian days,
			# 1.8e6 is year 216 in julian days and year 6787 in
			# modified julian days. It should be pretty safe to judge on
			# that for most use cases.
			if self.times[0] > 1.8e6:
				# julian days
				self.t_adj = 13.
			else:
				# modified julian days
				self.t_adj = -44.25
		super(ProxyModel, self).__init__(*args, **kwargs)

	def _lt_corr(self, t, tau, tmax=60.):

		"""Lifetime corrected values

		Corrects for a finite lifetime by summing over the last `tmax`
		days with an exponential decay given of lifetime(s) `taus`.
		"""
		bs = np.arange(self.dt, tmax + self.dt, self.dt)
		yp = np.zeros_like(t)
		tauexp = np.exp(-self.dt / tau)
		taufac = np.ones_like(tau)
		for b in bs:
			taufac *= tauexp
			yp += np.interp(t - self.lag - b / self.days_per_time_unit,
					self.times, self.values, left=0., right=0.) * taufac
		return yp * self.dt

	def _lt_corr_grad(self, t, tau, tmax=60.):

		"""Lifetime corrected gradient

		Corrects for a finite lifetime by summing over the last `tmax`
		days with an exponential decay given of lifetime(s) `taus`.
		"""
		bs = np.arange(self.dt, tmax + self.dt, self.dt)
		ypg = np.zeros_like(t)
		tauexp = np.exp(-self.dt / tau)
		taufac = np.ones_like(tau)
		for b in bs:
			taufac *= tauexp
			ypg += np.interp(t - self.lag - b / self.days_per_time_unit,
					self.times, self.values, left=0., right=0.) * taufac * b
		return ypg * self.dt / tau**2

	def get_value(self, t):
		t = np.atleast_1d(t)
		proxy_val = np.interp(t - self.lag,
				self.times, self.values, left=0., right=0.)
		if self.ltscan == 0:
			# no lifetime, nothing else to do
			return self.amp * proxy_val
		# annual variation of the proxy lifetime
		if self.sza_intp is not None:
			# using the solar zenith angle
			tau_cs = (self.taucos1 * np.cos(np.radians(self.sza_intp(t)))
					+ self.tausin1 * np.sin(np.radians(self.sza_intp(t))))
		elif self.fit_phase:
			# using time (cos) and phase (sin)
			tau_cs = (self.taucos1 * np.cos(1 * self.omega * (t + self.t_adj) + self.tausin1)
					+ self.taucos2 * np.cos(2 * self.omega * (t + self.t_adj) + self.tausin2))
		else:
			# using time
			tau_cs = (self.taucos1 * np.cos(1 * self.omega * (t + self.t_adj))
					+ self.tausin1 * np.sin(1 * self.omega * (t + self.t_adj))
					+ self.taucos2 * np.cos(2 * self.omega * (t + self.t_adj))
					+ self.tausin2 * np.sin(2 * self.omega * (t + self.t_adj)))
		tau_cs = np.maximum(0., tau_cs)  # clip to zero
		tau = self.tau0 + tau_cs
		if self.ltscan > 0:
			_ltscn = int(np.floor(self.ltscan))
		else:
			# infer the scan time from the maximal lifetime
			_ltscn = 3 * int(np.ceil(self.tau0 +
						np.sqrt(self.taucos1**2 + self.tausin1**2)))
		if np.all(tau > 0):
			proxy_val += self._lt_corr(t, tau, tmax=_ltscn)
		return self.amp * proxy_val

	def compute_gradient(self, t):
		t = np.atleast_1d(t)
		proxy_val = np.interp(t - self.lag,
				self.times, self.values, left=0., right=0.)
		proxy_val_grad0 = proxy_val.copy()
		# annual variation of the proxy lifetime
		if self.sza_intp is not None:
			# using the solar zenith angle
			dtau_cos1 = np.cos(np.radians(self.sza_intp(t)))
			dtau_sin1 = np.sin(np.radians(self.sza_intp(t)))
			dtau_cos2 = np.zeros_like(t)
			dtau_sin2 = np.zeros_like(t)
			tau_cs = self.taucos1 * dtau_cos1 + self.tausin1 * dtau_sin1
		elif self.fit_phase:
			# using time (cos) and phase (sin)
			dtau_cos1 = np.cos(1 * self.omega * (t + self.t_adj) + self.tausin1)
			dtau_sin1 = -self.taucos1 * np.sin(1 * self.omega * t + self.tausin1)
			dtau_cos2 = np.cos(2 * self.omega * (t + self.t_adj) + self.tausin2)
			dtau_sin2 = -self.taucos2 * np.sin(2 * self.omega * t + self.tausin2)
			tau_cs = self.taucos1 * dtau_cos1 + self.taucos2 * dtau_cos2
		else:
			# using time
			dtau_cos1 = np.cos(1 * self.omega * (t + self.t_adj))
			dtau_sin1 = np.sin(1 * self.omega * (t + self.t_adj))
			dtau_cos2 = np.cos(2 * self.omega * (t + self.t_adj))
			dtau_sin2 = np.sin(2 * self.omega * (t + self.t_adj))
			tau_cs = (self.taucos1 * dtau_cos1 + self.tausin1 * dtau_sin1 +
					self.taucos2 * dtau_cos2 + self.tausin2 * dtau_sin2)
		tau_cs = np.maximum(0., tau_cs)  # clip to zero
		tau = self.tau0 + tau_cs
		if self.ltscan > 0:
			_ltscn = int(np.floor(self.ltscan))
		else:
			# infer the scan time from the maximal lifetime
			_ltscn = 3 * int(np.ceil(self.tau0 +
						np.sqrt(self.taucos1**2 + self.tausin1**2)))
		if np.all(tau > 0):
			proxy_val += self._lt_corr(t, tau, tmax=_ltscn)
			proxy_val_grad0 += self._lt_corr_grad(t, tau, tmax=_ltscn)
		return np.array([proxy_val,
				# set the gradient wrt lag to zero for now
				np.zeros_like(t),
				self.amp * proxy_val_grad0,
				self.amp * proxy_val_grad0 * dtau_cos1,
				self.amp * proxy_val_grad0 * dtau_sin1,
				self.amp * proxy_val_grad0 * dtau_cos2,
				self.amp * proxy_val_grad0 * dtau_sin2,
				# set the gradient wrt lifetime scan to zero for now
				np.zeros_like(t)])

	def _log_prior_normal(self):
		l_prior = super(ProxyModel, self).log_prior()
		if not np.isfinite(l_prior):
			return -np.inf
		for n, p in self.get_parameter_dict().items():
			if n.startswith("tau"):
				# Gaussian prior for the lifetimes
				l_prior -= 0.5 * (p / self.lifetime_metric)**2
		return l_prior

	def _log_prior_exp(self):
		l_prior = super(ProxyModel, self).log_prior()
		if not np.isfinite(l_prior):
			return -np.inf
		for n, p in self.get_parameter_dict().items():
			if n.startswith("tau"):
				# exponential prior for the lifetimes
				l_prior -= np.abs(p / self.lifetime_metric)
		return l_prior

	def log_prior(self):
		_priors = {"exp": self._log_prior_exp,
				"normal": self._log_prior_normal}
		if self.lifetime_prior is None or self.lifetime_prior == "flat":
			return super(ProxyModel, self).log_prior()
		return _priors[self.lifetime_prior]()


class TraceGasModelSet(ModelSet):
	"""Combined model class for trace gases (and probably other data)

	Inherited from :class:`celerite.ModelSet`, provides `get_value()`
	and `compute_gradient()` methods.
	"""
	def get_value(self, t):
		t = np.atleast_1d(t)
		v = np.zeros_like(t)
		for m in self.models.values():
			v += m.get_value(t)
		return v

	def compute_gradient(self, t):
		t = np.atleast_1d(t)
		grad = []
		for m in self.models.values():
			grad.extend(list(m.compute_gradient(t)))
		return np.array(grad)


def setup_proxy_model_with_bounds(times, values,
		max_amp=1e10, max_days=100,
		**kwargs):
	# extract setup from `kwargs`
	center = kwargs.get("center", False)
	fit_phase = kwargs.get("fit_phase", False)
	lag = kwargs.get("lag", 0.)
	lt_metric = kwargs.get("lifetime_metric", 1)
	lt_prior = kwargs.get("lifetime_prior", "exp")
	lt_scan = kwargs.get("lifetime_scan", 60)
	positive = kwargs.get("positive", False)
	sza_intp = kwargs.get("sza_intp", None)
	time_format = kwargs.get("time_format", "jyear")

	return ProxyModel(times, values,
			center=center,
			sza_intp=sza_intp,
			fit_phase=fit_phase,
			lifetime_prior=lt_prior,
			lifetime_metric=lt_metric,
			days_per_time_unit=1 if time_format.endswith("d") else 365.25,
			amp=0.,
			lag=lag,
			tau0=0,
			taucos1=0, tausin1=0,
			taucos2=0, tausin2=0,
			ltscan=lt_scan,
			bounds=dict([
				("amp", [0, max_amp] if positive else [-max_amp, max_amp]),
				("lag", [0, max_days]),
				("tau0", [0, max_days]),
				("taucos1", [0, max_days] if fit_phase else [-max_days, max_days]),
				("tausin1", [-np.pi, np.pi] if fit_phase else [-max_days, max_days]),
				# semi-annual cycles for the life time
				("taucos2", [0, max_days] if fit_phase else [-max_days, max_days]),
				("tausin2", [-np.pi, np.pi] if fit_phase else [-max_days, max_days]),
				("ltscan", [0, 200])])
			)


def _default_proxy_config(tfmt="jyear"):

	from .load_data import load_dailymeanLya, load_dailymeanAE
	proxy_config = {}
	# Lyman-alpha
	plyat, plyadf = load_dailymeanLya(tfmt=tfmt)
	proxy_config.update({"Lya": {
		"times": plyat,
		"values": plyadf["Lya"],
		"center": False,
		"positive": False,
		"lifetime_scan": 0,
		}}
	)
	# AE index
	paet, paedf = load_dailymeanAE(name="GM", tfmt=tfmt)
	proxy_config.update({"GM": {
		"times": paet,
		"values": paedf["GM"],
		"center": False,
		"positive": True,
		"lifetime_scan": 60,
		}}
	)
	return proxy_config


def trace_gas_model(constant=True, freqs=None, proxy_config=None, **kwargs):
	"""Trace gas model setup

	Sets up the trace gas model for easy access. All parameters are optional,
	defaults to an offset, no harmonics, proxies are uncentered and unscaled
	Lyman-alpha and AE. AE with only positive amplitude and a seasonally
	varying lifetime.

	Parameters
	----------
	constant : bool, optional
		Whether or not to include a constant (offset) term, default is True.
	freqs : list, optional
		Frequencies of the harmonic terms in 1 / a^-1 (inverse years).
	proxy_config : dict, optional
		Proxy configuration if different from the standard setup.
	**kwargs : optional
		Additional keyword arguments, all of them are also passed on to
		the proxy setup. For now, supported are the following which are
		also passed along to the proxy setup with
		`setup_proxy_model_with_bounds()`:

		* fit_phase : bool
			fit amplitude and phase instead of sine and cosine
		* scale : float
			the factor by which the data is scaled, used to constrain
			the maximum and minimum amplitudes to be fitted.
		* time_format : string
			The `astropy.time.Time` format string to setup the time axis.
		* days_per_time_unit : float
			The number of days per time unit, used to normalize the frequencies
			for the harmonic terms. Use 365.25 if the times are in fractional years,
			1 if they are in days. Default: 365.25
		* max_amp : float
			Maximum magnitude of the coefficients, used to constrain the
			parameter search.
		* max_days : float
			Maximum magnitude of the lifetimes, used to constrain the
			parameter search.

	Returns
	-------
	model : :class:`TraceGasModelSet` (extends :class:`celerite.ModelSet`)
	"""
	fit_phase = kwargs.get("fit_phase", False)
	scale = kwargs.get("scale", 1e-6)
	tfmt = kwargs.get("time_format", "jyear")
	delta_t = kwargs.get("days_per_time_unit", 365.25)

	max_amp = kwargs.pop("max_amp", 1e10 * scale)
	max_days = kwargs.pop("max_days", 100)

	offset_model = []
	if constant:
		offset_model = [("offset",
				ConstantModel(value=0.,
						bounds={"value": [-max_amp, max_amp]}))]

	freqs = freqs or []
	harmonic_models = []
	for freq in freqs:
		if not fit_phase:
			harm = HarmonicModelCosineSine(freq=freq * delta_t / 365.25,
					cos=0, sin=0,
					bounds=dict([
						("cos", [-max_amp, max_amp]),
						("sin", [-max_amp, max_amp])])
			)
		else:
			harm = HarmonicModelAmpPhase(freq=freq * delta_t / 365.25,
					amp=0, phase=0,
					bounds=dict([
						("amp", [0, max_amp]),
						("phase", [-np.pi, np.pi])])
			)
		harm.freeze_parameter("freq")
		harmonic_models.append(("f{0:.0f}".format(freq), harm))

	proxy_config = proxy_config or _default_proxy_config(tfmt=tfmt)
	proxy_models = []
	for pn, conf in proxy_config.items():
		if "max_amp" not in conf:
			conf.update(dict(max_amp=max_amp))
		if "max_days" not in conf:
			conf.update(dict(max_days=max_days))
		kw = kwargs.copy()  # don't mess with the passed arguments
		kw.update(conf)
		proxy_models.append(
			(pn, setup_proxy_model_with_bounds(**kw))
		)

	return TraceGasModelSet(offset_model + harmonic_models + proxy_models)


1			# -- coding: utf-8 --
2			# vim:fileencoding=utf-8
3			#
4			# Copyright (c) 2017-2018 Stefan Bender
5			#
6			# This module is part of sciapy.
7			# sciapy is free software: you can redistribute it or modify
8			# it under the terms of the GNU General Public License as published
9			# by the Free Software Foundation, version 2.
10			# See accompanying LICENSE file or http://www.gnu.org/licenses/gpl-2.0.html.
11	1		"""SCIAMACHY regression models (celerite version)
12
13			Model classes for SCIAMACHY data regression fits using the
14			:mod:`celerite` [#]_ modelling protocol.
15
16			.. [#] https://celerite.readthedocs.io
17			"""
18	1		from __future__ import absolute_import, division, print_function
19
20	1		import numpy as np
21	1		from scipy.interpolate import interp1d
22
23	1		from celerite.modeling import Model, ModelSet, ConstantModel
24
25	1		__all__ = ["ConstantModel",
26			"HarmonicModelCosineSine", "HarmonicModelAmpPhase",
27			"ProxyModel", "TraceGasModelSet",
28			"setup_proxy_model_with_bounds", "trace_gas_model"]
29
30
31	1		class HarmonicModelCosineSine(Model):
32			"""Model for harmonic terms
33
34			Models harmonic terms using a cosine and sine part.
35			The total amplitude and phase can be inferred from that.
36
37			Parameters
38			----------
39			freq : float
40			The frequency in years^-1
41			cos : float
42			The amplitude of the cosine part
43			sin : float
44			The amplitude of the sine part
45			"""
46	1		parameter_names = ("freq", "cos", "sin")
47
48	1		def get_value(self, t):
49			t = np.atleast_1d(t)
50			return (self.cos * np.cos(self.freq * 2 * np.pi * t) +
51			self.sin * np.sin(self.freq * 2 * np.pi * t))
52
53	1		def get_amplitude(self):
54			return np.sqrt(self.cos2 + self.sin2)
55
56	1		def get_phase(self):
57			return np.arctan2(self.sin, self.cos)
58
59	1		def compute_gradient(self, t):
60			t = np.atleast_1d(t)
61			dcos = np.cos(self.freq * 2 * np.pi * t)
62			dsin = np.sin(self.freq * 2 * np.pi * t)
63			df = 2 * np.pi * t * (self.sin * dcos - self.cos * dsin)
64			return np.array([df, dcos, dsin])
65
66
67	1		class HarmonicModelAmpPhase(Model):
68			"""Model for harmonic terms
69
70			Models harmonic terms using a cosine and sine part.
71			The total amplitude and phase can be inferred from that.
72
73			Parameters
74			----------
75			freq : float
76			The frequency in years^-1
77			amp : float
78			The amplitude of the harmonic term
79			phase : float
80			The phase of the harmonic part
81			"""
82	1		parameter_names = ("freq", "amp", "phase")
83
84	1		def get_value(self, t):
85			t = np.atleast_1d(t)
86			return self.amp * np.cos(self.freq * 2 * np.pi * t + self.phase)
87
88	1		def get_amplitude(self):
89			return self.amp
90
91	1		def get_phase(self):
92			return self.phase
93
94	1		def compute_gradient(self, t):
95			t = np.atleast_1d(t)
96			damp = np.cos(self.freq * 2 * np.pi * t + self.phase)
97			dphi = -self.amp * np.sin(self.freq * 2 * np.pi * t + self.phase)
98			df = 2 * np.pi * t * dphi
99			return np.array([df, damp, dphi])
100
101
102	1		class ProxyModel(Model):
103			"""Model for proxy terms
104
105			Models proxy terms with a finite and (semi-)annually varying life time.
106
107			Parameters
108			----------
109			proxy_times : (N,) array_like
110			The data times of the proxy values
111			proxy_vals : (N,) array_like
112			The proxy values at `proxy_times`
113			amp : float
114			The amplitude of the proxy term
115			lag : float
116			The lag of the proxy value in years.
117			tau0 : float
118			The base life time of the proxy
119			taucos1 : float
120			The amplitude of the cosine part of the annual life time variation.
121			tausin1 : float
122			The amplitude of the sine part of the annual life time variation.
123			taucos2 : float
124			The amplitude of the cosine part of the semi-annual life time variation.
125			tausin2 : float
126			The amplitude of the sine part of the semi-annual life time variation.
127			ltscan : float
128			The number of days to sum the previous proxy values. If it is
129			negative, the value will be set to three times the maximal lifetime.
130			No lifetime adjustemets are calculated when set to zero.
131			center : bool, optional
132			Centers the proxy values by subtracting the overall mean. The mean is
133			calculated from the whole `proxy_vals` array and is stored in the
134			`mean` attribute.
135			Default: False
136			sza_intp : scipy.interpolate.interp1d() instance, optional
137			When not `None`, cos(sza) and sin(sza) are used instead
138			of the time to model the annual variation of the lifetime.
139			Semi-annual variations are not used in that case.
140			Default: None
141			fit_phase : bool, optional
142			Fit the phase shift directly instead of using sine and cosine
143			terms for the (semi-)annual lifetime variations. If True, the fitted
144			cosine parameter is the amplitude and the sine parameter the phase.
145			Default: False (= fit sine and cosine terms)
146			lifetime_prior : str, optional
147			The prior probability density for each coefficient of the lifetime.
148			Possible types are "flat" or `None` for a flat prior, "exp" for an
149			exponential density ~ :math:`\\text{exp}(-\|\\tau\| / \\text{metric})`,
150			and "normal" for a normal distribution
151			~ :math:`\\text{exp}(-\\tau^2 / (2 * \\text{metric}^2))`.
152			Default: None (= flat prior).
153			lifetime_metric : float, optional
154			The metric (scale) of the lifetime priors in days, see `prior`.
155			Default 1.
156			days_per_time_unit : float, optional
157			The number of days per time unit, used to normalize the lifetime
158			units. Use 365.25 if the times are in fractional years, or 1 if
159			they are in days.
160			Default: 365.25
161			"""
162	1		parameter_names = ("amp", "lag", "tau0",
163			"taucos1", "tausin1", "taucos2", "tausin2",
164			"ltscan")
165
166	1		def __init__(self, proxy_times, proxy_vals,
167			center=False,
168			sza_intp=None, fit_phase=False,
169			lifetime_prior=None, lifetime_metric=1.,
170			days_per_time_unit=365.25,
171			args, *kwargs):
172			self.mean = 0.
173			if center:
174			self.mean = np.nanmean(proxy_vals)
175			self.times = proxy_times
176			self.dt = 1.
177			self.values = proxy_vals - self.mean
178			self.sza_intp = sza_intp
179			self.fit_phase = fit_phase
180			self.days_per_time_unit = days_per_time_unit
181			self.omega = 2 * np.pi * days_per_time_unit / 365.25
182			self.lifetime_prior = lifetime_prior
183			self.lifetime_metric = lifetime_metric
184			# Makes "(m)jd" and "jyear" compatible for the lifetime
185			# seasonal variation. The julian epoch (the default)
186			# is slightly offset with respect to (modified) julian days.
187			self.t_adj = 0.
188			if self.days_per_time_unit == 1:
189			# discriminate between julian days and modified julian days,
190			# 1.8e6 is year 216 in julian days and year 6787 in
191			# modified julian days. It should be pretty safe to judge on
192			# that for most use cases.
193			if self.times[0] > 1.8e6:
194			# julian days
195			self.t_adj = 13.
196			else:
197			# modified julian days
198			self.t_adj = -44.25
199			super(ProxyModel, self).__init__(args, *kwargs)
200
201	1	View Code Duplication	def _lt_corr(self, t, tau, tmax=60.):
			0 ignored issues – show Duplication introduced 2020-11-04 16:18 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
202			"""Lifetime corrected values
203
204			Corrects for a finite lifetime by summing over the last `tmax`
205			days with an exponential decay given of lifetime(s) `taus`.
206			"""
207			bs = np.arange(self.dt, tmax + self.dt, self.dt)
208			yp = np.zeros_like(t)
209			tauexp = np.exp(-self.dt / tau)
210			taufac = np.ones_like(tau)
211			for b in bs:
212			taufac *= tauexp
213			yp += np.interp(t - self.lag - b / self.days_per_time_unit,
214			self.times, self.values, left=0., right=0.) * taufac
215			return yp * self.dt
216
217	1	View Code Duplication	def _lt_corr_grad(self, t, tau, tmax=60.):
			0 ignored issues – show Duplication introduced 2020-11-04 16:18 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
218			"""Lifetime corrected gradient
219
220			Corrects for a finite lifetime by summing over the last `tmax`
221			days with an exponential decay given of lifetime(s) `taus`.
222			"""
223			bs = np.arange(self.dt, tmax + self.dt, self.dt)
224			ypg = np.zeros_like(t)
225			tauexp = np.exp(-self.dt / tau)
226			taufac = np.ones_like(tau)
227			for b in bs:
228			taufac *= tauexp
229			ypg += np.interp(t - self.lag - b / self.days_per_time_unit,
230			self.times, self.values, left=0., right=0.) * taufac * b
231			return ypg * self.dt / tau**2
232
233	1		def get_value(self, t):
234			t = np.atleast_1d(t)
235			proxy_val = np.interp(t - self.lag,
236			self.times, self.values, left=0., right=0.)
237			if self.ltscan == 0:
238			# no lifetime, nothing else to do
239			return self.amp * proxy_val
240			# annual variation of the proxy lifetime
241			if self.sza_intp is not None:
242			# using the solar zenith angle
243			tau_cs = (self.taucos1 * np.cos(np.radians(self.sza_intp(t)))
244			+ self.tausin1 * np.sin(np.radians(self.sza_intp(t))))
245			elif self.fit_phase:
246			# using time (cos) and phase (sin)
247			tau_cs = (self.taucos1 * np.cos(1 * self.omega * (t + self.t_adj) + self.tausin1)
248			+ self.taucos2 * np.cos(2 * self.omega * (t + self.t_adj) + self.tausin2))
249			else:
250			# using time
251			tau_cs = (self.taucos1 * np.cos(1 * self.omega * (t + self.t_adj))
252			+ self.tausin1 * np.sin(1 * self.omega * (t + self.t_adj))
253			+ self.taucos2 * np.cos(2 * self.omega * (t + self.t_adj))
254			+ self.tausin2 * np.sin(2 * self.omega * (t + self.t_adj)))
255			tau_cs = np.maximum(0., tau_cs) # clip to zero
256			tau = self.tau0 + tau_cs
257			if self.ltscan > 0:
258			_ltscn = int(np.floor(self.ltscan))
259			else:
260			# infer the scan time from the maximal lifetime
261			_ltscn = 3 * int(np.ceil(self.tau0 +
262			np.sqrt(self.taucos12 + self.tausin12)))
263			if np.all(tau > 0):
264			proxy_val += self._lt_corr(t, tau, tmax=_ltscn)
265			return self.amp * proxy_val
266
267	1		def compute_gradient(self, t):
268			t = np.atleast_1d(t)
269			proxy_val = np.interp(t - self.lag,
270			self.times, self.values, left=0., right=0.)
271			proxy_val_grad0 = proxy_val.copy()
272			# annual variation of the proxy lifetime
273			if self.sza_intp is not None:
274			# using the solar zenith angle
275			dtau_cos1 = np.cos(np.radians(self.sza_intp(t)))
276			dtau_sin1 = np.sin(np.radians(self.sza_intp(t)))
277			dtau_cos2 = np.zeros_like(t)
278			dtau_sin2 = np.zeros_like(t)
279			tau_cs = self.taucos1 * dtau_cos1 + self.tausin1 * dtau_sin1
280			elif self.fit_phase:
281			# using time (cos) and phase (sin)
282			dtau_cos1 = np.cos(1 * self.omega * (t + self.t_adj) + self.tausin1)
283			dtau_sin1 = -self.taucos1 * np.sin(1 * self.omega * t + self.tausin1)
284			dtau_cos2 = np.cos(2 * self.omega * (t + self.t_adj) + self.tausin2)
285			dtau_sin2 = -self.taucos2 * np.sin(2 * self.omega * t + self.tausin2)
286			tau_cs = self.taucos1 * dtau_cos1 + self.taucos2 * dtau_cos2
287			else:
288			# using time
289			dtau_cos1 = np.cos(1 * self.omega * (t + self.t_adj))
290			dtau_sin1 = np.sin(1 * self.omega * (t + self.t_adj))
291			dtau_cos2 = np.cos(2 * self.omega * (t + self.t_adj))
292			dtau_sin2 = np.sin(2 * self.omega * (t + self.t_adj))
293			tau_cs = (self.taucos1 * dtau_cos1 + self.tausin1 * dtau_sin1 +
294			self.taucos2 * dtau_cos2 + self.tausin2 * dtau_sin2)
295			tau_cs = np.maximum(0., tau_cs) # clip to zero
296			tau = self.tau0 + tau_cs
297			if self.ltscan > 0:
298			_ltscn = int(np.floor(self.ltscan))
299			else:
300			# infer the scan time from the maximal lifetime
301			_ltscn = 3 * int(np.ceil(self.tau0 +
302			np.sqrt(self.taucos12 + self.tausin12)))
303			if np.all(tau > 0):
304			proxy_val += self._lt_corr(t, tau, tmax=_ltscn)
305			proxy_val_grad0 += self._lt_corr_grad(t, tau, tmax=_ltscn)
306			return np.array([proxy_val,
307			# set the gradient wrt lag to zero for now
308			np.zeros_like(t),
309			self.amp * proxy_val_grad0,
310			self.amp * proxy_val_grad0 * dtau_cos1,
311			self.amp * proxy_val_grad0 * dtau_sin1,
312			self.amp * proxy_val_grad0 * dtau_cos2,
313			self.amp * proxy_val_grad0 * dtau_sin2,
314			# set the gradient wrt lifetime scan to zero for now
315			np.zeros_like(t)])
316
317	1		def _log_prior_normal(self):
318			l_prior = super(ProxyModel, self).log_prior()
319			if not np.isfinite(l_prior):
320			return -np.inf
321			for n, p in self.get_parameter_dict().items():
322			if n.startswith("tau"):
323			# Gaussian prior for the lifetimes
324			l_prior -= 0.5 * (p / self.lifetime_metric)**2
325			return l_prior
326
327	1		def _log_prior_exp(self):
328			l_prior = super(ProxyModel, self).log_prior()
329			if not np.isfinite(l_prior):
330			return -np.inf
331			for n, p in self.get_parameter_dict().items():
332			if n.startswith("tau"):
333			# exponential prior for the lifetimes
334			l_prior -= np.abs(p / self.lifetime_metric)
335			return l_prior
336
337	1		def log_prior(self):
338			_priors = {"exp": self._log_prior_exp,
339			"normal": self._log_prior_normal}
340			if self.lifetime_prior is None or self.lifetime_prior == "flat":
341			return super(ProxyModel, self).log_prior()
342			return _priors[self.lifetime_prior]()
343
344
345	1		class TraceGasModelSet(ModelSet):
346			"""Combined model class for trace gases (and probably other data)
347
348			Inherited from :class:`celerite.ModelSet`, provides `get_value()`
349			and `compute_gradient()` methods.
350			"""
351	1		def get_value(self, t):
352			t = np.atleast_1d(t)
353			v = np.zeros_like(t)
354			for m in self.models.values():
355			v += m.get_value(t)
356			return v
357
358	1		def compute_gradient(self, t):
359			t = np.atleast_1d(t)
360			grad = []
361			for m in self.models.values():
362			grad.extend(list(m.compute_gradient(t)))
363			return np.array(grad)
364
365
366	1		def setup_proxy_model_with_bounds(times, values,
367			max_amp=1e10, max_days=100,
368			**kwargs):
369			# extract setup from `kwargs`
370			center = kwargs.get("center", False)
371			fit_phase = kwargs.get("fit_phase", False)
372			lag = kwargs.get("lag", 0.)
373			lt_metric = kwargs.get("lifetime_metric", 1)
374			lt_prior = kwargs.get("lifetime_prior", "exp")
375			lt_scan = kwargs.get("lifetime_scan", 60)
376			positive = kwargs.get("positive", False)
377			sza_intp = kwargs.get("sza_intp", None)
378			time_format = kwargs.get("time_format", "jyear")
379
380			return ProxyModel(times, values,
381			center=center,
382			sza_intp=sza_intp,
383			fit_phase=fit_phase,
384			lifetime_prior=lt_prior,
385			lifetime_metric=lt_metric,
386			days_per_time_unit=1 if time_format.endswith("d") else 365.25,
387			amp=0.,
388			lag=lag,
389			tau0=0,
390			taucos1=0, tausin1=0,
391			taucos2=0, tausin2=0,
392			ltscan=lt_scan,
393			bounds=dict([
394			("amp", [0, max_amp] if positive else [-max_amp, max_amp]),
395			("lag", [0, max_days]),
396			("tau0", [0, max_days]),
397			("taucos1", [0, max_days] if fit_phase else [-max_days, max_days]),
398			("tausin1", [-np.pi, np.pi] if fit_phase else [-max_days, max_days]),
399			# semi-annual cycles for the life time
400			("taucos2", [0, max_days] if fit_phase else [-max_days, max_days]),
401			("tausin2", [-np.pi, np.pi] if fit_phase else [-max_days, max_days]),
402			("ltscan", [0, 200])])
403			)
404
405
406	1	View Code Duplication	def _default_proxy_config(tfmt="jyear"):
			0 ignored issues – show Duplication introduced 2022-04-06 09:56 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
407			from .load_data import load_dailymeanLya, load_dailymeanAE
408			proxy_config = {}
409			# Lyman-alpha
410			plyat, plyadf = load_dailymeanLya(tfmt=tfmt)
411			proxy_config.update({"Lya": {
412			"times": plyat,
413			"values": plyadf["Lya"],
414			"center": False,
415			"positive": False,
416			"lifetime_scan": 0,
417			}}
418			)
419			# AE index
420			paet, paedf = load_dailymeanAE(name="GM", tfmt=tfmt)
421			proxy_config.update({"GM": {
422			"times": paet,
423			"values": paedf["GM"],
424			"center": False,
425			"positive": True,
426			"lifetime_scan": 60,
427			}}
428			)
429			return proxy_config
430
431
432	1		def trace_gas_model(constant=True, freqs=None, proxy_config=None, **kwargs):
433			"""Trace gas model setup
434
435			Sets up the trace gas model for easy access. All parameters are optional,
436			defaults to an offset, no harmonics, proxies are uncentered and unscaled
437			Lyman-alpha and AE. AE with only positive amplitude and a seasonally
438			varying lifetime.
439
440			Parameters
441			----------
442			constant : bool, optional
443			Whether or not to include a constant (offset) term, default is True.
444			freqs : list, optional
445			Frequencies of the harmonic terms in 1 / a^-1 (inverse years).
446			proxy_config : dict, optional
447			Proxy configuration if different from the standard setup.
448			**kwargs : optional
449			Additional keyword arguments, all of them are also passed on to
450			the proxy setup. For now, supported are the following which are
451			also passed along to the proxy setup with
452			`setup_proxy_model_with_bounds()`:
453
454			* fit_phase : bool
455			fit amplitude and phase instead of sine and cosine
456			* scale : float
457			the factor by which the data is scaled, used to constrain
458			the maximum and minimum amplitudes to be fitted.
459			* time_format : string
460			The `astropy.time.Time` format string to setup the time axis.
461			* days_per_time_unit : float
462			The number of days per time unit, used to normalize the frequencies
463			for the harmonic terms. Use 365.25 if the times are in fractional years,
464			1 if they are in days. Default: 365.25
465			* max_amp : float
466			Maximum magnitude of the coefficients, used to constrain the
467			parameter search.
468			* max_days : float
469			Maximum magnitude of the lifetimes, used to constrain the
470			parameter search.
471
472			Returns
473			-------
474			model : :class:`TraceGasModelSet` (extends :class:`celerite.ModelSet`)
475			"""
476			fit_phase = kwargs.get("fit_phase", False)
477			scale = kwargs.get("scale", 1e-6)
478			tfmt = kwargs.get("time_format", "jyear")
479			delta_t = kwargs.get("days_per_time_unit", 365.25)
480
481			max_amp = kwargs.pop("max_amp", 1e10 * scale)
482			max_days = kwargs.pop("max_days", 100)
483
484			offset_model = []
485			if constant:
486			offset_model = [("offset",
487			ConstantModel(value=0.,
488			bounds={"value": [-max_amp, max_amp]}))]
489
490			freqs = freqs or []
491			harmonic_models = []
492			for freq in freqs:
493			if not fit_phase:
494			harm = HarmonicModelCosineSine(freq=freq * delta_t / 365.25,
495			cos=0, sin=0,
496			bounds=dict([
497			("cos", [-max_amp, max_amp]),
498			("sin", [-max_amp, max_amp])])
499			)
500			else:
501			harm = HarmonicModelAmpPhase(freq=freq * delta_t / 365.25,
502			amp=0, phase=0,
503			bounds=dict([
504			("amp", [0, max_amp]),
505			("phase", [-np.pi, np.pi])])
506			)
507			harm.freeze_parameter("freq")
508			harmonic_models.append(("f{0:.0f}".format(freq), harm))
509
510			proxy_config = proxy_config or _default_proxy_config(tfmt=tfmt)
511			proxy_models = []
512			for pn, conf in proxy_config.items():
513			if "max_amp" not in conf:
514			conf.update(dict(max_amp=max_amp))
515			if "max_days" not in conf:
516			conf.update(dict(max_days=max_days))
517			kw = kwargs.copy() # don't mess with the passed arguments
518			kw.update(conf)
519			proxy_models.append(
520			(pn, setup_proxy_model_with_bounds(**kw))
521			)
522
523			return TraceGasModelSet(offset_model + harmonic_models + proxy_models)
524

st-bender / sciapy

Push — master ( af9f8c...efc6a8 )

sciapy.regress.models_cel.ProxyModel.log_prior() A

Complexity

Size

Duplication

Code Coverage

Importance

Duplication Side-by-Side

Filter issues like