ocpy package¶
Submodules¶
ocpy.custom_types module¶
ocpy.data module¶
- class ocpy.data.Data(minimum_time, minimum_time_error=None, weights=None, minimum_type=None, labels=None)[source]¶
Bases:
DataModelContainer for eclipse minimum timing data.
The Data class stores and manages observational times of minima (e.g., eclipse timings of binary stars) together with optional metadata such as timing uncertainties, observational weights, minimum type (primary/secondary), and labels.
Internally the data are stored in a
pandas.DataFramewith standardized column names. The class provides utilities for safely manipulating the dataset, including:Filling or computing timing uncertainties and weights
Loading data from files
Computing O–C (Observed minus Calculated) values
Grouping or merging datasets
The class is designed to behave like a lightweight table object while preserving domain-specific semantics required for O–C analysis of eclipsing binaries.
- Parameters:
minimum_time (
list-like) – Times of observed minima (typically in Julian Date or BJD). This field is required.minimum_time_error (
list-like, optional) – Uncertainties associated with each minimum time. If provided, the length must matchminimum_time.weights (
list-like, optional) – Weights assigned to each observation.minimum_type (
BinarySeq, optional) – Indicator of minimum type (e.g., primary or secondary eclipse). Accepted representations may include integers, strings (e.g.,"I","II","primary","secondary"), or other binary encodings.labels (
list-like, optional) – Optional labels or identifiers for each observation.
- Raises:
ValueError – If
minimum_timeisNone.LengthCheckError – If provided sequences do not match the length of
minimum_time.
Notes
All input sequences are normalized internally so that their lengths match the length of
minimum_time. Missing values are automatically expanded or filled using utilities fromocpy.utils.The underlying storage is a
pandas.DataFramewith the following standard columns:minimum_timeminimum_time_errorweightsminimum_typelabels
Most modification methods return a new `Data` instance rather than mutating the existing object.
Examples
Create a dataset with minimum times:
>>> from ocpy import Data >>> d = Data(minimum_time=[2450000.1, 2450001.3, 2450002.5])
With uncertainties and types:
>>> d = Data( ... minimum_time=[2450000.1, 2450001.3], ... minimum_time_error=[0.0002, 0.0003], ... minimum_type=["I", "II"] ... )
Access the underlying table:
>>> d.data
- classmethod from_file(file, columns=None)[source]¶
Create a Data object from a tabular file.
This method reads timing data from a CSV or Excel file and converts it into a Data instance. Column names in the file can optionally be mapped to the standardized column names used by the Data class.
- Parameters:
file (
strorpathlib.Path) –Path to the input file. Supported formats are:
.csv.xls.xlsx
columns (
dictofstrtostr, optional) –Mapping between the standardized Data column names and the column names present in the file.
Two mapping styles are accepted:
Standard → file column name
Example:
>>> columns = {"minimum_time": "BJD", "minimum_time_error": "err"}
File column name → standard
Example:
>>> columns = {"BJD": "minimum_time", "err": "minimum_time_error"}
The standard column names recognized by the Data class are:
minimum_timeminimum_time_errorweightsminimum_typelabels
- Returns:
A new Data instance containing the imported observations.
- Return type:
- Raises:
ValueError – If the file format is not supported.
ValueError – If the file does not contain a column corresponding to
minimum_time.
Notes
The
minimum_timecolumn is required to construct a valid dataset. All other columns are optional and will be set toNoneif not present in the file.Internally, the file is loaded using
pandas.read_csv()orpandas.read_excel().Examples
Load a dataset from a CSV file:
>>> d = Data.from_file("minima.csv")
Load a file with custom column names:
>>> d = Data.from_file( ... "observations.csv", ... columns={"BJD": "minimum_time", "err": "minimum_time_error"} ... )
Load an Excel file:
>>> d = Data.from_file("minima.xlsx")
- fill_errors(errors, override=False)[source]¶
Fill or assign timing uncertainties for the dataset.
This method returns a new Data object in which the
minimum_time_errorcolumn is populated using the provided values. Existing values can optionally be preserved or replaced.- Parameters:
errors (
list,tuple,numpy.ndarray, orfloat) –Uncertainty values to assign to
minimum_time_error.If a scalar is provided, it will be applied to all rows.
If an array-like object is provided, its length must match the number of observations in the dataset.
override (
bool, defaultFalse) –If
True, all existing values in theminimum_time_errorcolumn are replaced.If
False, only entries that are currentlyNaNare filled, leaving existing values unchanged.
- Returns:
A new Data instance with updated
minimum_time_errorvalues.- Return type:
- Raises:
LengthCheckError – If
errorsis array-like and its length does not match the number of rows in the dataset.
Notes
This method does not modify the original object. Instead, it returns a new Data instance with the updated values.
Internally, column assignment is handled by the private
_assign_or_fill()method to ensure consistent behavior across similar operations.Examples
Assign a constant uncertainty to all observations:
>>> d2 = d.fill_errors(0.0002)
Fill only missing uncertainties:
>>> d2 = d.fill_errors([0.0002, 0.0003, 0.00025])
Replace all existing uncertainties:
>>> d2 = d.fill_errors(0.0002, override=True)
- fill_weights(weights, override=False)[source]¶
Fill or assign observational weights for the dataset.
This method returns a new Data object in which the
weightscolumn is populated using the provided values. Existing weights can optionally be preserved or replaced.- Parameters:
weights (
list,tuple,numpy.ndarray, orfloat) –Weight values to assign to the
weightscolumn.If a scalar is provided, it will be applied to all rows.
If an array-like object is provided, its length must match the number of observations in the dataset.
override (
bool, defaultFalse) –If
True, all existing values in theweightscolumn are replaced.If
False, only entries that are currentlyNaNare filled, leaving existing values unchanged.
- Returns:
A new Data instance with updated
weightsvalues.- Return type:
- Raises:
LengthCheckError – If
weightsis array-like and its length does not match the number of rows in the dataset.
Notes
This method does not modify the original object. Instead, it returns a new Data instance with the updated values.
Internally, column assignment is handled by the private
_assign_or_fill()method to ensure consistent behavior across similar operations.Examples
Assign a constant weight to all observations:
>>> d2 = d.fill_weights(1.0)
Fill only missing weights:
>>> d2 = d.fill_weights([1.0, 0.5, 2.0])
Replace all existing weights:
>>> d2 = d.fill_weights(1.0, override=True)
- calculate_weights(method=None, override=True)[source]¶
Calculate observational weights based on timing uncertainties.
This method computes weights for each observation, typically using the inverse-variance method, and returns a new Data instance with the updated
weightscolumn.- Parameters:
method (
callable, optional) –A custom function to compute weights from the
minimum_time_errorseries. It must accept apandas.Seriesof errors and return apandas.Seriesof weights.If
None(default), the inverse-variance method is used:\[w_i = \frac{1}{\sigma_i^2}\]where \(\sigma_i\) is the timing uncertainty for the i-th observation.
override (
bool, defaultTrue) – IfTrue, existingweightsvalues are replaced. IfFalse, only missing entries (NaN) are filled.
- Returns:
A new Data instance with updated weights.
- Return type:
- Raises:
ValueError – If
minimum_time_errorcontainsNaNvalues.ValueError – If
minimum_time_errorcontains zero values, which would cause division by zero in the default method.TypeError – If
methodis provided but is not callable.
Notes
The default inverse-variance weighting gives higher weight to observations with smaller uncertainties.
This method does not modify the original Data instance; it returns a new instance with updated weights.
Internally, column assignment uses
_assign_or_fill()to respect theoverrideflag.
Examples
Compute default inverse-variance weights:
>>> d2 = d.calculate_weights()
Compute weights with a custom method:
>>> def custom_weights(errors): ... return 1 / errors >>> d2 = d.calculate_weights(method=custom_weights, override=True)
Fill only missing weights without overwriting existing ones:
>>> d2 = d.calculate_weights(override=False)
- calculate_oc(reference_minimum, reference_period, model_type='lmfit')[source]¶
Compute Observed minus Calculated (O–C) values for the dataset.
This method calculates the O–C values for each observed minimum based on a reference minimum time and period. The O–C values quantify the difference between observed and predicted timings, which is fundamental for analyzing period variations in eclipsing binaries.
- reference_minimumfloat
The reference time of minimum (e.g., initial epoch) used to compute predicted minima.
- reference_periodfloat
The reference orbital period of the system. This is used to compute the expected timing of each cycle.
- model_typestr, default=’lmfit’
Specifies the type of O–C model to return. Supported options:
'lmfit'or'lmfit_model'– returns anOCLMFitobject.'pymc'or'pymc_model'– returns anOCPyMCobject.Any other string – returns a generic
OCobject.
- OC
An instance of the appropriate O–C model class (
OC,OCLMFit, orOCPyMC) containing:minimum_time– observed minimacycle– computed cycle numbers (integer or half-integer for secondary minima)oc– O–C valuesAdditional columns from the original Data (errors, weights, labels, minimum_type)
- ValueError
If the
minimum_timecolumn is missing.
Cycle calculation: The phase of each observation is computed as:
\[ext{phase} =\]
rac{t - ext{reference_minimum}}{ ext{reference_period}}
The cycle number is the nearest integer to the phase.
Secondary minima: If
minimum_typeis present and indicates a secondary eclipse (e.g., “II”, “secondary”, 2), the cycle is adjusted to half-integer values:\[ext{cycle}_{ ext{sec}} = ext{round}( ext{phase} - 0.5) + 0.5\]O–C computation: The O–C value for each observation is:
\[ext{O–C} = t_{ ext{obs}} - ( ext{reference_minimum} + ext{cycle} imes ext{reference_period} )\]The method does not modify the original `Data`. The returned object contains a copy of the original data along with the computed
cycleandocarrays.The
model_typedetermines which class is instantiated for further modeling of the O–C diagram.
Compute O–C values using the default LMFit model:
>>> oc_model = d.calculate_oc(reference_minimum=2450000.0, reference_period=1.2345)
Compute O–C values using a PyMC model:
>>> oc_model = d.calculate_oc( ... reference_minimum=2450000.0, ... reference_period=1.2345, ... model_type="pymc" ... )
Access the computed O–C values:
>>> oc_model.oc [0.0001, -0.0002, 0.0003, ...]
- merge(data)[source]¶
Merge the current dataset with another Data object.
This method concatenates the rows of the current Data instance with those of another Data object, returning a new Data instance. Column alignment is based on column names.
- Parameters:
data (
Data) – Another Data instance to merge with the current dataset.- Returns:
A new Data instance containing all rows from both datasets.
- Return type:
Notes
The original datasets are not modified.
Missing columns in either dataset will result in
NaNvalues in the merged dataset, following pandas’ concatenation rules.Indexes are reset in the merged dataset for consistency.
Examples
>>> d1 = Data(minimum_time=[2450000.1, 2450001.2]) >>> d2 = Data(minimum_time=[2450002.3, 2450003.4]) >>> d_merged = d1.merge(d2) >>> len(d_merged) 4
- group_by(column)[source]¶
Split the dataset into groups based on a column.
This method groups the Data object by the values in a specified column and returns a list of new Data instances, each containing one group of rows.
- Parameters:
column (
str) – Name of the column to group by.- Returns:
A list of Data objects, each corresponding to one group. If the column is missing or contains only NaN values, a list with a single copy of the original dataset is returned.
- Return type:
listofData
Notes
Grouping is performed using
pandas.DataFrame.groupby().The original Data object is not modified; each group is a deep copy.
NaN values are treated as a separate group unless
dropna=Truein the internal pandas grouping.
Examples
Group a dataset by the
minimum_typecolumn:>>> groups = d.group_by("minimum_type") >>> len(groups) 2 # e.g., one group for primary, one for secondary minima
Access the first group:
>>> groups[0].data minimum_time minimum_type 2450000.1 I 2450002.3 I
If the grouping column does not exist:
>>> groups = d.group_by("nonexistent_column") >>> len(groups) 1 # returns a single copy of the original Data
ocpy.errors module¶
ocpy.model_data module¶
- class ocpy.model_data.DataModel(minimum_time, minimum_time_error=None, weights=None, minimum_type=None, labels=None, ecorr=None, oc=None)[source]¶
Bases:
ABC- abstractmethod calculate_weights(method=None, override=True)[source]¶
Calculates weights using errors
ocpy.model_oc module¶
- class ocpy.model_oc.ParameterModel[source]¶
Bases:
ABC- value: float¶
- min: float¶
- max: float¶
- std: float¶
- fixed: bool = False¶
ocpy.newtonian module¶
- class ocpy.newtonian.NewtonianModel(*, integrator='ias15', dt=0.01, integrator_params=None, units=None, reference_time=0.0, t_start=None, t_end=None, stop_at_exact_time=True, escape_radius=100000.0, min_distance=0.0001, precision_integration_steps=0, integration_grid=None, central_mass=1.0, bodies=None, orbit_type='heliocentric', orbit_output_type='heliocentric', T0_ref=0.0, P_ref=1.0, compute_xyz=True, compute_orbital=True, name=None)[source]¶
Bases:
ModelComponent- name = 'newtonian'¶
- class ocpy.newtonian.NewtonianOp(model)[source]¶
Bases:
Op- make_node(x, *args)[source]¶
Construct an Apply node that represent the application of this operation to the given inputs.
This must be implemented by sub-classes.
- Returns:
node – The constructed Apply node.
- Return type:
Apply
- perform(node, inputs, outputs)[source]¶
Calculate the function on the inputs and put the variables in the output storage.
- Parameters:
node – The symbolic Apply node that represents this computation.
inputs – Immutable sequence of non-symbolic/numeric inputs. These are the values of each Variable in
node.inputs.output_storage – List of mutable single-element lists (do not change the length of these lists). Each sub-list corresponds to value of each Variable in
node.outputs. The primary purpose of this method is to set the values of these sub-lists.
Notes
The output_storage list might contain data. If an element of output_storage is not
None, it has to be of the right type, for instance, for a TensorVariable, it has to be a NumPyndarraywith the right number of dimensions and the correct dtype. Its shape and stride pattern can be arbitrary. It is not guaranteed that such pre-set values were produced by a previous call to thisOp.perform(); they could’ve been allocated by another Op’s perform method. An Op is free to reuse output_storage as it sees fit, or to discard it and allocate new memory.
- grad(inputs, output_gradients)[source]¶
Construct a graph for the gradient with respect to each input variable.
Each returned Variable represents the gradient with respect to that input computed based on the symbolic gradients with respect to each output. If the output is not differentiable with respect to an input, then this method should return an instance of type NullType for that input.
Using the reverse-mode AD characterization given in [1], for a \(C = f(A, B)\) representing the function implemented by the Op and its two arguments \(A\) and \(B\), given by the Variables in inputs, the values returned by Op.grad represent the quantities \(\bar{A} \equiv \frac{\partial S_O}{A}\) and \(\bar{B}\), for some scalar output term \(S_O\) of \(C\) in
\[\operatorname{Tr}\left(\bar{C}^\top dC\right) = \operatorname{Tr}\left(\bar{A}^\top dA\right) + \operatorname{Tr}\left(\bar{B}^\top dB\right)\]- Parameters:
inputs – The input variables.
output_grads – The gradients of the output variables.
- Returns:
The gradients with respect to each Variable in inputs.
- Return type:
grads
References
ocpy.oc module¶
- class ocpy.oc.Parameter(value=None, min=None, max=None, std=None, fixed=False, distribution='truncatednormal')[source]¶
Bases:
ParameterModelRepresents a model parameter with optional bounds, uncertainty, and distribution.
- Parameters:
value (
float, optional) – The current or initial value of the parameter.min (
float, optional) – Lower bound for the parameter, if applicable.max (
float, optional) – Upper bound for the parameter, if applicable.std (
float, optional) – Standard deviation of the parameter, used for uncertainty.fixed (
bool, defaultFalse) – Whether the parameter is fixed (not free to vary) during fitting.distribution (
str, default'truncatednormal') – The assumed probability distribution for Bayesian sampling (e.g., ‘truncatednormal’, ‘uniform’, ‘normal’).
Notes
This class inherits from
ParameterModel.Can be used to define parameters for any
ModelComponent.Supports integration with fitting routines and Bayesian inference.
Examples
Define a free parameter with initial value 1.0:
>>> p = Parameter(value=1.0)
Define a bounded parameter with uncertainty:
>>> p = Parameter(value=0.5, min=0.0, max=1.0, std=0.1)
- value: float | None = None¶
- min: float | None = None¶
- max: float | None = None¶
- std: float | None = None¶
- fixed: bool | None = False¶
- distribution: str = 'truncatednormal'¶
- class ocpy.oc.ModelComponent(args=None)[source]¶
Bases:
ModelComponentModelBase class for a mathematical model component with named parameters.
This class provides the foundation for any model component used in O–C analysis or other curve-fitting tasks. It manages parameters, allows updating from inference data, and supports switching between numerical backends such as NumPy or PyMC.
- params¶
A dictionary of parameter names to
Parameterobjects. These parameters define the model and can be updated during fitting.- Type:
dictofstr -> Parameter
- math_class¶
Numerical backend used for mathematical operations (default: numpy).
- Type:
module
- _atan2¶
Function used for two-argument arctangent, adjusted to the current math backend.
- Type:
callable
- update_from_idata(inference_data, group='posterior', stat='median')[source]¶
Update parameters from a PyMC
InferenceDataobject.
Notes
Intended to be subclassed by concrete model components such as
Linear,Sinusoidal, orKeplerian.Supports both classical and Bayesian fitting workflows.
Provides compatibility with different computational backends, including NumPy and PyMC (if available).
Parameter values can be updated manually or via inference results.
Examples
Create a simple linear model component:
>>> linear = Linear(a=1.0, b=0.0) >>> linear.params {'a': Parameter(value=1.0), 'b': Parameter(value=0.0)}
Switch to a PyMC math backend (if PyMC installed):
>>> linear.set_math(pm.math)
Update parameters manually:
>>> linear.update_parameters({'a': 2.0}) >>> linear.params['a'].value 2.0
- math_class = <module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/oc-py/envs/latest/lib/python3.12/site-packages/numpy/__init__.py'>¶
- set_math(mathmod)[source]¶
Set the mathematical backend for model computations.
This method allows switching the numerical backend used in all mathematical operations of the component, including functions like sin, cos, sqrt, and arctan2. Supports NumPy by default and PyMC if available.
- Parameters:
mathmod (
module) – The module providing mathematical functions (e.g., numpy, pymc.math, or pytensor.tensor).- Returns:
self – Returns the instance itself to allow method chaining.
- Return type:
Notes
The _atan2 function is updated to match the chosen backend.
If PyMC is available and mathmod corresponds to PyMC’s math module, _atan2 uses pm.math.arctan2 or pytensor.tensor.arctan2.
If the backend does not provide arctan2, it falls back to numpy.arctan2.
This is useful for switching between classical numerical computation and probabilistic programming frameworks.
Examples
Use NumPy as the backend (default):
>>> component.set_math(np)
Use PyMC backend if installed:
>>> import pymc as pm >>> component.set_math(pm.math)
- model_function()[source]¶
Return the underlying model function of the component.
The model function defines the mathematical relationship between the independent variable(s) and the output, based on the component’s parameters. It is intended to be used for evaluation, fitting, or plotting.
- Returns:
The function implementing the model. The callable typically has the signature:
f(x, param1, param2, ...)where x is the independent variable and the remaining arguments are the component’s parameters.
- Return type:
callable
Notes
Subclasses must define the model_func method, which this property returns.
Useful for passing to optimization routines or plotting functions that require a callable.
Examples
Retrieve and call the model function of a linear component:
>>> linear = Linear(a=2.0, b=1.0) >>> f = linear.model_function() >>> f(3, linear.params['a'].value, linear.params['b'].value) 7.0
- update_parameters(params_dict)[source]¶
Update the values of existing parameters from a dictionary.
- Parameters:
params_dict (
dictofstr -> float) – Dictionary mapping parameter names to their new numeric values. Only parameters that already exist in the component will be updated; any unknown keys are ignored.- Returns:
self – Returns the instance itself to allow method chaining.
- Return type:
Notes
This method does not create new parameters; it only updates values of existing ones.
All values are cast to float before assignment.
Useful for manually adjusting parameters between fitting steps or after loading results from inference.
Examples
Update a linear model component’s parameters:
>>> linear = Linear(a=1.0, b=0.0) >>> linear.update_parameters({'a': 2.5, 'b': -1.0}) >>> linear.params['a'].value 2.5 >>> linear.params['b'].value -1.0
Attempting to update a non-existent parameter has no effect:
>>> linear.update_parameters({'c': 10.0}) >>> 'c' in linear.params False
- update_from_idata(inference_data, group='posterior', stat='median')[source]¶
Update model parameters from a PyMC InferenceData object.
This method extracts parameter values from a Bayesian posterior (or other group) and updates the component’s parameters accordingly. It automatically matches variables that are prefixed with the component’s name.
- Parameters:
inference_data (
arviz.InferenceData) – The InferenceData object containing posterior samples from a PyMC model.group (
str, default'posterior') – The group within the InferenceData to use (e.g., ‘posterior’, ‘prior’).stat (
str, default'median') – Statistic to use for updating parameters. Options include: - ‘median’: use the median of the posterior samples. - ‘mean’: use the mean of the posterior samples. - Any other value: uses the first sample as a fallback.
- Returns:
self – The component instance with updated parameter values.
- Return type:
Notes
Only parameters present in self.params and matching the component’s prefix are updated.
Useful for propagating posterior estimates back into model components for further computation or plotting.
Requires PyMC and ArviZ to be installed.
Examples
Suppose idata contains posterior samples for a linear component named ‘linear’:
>>> linear = Linear(a=1.0, b=0.0) >>> linear.update_from_idata(idata, group='posterior', stat='median') >>> linear.params['a'].value # updated from posterior median 2.3 >>> linear.params['b'].value -0.1
- class ocpy.oc.Linear(a=1.0, b=0.0, *, name=None)[source]¶
Bases:
ModelComponentLinear model component of the form f(x) = a * x + b.
Represents a simple linear relationship between an independent variable x and the output, parameterized by slope a and intercept b.
- Parameters:
a (
floatorParameter, default1.0) – Slope of the linear model. Can be a numeric value or a Parameter object.b (
floatorParameter, default0.0) – Intercept of the linear model. Can be a numeric value or a Parameter object.name (
str, optional) – Optional name for the component. Defaults to “linear”.
- params¶
Dictionary of parameters: {‘a’: Parameter, ‘b’: Parameter}.
- Type:
dict
- name¶
Name of the component (default: “linear”).
- Type:
str
Examples
Create a linear component with default parameters:
>>> linear = Linear() >>> linear.model_func(3, linear.params['a'].value, linear.params['b'].value) 3.0
Create a linear component with custom slope and intercept:
>>> linear = Linear(a=2.0, b=1.0) >>> linear.model_func(3, linear.params['a'].value, linear.params['b'].value) 7.0
- name = 'linear'¶
- model_func(x, a, b)[source]¶
Evaluate the linear model f(x) = a * x + b.
- Parameters:
x (
floatorarray-like) – Independent variable(s) at which to evaluate the model.a (
float) – Slope of the linear function.b (
float) – Intercept of the linear function.
- Returns:
The computed value(s) of the linear function at x.
- Return type:
floatornp.ndarray
Examples
Evaluate a linear component at a single point:
>>> linear = Linear(a=2.0, b=1.0) >>> linear.model_func(3, linear.params['a'].value, linear.params['b'].value) 7.0
Evaluate at multiple points:
>>> import numpy as np >>> x = np.array([0, 1, 2, 3]) >>> linear.model_func(x, linear.params['a'].value, linear.params['b'].value) array([1., 3., 5., 7.])
- class ocpy.oc.Quadratic(q=0.0, *, name=None)[source]¶
Bases:
ModelComponentQuadratic model component of the form f(x) = q * x^2.
Represents a simple quadratic relationship between an independent variable x and the output, parameterized by a single coefficient q.
- Parameters:
q (
floatorParameter, default0.0) – Quadratic coefficient. Can be a numeric value or a Parameter object.name (
str, optional) – Optional name for the component. Defaults to “quadratic”.
- params¶
Dictionary of parameters: {‘q’: Parameter}.
- Type:
dict
- name¶
Name of the component (default: “quadratic”).
- Type:
str
Examples
Create a quadratic component with default parameter:
>>> quad = Quadratic() >>> quad.model_func(2, quad.params['q'].value) 0.0
Create a quadratic component with custom coefficient:
>>> quad = Quadratic(q=3.0) >>> quad.model_func(2, quad.params['q'].value) 12.0
- name = 'quadratic'¶
- model_func(x, q)[source]¶
Evaluate the quadratic model f(x) = q * x^2.
- Parameters:
x (
floatorarray-like) – Independent variable(s) at which to evaluate the model.q (
float) – Coefficient of the quadratic term.
- Returns:
The computed value(s) of the quadratic function at x.
- Return type:
floatornp.ndarray
Examples
Evaluate a quadratic component at a single point:
>>> quad = Quadratic(q=2.0) >>> quad.model_func(3, quad.params['q'].value) 18.0
Evaluate at multiple points:
>>> import numpy as np >>> x = np.array([0, 1, 2, 3]) >>> quad.model_func(x, quad.params['q'].value) array([ 0., 2., 8., 18.])
- class ocpy.oc.Sinusoidal(*, amp=None, P=None, name=None)[source]¶
Bases:
ModelComponentSinusoidal model component of the form f(x) = amp * sin(2π * x / P).
Represents a simple periodic function with amplitude amp and period P.
- Parameters:
amp (
floatorParameter, optional) – Amplitude of the sinusoidal function. Can be a numeric value or a Parameter object.P (
floatorParameter, optional) – Period of the sinusoidal function. Can be a numeric value or a Parameter object.name (
str, optional) – Optional name for the component. Defaults to “sinusoidal”.
- params¶
Dictionary of parameters: {‘amp’: Parameter, ‘P’: Parameter}.
- Type:
dict
- name¶
Name of the component (default: “sinusoidal”).
- Type:
str
- model_func(x, amp, P)[source]¶
Evaluate the sinusoidal model for input x using parameters amp and P.
Examples
Create a sinusoidal component with amplitude 1.0 and period 2.0:
>>> sinus = Sinusoidal(amp=1.0, P=2.0) >>> sinus.model_func(0, sinus.params['amp'].value, sinus.params['P'].value) 0.0
Evaluate at multiple points:
>>> import numpy as np >>> x = np.array([0, 0.5, 1.0, 1.5]) >>> sinus.model_func(x, sinus.params['amp'].value, sinus.params['P'].value) array([0. , 1. , 0. , -1.])
- name = 'sinusoidal'¶
- model_func(x, amp, P)[source]¶
Evaluate the sinusoidal model f(x) = amp * sin(2π * x / P).
- Parameters:
x (
floatorarray-like) – Independent variable(s) at which to evaluate the model.amp (
float) – Amplitude of the sinusoidal function.P (
float) – Period of the sinusoidal function.
- Returns:
The computed value(s) of the sinusoidal function at x.
- Return type:
floatornp.ndarray
Examples
Evaluate at a single point:
>>> sinus = Sinusoidal(amp=2.0, P=4.0) >>> sinus.model_func(1, sinus.params['amp'].value, sinus.params['P'].value) 2.0
Evaluate at multiple points:
>>> import numpy as np >>> x = np.array([0, 1, 2, 3, 4]) >>> sinus.model_func(x, sinus.params['amp'].value, sinus.params['P'].value) array([ 0., 2., 0., -2., -0.])
- class ocpy.oc.Keplerian(*, amp=None, e=0.0, omega=0.0, P=None, T0=None, name=None)[source]¶
Bases:
ModelComponentKeplerian orbital model component for radial velocity variations.
Represents a Keplerian orbit with parameters corresponding to the orbital amplitude, eccentricity, argument of periastron, orbital period, and time of periastron passage.
The model computes radial velocity or timing variations as a function of time based on classical Keplerian motion.
- Parameters:
amp (
floatorParameter, optional) – Amplitude of the Keplerian signal.e (
floatorParameter, default0.0) – Orbital eccentricity. Must be between 0 and 1.omega (
floatorParameter, default0.0) – Argument of periastron in degrees.P (
floatorParameter, optional) – Orbital period.T0 (
floatorParameter, optional) – Time of periastron passage.name (
str, optional) – Optional name for the component. Defaults to “keplerian”.
- params¶
Dictionary of parameters: {‘amp’, ‘e’, ‘omega’, ‘P’, ‘T0’}.
- Type:
dict
- name¶
Name of the component (default: “keplerian”).
- Type:
str
Notes
The input angles are in degrees and internally converted to radians.
Uses Newton-Raphson iteration to solve Kepler’s equation.
Suitable for modeling radial velocity or timing variations in binary stars or exoplanets.
Examples
Create a Keplerian component with default parameters:
>>> kepler = Keplerian() >>> kepler.params['e'].value 0.0
Custom parameters:
>>> kepler = Keplerian(amp=2.0, e=0.5, omega=45.0, P=10.0, T0=0.0) >>> kepler.params['amp'].value 2.0 >>> kepler.params['P'].value 10.0
- name = 'keplerian'¶
- model_func(x, amp, e, omega, P, T0)[source]¶
Evaluate the Keplerian model at given times.
Computes the Keplerian signal (e.g., radial velocity or timing variations) based on orbital parameters: amplitude, eccentricity, argument of periastron, period, and time of periastron passage.
- Parameters:
x (
floatorarray-like) – Times at which to evaluate the model.amp (
float) – Amplitude of the Keplerian signal.e (
float) – Orbital eccentricity (0 ≤ e < 1).omega (
float) – Argument of periastron in degrees.P (
float) – Orbital period.T0 (
float) – Time of periastron passage.
- Returns:
Keplerian model values at the specified times.
- Return type:
floatornp.ndarray
Notes
The argument of periastron omega is internally converted to radians.
Uses _kepler_solve to determine the eccentric anomaly.
Handles both scalar and array inputs for x.
Suitable for modeling timing variations in eclipsing binaries or radial velocity signals.
Examples
Evaluate at a single time:
>>> kepler = Keplerian(amp=1.0, e=0.1, omega=30.0, P=10.0, T0=0.0) >>> kepler.model_func(5, 1.0, 0.1, 30.0, 10.0, 0.0) 0.85 # Example value (approximate)
Evaluate at multiple times:
>>> import numpy as np >>> times = np.linspace(0, 10, 5) >>> kepler.model_func(times, 1.0, 0.1, 30.0, 10.0, 0.0) array([0.0, 0.52, 0.85, 0.85, 0.52])
- class ocpy.oc.KeplerianOld(*, amp=None, e=0.0, omega=0.0, P=None, T0=None, name=None)[source]¶
Bases:
ModelComponent- name = 'keplerian'¶
- class ocpy.oc.OC(oc, minimum_time=None, minimum_time_error=None, weights=None, minimum_type=None, labels=None, cycle=None)[source]¶
Bases:
OCModelObserved-minus-Calculated (O–C) data container and analysis class.
Represents timing variations of observed events (e.g., eclipses, pulsations) relative to a reference ephemeris. Provides methods for binning, merging, computing O–C values, and fitting models to timing variations.
- data¶
Internal DataFrame storing columns: ‘minimum_time’, ‘minimum_time_error’, ‘weights’, ‘minimum_type’, ‘labels’, ‘cycle’, ‘oc’.
- Type:
pd.DataFrame
- bin(bin_count=1, bin_method=None, bin_error_method=None, bin_style=None)[source]¶
Bin O–C data using weighted averages.
- calculate_oc(reference_minimum, reference_period, model_type='lmfit')[source]¶
Compute O–C values relative to a reference ephemeris.
- plot(model=None, ax=None, res_ax=None, res=True, ...)[source]¶
Visualize O–C data and optional model fits.
Notes
Supports both scalar and array-like input for all columns.
cycle is automatically computed or can be provided explicitly.
Binning can be performed with custom reducers and bin styles.
Integration with ModelComponent subclasses enables model fitting.
Examples
Create an OC instance from raw data:
>>> oc = OC(oc=[0.1, -0.05, 0.0], ... minimum_time=[2450000.0, 2450001.0, 2450002.0], ... weights=[1.0, 1.0, 1.0])
Access O–C column:
>>> oc['oc'] 0 0.10 1 -0.05 2 0.00 Name: oc, dtype: float64
Bin O–C data:
>>> binned = oc.bin(bin_count=2)
Merge with another OC:
>>> merged = oc.merge(other_oc)
- classmethod from_file(file, columns=None)[source]¶
Load O–C data from a CSV or Excel file into an OC instance.
- Parameters:
file (
strorPath) – Path to the input file. Supported formats: .csv, .xls, .xlsx.columns (
dict, optional) – Mapping of file column names to OC attribute names. Keys are original column names in the file; values are the corresponding attribute names (minimum_time, oc, weights, cycle, etc.).
- Returns:
An OC instance containing data loaded from the file.
- Return type:
- Raises:
ValueError – If the file type is unsupported.
Notes
Any columns not present in the file are filled with None.
Column renaming allows flexibility for files with nonstandard column names.
The method automatically handles array-like columns and converts them for internal storage in a pandas.DataFrame.
Examples
Load a CSV file with standard columns:
>>> oc = OC.from_file("observations.csv")
Load an Excel file with custom column names:
>>> oc = OC.from_file("observations.xlsx", ... columns={"Time": "minimum_time", ... "O-C": "oc", ... "Weight": "weights"})
- bin(bin_count=1, bin_method=None, bin_error_method=None, bin_style=None)[source]¶
Bin the O–C data along the cycle axis and compute weighted averages and errors.
- Parameters:
bin_count (
int, optional) – Number of bins to create. Default is 1.bin_method (
callable, optional) – Function to compute binned values from array and weights: func(array: np.ndarray, weights: np.ndarray) -> float. Default is weighted mean.bin_error_method (
callable, optional) – Function to compute binned errors from weights: func(weights: np.ndarray) -> float. Default is 1 / sqrt(sum(weights)).bin_style (
callable, optional) – Function to define custom bin edges: func(df: pd.DataFrame, bin_count: int) -> np.ndarray of shape (N, 2).
- Returns:
New OC instance containing binned data.
- Return type:
- Raises:
ValueError – If cycle, oc, or weights columns are missing or contain NaNs. If bin_count is not positive.
Notes
Uses cycle as the x-axis and oc as the y-axis.
By default, bins are equal-width and weighted averages are used.
Supports custom binning logic via bin_style.
- Binned rows have:
labels set to “Binned”
minimum_time and weights set to NaN
minimum_type set to None
Examples
Bin the data into 5 equal-width bins:
>>> binned_oc = oc.bin(bin_count=5)
Bin using a custom error function and custom bin edges:
>>> def custom_bin_style(df, n_bins): ... return np.array([[0, 2], [2, 5], [5, 10]]) >>> binned_oc = oc.bin(bin_count=3, bin_error_method=lambda w: np.sqrt(np.sum(w)), bin_style=custom_bin_style)
- merge(oc)[source]¶
Merge another OC instance into the current one.
- Parameters:
oc (
OC) – Another OC instance whose data will be appended.- Returns:
A new OC instance containing the combined data of both instances.
- Return type:
Notes
The merge is performed row-wise, preserving all columns.
Indexes are reset in the resulting DataFrame.
The original OC instances are not modified; the merge returns a new instance.
Examples
Merge two OC datasets:
>>> oc1 = OC.from_file("observations1.csv") >>> oc2 = OC.from_file("observations2.csv") >>> merged = oc1.merge(oc2) >>> len(merged) == len(oc1) + len(oc2) True
- calculate_oc(reference_minimum, reference_period, model_type='lmfit')[source]¶
Compute O–C (Observed minus Calculated) values based on a reference ephemeris.
- Parameters:
reference_minimum (
float) – Reference time of a primary minimum (e.g., first observed eclipse).reference_period (
float) – Reference period for the cycle calculation.model_type (
str, optional) – Type of model to return: - “lmfit” : returns an OCLMFit instance if available. - other : returns an OC instance. Default is “lmfit”.
- Returns:
New OC (or OCLMFit) instance with calculated cycle and oc columns.
- Return type:
- Raises:
ValueError – If the minimum_time column is missing. If minimum_type contains unrecognized values for secondary/primary distinction.
Notes
The method computes the cycle number as (t - reference_minimum) / reference_period, rounded to the nearest integer.
If minimum_type is present, secondary minima (e.g., “II”, “sec”) are offset by 0.5 cycles.
Calculated O–C values are added as the oc column; cycle numbers as cycle.
Other columns (minimum_time_error, weights, labels, etc.) are preserved.
Examples
Compute O–C for an OC dataset:
>>> oc = OC.from_file("observations.csv") >>> oc_calc = oc.calculate_oc(reference_minimum=2450000.0, reference_period=1.2345) >>> oc_calc.data[['cycle', 'oc']].head() cycle oc 0 0.0 0.012 1 1.0 -0.005 2 2.0 0.003
- residue(coefficients)[source]¶
Compute the residuals of the O–C data relative to a fitted model.
- Parameters:
coefficients (
ModelResult) – Fitted model object (from lmfit) containing the optimized parameters used to compute model-predicted O–C values.- Returns:
New OC instance containing the residuals (observed - model) in the oc column. Other columns are preserved from the original instance.
- Return type:
Notes
This method does not modify the original OC instance; a new instance is returned.
Residuals are computed using the model defined by coefficients applied to the cycle column.
Typically used after OC.fit or similar model fitting methods.
Examples
Compute residuals after fitting a linear model:
>>> result = oc.fit_linear(a=1.0, b=0.0) # returns a ModelResult >>> oc_resid = oc.residue(result) >>> oc_resid.data[['cycle', 'oc']].head() cycle oc 0 0.0 0.002 1 1.0 -0.001 2 2.0 0.003
- fit(functions)[source]¶
Fit one or more model components to the O–C data.
- Parameters:
functions (
ModelComponentModelorlistofModelComponentModel) – Single or multiple model components (e.g., Linear, Quadratic, Sinusoidal, Keplerian) to fit to the oc data. Each component defines its own functional form and parameters.- Returns:
Fitted model result (from lmfit) containing optimized parameters, fit statistics, and the model-predicted O–C values.
- Return type:
ModelResult- Raises:
ValueError – If the required columns (cycle, oc, or weights) are missing or contain NaNs.
TypeError – If functions is not a ModelComponentModel or list thereof.
Notes
Supports combining multiple model components into a single composite fit.
Uses weighted fitting based on the weights column.
The resulting ModelResult can be used to evaluate residuals or plot the fitted model.
Internally, converts functions to a composite lmfit model if multiple components are provided.
Examples
Fit a linear and sinusoidal model to O–C data:
>>> linear = Linear(a=0.0, b=0.0) >>> sinus = Sinusoidal(amp=0.01, P=10.0) >>> result = oc.fit([linear, sinus]) >>> result.best_values {'a': 0.0012, 'b': -0.0023, 'amp': 0.0098, 'P': 9.97}
- fit_keplerian(*, amp=None, e=None, omega=None, P=None, T=None)[source]¶
Fit a Keplerian (elliptical orbit) model to the O–C data.
- Parameters:
amp (
ParameterModel, optional) – Amplitude of the Keplerian signal.e (
ParameterModel, optional) – Orbital eccentricity.omega (
ParameterModel, optional) – Argument of periastron in degrees.P (
ParameterModel, optional) – Orbital period.T (
ParameterModel, optional) – Time of periastron passage.
- Returns:
Fitted Keplerian model component with optimized parameters.
- Return type:
ModelComponentModel
Notes
Uses the Keplerian formula to model periodic variations in O–C data.
Supports optional ParameterModel objects for each orbital parameter; unspecified parameters can be estimated from the data.
The returned model component can be used to evaluate predicted O–C values or for further combination with other model components.
Examples
Fit a Keplerian model with initial guesses:
>>> kepler_model = oc.fit_keplerian( ... amp=Parameter(value=0.01), ... e=Parameter(value=0.1), ... omega=Parameter(value=90), ... P=Parameter(value=1000), ... T=Parameter(value=2450000) ... ) >>> kepler_model.params['amp'].value 0.0098
- fit_lite(*, amp=None, e=None, omega=None, P=None, T=None)[source]¶
Fit a simplified Keplerian model to the O–C data.
- Parameters:
amp (
ParameterModel, optional) – Amplitude of the Keplerian signal.e (
ParameterModel, optional) – Orbital eccentricity.omega (
ParameterModel, optional) – Argument of periastron in degrees.P (
ParameterModel, optional) – Orbital period.T (
ParameterModel, optional) – Time of periastron passage.
- Returns:
Fitted simplified Keplerian model component with optimized parameters.
- Return type:
ModelComponentModel
Notes
This method provides a lighter or approximate version of a full Keplerian fit.
Useful for initial estimates or for datasets where full Keplerian fitting is unnecessary.
Supports optional ParameterModel objects for each orbital parameter; unspecified parameters can be estimated from the data.
The returned model component can be evaluated for predicted O–C values or combined with other model components.
Examples
Fit a simplified Keplerian model:
>>> kepler_lite = oc.fit_lite( ... amp=Parameter(value=0.01), ... e=Parameter(value=0.05), ... P=Parameter(value=1000), ... ) >>> kepler_lite.params['P'].value 1000.0
- fit_linear(*, a=None, b=None)[source]¶
Fit a linear model to the O–C data.
- Parameters:
a (
ParameterModel, optional) – Slope of the linear model. If not provided, it will be estimated from the data.b (
ParameterModel, optional) – Intercept of the linear model. If not provided, it will be estimated from the data.
- Returns:
Fitted linear model component with optimized parameters.
- Return type:
ModelComponentModel
Notes
The linear model has the form oc = a * cycle + b.
Optional ParameterModel objects can be used to fix or initialize parameters.
The returned model component can be used to compute predicted O–C values or combined with other model components for composite fitting.
Examples
Fit a linear trend to O–C data:
>>> linear_model = oc.fit_linear(a=Parameter(value=0.001), b=Parameter(value=-0.002)) >>> linear_model.params['a'].value 0.001 >>> linear_model.params['b'].value -0.002
- fit_quadratic(*, q=None)[source]¶
Fit a quadratic model to the O–C data.
- Parameters:
q (
ParameterModel, optional) – Quadratic coefficient. If not provided, it will be estimated from the data.- Returns:
Fitted quadratic model component with optimized parameter(s).
- Return type:
ModelComponentModel
Notes
The quadratic model has the form oc = q * cycle^2.
Optional ParameterModel can be used to fix or initialize the quadratic coefficient.
The returned model component can be used to compute predicted O–C values or combined with other model components for composite fitting.
Examples
Fit a quadratic trend to O–C data:
>>> quad_model = oc.fit_quadratic(q=Parameter(value=0.0001)) >>> quad_model.params['q'].value 0.0001
- fit_sinusoidal(*, amp=None, P=None)[source]¶
Fit a sinusoidal model to the O–C data.
- Parameters:
amp (
ParameterModel, optional) – Amplitude of the sinusoidal signal. If not provided, it will be estimated from the data.P (
ParameterModel, optional) – Period of the sinusoidal signal. If not provided, it will be estimated from the data.
- Returns:
Fitted sinusoidal model component with optimized parameters.
- Return type:
ModelComponentModel
Notes
The sinusoidal model has the form oc = amp * sin(2 * pi * cycle / P).
Optional ParameterModel objects can be used to fix or initialize parameters.
The returned model component can be evaluated for predicted O–C values or combined with other model components in a composite fit.
Examples
Fit a sinusoidal variation to O–C data:
>>> sinus_model = oc.fit_sinusoidal(amp=Parameter(value=0.01), P=Parameter(value=1000)) >>> sinus_model.params['amp'].value 0.01 >>> sinus_model.params['P'].value 1000
- fit_parabola(*, q=None, a=None, b=None)[source]¶
Fit a parabolic model to the O–C data.
- Parameters:
q (
ParameterModel, optional) – Quadratic coefficient. If not provided, it will be estimated from the data.a (
ParameterModel, optional) – Linear coefficient. If not provided, it will be estimated from the data.b (
ParameterModel, optional) – Constant term (intercept). If not provided, it will be estimated from the data.
- Returns:
Fitted parabolic model component with optimized parameters.
- Return type:
ModelComponentModel
Notes
The parabolic model has the form oc = q * cycle^2 + a * cycle + b.
Optional ParameterModel objects can be used to fix or initialize coefficients.
The returned model component can be used to compute predicted O–C values or combined with other model components for composite fitting.
Examples
Fit a parabolic trend to O–C data:
>>> parab_model = oc.fit_parabola(q=Parameter(value=0.0001), ... a=Parameter(value=0.001), ... b=Parameter(value=-0.002)) >>> parab_model.params['q'].value 0.0001 >>> parab_model.params['a'].value 0.001 >>> parab_model.params['b'].value -0.002
- plot(model=None, *, ax=None, res_ax=None, res=True, title=None, x_col='cycle', y_col='oc', fig_size=(10, 7), plot_kwargs=None, extension_factor=0.05, model_components=None)[source]¶
Plot the O–C data along with optional model predictions and residuals.
- Parameters:
model (
InferenceData,ModelResult, orlistofModelComponent, optional) – Model or model components to overlay on the O–C plot.ax (
matplotlib.axes.Axes, optional) – Axes object for the main plot. If None, a new figure and axes are created.res_ax (
matplotlib.axes.Axes, optional) – Axes object for residuals plot. Only used if res=True.res (
bool, defaultTrue) – If True, show residuals between the data and model on a separate subplot.title (
str, optional) – Title for the plot.x_col (
str, default"cycle") – Column name in the data to use for the x-axis.y_col (
str, default"oc") – Column name in the data to use for the y-axis.fig_size (
tuple, default(10,7)) – Figure size if a new figure is created.plot_kwargs (
dict, optional) – Additional keyword arguments passed to the plotting function (e.g., marker, linestyle, color).extension_factor (
float, default0.05) – Fractional extension for axis limits to improve visual spacing.model_components (
list, optional) – List of individual model components to display separately on the plot.
- Returns:
If res=False, returns the main Axes object. If res=True, returns a dictionary with keys ‘main_ax’ and ‘res_ax’ containing the corresponding Axes objects.
- Return type:
matplotlib.axes.Axesordict
Notes
This method uses the Plot.plot function from the visualization module.
Supports overlaying multiple model components or posterior predictive samples.
Automatically handles binning or weighting if present in the data.
Residuals plot shows data minus model values for quick assessment of fit quality.
Examples
Plot O–C data with a fitted sinusoidal model:
>>> sinus_model = oc.fit_sinusoidal(amp=Parameter(value=0.01), P=Parameter(value=1000)) >>> oc.plot(model=sinus_model, title="O–C with Sinusoidal Fit")
Plot data with residuals on a separate subplot:
>>> oc.plot(model=sinus_model, res=True, fig_size=(12, 8))
ocpy.oc_lmfit module¶
- class ocpy.oc_lmfit.OCLMFit(oc, minimum_time=None, minimum_time_error=None, weights=None, minimum_type=None, labels=None, cycle=None)[source]¶
Bases:
OCO–C data handler with fitting capabilities using lmfit.
This class extends OC by providing methods to fit O–C data to linear, quadratic, sinusoidal, Keplerian, and combined models using lmfit.Model. It supports parameter constraints, fixed parameters, and custom initial values. Fitted models can be used to compute residuals or generate new O–C predictions.
- math¶
Mathematical module used in component evaluations. Defaults to numpy.
- Type:
module
Notes
Each fit method (fit_linear, fit_quadratic, fit_sinusoidal, etc.) returns an lmfit.ModelResult with the best-fit parameters.
Parameters can be passed as Parameter instances, floats, or None. None defaults are used.
Residuals can be computed using residue() to generate a new OCLMFit object.
Designed for O–C analysis of periodic phenomena in astronomy, such as eclipsing binaries.
Examples
Fit a linear model to the O–C data:
>>> oc_fit = OCLMFit(oc, minimum_time=times, cycle=cycles, weights=weights) >>> result = oc_fit.fit_linear(a=0.0, b=0.0)
Fit a Keplerian model with custom initial parameters:
>>> result = oc_fit.fit_keplerian(amp=0.01, e=0.1, P=1000.0, T0=2450000.0)
Compute residuals after fitting:
>>> oc_resid = oc_fit.residue(result)
- math = <module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/oc-py/envs/latest/lib/python3.12/site-packages/numpy/__init__.py'>¶
- fit(model_components, *, nan_policy='raise', method='leastsq', **kwargs)[source]¶
Fit the O–C data using one or more model components via lmfit.
This method combines one or more ModelComponentModel instances into a single composite lmfit model, applies parameter constraints, and fits the O–C data (cycle vs oc) according to the provided weights.
- Parameters:
model_components (
listofModelComponentModel) – A list of model components (e.g., Linear, Quadratic, Keplerian, Sinusoidal) to fit simultaneously.nan_policy (
{'raise', 'omit', 'propagate'}, optional) – How to handle NaN values in the data. Default is ‘raise’.method (
str, optional) – The fitting algorithm to use. Passed directly to lmfit.Model.fit. Default is ‘leastsq’.**kwargs – Additional keyword arguments forwarded to lmfit.Model.fit.
- Returns:
The lmfit.ModelResult object containing best-fit parameters, uncertainties, and fit statistics.
- Return type:
ModelResult- Raises:
ValueError – If the weights column contains NaN values.
TypeError – If any model component is not compatible with lmfit.
Notes
Each component can have initial parameter values, bounds (min, max), and fixed status.
The method automatically generates unique prefixes for each component to avoid parameter name collisions.
Fitted parameters are applied to the respective ModelComponentModel instances internally.
Supports multiple components of the same type with unique prefixes.
Examples
Fit a linear and sinusoidal model simultaneously:
>>> oc_fit = OCLMFit(oc, minimum_time=times, cycle=cycles, weights=weights) >>> linear_comp = Linear(a=0.0, b=0.0) >>> sinusoidal_comp = Sinusoidal(amp=0.01, P=1000) >>> result = oc_fit.fit([linear_comp, sinusoidal_comp], method='leastsq')
Access fitted parameters:
>>> result.params['linear_a'].value >>> result.params['sinusoidal_amp'].value
- residue(coefficients, *, x_col='cycle', y_col='oc')[source]¶
Compute the residual O–C values after subtracting a fitted model.
This method evaluates the fitted model at the observed cycles and subtracts the model values from the observed O–C data, returning a new OCLMFit instance containing only the residuals.
- Parameters:
coefficients (
ModelResult) – The fitted lmfit.ModelResult object containing the best-fit parameters for the model.x_col (
str, optional) – Name of the column in self.data to use as the independent variable. Default is ‘cycle’.y_col (
str, optional) – Name of the column in self.data to use as the dependent variable (O–C values). Default is ‘oc’.
- Returns:
A new OCLMFit instance with the same metadata as the original object, but with the oc column replaced by residuals (observed minus fitted values).
- Return type:
- Raises:
ValueError – If x_col or y_col is not present in the data.
Notes
The returned object preserves all other columns (minimum_time, weights, labels, etc.) from the original dataset.
Useful for iterative fitting, diagnostics, or residual analysis.
Examples
>>> oc_fit = OCLMFit(oc, minimum_time=times, cycle=cycles, weights=weights) >>> result = oc_fit.fit_linear(a=0.0, b=0.0) >>> residuals = oc_fit.residue(result) >>> residuals.data.head()
- fit_linear(*, a=None, b=None, **kwargs)[source]¶
Fit a linear model to the O–C data using least-squares optimization.
The model has the form:
y = a * x + b
where x is the cycle number and y is the O–C value. Both slope (a) and intercept (b) can be provided as fixed values, free parameters, or left to be estimated from the data.
- Parameters:
a (
Parameter,float, orNone, optional) – Initial guess or fixed value for the slope. If None, defaults to 0.0.b (
Parameter,float, orNone, optional) – Initial guess or fixed value for the intercept. If None, defaults to 0.0.**kwargs – Additional keyword arguments are passed to the underlying fit method, e.g., method, nan_policy, etc.
- Returns:
An lmfit.ModelResult object containing the optimized parameters, fit statistics, and methods to evaluate the model or extract residuals.
- Return type:
ModelResult- Raises:
ValueError – If the O–C data contains NaN values in the weights column, which are required for weighted fitting.
Notes
This method internally wraps the a and b values as Parameter objects if they are not already, allowing fixed or free fitting.
Useful for detecting long-term trends in O–C diagrams.
Examples
>>> oc_fit = OCLMFit(oc, minimum_time=times, cycle=cycles, weights=weights) >>> result = oc_fit.fit_linear(a=0.0, b=0.0) >>> print(result.best_values) {'a': 0.0012, 'b': -0.03} >>> residuals = oc_fit.residue(result)
- fit_quadratic(*, q=None, **kwargs)[source]¶
Fit a quadratic model to the O–C data using least-squares optimization.
The model has the form:
y = q * x^2
where x is the cycle number and y is the O–C value. The quadratic coefficient q can be provided as a fixed value, a free parameter, or left to be estimated from the data.
- Parameters:
q (
Parameter,float, orNone, optional) – Initial guess or fixed value for the quadratic coefficient. If None, defaults to 0.0.**kwargs – Additional keyword arguments are passed to the underlying fit method, e.g., method, nan_policy, etc.
- Returns:
An lmfit.ModelResult object containing the optimized parameter q, fit statistics, and methods to evaluate the model or extract residuals.
- Return type:
ModelResult- Raises:
ValueError – If the O–C data contains NaN values in the weights column, which are required for weighted fitting.
Notes
This method internally wraps q as a Parameter object if it is not already, allowing fixed or free fitting.
Useful for detecting parabolic trends in O–C diagrams, e.g., due to period changes.
Examples
>>> oc_fit = OCLMFit(oc, minimum_time=times, cycle=cycles, weights=weights) >>> result = oc_fit.fit_quadratic(q=1e-6) >>> print(result.best_values) {'q': 2.3e-6} >>> residuals = oc_fit.residue(result)
- fit_parabola(*, q=None, a=None, b=None, **kwargs)[source]¶
Fit a combined quadratic and linear model to the O–C data using least-squares optimization.
The model has the form:
y = q * x^2 + a * x + b
where x is the cycle number and y is the O–C value. Each parameter (q, a, b) can be provided as a fixed value, a free parameter, or left to be estimated from the data.
- Parameters:
q (
Parameter,float, orNone, optional) – Initial guess or fixed value for the quadratic coefficient. Defaults to 0.0 if None.a (
Parameter,float, orNone, optional) – Initial guess or fixed value for the linear coefficient. Defaults to 0.0 if None.b (
Parameter,float, orNone, optional) – Initial guess or fixed value for the constant term. Defaults to 0.0 if None.**kwargs – Additional keyword arguments passed to the underlying fit method, e.g., method, nan_policy, weights, etc.
- Returns:
An lmfit.ModelResult object containing optimized parameters q, a, b, fit statistics, and methods to evaluate the model or extract residuals.
- Return type:
ModelResult- Raises:
ValueError – If the O–C data contains NaN values in the weights column.
Notes
Internally wraps all input values as Parameter objects for consistent handling.
Useful for modeling O–C diagrams that exhibit both parabolic trends (period changes) and linear trends (systematic offsets).
Examples
>>> oc_fit = OCLMFit(oc, minimum_time=times, cycle=cycles, weights=weights) >>> result = oc_fit.fit_parabola(q=1e-6, a=1e-4, b=0.0) >>> print(result.best_values) {'q': 2.3e-6, 'a': 1.1e-4, 'b': 0.0} >>> residuals = oc_fit.residue(result)
- fit_lite(*, amp=None, e=None, omega=None, P=None, T0=None, **kwargs)[source]¶
Fit a simplified Keplerian model to O–C data using least-squares optimization.
This “lite” Keplerian model describes a sinusoidal O–C variation approximating the orbital motion of a binary system or planet, with the form:
y = amp * ( (1-e^2)/(1+e*cos(true_anomaly)) * sin(true_anomaly + omega)/denom + e*sin(omega)/denom )
where true_anomaly is calculated from the mean anomaly using Kepler’s equation.
- Parameters:
amp (
Parameter,float, orNone, optional) – Semi-amplitude of the O–C variation. Defaults to 1e-3 if None.e (
Parameter,float, orNone, optional) – Orbital eccentricity (0 ≤ e < 0.95). Defaults to 0.0 if None.omega (
Parameter,float, orNone, optional) – Argument of periastron in degrees. Defaults to 90.0 if None.P (
Parameter,float, orNone, optional) – Orbital period. Defaults to 3000.0 if None.T0 (
Parameter,float, orNone, optional) – Reference time of periastron passage. Defaults to 0.0 if None.**kwargs – Additional keyword arguments passed to the underlying fit method, such as method, nan_policy, or weights.
- Returns:
An lmfit.ModelResult object containing optimized parameters amp, e, omega, P, T0, fit statistics, and methods for evaluation or residual extraction.
- Return type:
ModelResult- Raises:
ValueError – If the O–C data contains NaN values in the weights column.
Notes
Parameters are internally wrapped as Parameter objects if they are not already.
Useful for rapid fitting of periodic variations in O–C diagrams without modeling additional complexities of full Keplerian motion.
Examples
>>> oc_fit = OCLMFit(oc, minimum_time=times, cycle=cycles, weights=weights) >>> result = oc_fit.fit_lite(amp=0.01, e=0.1, P=5000.0) >>> print(result.best_values) {'amp': 0.012, 'e': 0.098, 'omega': 90.0, 'P': 4990.2, 'T0': 0.0} >>> residuals = oc_fit.residue(result)
- fit_keplerian(*, amp=None, e=None, omega=None, P=None, T0=None, **kwargs)[source]¶
Fit a Keplerian model to O–C data by delegating to the fit_lite method.
This method is a convenient alias for fitting standard Keplerian variations in O–C diagrams, modeling periodic deviations due to orbital motion. Internally, it calls fit_lite with the given parameters.
- Parameters:
amp (
Parameter,float, orNone, optional) – Semi-amplitude of the O–C variation. Defaults to 1e-3 if None.e (
Parameter,float, orNone, optional) – Orbital eccentricity (0 ≤ e < 0.95). Defaults to 0.0 if None.omega (
Parameter,float, orNone, optional) – Argument of periastron in degrees. Defaults to 90.0 if None.P (
Parameter,float, orNone, optional) – Orbital period. Defaults to 3000.0 if None.T0 (
Parameter,float, orNone, optional) – Reference time of periastron passage. Defaults to 0.0 if None.**kwargs – Additional keyword arguments passed to fit_lite and ultimately to OCLMFit.fit, such as method, nan_policy, or weights.
- Returns:
An lmfit.ModelResult object containing optimized parameters amp, e, omega, P, T0, fit statistics, and methods for evaluation or residual extraction.
- Return type:
ModelResult- Raises:
ValueError – If the O–C data contains NaN values in the weights column (raised in fit_lite).
Notes
This method provides semantic clarity for users who want a Keplerian fit, but it is functionally identical to fit_lite.
Parameters are internally wrapped as Parameter objects if they are not already.
Examples
>>> oc_fit = OCLMFit(oc, minimum_time=times, cycle=cycles, weights=weights) >>> result = oc_fit.fit_keplerian(amp=0.01, e=0.1, P=5000.0) >>> print(result.best_values) {'amp': 0.012, 'e': 0.098, 'omega': 90.0, 'P': 4990.2, 'T0': 0.0} >>> residuals = oc_fit.residue(result)
- fit_sinusoidal(*, amp=None, P=None, **kwargs)[source]¶
Fit a sinusoidal model to O–C data.
This method fits a simple sinusoidal variation to the observed O–C values, which can model periodic effects such as light-time effects in binary systems.
- Parameters:
amp (
Parameter,float, orNone, optional) – Amplitude of the sinusoidal O–C variation. Defaults to 1e-3 if None.P (
Parameter,float, orNone, optional) – Period of the sinusoidal variation. Defaults to 3000.0 if None.**kwargs – Additional keyword arguments passed to OCLMFit.fit, such as method, nan_policy, or weights.
- Returns:
An lmfit.ModelResult object containing optimized parameters amp and P, fit statistics, and methods for evaluation or residual extraction.
- Return type:
ModelResult- Raises:
ValueError – If the O–C data contains NaN values in the weights column.
Notes
This method internally wraps numerical values into Parameter objects if they are not already.
Useful for detecting or modeling sinusoidal periodicities in O–C diagrams.
Examples
>>> oc_fit = OCLMFit(oc, minimum_time=times, cycle=cycles, weights=weights) >>> result = oc_fit.fit_sinusoidal(amp=0.01, P=5000.0) >>> print(result.best_values) {'amp': 0.011, 'P': 4985.0} >>> residuals = oc_fit.residue(result)
ocpy.oc_pymc module¶
- class ocpy.oc_pymc.OCPyMC(oc, minimum_time=None, minimum_time_error=None, weights=None, minimum_type=None, labels=None, cycle=None)[source]¶
Bases:
OCA probabilistic O–C (Observed minus Calculated) modeling class using PyMC.
This class extends the OC base class to provide Bayesian inference for O–C data, allowing modeling with various components such as linear, quadratic, sinusoidal, and Keplerian functions. It leverages PyMC and ArviZ for probabilistic sampling, posterior analysis, and visualization.
- Parameters:
OC (Inherited from) – Accepts the same initialization parameters as OC, such as cycle counts, O–C values, errors, weights, labels, and minimum types.
- math¶
Set to pymc.math to allow model components to use PyMC tensor operations.
- Type:
module
- fit(model_components, \*, draws, tune, chains, cores, target_accept, random_seed, progressbar, return_model, \*\*kwargs)[source]¶
Fit the provided model components to the O–C data using Bayesian sampling. Returns an InferenceData object with posterior samples.
- clean(inference_data, drop_chains=0, filter_outliers=True, iqr_multiplier=4.0)[source]¶
Post-process InferenceData to drop chains or filter extreme outliers.
- residue(inference_data, \*, x_col="cycle", y_col="oc")[source]¶
Return a new OCPyMC instance containing residuals between observed and model-predicted O–C values.
- fit_linear(*, a=None, b=None, cores=None, \*\*kwargs)[source]¶
Fit a linear component (a * x + b) to the data.
- fit_quadratic(*, q=None, cores=None, \*\*kwargs)[source]¶
Fit a quadratic component (q * x^2) to the data.
- fit_parabola(*, q=None, a=None, b=None, cores=None, \*\*kwargs)[source]¶
Fit a combination of quadratic and linear components to the data.
- fit_sinusoidal(*, amp=None, P=None, cores=None, \*\*kwargs)[source]¶
Fit a sinusoidal component (amp * sin(2π x / P)) to the data.
- fit_keplerian(*, amp=None, e=None, omega=None, P=None, T0=None, name=None, cores=None, \*\*kwargs)[source]¶
Fit a Keplerian (orbital) component to the data.
- fit_lite(**kwargs)[source]¶
Alias for fit_keplerian; provides a simplified interface for fitting a single Keplerian component.
- corner(inference_data, cornerstyle="corner", units=None, \*\*kwargs)[source]¶
Plot a corner plot of posterior distributions using Plot.plot_corner.
- trace(inference_data, \*\*kwargs)[source]¶
Plot trace plots for posterior chains using Plot.plot_trace.
Notes
This class is designed for O–C analyses in eclipsing binaries, exoplanet timing, and other periodic astrophysical phenomena.
Supports both deterministic and probabilistic parameters, including truncated distributions for constrained parameters.
Works seamlessly with multiple model components, automatically managing parameter naming and aggregation.
- math = <module 'pymc.math' from '/home/docs/checkouts/readthedocs.org/user_builds/oc-py/envs/latest/lib/python3.12/site-packages/pymc/math.py'>¶
- fit(model_components, *, draws=2000, tune=2000, chains=4, cores=None, target_accept=None, random_seed=None, progressbar=True, return_model=False, **kwargs)[source]¶
Fit a model to the observed O–C (observed minus calculated) data using Bayesian inference.
This method constructs a PyMC probabilistic model using the specified model_components, applies priors based on their parameters, and samples from the posterior distribution.
- Parameters:
model_components (
listofModelComponent) – List of model components to fit (e.g., Linear, Quadratic, Sinusoidal, Keplerian). Each component must define a model_func and optional params.draws (
int, optional) – Number of posterior samples to draw after tuning (default: 2000).tune (
int, optional) – Number of tuning (burn-in) steps (default: 2000).chains (
int, optional) – Number of Markov chains to run (default: 4).cores (
intorNone, optional) – Number of CPU cores to use. Defaults to number of chains if None.target_accept (
floatorNone, optional) – Target acceptance probability for NUTS sampler (default: None).random_seed (
intorNone, optional) – Random seed for reproducibility.progressbar (
bool, optional) – Whether to show the sampling progress bar (default: True).return_model (
bool, optional) – If True, return the PyMC Model object without sampling.**kwargs – Additional keyword arguments passed to pm.sample() or sampler step methods.
- Returns:
If return_model is False, returns an ArviZ InferenceData object containing posterior samples, posterior predictive samples, and deterministic variables for model evaluation. If return_model is True, returns the PyMC Model object without sampling.
- Return type:
arviz.InferenceDataorpm.Model
Notes
All model components’ parameters are converted to Parameter objects using _to_param.
Observational uncertainties are taken from minimum_time_error; NaNs will raise a ValueError.
Deterministic variables y_model and, if possible, y_model_dense are created for visualization.
Components marked as _expensive bypass dense evaluation to avoid heavy computations.
- clean(inference_data, drop_chains=0, filter_outliers=True, iqr_multiplier=4.0)[source]¶
Clean an ArviZ InferenceData object by optionally removing chains and filtering outliers.
This method allows post-processing of MCMC samples to improve data quality for visualization or further analysis. Chains with extreme deviations can be dropped, and outliers in the posterior distributions can be filtered using the interquartile range (IQR) method.
- Parameters:
inference_data (
arviz.InferenceData) – The posterior samples to clean, typically returned by OCPyMC.fit().drop_chains (
int, optional) – Number of chains to drop based on deviation from the median across chains (default: 0, i.e., keep all chains).filter_outliers (
bool, optional) – Whether to remove outliers from the posterior samples (default: True).iqr_multiplier (
float, optional) – Multiplier for the interquartile range (IQR) to define outlier bounds (default: 4.0). Samples outside [Q1 - multiplier*IQR, Q3 + multiplier*IQR] are removed.
- Returns:
A new InferenceData object with cleaned posterior samples, while preserving deterministic variables and model component metadata (_model_components and _model_prefixes).
- Return type:
arviz.InferenceData
Notes
Outlier filtering is applied independently to each parameter in the posterior (excluding deterministic variables like y_model, y_model_dense, y_obs, dense_x).
Chains are ranked by their deviation from the overall median, and the drop_chains largest deviations are removed.
The cleaned dataset preserves the original structure of InferenceData, including posterior, prior, and observed_data groups.
Any errors encountered during filtering are issued as warnings without stopping execution.
- residue(inference_data, *, x_col='cycle', y_col='oc')[source]¶
Compute the residuals between observed O–C values and the model predictions.
This method generates a new OCPyMC instance where the O–C values are replaced by the difference between the observed values and the median model predictions from a fitted PyMC InferenceData object.
- Parameters:
inference_data (
arviz.InferenceData) – The posterior samples returned by OCPyMC.fit() containing the fitted model.x_col (
str, optional) – Name of the column in self.data containing the independent variable (cycle), used to evaluate the model (default: “cycle”).y_col (
str, optional) – Name of the column in self.data containing the observed O–C values (default: “oc”).
- Returns:
A new OCPyMC instance with the same metadata as the original, but with the oc values replaced by the residuals (observed minus fitted).
- Return type:
Notes
The model predictions are computed as the median across all chains and draws in the posterior (inference_data.posterior[“y_model”]).
Only the oc values are updated; all other attributes (weights, labels, minimum times, etc.) are preserved.
This method is useful for examining the remaining structure or noise after fitting a model.
- fit_linear(*, a=None, b=None, cores=None, **kwargs)[source]¶
Fit a linear O–C model (O–C = a * cycle + b) using PyMC Bayesian inference.
- The linear model is represented as:
O–C = a * x + b
where x is the cycle number.
- Parameters:
a (
floatorParameter, optional) – Initial value or Parameter object for the slope of the line. If None, defaults to 0.0.b (
floatorParameter, optional) – Initial value or Parameter object for the intercept of the line. If None, defaults to 0.0.cores (
int, optional) – Number of CPU cores to use for sampling. If None, the number of cores defaults to the number of chains.**kwargs – Additional keyword arguments passed to OCPyMC.fit(), such as draws, tune, chains, target_accept, etc.
- Returns:
Posterior samples for the linear model parameters and deterministic predictions for the O–C values.
- Return type:
arviz.InferenceData
Notes
The slope (a) and intercept (b) are treated as PyMC random variables unless provided as fixed Parameter objects with fixed=True.
Use residue() on the returned model to compute residuals after fitting.
This method is a convenience wrapper that internally creates a Linear component and calls the general fit() method.
- fit_quadratic(*, q=None, cores=None, **kwargs)[source]¶
Fit a quadratic O–C model (O–C = q * cycle^2) using PyMC Bayesian inference.
- The quadratic model is represented as:
O–C = q * x^2
where x is the cycle number.
- Parameters:
q (
floatorParameter, optional) – Initial value or Parameter object for the quadratic coefficient. If None, defaults to 0.0.cores (
int, optional) – Number of CPU cores to use for sampling. If None, the number of cores defaults to the number of chains.**kwargs – Additional keyword arguments passed to OCPyMC.fit(), such as draws, tune, chains, target_accept, etc.
- Returns:
Posterior samples for the quadratic coefficient and deterministic predictions for the O–C values.
- Return type:
arviz.InferenceData
Notes
The coefficient q is treated as a PyMC random variable unless provided as a fixed Parameter with fixed=True.
Use residue() on the returned model to compute residuals after fitting.
This method is a convenience wrapper that internally creates a Quadratic component and calls the general fit() method.
- fit_sinusoidal(*, amp=None, P=None, cores=None, **kwargs)[source]¶
Fit a sinusoidal O–C model (O–C = amp * sin(2π * cycle / P)) using PyMC Bayesian inference.
- The sinusoidal model is represented as:
O–C = amp * sin(2π * x / P)
where x is the cycle number.
- Parameters:
amp (
floatorParameter, optional) – Initial value or Parameter object for the sinusoidal amplitude. If None, defaults to 1e-3.P (
floatorParameter, optional) – Initial value or Parameter object for the period of the sinusoid. If None, defaults to 1000.0.cores (
int, optional) – Number of CPU cores to use for sampling. If None, the number of cores defaults to the number of chains.**kwargs – Additional keyword arguments passed to OCPyMC.fit(), such as draws, tune, chains, target_accept, etc.
- Returns:
Posterior samples for the sinusoidal parameters and deterministic predictions for the O–C values.
- Return type:
arviz.InferenceData
Notes
The parameters amp and P are treated as PyMC random variables unless provided as fixed Parameter objects with fixed=True.
Use residue() on the returned model to compute residuals after fitting.
This method is a convenience wrapper that internally creates a Sinusoidal component and calls the general fit() method.
- fit_keplerian(*, amp=None, e=None, omega=None, P=None, T0=None, name=None, cores=None, **kwargs)[source]¶
Fit a Keplerian O–C model using PyMC Bayesian inference.
The Keplerian model represents the O–C variations as caused by a companion in an orbit, parameterized by:
amp : amplitude of the timing variation
e : orbital eccentricity
omega: argument of periastron (degrees)
P : orbital period
T0 : time of periastron passage
- Parameters:
amp (
floatorParameter, optional) – Initial value or Parameter object for the amplitude of the O–C signal. Defaults to 0.001 if None.e (
floatorParameter, optional) – Initial value or Parameter object for orbital eccentricity. Defaults to 0.1 if None.omega (
floatorParameter, optional) – Initial value or Parameter object for the argument of periastron (degrees). Defaults to 90.0 if None.P (
floatorParameter, optional) – Initial value or Parameter object for the orbital period. Defaults to 1000.0 if None.T0 (
floatorParameter, optional) – Initial value or Parameter object for the time of periastron passage. Defaults to 0.0 if None.name (
str, optional) – Optional name for the Keplerian component. If None, defaults to ‘keplerian1’.cores (
int, optional) – Number of CPU cores to use for sampling. If None, defaults to the number of chains.**kwargs – Additional keyword arguments passed to OCPyMC.fit(), e.g., draws, tune, chains, target_accept, etc.
- Returns:
Posterior samples for all Keplerian parameters, along with deterministic predictions for the O–C curve.
- Return type:
arviz.InferenceData
Notes
The parameters are treated as PyMC random variables unless provided as fixed Parameter objects with fixed=True.
Use residue() on the returned model to compute residuals after fitting.
This method internally creates a Keplerian component and calls the general fit() method for Bayesian inference.
- fit_lite(**kwargs)[source]¶
Fit a simplified Keplerian O–C model using PyMC Bayesian inference.
This is a convenience wrapper around fit_keplerian with default Keplerian parameters. It is designed for quick fitting without specifying all orbital parameters explicitly.
- Parameters:
**kwargs – Keyword arguments passed directly to fit_keplerian, including: - amp : amplitude of the O–C signal - e : orbital eccentricity - omega: argument of periastron (degrees) - P : orbital period - T0 : time of periastron passage - cores: number of CPU cores to use - draws, tune, chains, target_accept, etc.
- Returns:
Posterior samples for the Keplerian parameters and deterministic predictions of the O–C curve.
- Return type:
arviz.InferenceData
Notes
- This method internally calls fit_keplerian with default parameter values:
amp=0.001, e=0.1, omega=90.0, P=1000.0, T0=0.0
Use residue() on the returned model to compute residuals after fitting.
Suitable for quick exploratory fits or as an initial guess for more complex models.
- fit_parabola(*, q=None, a=None, b=None, cores=None, **kwargs)[source]¶
Fit a combined quadratic and linear O–C model using PyMC Bayesian inference.
This method models the O–C variations as a sum of a quadratic component (representing secular period change) and a linear component (representing a constant period offset).
- Parameters:
q (
floatorParameter, optional) – Quadratic coefficient for the parabolic term. Default is 0.0.a (
floatorParameter, optional) – Linear coefficient for the linear term. Default is 0.0.b (
floatorParameter, optional) – Constant offset term. Default is 0.0.cores (
int, optional) – Number of CPU cores to use for PyMC sampling. Defaults to the number of chains.**kwargs – Additional keyword arguments passed to the underlying fit method, such as draws, tune, chains, target_accept, and random_seed.
- Returns:
Posterior samples for the quadratic and linear coefficients, along with deterministic predictions of the O–C curve.
- Return type:
arviz.InferenceData
Notes
The quadratic component captures long-term period changes.
The linear component models any overall offset in the period.
Use the residue() method on the returned object to compute residuals of the fit.
Recommended as a first approach when O–C variations show a parabolic trend.
- corner(inference_data, cornerstyle='corner', units=None, **kwargs)[source]¶
Generate a corner (pairwise posterior) plot from PyMC inference results.
This method visualizes the posterior distributions of model parameters, showing marginal distributions along the diagonal and pairwise correlations in the off-diagonal plots. It can use either the corner library or ArviZ built-in plotting depending on the cornerstyle argument.
- Parameters:
inference_data (
arviz.InferenceData) – The PyMC sampling results, typically returned by fit(), containing posterior samples.cornerstyle (
{"corner", "arviz"}, default"corner") – Choose the plotting backend: - “corner”: use the corner.py library for plotting. - “arviz”: use ArviZ’s built-in pairplot functionality.units (
dict, optional) – A mapping from parameter names to strings describing their units. These units will be displayed in axis labels.**kwargs – Additional keyword arguments passed to the underlying plotting function.
- Returns:
matplotlib.figure.Figure– The figure object for the ‘corner’ style.arviz.plot_pair axesorFigure– The ArviZ axes object when cornerstyle=’arviz’.
- Return type:
Figure | plot_pair
Notes
Useful for visual inspection of posterior distributions and parameter correlations.
Recommended to use after fit() and optionally after clean() for outlier removal.
Units help clarify parameter scales in the plot.
- trace(inference_data, **kwargs)[source]¶
Generate trace plots for PyMC posterior samples.
Trace plots visualize the sampled parameter values across chains, allowing assessment of convergence, mixing, and sampling behavior. Each subplot shows the parameter values over draws for one chain.
- Parameters:
inference_data (
arviz.InferenceData) – The PyMC sampling results, typically returned by fit(), containing posterior samples.**kwargs – Additional keyword arguments passed to the underlying plotting function.
- Returns:
Array of matplotlib axes objects containing the trace plots.
- Return type:
matplotlib.axes.Axes
Notes
Useful for evaluating MCMC convergence and diagnosing sampling issues.
Typically used after fit() and optionally after clean() to inspect cleaned chains.
Can be combined with corner plots for a complete posterior overview.
ocpy.utils module¶
ocpy.visualization module¶
- class ocpy.visualization.Plot[source]¶
Bases:
object- static plot_data(data, *, ax=None, x_col='cycle', y_col='oc', plot_kwargs=None)[source]¶
Plot the raw O−C data with optional error bars and labeling.
- Parameters:
data (
OC) – The observational O−C dataset. Must have a data attribute (pandas DataFrame) containing at least columns for x_col and y_col. Optionally, it can include ‘minimum_time_error’ for y-error bars and ‘labels’ for grouping data by category.ax (
matplotlib.axes.Axes, optional) – Axes object on which to plot. If None, a new figure and axes are created.x_col (
str, default"cycle") – Name of the column to use for the x-axis.y_col (
str, default"oc") – Name of the column to use for the y-axis.plot_kwargs (
dict, optional) – Additional keyword arguments passed to matplotlib.pyplot.errorbar.
- Returns:
The axes containing the plotted data.
- Return type:
matplotlib.axes.Axes
Notes
If ‘labels’ exist in the DataFrame, points are color-coded per unique label.
Unlabeled points are plotted in gray.
If ‘minimum_time_error’ exists, it is used as y-error bars.
Axes are automatically labeled and a grid is added.
- classmethod plot_model_pymc(inference_data, data, *, ax=None, x_col='cycle', n_points=800, sum_kwargs=None, comp_kwargs=None, plot_kwargs=None, plot_band=True, extension_factor=0.1, model_components=None)[source]¶
Plot a model fit to O−C data using a PyMC inference result, with optional uncertainty bands and component decomposition.
- Parameters:
inference_data (
arviz.InferenceData) – Posterior samples from a PyMC model. Expected to contain parameter variables (2D arrays with shape [chain, draw]), and optionally ‘y_model’ or ‘y_model_dense’ for precomputed model fits.data (
OCPyMC) – Observational O−C dataset to plot against. Must have a data attribute (pandas DataFrame) containing at least the x_col column.ax (
matplotlib.axes.Axes, optional) – Axes object on which to plot. If None, uses current axes or creates a new figure.x_col (
str, default"cycle") – Column in data to use as the x-axis.n_points (
int, default800) – Number of points to use when plotting continuous model curves.sum_kwargs (
dict, optional) – Additional keyword arguments for the summed model line.comp_kwargs (
dict, optional) – Additional keyword arguments for individual model component lines.plot_kwargs (
dict, optional) – Additional keyword arguments for the data points plot (markers, color, etc.).plot_band (
bool, defaultTrue) – Whether to display a 1σ uncertainty band from posterior samples.extension_factor (
float, default0.1) – Fractional extension beyond the data range for plotting the fit curve.model_components (
list, optional) – List of model component objects (Linear, Quadratic, Sinusoidal, Keplerian, etc.) to reconstruct and plot individual contributions.
- Returns:
The axes containing the plotted data, model fit, components, and optional uncertainty band.
- Return type:
matplotlib.axes.Axes
Notes
If ‘y_model_dense’ and ‘dense_x’ exist in inference_data, they are used for a smooth fit.
Otherwise, the method reconstructs the model from posterior medians of scalar parameters.
If model_components are provided, each component is plotted individually.
Uncertainty bands are computed from a subset of posterior samples (default 200 draws).
Automatically handles multiple components, labeling, and plotting the sum of components.
- classmethod plot_model_lmfit(result, data, *, ax=None, x_col='cycle', n_points=500, plot_kwargs=None, extension_factor=0.1)[source]¶
Plot a model fit to O−C data using an lmfit.ModelResult object, with optional uncertainty bands.
- Parameters:
result (
lmfit.model.ModelResult) – The fitted model result returned by lmfit.Model.fit. Expected to provide eval(x=…) and optionally eval_uncertainty(x=…, sigma=1).data (
OCLMFit) – Observational O−C dataset to plot against. Must have a data attribute (pandas DataFrame) containing at least the x_col column.ax (
matplotlib.axes.Axes, optional) – Axes object on which to plot. If None, creates a new figure and axes.x_col (
str, default"cycle") – Column in data to use as the x-axis.n_points (
int, default500) – Number of points to evaluate for a smooth model curve.plot_kwargs (
dict, optional) – Additional keyword arguments for the fit line (color, linewidth, label, etc.).extension_factor (
float, default0.1) – Fractional extension beyond the data range for plotting the fit curve.
- Returns:
The axes containing the plotted data, model fit, and optional uncertainty band.
- Return type:
matplotlib.axes.Axes
Notes
The method evaluates the fitted model on a dense set of points across the data range, optionally extended by extension_factor.
If result.eval_uncertainty is available, a 1σ uncertainty band is plotted around the fit.
Data points are plotted separately using the plot_data method (if called externally).
- classmethod plot_model_components(model_components, xline, *, ax=None, sum_kwargs=None, comp_kwargs=None, uncertainty_band=None)[source]¶
Plot individual model components and their sum over a specified x-range, with optional uncertainty bands.
- Parameters:
model_components (
list) – List of model component objects. Each component must have a model_func(x, **params) method and a params attribute (dict of Parameter objects or numeric values).xline (
np.ndarray) – Array of x-values to evaluate the component models.ax (
matplotlib.axes.Axes, optional) – Axes object on which to plot. If None, a new figure and axes are created.sum_kwargs (
dict, optional) – Keyword arguments for the summed model curve. Default color is red, linewidth 2.6, alpha 0.95.comp_kwargs (
dict, optional) – Keyword arguments for individual component curves. Default linewidth 1.5, alpha 0.9, linestyle ‘–‘.uncertainty_band (
tuple, optional) – Tuple (x_band, y_low, y_high) representing an uncertainty envelope around the sum of components. If provided, plotted as a filled area behind the curves.
- Returns:
The axes containing the component curves, sum curve, and optional uncertainty band.
- Return type:
matplotlib.axes.Axes
Notes
Each component is evaluated using its model_func and current parameter values.
The sum curve is drawn on top of the components, optionally with an uncertainty band.
Components without parameters (or with missing parameters) will raise a KeyError.
Useful for visualizing contributions of multiple additive model components in O−C analysis or similar time-series modeling contexts.
- classmethod plot(data, model=None, *, ax=None, res_ax=None, res=True, title=None, x_col='cycle', y_col='oc', fig_size=(10, 7), plot_kwargs=None, extension_factor=0.1, model_components=None)[source]¶
Plot data with optional model fit and residuals.
This is a high-level plotting function that can: - Display raw O−C data points (with optional labels and error bars), - Overlay model fits from PyMC (InferenceData), lmfit (ModelResult), or a list of model components, - Display residuals below the main plot if requested.
- Parameters:
data (
OC) – The data object containing O−C measurements. Must have data attribute (pandas DataFrame) with at least columns specified by x_col and y_col.model (
Union[InferenceData,ModelResult,List[ModelComponent]], optional) – Model to overlay on the data: - PyMC model (InferenceData from arviz) with posterior samples, - lmfit result (ModelResult) with .eval() method, - List of component objects with .model_func and .params. If None, only the raw data is plotted.ax (
matplotlib.axes.Axes, optional) – Axes for the main data/fit plot. If None, a new figure is created.res_ax (
matplotlib.axes.Axes, optional) – Axes for residuals plot. If None and res=True, a new subplot is created.res (
bool, defaultTrue) – Whether to plot residuals beneath the main plot.title (
str, optional) – Title for the main plot.x_col (
str, default"cycle") – Column in data.data used for x-axis.y_col (
str, default"oc") – Column in data.data used for y-axis.fig_size (
tuple, default(10,7)) – Figure size in inches.plot_kwargs (
dict, optional) – Keyword arguments passed to the main data plot (color, markers, alpha, etc.).extension_factor (
float, default0.1) – Fractional extension beyond the data range for plotting model fits.model_components (
list, optional) – If provided, used for plotting PyMC components in plot_model_pymc.
- Returns:
If res=True, returns a tuple (ax, res_ax) for main and residual plots. Otherwise, returns only ax.
- Return type:
matplotlib.axes.Axesortuple(matplotlib.axes.Axes,matplotlib.axes.Axes)
Notes
Automatically handles labeled and unlabeled data points.
Residuals are computed as y - y_model if model is provided.
Supports both PyMC posterior models (median and uncertainty bands) and lmfit fits.
Useful for O−C analysis in astronomy or any time-series with additive models.
- static plot_corner(inference_data, var_names=None, cornerstyle='corner', units=None, **kwargs)[source]¶
Generate a corner plot (pairwise parameter plot) from PyMC/ArviZ inference data.
- Parameters:
inference_data (
arviz.InferenceData) – The posterior samples from a Bayesian model.var_names (
listofstr, optional) – Names of parameters to include. If None, automatically selects all posterior parameters with 2 dimensions (chain, draw) excluding model output.cornerstyle (
{'corner', 'arviz'}, default'corner') – Which library/style to use for the corner plot: - ‘corner’ : uses the corner package. - ‘arviz’ : uses ArviZ’s plot_pair.units (
dict, optional) – Mapping from parameter name to unit string, e.g., {‘P’: ‘days’}. Units are appended to axis labels.**kwargs – Additional keyword arguments passed to corner.corner or arviz.plot_pair.
- Returns:
matplotlib.figure.Figure– The figure object for the ‘corner’ style.arviz.plot_pair axesorFigure– The ArviZ axes object when cornerstyle=’arviz’.
- Return type:
Figure | plot_pair
Notes
Automatically ignores fixed parameters with negligible variation.
Adds median values as “truths” in the plot if using ‘corner’ and no truths provided.
If units are provided, labels are formatted with LaTeX, e.g., r”$P$ [days]”.
Supports up to 200 subplots in ArviZ style using plot.max_subplots context.
- Raises:
ImportError – If corner library is required but not installed.
ValueError – If no suitable parameters are found for plotting.
- static plot_trace(inference_data, var_names=None, **kwargs)[source]¶
Generate trace plots for posterior samples from a PyMC InferenceData object.
Trace plots show the sampled parameter values across chains and draws, allowing evaluation of convergence, mixing, and sampling behavior.
- Parameters:
inference_data (
arviz.InferenceData) – The posterior sampling results, typically returned by a PyMC fit.var_names (
listofstr, optional) – List of parameter names to include in the trace plot. If None, all variable parameters with variation are included.**kwargs – Additional keyword arguments passed to arviz.plot_trace.
- Returns:
Array of matplotlib axes objects containing the trace plots.
- Return type:
matplotlib.axes.Axes
Notes
Automatically excludes fixed parameters (near-zero variance) unless all are fixed.
Figures are automatically tightened with tight_layout.
Designed to complement plot_corner for full posterior visualization.