sciapy.regress.models_cel.HarmonicModelAmpPhase.compute_gradient() - Code Metrics - Inspection of "regress: Default proxy loader fix" - st-bender/sciapy - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Push — master ( 7bc10c...b988a0 )

by Stefan

created 2019-02-20 16:54 UTC

HarmonicModelAmpPhase.compute_gradient() A

↳ Parent: sciapy.regress.models_cel

Complexity

Conditions

Size

Total Lines	5
Code Lines	5

Duplication

Lines	5
Ratio	100 %

Code Coverage

Tests	1
CRAP Score	1.512

Importance

Changes

Metric	Value
cc	1
eloc	5
nop	2
dl	5
loc	5
ccs	1
cts	5
cp	0.2
crap	1.512
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` [#]_ modeling 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", "CeleriteModelSet"]

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):
		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):
		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):
		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):
		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.intp = interp1d(proxy_times, proxy_vals - self.mean,
				bounds_error=False)
		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
		super(ProxyModel, self).__init__(*args, **kwargs)

	def get_value(self, t):
		proxy_val = self.intp(t - self.lag)
		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.tausin1)
					+ self.taucos2 * np.cos(2 * self.omega * t + self.tausin2))
		else:
			# using time
			tau_cs = (self.taucos1 * np.cos(1 * self.omega * t)
					+ self.tausin1 * np.sin(1 * self.omega * t)
					+ self.taucos2 * np.cos(2 * self.omega * t)
					+ self.tausin2 * np.sin(2 * self.omega * t))
		tau_cs[tau_cs < 0] = 0.  # 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):
			bs = np.arange(1, _ltscn + 1, 1.)[None, :]
			taufacs = np.exp(-bs / tau[:, None])
			proxy_val += np.sum(
					self.intp(t[:, None] - self.lag -
						bs / self.days_per_time_unit) * taufacs,
					axis=1)
		return self.amp * proxy_val

	def compute_gradient(self, t):
		proxy_val = self.intp(t - self.lag)
		proxy_val_grad0 = self.intp(t - self.lag)
		# 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.tausin1)
			dtau_sin1 = -self.taucos1 * np.sin(1 * self.omega * t + self.tausin1)
			dtau_cos2 = np.cos(2 * self.omega * t + 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)
			dtau_sin1 = np.sin(1 * self.omega * t)
			dtau_cos2 = np.cos(2 * self.omega * t)
			dtau_sin2 = np.sin(2 * self.omega * t)
			tau_cs = (self.taucos1 * dtau_cos1 + self.tausin1 * dtau_sin1 +
					self.taucos2 * dtau_cos2 + self.tausin2 * dtau_sin2)
		tau_cs[tau_cs < 0] = 0.  # 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):
			bs = np.arange(1, _ltscn + 1, 1.)[None, :]
			taufacs = np.exp(-bs / tau[:, None])
			proxy_ts = self.intp(t[:, None] - self.lag -
					bs / self.days_per_time_unit) * taufacs
			proxy_val += np.sum(proxy_ts, axis=1)
			proxy_val_grad0 += np.sum(proxy_ts * bs / tau[:, None]**2, axis=1)
		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 CeleriteModelSet(ModelSet):


	def get_value(self, t):
		v = np.zeros_like(t)
		for m in self.models.values():
			v += m.get_value(t)
		return v

	def compute_gradient(self, 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])])
			)


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` [#]_ modeling 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", "CeleriteModelSet"]
28
29	1	View Code Duplication	class HarmonicModelCosineSine(Model):
			0 ignored issues – show Duplication introduced 2018-07-02 09:55 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
30			"""Model for harmonic terms
31
32			Models harmonic terms using a cosine and sine part.
33			The total amplitude and phase can be inferred from that.
34
35			Parameters
36			----------
37			freq : float
38			The frequency in years^-1
39			cos : float
40			The amplitude of the cosine part
41			sin : float
42			The amplitude of the sine part
43			"""
44	1		parameter_names = ("freq", "cos", "sin")
45
46	1		def get_value(self, t):
47			return (self.cos * np.cos(self.freq * 2 * np.pi * t) +
48			self.sin * np.sin(self.freq * 2 * np.pi * t))
49
50	1		def get_amplitude(self):
51			return np.sqrt(self.cos2 + self.sin2)
52
53	1		def get_phase(self):
54			return np.arctan2(self.sin, self.cos)
55
56	1		def compute_gradient(self, t):
57			dcos = np.cos(self.freq * 2 * np.pi * t)
58			dsin = np.sin(self.freq * 2 * np.pi * t)
59			df = 2 * np.pi * t * (self.sin * dcos - self.cos * dsin)
60			return np.array([df, dcos, dsin])
61
62
63	1	View Code Duplication	class HarmonicModelAmpPhase(Model):
			0 ignored issues – show Duplication introduced 2018-07-02 09:55 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
64			"""Model for harmonic terms
65
66			Models harmonic terms using a cosine and sine part.
67			The total amplitude and phase can be inferred from that.
68
69			Parameters
70			----------
71			freq : float
72			The frequency in years^-1
73			amp : float
74			The amplitude of the harmonic term
75			phase : float
76			The phase of the harmonic part
77			"""
78	1		parameter_names = ("freq", "amp", "phase")
79
80	1		def get_value(self, t):
81			return self.amp * np.cos(self.freq * 2 * np.pi * t + self.phase)
82
83	1		def get_amplitude(self):
84			return self.amp
85
86	1		def get_phase(self):
87			return self.phase
88
89	1		def compute_gradient(self, t):
90			damp = np.cos(self.freq * 2 * np.pi * t + self.phase)
91			dphi = -self.amp * np.sin(self.freq * 2 * np.pi * t + self.phase)
92			df = 2 * np.pi * t * dphi
93			return np.array([df, damp, dphi])
94
95
96	1	View Code Duplication	class ProxyModel(Model):
			0 ignored issues – show Duplication introduced 2018-07-02 09:55 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
97			"""Model for proxy terms
98
99			Models proxy terms with a finite and (semi-)annually varying life time.
100
101			Parameters
102			----------
103			proxy_times : (N,) array_like
104			The data times of the proxy values
105			proxy_vals : (N,) array_like
106			The proxy values at `proxy_times`
107			amp : float
108			The amplitude of the proxy term
109			lag : float
110			The lag of the proxy value in years.
111			tau0 : float
112			The base life time of the proxy
113			taucos1 : float
114			The amplitude of the cosine part of the annual life time variation.
115			tausin1 : float
116			The amplitude of the sine part of the annual life time variation.
117			taucos2 : float
118			The amplitude of the cosine part of the semi-annual life time variation.
119			tausin2 : float
120			The amplitude of the sine part of the semi-annual life time variation.
121			ltscan : float
122			The number of days to sum the previous proxy values. If it is
123			negative, the value will be set to three times the maximal lifetime.
124			No lifetime adjustemets are calculated when set to zero.
125			center : bool, optional
126			Centers the proxy values by subtracting the overall mean. The mean is
127			calculated from the whole `proxy_vals` array and is stored in the
128			`mean` attribute.
129			Default: False
130			sza_intp : scipy.interpolate.interp1d() instance, optional
131			When not `None`, cos(sza) and sin(sza) are used instead
132			of the time to model the annual variation of the lifetime.
133			Semi-annual variations are not used in that case.
134			Default: None
135			fit_phase : bool, optional
136			Fit the phase shift directly instead of using sine and cosine
137			terms for the (semi-)annual lifetime variations. If True, the fitted
138			cosine parameter is the amplitude and the sine parameter the phase.
139			Default: False (= fit sine and cosine terms)
140			lifetime_prior : str, optional
141			The prior probability density for each coefficient of the lifetime.
142			Possible types are "flat" or `None` for a flat prior, "exp" for an
143			exponential density ~ :math:`\\text{exp}(-\|\\tau\| / \\text{metric})`,
144			and "normal" for a normal distribution
145			~ :math:`\\text{exp}(-\\tau^2 / (2 * \\text{metric}^2))`.
146			Default: None (= flat prior).
147			lifetime_metric : float, optional
148			The metric (scale) of the lifetime priors in days, see `prior`.
149			Default 1.
150			days_per_time_unit : float, optional
151			The number of days per time unit, used to normalize the lifetime
152			units. Use 365.25 if the times are in fractional years, or 1 if
153			they are in days.
154			Default: 365.25
155			"""
156	1		parameter_names = ("amp", "lag", "tau0",
157			"taucos1", "tausin1", "taucos2", "tausin2",
158			"ltscan")
159
160	1		def __init__(self, proxy_times, proxy_vals,
161			center=False,
162			sza_intp=None, fit_phase=False,
163			lifetime_prior=None, lifetime_metric=1.,
164			days_per_time_unit=365.25,
165			args, *kwargs):
166			self.mean = 0.
167			if center:
168			self.mean = np.nanmean(proxy_vals)
169			self.intp = interp1d(proxy_times, proxy_vals - self.mean,
170			bounds_error=False)
171			self.sza_intp = sza_intp
172			self.fit_phase = fit_phase
173			self.days_per_time_unit = days_per_time_unit
174			self.omega = 2 * np.pi * days_per_time_unit / 365.25
175			self.lifetime_prior = lifetime_prior
176			self.lifetime_metric = lifetime_metric
177			super(ProxyModel, self).__init__(args, *kwargs)
178
179	1		def get_value(self, t):
180			proxy_val = self.intp(t - self.lag)
181			if self.ltscan == 0:
182			# no lifetime, nothing else to do
183			return self.amp * proxy_val
184			# annual variation of the proxy lifetime
185			if self.sza_intp is not None:
186			# using the solar zenith angle
187			tau_cs = (self.taucos1 * np.cos(np.radians(self.sza_intp(t)))
188			+ self.tausin1 * np.sin(np.radians(self.sza_intp(t))))
189			elif self.fit_phase:
190			# using time (cos) and phase (sin)
191			tau_cs = (self.taucos1 * np.cos(1 * self.omega * t + self.tausin1)
192			+ self.taucos2 * np.cos(2 * self.omega * t + self.tausin2))
193			else:
194			# using time
195			tau_cs = (self.taucos1 * np.cos(1 * self.omega * t)
196			+ self.tausin1 * np.sin(1 * self.omega * t)
197			+ self.taucos2 * np.cos(2 * self.omega * t)
198			+ self.tausin2 * np.sin(2 * self.omega * t))
199			tau_cs[tau_cs < 0] = 0. # clip to zero
200			tau = self.tau0 + tau_cs
201			if self.ltscan > 0:
202			_ltscn = int(np.floor(self.ltscan))
203			else:
204			# infer the scan time from the maximal lifetime
205			_ltscn = 3 * int(np.ceil(self.tau0 +
206			np.sqrt(self.taucos12 + self.tausin12)))
207			if np.all(tau > 0):
208			bs = np.arange(1, _ltscn + 1, 1.)[None, :]
209			taufacs = np.exp(-bs / tau[:, None])
210			proxy_val += np.sum(
211			self.intp(t[:, None] - self.lag -
212			bs / self.days_per_time_unit) * taufacs,
213			axis=1)
214			return self.amp * proxy_val
215
216	1		def compute_gradient(self, t):
217			proxy_val = self.intp(t - self.lag)
218			proxy_val_grad0 = self.intp(t - self.lag)
219			# annual variation of the proxy lifetime
220			if self.sza_intp is not None:
221			# using the solar zenith angle
222			dtau_cos1 = np.cos(np.radians(self.sza_intp(t)))
223			dtau_sin1 = np.sin(np.radians(self.sza_intp(t)))
224			dtau_cos2 = np.zeros_like(t)
225			dtau_sin2 = np.zeros_like(t)
226			tau_cs = self.taucos1 * dtau_cos1 + self.tausin1 * dtau_sin1
227			elif self.fit_phase:
228			# using time (cos) and phase (sin)
229			dtau_cos1 = np.cos(1 * self.omega * t + self.tausin1)
230			dtau_sin1 = -self.taucos1 * np.sin(1 * self.omega * t + self.tausin1)
231			dtau_cos2 = np.cos(2 * self.omega * t + self.tausin2)
232			dtau_sin2 = -self.taucos2 * np.sin(2 * self.omega * t + self.tausin2)
233			tau_cs = self.taucos1 * dtau_cos1 + self.taucos2 * dtau_cos2
234			else:
235			# using time
236			dtau_cos1 = np.cos(1 * self.omega * t)
237			dtau_sin1 = np.sin(1 * self.omega * t)
238			dtau_cos2 = np.cos(2 * self.omega * t)
239			dtau_sin2 = np.sin(2 * self.omega * t)
240			tau_cs = (self.taucos1 * dtau_cos1 + self.tausin1 * dtau_sin1 +
241			self.taucos2 * dtau_cos2 + self.tausin2 * dtau_sin2)
242			tau_cs[tau_cs < 0] = 0. # clip to zero
243			tau = self.tau0 + tau_cs
244			if self.ltscan > 0:
245			_ltscn = int(np.floor(self.ltscan))
246			else:
247			# infer the scan time from the maximal lifetime
248			_ltscn = 3 * int(np.ceil(self.tau0 +
249			np.sqrt(self.taucos12 + self.tausin12)))
250			if np.all(tau > 0):
251			bs = np.arange(1, _ltscn + 1, 1.)[None, :]
252			taufacs = np.exp(-bs / tau[:, None])
253			proxy_ts = self.intp(t[:, None] - self.lag -
254			bs / self.days_per_time_unit) * taufacs
255			proxy_val += np.sum(proxy_ts, axis=1)
256			proxy_val_grad0 += np.sum(proxy_ts * bs / tau[:, None]**2, axis=1)
257			return np.array([proxy_val,
258			# set the gradient wrt lag to zero for now
259			np.zeros_like(t),
260			self.amp * proxy_val_grad0,
261			self.amp * proxy_val_grad0 * dtau_cos1,
262			self.amp * proxy_val_grad0 * dtau_sin1,
263			self.amp * proxy_val_grad0 * dtau_cos2,
264			self.amp * proxy_val_grad0 * dtau_sin2,
265			# set the gradient wrt lifetime scan to zero for now
266			np.zeros_like(t)])
267
268	1		def _log_prior_normal(self):
269			l_prior = super(ProxyModel, self).log_prior()
270			if not np.isfinite(l_prior):
271			return -np.inf
272			for n, p in self.get_parameter_dict().items():
273			if n.startswith("tau"):
274			# Gaussian prior for the lifetimes
275			l_prior -= 0.5 * (p / self.lifetime_metric)**2
276			return l_prior
277
278	1		def _log_prior_exp(self):
279			l_prior = super(ProxyModel, self).log_prior()
280			if not np.isfinite(l_prior):
281			return -np.inf
282			for n, p in self.get_parameter_dict().items():
283			if n.startswith("tau"):
284			# exponential prior for the lifetimes
285			l_prior -= np.abs(p / self.lifetime_metric)
286			return l_prior
287
288	1		def log_prior(self):
289			_priors = {"exp": self._log_prior_exp,
290			"normal": self._log_prior_normal}
291			if self.lifetime_prior is None or self.lifetime_prior == "flat":
292			return super(ProxyModel, self).log_prior()
293			return _priors[self.lifetime_prior]()
294
295
296	1	View Code Duplication	class CeleriteModelSet(ModelSet):
			0 ignored issues – show Duplication introduced 2018-07-02 09:55 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
297
298	1		def get_value(self, t):
299			v = np.zeros_like(t)
300			for m in self.models.values():
301			v += m.get_value(t)
302			return v
303
304	1		def compute_gradient(self, t):
305			grad = []
306			for m in self.models.values():
307			grad.extend(list(m.compute_gradient(t)))
308			return np.array(grad)
309
310
311	1	View Code Duplication	def _setup_proxy_model_with_bounds(times, values,
			0 ignored issues – show Duplication introduced 2019-02-20 16:59 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
312			max_amp=1e10, max_days=100,
313			**kwargs):
314			# extract setup from `kwargs`
315			center = kwargs.get("center", False)
316			fit_phase = kwargs.get("fit_phase", False)
317			lag = kwargs.get("lag", 0.)
318			lt_metric = kwargs.get("lifetime_metric", 1)
319			lt_prior = kwargs.get("lifetime_prior", "exp")
320			lt_scan = kwargs.get("lifetime_scan", 60)
321			positive = kwargs.get("positive", False)
322			sza_intp = kwargs.get("sza_intp", None)
323			time_format = kwargs.get("time_format", "jyear")
324
325			return ProxyModel(times, values,
326			center=center,
327			sza_intp=sza_intp,
328			fit_phase=fit_phase,
329			lifetime_prior=lt_prior,
330			lifetime_metric=lt_metric,
331			days_per_time_unit=1 if time_format.endswith("d") else 365.25,
332			amp=0.,
333			lag=lag,
334			tau0=0,
335			taucos1=0, tausin1=0,
336			taucos2=0, tausin2=0,
337			ltscan=lt_scan,
338			bounds=dict([
339			("amp", [0, max_amp] if positive else [-max_amp, max_amp]),
340			("lag", [0, max_days]),
341			("tau0", [0, max_days]),
342			("taucos1", [0, max_days] if fit_phase else [-max_days, max_days]),
343			("tausin1", [-np.pi, np.pi] if fit_phase else [-max_days, max_days]),
344			# semi-annual cycles for the life time
345			("taucos2", [0, max_days] if fit_phase else [-max_days, max_days]),
346			("tausin2", [-np.pi, np.pi] if fit_phase else [-max_days, max_days]),
347			("ltscan", [0, 200])])
348			)
349

st-bender / sciapy

Push — master ( 7bc10c...b988a0 )

HarmonicModelAmpPhase.compute_gradient() A

Complexity

Size

Duplication

Code Coverage

Importance

Duplication Side-by-Side

Filter issues like