Option to use CBVs and better priors

src/psfmachine/machine.py

@@ -1,6 +1,7 @@
 """
 Defines the main Machine object that fit a mean PRF model to sources
 """
+import logging
 import numpy as np
 import pandas as pd
 from scipy import sparse
@@ -19,6 +20,7 @@
 from .aperture import optimize_aperture, compute_FLFRCSAP, compute_CROWDSAP
 from .perturbation import PerturbationMatrix3D
 
+log = logging.getLogger(__name__)
 __all__ = ["Machine"]
 
 
@@ -610,12 +612,14 @@ def build_time_model(
         self,
         plot=False,
         bin_method="bin",
+        poly_order=3,
         segments=False,
         focus=False,
         focus_exptime=50,
         pca_ncomponents=0,
         pca_smooth_time_scale=0,
         positions=False,
+        other_vectors=None,
     ):
         """
         Builds a time model that moves the PRF model to account for the scene movement
@@ -632,6 +636,8 @@ def build_time_model(
             Plot a diagnostic figure.
         bin_method: string
             Type of bin method, options are "bin" and "downsample".
+        poly_order: int
+            Degree of the time polynomial used for the time model. Default is 3.
         segments : boolean
             If `True` will split the light curve into segments to fit different time
             models with a common pixel normalization. If `False` will fit the full
@@ -650,6 +656,11 @@ def build_time_model(
         positions : boolean or string
             If one of strings `"poscorr", "centroid"` then the perturbation matrix will
             add other vectors accounting for position shift.
+        other_vectors : list or numpy.ndarray
+            Set of other components used to include  in the perturbed model.
+            See `psfmachine.perturbation.PerturbationMatrix` object for details.
+            Posible use case are using Kepler CBVs or engeneering data. Shape has to be
+            (ncomponents, ntimes). Default is `None`.
         """
         # create the time and space basis
         _whitened_time = (self.time - self.time.mean()) / (
@@ -667,6 +678,8 @@ def build_time_model(
 
         # add other vectors if asked, centroids or poscorrs
         if positions:
+            if other_vectors is not None:
+                raise ValueError("When using `positions` do not provide `other_vector`")
             if positions == "poscorr" and hasattr(self, "pos_corr1"):
                 mpc1 = np.nanmedian(self.pos_corr1, axis=0)
                 mpc2 = np.nanmedian(self.pos_corr2, axis=0)
@@ -679,10 +692,15 @@ def build_time_model(
                 )
 
             # smooth the vectors
+            if not segments:
+                log.warning(
+                    "Segments will still be used to smooth the position vectors."
+                    "See https://github.com/SSDataLab/psfmachine/pull/63 for details."
+                )
             mpc1_smooth, mpc2_smooth = bspline_smooth(
                 [mpc1, mpc2],
                 x=self.time,
-                do_segments=segments,
+                do_segments=True,
                 n_knots=50,
             )
             # normalize components
@@ -694,14 +712,16 @@ def build_time_model(
             )
             # combine them as the first order
             other_vectors = [mpc1_smooth, mpc2_smooth, mpc1_smooth * mpc2_smooth]
-        else:
-            other_vectors = None
+        if other_vectors is not None:
+            if not isinstance(other_vectors, (list, np.ndarray)):
+                raise ValueError("`other vector` is not a list of arrays or a ndarray")
 
         # create a 3D perturbation matrix
         P = PerturbationMatrix3D(
             time=_whitened_time,
             dx=dx,
             dy=dy,
+            poly_order=poly_order,
             segments=segments,
             focus=focus,
             other_vectors=other_vectors,
@@ -879,7 +899,7 @@ def plot_time_model(self):
             vmin=0.9,
             vmax=1.1,
         )
-        ax[1].set(xlabel="Binned Time Index", title="Perturbed Mode")
+        ax[1].set(xlabel="Binned Time Index", title="Perturbed Model")
 
         cbar = fig2.colorbar(
             im, ax=ax[:2], shrink=0.7, orientation="horizontal", location="bottom"

src/psfmachine/perturbation.py

@@ -85,11 +85,11 @@ def __repr__(self):
 
     @property
     def prior_mu(self):
-        return np.zeros(self.shape[1])
+        return np.ones(self.shape[1])
 
     @property
     def prior_sigma(self):
-        return np.ones(self.shape[1]) * 1e4
+        return np.ones(self.shape[1]) * 0.5
 
     @property
     def breaks(self):

tests/test_perturbation.py

@@ -86,10 +86,8 @@ def test_perturbation_matrix3d():
     dx, dy = dx.ravel(), dy.ravel()
 
     # ntime x npixels
-    flux = (
-        np.random.normal(0, 0.01, size=(10, 169))
-        + 2 * dx[None, :] ** 2 * 3 * dy[None, :]
-    )
+    flux = np.random.normal(1, 0.01, size=(10, 169)) + dx[None, :] * dy[None, :]
+    # the perturbation model assumes the perturbation is around 1
     flux_err = np.ones((10, 169)) * 0.01
 
     p3 = PerturbationMatrix3D(
@@ -117,10 +115,8 @@ def test_perturbation_matrix3d():
     assert chi < 1.5
 
     time = np.arange(0, 100, 1)
-    flux = (
-        np.random.normal(0, 0.01, size=(100, 169))
-        + 2 * dx[None, :] ** 2 * 3 * dy[None, :]
-    )
+    flux = np.random.normal(1, 0.01, size=(100, 169)) + dx[None, :] * dy[None, :]
+    # the perturbation model assumes the perturbation is around 1
     flux_err = np.ones((100, 169)) * 0.01
 
     for bin_method in ["downsample", "bin"]: