OxIonics / ionics_fits

Small python fitting library with an emphasis on Atomic Molecular and Optical Physics
Apache License 2.0
1 stars 0 forks source link

Phi not fixed in `Sinusoid` model #188

Closed amycxh closed 3 days ago

amycxh commented 5 days ago

The following should fix the parameter phi to 0 in the fit, and float x0 instead, but phi is not fixed.

model=Sinusoid()
model.parameters["phi"].fixed_to = 0.0
model_phase.parameters["x0"].fixed_to = None

Details

I have a simple experiment to fit to fake sinusoidal data using the Sinusoid model. I provide the phase in turns, and want to associate the fitted phase offset with the phase parameter, so I fix omega=2*np.pi, fix phi (in radians) to zero, and float x0 (in turns).

However, whether I set my modelled phi_offset_turns to 0, 0.25 or 0.5, I always get a fitted x0 = 0. On further inspection, printing the other fitted values, it seems phi is in fact not fixed to 0 as requested.

ionics_fragments.ndscan_extensions.analysis_exp:phase_analysis succeeded: x0=0.0 turns, a=0.500000(6), y0=0.500000(6), phi=-1.6(1) (significance 1.00)

ionics_fragments.ndscan_extensions.analysis_exp:phase_analysis succeeded: x0=0.0 turns, a=0.50(1), y0=0.500(10), phi=-3.14(8) (significance 1.00)

@hartytp suggested this was a quirk of the PeriodicModelParameter when the fixed_to value is on a boundary, but the default range for phi as define in the Sinusoid model is {-pi,pi} so 0 is not on the boundary: https://github.com/OxIonics/ionics_fits/blob/48968cf7742fd97009cb48615f6c5961b53c156d/ionics_fits/models/sinusoid.py#L40C1-L44C11

Example to reproduce

import numpy as np

import ndscan.experiment as nd

from ionics_fits.models.sinusoid import Sinusoid
from ionics_fragments.ndscan_extensions import AnalysisDescription, AnalysisExpFragment
from ionics_fragments.ndscan_extensions.utils import make_fragment_scan_exp

class MyFragment(AnalysisExpFragment):
    def build_fragment(self, phi_offset_turns=0.5):
        super().build_fragment()
        self.setattr_param(
            "phase",
            nd.FloatParam,
            "phase",
            default=0.0,
            min=0.0,
            max=1.0,
            unit="turns",
            scale=1.0,
        )
        self.setattr_result(
            "p",
            nd.FloatChannel,
            description="Bright fraction",
            min=0.0,
            max=1.0,
        )
        self.setattr_result(
            "p_err", min=0.0, display_hints={"error_bar_for": self.p.path}
        )
        model_phase = Sinusoid()
        model_phase.parameters["omega"].fixed_to = 2 * np.pi
        model_phase.parameters["a"].fixed_to = 0.5
        model_phase.parameters["y0"].fixed_to = 0.5
        model_phase.parameters["phi"].fixed_to = 0.0
        model_phase.parameters["x0"].fixed_to = None

        self.setattr_analysis(
            AnalysisDescription(
                model=Sinusoid(),
                x=self.phase,
                y=[self.p],
                fitter_args={"sigma": [self.p_err]},
                results={"x0": self.phase},#, "a": nd.FloatChannel("a"), "y0": nd.FloatChannel("y0"), "phi": nd.FloatChannel("phi")},
                log=True,
            )
        )
        self.phi_offset_turns = phi_offset_turns

    def run_once(self):
        N = 50
        p =  0.5 * np.sin(2 * np.pi * (self.phase.get() - self.phi_offset_turns)) + 0.5
        self.p.push(p)
        self.p_err.push(np.sqrt(p * (1-p) / N) + 1e-5)

ScanMyFragment = make_fragment_scan_exp(MyFragment)
hartytp commented 3 days ago

@amycxh spot the typo!

self.setattr_analysis(
            AnalysisDescription(
                model=Sinusoid(),
                x=self.phase,
                y=[self.p],
                fitter_args={"sigma": [self.p_err]},
                results={"x0": self.phase},#, "a": nd.FloatChannel("a"), "y0": nd.FloatChannel("y0"), "phi": nd.FloatChannel("phi")},
                log=True,
            )
        )

You're not passing the constrained model into the AnalysisDescription!

FWIW I did write a test case for this, which passes (see below). I'm going to close this on the assumption it's just a typo, but do reopen if you think that's not the case

import numpy as np

from ionics_fits.normal import NormalFitter
from ionics_fits.models.sinusoid import Sinusoid
from ..common import is_close

def test_188(plot_failures):
    """Test that fixing the phase parameter of a Sinusoid model works correctly"""
    model = Sinusoid()
    model.parameters["phi"].fixed_to = 0.0
    model.parameters["x0"].fixed_to = None

    params = {
        "phi": 0.,
        "omega": 2 * np.pi,
        "a": 1,
        "y0": 0.5,
        "x0": 0.25,
    }

    x = np.linspace(0, 1)
    y = model(x, **params)
    fit = NormalFitter(x=x, y=y, model=model)

    for param in params.keys():
        assert is_close(params[param], fit.values[param], 1e-3)
amycxh commented 3 days ago

🤦‍♀️ whoops! Thank you 😅