models.py

import logging
import math
import random
from enum import Enum
from pathlib import Path
from typing import Dict, List, Optional

import fsspec
import numpy as np
import pandas as pd
from pydantic import BaseModel, Field, field_validator, model_validator

logger = logging.getLogger(__name__)


def generate_reagent_field(**kwargs):
    """
    Create a Pydantic Field with reagent-specific metadata.
    """
    cost = kwargs.get("cost", 0.0)
    cost_ginkgo = kwargs.get("cost_ginkgo", 0.0)
    concentration = kwargs.get("concentration", 0.0)
    category = kwargs.get("category", "")
    components = kwargs.get("components", "")
    description = kwargs.get("description", "")
    out_of_stock = kwargs.get("out_of_stock", False)
    default = kwargs.get("default", 0.0)

    # Create the field description for Pydantic
    out_of_stock_text = " Out of stock" if out_of_stock else ""
    field_description = f"{description}, Components: {components}, Cost: {cost} USD/mL, Concentration: {concentration}, Category: {category}. {out_of_stock_text}"

    # Return a properly configured Pydantic Field
    f = Field(
        default=default,
        description=field_description,
        cost=cost,
        cost_ginkgo=cost_ginkgo,
        concentration=concentration,
        category=category,
        components=components,
        out_of_stock=out_of_stock,
        ge=0.0,
    )
    return f

class ReagentList(BaseModel):
    # Reagents
    # All values are in nanoliters with increments of 25nL

    # Control flags (excluded from serialization so they're not treated as reagents)
    skip_validation: bool = Field(default=False, exclude=True)

    # Core Salt & Buffer Components
    potassium_glutamate: float = generate_reagent_field(
        concentration="0.85M",
        category="salts",
        components="L-Glutamic acid potassium salt monohydrate",
        description="Potassium glutamate.",
        cost=0.0700,
        cost_ginkgo=0.0710,
        default=0.0,
        out_of_stock=False,
    )

    magnesium_glutamate: float = generate_reagent_field(
        concentration="0.5M",
        category="salts",
        components="L-Glutamic acid hemimagnesium salt tetrahydrate",
        description="Magnesium glutamate.",
        cost=0.065,
        cost_ginkgo=0.0627,
        default=0.0,
        out_of_stock=False,
    )

    # Buffers
    base_buffer: float = generate_reagent_field(
        concentration="fixed",
        category="buffers",
        components="1.5M KGlu, 0.026M MgGlu, 0.3M HEPES, 0.01M 17AA mix, 0.01M Cys, 0.01M Tyr",
        description="Base buffer volume is fixed at 2000nL, do not modify this value.",
        cost=0.22,
        cost_ginkgo=0.22,
        default=2000.0,
    ) # Base buffer volume is fixed at 2000nL

    hepes_koh: float = generate_reagent_field(
        concentration="1.0M",
        category="buffers",
        components="HEPES, potassium hydroxide",
        description="HEPES-KOH pH7.5.",
        cost=0.13,
        cost_ginkgo=0.254,
        default=0.0,
        out_of_stock=False,
    )  # HEPES-KOH  pH7.5, stock concentration 1.0M *

    # Core Transcription Components
    atp: float = generate_reagent_field(
        concentration="0.1M",
        category="transcription components",
        components="Adenosine 5′-triphosphate disodium salt hydrate, potassium hydroxide",
        description="ATP.",
        cost=1.539,
        cost_ginkgo=1.5388,
        default=0.0,
    )  # ATP, stock concentration 0.1M *
    ctp: float = generate_reagent_field(
        concentration="0.1M",
        category="transcription components",
        components="Cytidine 5′-triphosphate disodium salt hydrate, potassium hydroxide",
        description="CTP.",
        cost=27.9374,
        cost_ginkgo=27.9374,
        default=0.0,
    )  # CTP, stock concentration 0.1M *
    gtp: float = generate_reagent_field(
        concentration="0.1M",
        category="transcription components",
        components="Guanosine 5′-triphosphate sodium salt hydrate, potassium hydroxide",
        description="GTP.",
        cost=39.3431,
        cost_ginkgo=39.3431,
        default=0.0,
    )  # GTP, stock concentration 0.1M *
    utp: float = generate_reagent_field(
        concentration="0.1M",
        category="transcription components",
        components="Uridine 5′-triphosphate trisodium salt hydrate, potassium hydroxide",
        description="UTP.",
        cost=37.8462,
        cost_ginkgo=37.8462,
        default=0.0,
        out_of_stock=False,
    )  # UTP, stock concentration 0.1M *

    # Core Translation Components
    amino_acid_mix: float = generate_reagent_field(
        concentration="0.050M",
        category="translation components",
        components="18 amino acids in glacial acetic acid",
        description="18 amino acid mix. Does not include cysteine and glutamine.",
        cost=0.09325,
        cost_ginkgo=0.1013,
        default=0.0,
        out_of_stock=True,
    )  # 18 amino acid mix, stock concentration 0.05M
    amino_acid_mix_17: float = generate_reagent_field(
        concentration="0.050M",
        category="translation components",
        components="17 amino acids in glacial acetic acid",
        description="17 amino acid mix. Does not include tyrosine, cysteine, or glutamine.",
        cost=0.08275,
        cost_ginkgo=0.0938,
        default=0.0,
    )  # 17 amino acid mix, stock concentration 0.05M
    tyrosine: float = generate_reagent_field(
        concentration="0.050M",
        category="translation components",
        components="Tyrosine, Potassium hydroxide",
        description="Soluble tyrosine. pH 12.0",
        cost=0.0105,
        cost_ginkgo=0.0075,
        default=0.0,
    )  # Soluble tyrosine, stock concentration 0.05M
    cysteine: float = generate_reagent_field(
        concentration="0.2M",
        category="translation components",
        components="L-Cysteine",
        description="Cysteine.",
        cost=0.024,
        cost_ginkgo=0.0230,
        default=0.0,
    )  # Cysteine, stock concentration 0.2M

    # Core ATP Regeneration Components (E. coli CFPS)
    three_phosphoglycerate: float = generate_reagent_field(
        concentration="1M",
        category="ATP regeneration components",
        components="D-(−)-3-Phosphoglyceric acid disodium salt, potassium hydroxide",
        description="Three phosphoglycerate.",
        cost=77.29,
        cost_ginkgo=77.2867,
        default=0.0,
        out_of_stock=True,
    )  # 3PGA
    sodium_pyruvate: float = generate_reagent_field(
        concentration="100g/L",
        category="ATP regeneration components",
        components="Sodium pyruvate",
        description="Sodium pyruvate.",
        cost=0.16358,
        cost_ginkgo=0.0936,
        default=0.0,
    )  # Sodium pyruvate
    maltodextrin: float = generate_reagent_field(
        concentration="300g/L",
        category="ATP regeneration components",
        components="Maltodextrin",
        description="Maltodextrin.",
        cost=0.11571,
        cost_ginkgo=0.11571,
        default=0.0,
        out_of_stock=True,
    )  # Maltodextrin
    maltodextrin_17: float = generate_reagent_field(
        concentration="300g/L",
        category="ATP regeneration components",
        components="Maltodextrin 17",
        description="Maltodextrin 17.",
        cost=0.1224,
        cost_ginkgo=0.1224,
        default=0.0,
    )  # Maltodextrin 17: 17 refers to the chain length of the maltodextrin
    glucose: float = generate_reagent_field(
        concentration="200g/L",
        category="ATP regeneration components",
        components="Glucose",
        description="Glucose.",
        cost=0.00222,
        cost_ginkgo=0.002784,
        default=0.0,
    )  # Glucose
    maltose: float = generate_reagent_field(
        concentration="50g/L",
        category="ATP regeneration components",
        components="D-(+)-Maltose monohydrate",
        description="Maltose.",
        cost=0.01995,
        cost_ginkgo=0.01995,
        default=0.0,
    )  # Maltose
    ribose: float = generate_reagent_field(
        concentration="100g/L",
        category="ATP regeneration components",
        components="D-(-)-Ribose",
        description="Ribose.",
        cost=0.17984,
        cost_ginkgo=0.183,
        default=0.0,
    )  # Ribose
    pep_mono: float = generate_reagent_field(
        concentration="0.1M",
        category="ATP regeneration components",
        components="Phosphoenolpyruvic_acid_monocyclohexylammonium_salt",
        description="PEP mono",
        cost=2.2019,
        cost_ginkgo=2.2019,
        default=0.0,
    )  # PEP mono
    pep_tri: float = generate_reagent_field(
        concentration="0.5M",
        category="ATP regeneration components",
        components="Phosphoenolpyruvic_acid_tricyclohexylammonium_salt",
        description="PEP tricyclo",
        cost=5.43,
        cost_ginkgo=5.43,
        default=0.0,
        out_of_stock=True,
    )  # PEP tricyclo

    # Co factors
    folinic_acid: float = generate_reagent_field(
        concentration="10g/L",
        category="co factors",
        components="Folinic acid, potassium hydroxide",
        description="Folinic acid.",
        cost=8.9775,
        cost_ginkgo=4.78,
        default=0.0,
    )  # Folinic acid
    nad: float = generate_reagent_field(
        concentration="0.1M",
        category="co factors",
        components="β-Nicotinamide adenine dinucleotide, potassium hydroxide",
        description="NAD.",
        cost=33.348,
        cost_ginkgo=7.0992,
        default=0.0,
    )  # NAD
    nicotinamide: float = generate_reagent_field(
        concentration="0.1M",
        category="co factors",
        components="nicotinamide",
        description="nicotinamide",
        cost=0.002,
        cost_ginkgo=0.001478,
        default=0.0,
    )  # nicotinamide
    coa: float = generate_reagent_field(
        concentration="0.05M",
        category="co factors",
        components="Coenzyme A hydrate, potassium hydroxide",
        description="Coenzyme A.",
        cost=104.0,
        cost_ginkgo=157.3437,
        default=0.0,
        out_of_stock=True,
    )  # Coenzyme A -
    camp: float = generate_reagent_field(
        concentration="0.2M",
        category="co factors",
        components="Adenosine 3′,5′-cyclic monophosphate, potassium hydroxide",
        description="cAMP.",
        cost=7.48,
        cost_ginkgo=7.4797,
        default=0.0,
    )  # cAMP
    spermidine: float = generate_reagent_field(
        concentration="0.250M",
        category="co-factors",
        components="Spermidine",
        description="Spermidine.",
        cost=1.045,
        cost_ginkgo=1.1835,
        default=0.0,
    )  # Spermidine

    brij_35: float = generate_reagent_field(
        concentration="10% volume",
        category="detergents",
        components="Brij-35",
        description="Brij-35.",
        cost=0.0266,
        cost_ginkgo=0.0266,
        default=0.0,
    )  # Brij-35
    dmso: float = generate_reagent_field(
        concentration="60% volume",
        category="solvents",
        components="Dimethyl sulfoxide",
        description="DMSO.",
        cost=0.2028,
        cost_ginkgo=0.2028,
        default=0.0,
    )  # DMSO
    pyruvate_oxidase: float = generate_reagent_field(
        concentration="10.000g/L",
        category="other",
        components="pyruvate oxidase, potassium phosphate monobasic, potassium phosphate dibasic",
        description="Pyruvate oxidase.",
        cost=48.15,
        cost_ginkgo=48.15,
        default=0.0,
        out_of_stock=True,
    )  # Pyruvate oxidase
    catalase: float = generate_reagent_field(
        concentration="50000.000U/mL",
        category="other",
        components="catalase",
        description="Catalase (>50000U/ml)",
        cost=5.0,
        cost_ginkgo=5.0,
        default=0.0,
    )  # Catalase
    sodium_hexametaphosphate: float = generate_reagent_field(
        concentration="0.150M",
        category="other",
        components="sodium hexametaphosphate",
        description="Sodium hexametaphosphate.",
        cost=0.0044,
        cost_ginkgo=0.0044,
        default=0.0,
    )  # Sodium hexametaphosphate
    potassium_phosphate_monobasic: float = generate_reagent_field(
        concentration="0.500M",
        category="other",
        components="potassium phosphate monbasic",
        description="Potassium phosphate monobasic.",
        cost=0.01,
        cost_ginkgo=0.0041,
        default=0.0,
    )  # Potassium phosphate monobasic
    potassium_phosphate_dibasic: float = generate_reagent_field(
        concentration="0.500M",
        category="other",
        components="potassium phosphate dibasic",
        description="Potassium phosphate dibasic.",
        cost=0.035,
        cost_ginkgo=0.0144,
        default=0.0,
    )  # Potassium phosphate dibasic
    kpo_monobasic_mix: float = generate_reagent_field(
        concentration="0.500M",
        category="other",
        components="potassium phosphate monobasic, potassium phosphate dibasic",
        description="0.5M stock concentration, 1.6x monobasic:dibasic molar ratio, pH 7.0",
        cost=0.0177,
        cost_ginkgo=0.0081,
        default=0.0,
    )
    kpo_dibasic_mix: float = generate_reagent_field(
        concentration="0.500M",
        category="other",
        components="potassium phosphate monobasic, potassium phosphate dibasic",
        description="0.5M stock concentration, 1.6x dibasic:monobasic molar ratio, pH 7.0",
        cost=0.0238,
        cost_ginkgo=0.0104,
        default=0.0,
    )
    dilithium_acetyl_phosphate: float = generate_reagent_field(
        concentration="0.050M",
        category="other",
        components="dilithium acetyl phosphate, potassium hydroxide",
        description="Potassium phosphate dibasic.",
        cost=0.2730,
        cost_ginkgo=0.2730,
        default=0.0,
    )
    amp: float = generate_reagent_field(
        concentration="0.10M",
        category="other",
        components="Adenosine 5’-monophosphate disodium salt, acetic acid",
        description="Adenosine 5’-monophosphate disodium salt",
        cost=0.699,
        cost_ginkgo=0.6994,
        default=0.0,
    )
    cmp: float = generate_reagent_field(
        concentration="0.10M",
        category="other",
        components="Cytidine 5’-monophosphate disodium salt, acetic acid",
        description="Cytidine 5’-monophosphate disodium salt",
        cost=1.27,
        cost_ginkgo=1.2704,
        default=0.0,
    )
    gmp: float = generate_reagent_field(
        concentration="0.10M",
        category="other",
        components="Guanosine 5’-monophosphate disodium salt, acetic acid",
        description="Guanosine 5’-monophosphate disodium salt",
        cost=0.533,
        cost_ginkgo=0.5334,
        default=0.0,
    )
    ump: float = generate_reagent_field(
        concentration="0.10M",
        category="other",
        components="Uridine 5’-monophosphate disodium salt, acetic acid",
        description="Uridine 5’-monophosphate disodium salt",
        cost=0.839,
        cost_ginkgo=0.8394,
        default=0.0,
    )
    nmp_mix: float = generate_reagent_field(
        concentration="0.10M",
        category="other",
        components=(
            "Adenosine 5’-monophosphate disodium salt, "
            "Cytidine 5’-monophosphate disodium salt, "
            "Guanosine 5’-monophosphate disodium salt, "
            "Uridine 5’-monophosphate disodium salt, "
            "acetic acid"
        ),
        description="Mix of 0.115M AMP, 0.100M CMP, 0.100M GMP, 0.100M UMP",
        cost=3.476,
        cost_ginkgo=3.476,
        default=0.0,
    )
    adenosine: float = generate_reagent_field(
        concentration="0.025M",
        category="other",
        components="adenosine, potassium hydroxide",
        description="Adenosine",
        cost=0.002,
        cost_ginkgo=0.0202,
        default=0.0,
    )
    cytidine: float = generate_reagent_field(
        concentration="0.025M",
        category="other",
        components="cytidine",
        description="Cytidine",
        cost=0.0465,
        cost_ginkgo=0.0738,
        default=0.0,
    )
    guanosine: float = generate_reagent_field(
        concentration="0.025M",
        category="other",
        components="guanosine, potassium hydroxide",
        description="Guanosine",
        cost=0.01225,
        cost_ginkgo=0.0303,
        default=0.0,
    )
    uridine: float = generate_reagent_field(
        concentration="0.025M",
        category="other",
        components="uridine, potassium hydroxide",
        description="Uridine",
        cost=0.00275,
        cost_ginkgo=0.0268,
        default=0.0,
    )
    thymidine: float = generate_reagent_field(
        concentration="0.025M",
        category="other",
        components="thymidine, potassium hydroxide",
        description="Thymidine",
        cost=0.1080,
        cost_ginkgo=0.1080,
        default=0.0,
    )
    adenine: float = generate_reagent_field(
        concentration="0.025M",
        category="other",
        components="adenine, potassium hydroxide",
        description="Adenine",
        cost=0.0095,
        cost_ginkgo=0.0095,
        default=0.0,
    )
    guanine: float = generate_reagent_field(
        concentration="0.025M",
        category="other",
        components="guanine, potassium hydroxide",
        description="Guanine",
        cost=0.00325,
        cost_ginkgo=0.0034,
        default=0.0,
    )
    cytosine: float = generate_reagent_field(
        concentration="0.025M",
        category="other",
        components="cytosine, acetic acid",
        description="Cytosine",
        cost=0.034,
        cost_ginkgo=0.0339,
        default=0.0,
    )
    uracil: float = generate_reagent_field(
        concentration="0.025M",
        category="other",
        components="uracil, potassium hydroxide",
        description="Uracil",
        cost=0.00175,
        cost_ginkgo=0.00325,
        default=0.0,
    )
    thymine: float = generate_reagent_field(
        concentration="0.025M",
        category="other",
        components="thymine, potassium hydroxide",
        description="Thymine",
        cost=0.0069,
        cost_ginkgo=0.0069,
        default=0.0,
    )
    # Other
    lysate: float = generate_reagent_field(
        concentration="fixed",
        category="other",
        components="Cell lysate",
        description="Lysate volume is fixed at 5000nL, do not modify this value",
        cost=1.2,
        cost_ginkgo=20.0,
        default=5000.0,
    )  # Lysate volume is fixed at 5000nL
    nuclease_free_water: float = generate_reagent_field(
        concentration="n/a",
        category="other",
        components="Nuclease free water",
        description="Nuclease free water, use this to fill the rest of the volume up to 20000nL",
        cost=0.0,
        cost_ginkgo=0.0703,
        default=0.0,
    )  # Nuclease free water, use this to fill the rest of the volume up to 20000nL
    # DNA
    dna: float = generate_reagent_field(
        concentration="100nM",
        category="dna",
        components="plasmid DNA coding for HiBit tagged erythropoietin",
        description="DNA volume is fixed at 2000nL. Do not modify this value",
        cost=1.933,
        cost_ginkgo=6.418,
        default=2000.0,
    )  # DNA volume is fixed at 2000nL
    master_mix: float = generate_reagent_field(
        concentration="fixed",
        category="other",
        components="Master mix",
        description="Master mix volume is fixed at 0 for experimental samples. Do not modify this value",
        cost=0.0,
        cost_ginkgo=0.0,
        default=0.0,
    )  # Master mix volume is fixed at 2000nL

    @model_validator(mode="after")
    def validate_unavailable_reagents(self) -> bool:
        unavailable_reagents = [
            reagent for reagent, info in self.model_fields.items()
            if not info.exclude and info.json_schema_extra.get("out_of_stock", False)
        ]
        if not self.skip_validation:
            for reagent in unavailable_reagents:
                assert getattr(self, reagent) == 0, (
                    f"{reagent} is out of stock and must be set to 0"
                )
        return self

    @model_validator(mode="after")
    def validate_25nL_increments(self) -> bool:
        """
        Due to equipment limitations, we can only use 25nL increments.
        """
        for key, value in self.model_dump().items():
            if value != 0:
                assert int(value) % 25 == 0, (
                    f"{key} must be a multiple of 25, current value: {value}"
                )
        return self

    @model_validator(mode="after")
    def validate_non_zero_default_values(self) -> bool:
        """
        Values for non-zero reagents should match exactly what is defined in the model.
        """
        if self.skip_validation:
            return self
        for key, value in self.model_dump().items():
            default_value = self.model_fields[key].default
            if default_value > 0:
                assert value == default_value, (
                    f"{key} must be {default_value}, current value: {value}"
                )
        return self

    @property
    def final_volume(self) -> float:
        """
        Final volume of the reagent list.
        """
        return sum(self.model_dump().values())


class SampleType(str, Enum):
    POSITIVE_CONTROL = "positive_control"
    POSITIVE_CONTROL_MIXED = "positive_control_mixed"
    EXPERIMENTAL = "experimental"
    EMPTY = "empty"
    NEGATIVE = "negative"
    STANDARD = "standard"

    def __str__(self) -> str:
        return self.value


class Sample(BaseModel):
    sample_id: str  # unique identifier for the sample
    sample_type: SampleType  # type of the sample

    reagent_list: Optional[ReagentList] = None  # Reagent list
    metadata: Optional[dict] = None  # Metadata for the sample
    skip_validation: bool = False  # Skip validation for the sample

    @model_validator(mode="after")
    def validate_reagent_list(self) -> bool:
        """
        Reagent list is required for experimental samples.
        """
        if self.sample_type == SampleType.EXPERIMENTAL:
            assert self.reagent_list is not None, (
                "Reagent list is required for experimental samples"
            )
        elif self.sample_type not in {SampleType.POSITIVE_CONTROL_MIXED}:
            assert self.reagent_list is None, (
                "Reagent list is not allowed for positive, negative, or hibit standard samples"
            )
        return self

    @model_validator(mode="after")
    def validate_total_volume(
        self,
    ) -> bool:  # Not used in validator, to make it easy to
        if not self.skip_validation:
            if self.sample_type == SampleType.EXPERIMENTAL:
                total_volume = 20000  # 20 microliters
                assert self.reagent_list.final_volume == total_volume, (
                    f"Total volume of the reagent list must be {total_volume}nL, current {self.reagent_list.final_volume}nL"
                )
        return self

    @property
    def sample_cost(self) -> float:
        """
        Calculate the cost of the sample.
        """
        if self.reagent_list is None:
            if self.sample_type == SampleType.POSITIVE_CONTROL_MIXED or self.sample_type == SampleType.POSITIVE_CONTROL:
                return 0.02795 # Calculated from the cost of the reagents in the positive control mix
            else:
                return 0.0

        cost_per_sample = 0.0
        for reagent_name, reagent_value in self.reagent_list.model_dump().items():
            reagent_field = ReagentList.model_fields[reagent_name]
            cost = reagent_field.json_schema_extra["cost"]
            # cost is in $/mL, but reagent_value is in nL
            # 1 mL = 1,000,000 nL, so $/nL = cost / 1_000_000
            cost_per_sample += (cost / 1_000_000) * reagent_value
        return cost_per_sample



class BasePherastarResult(BaseModel):
    @staticmethod
    def parse_pherastar_well_luminosity(file_path: str) -> dict:
        """
        Parse a PHERAstar CSV and return a dict mapping well names (e.g., 'A10') to their luminescence values.
        """
        well_luminosity = {}
        with fsspec.open(file_path, "r", encoding="latin1") as f:
            lines = f.readlines()
        # Find the start of the measurement data
        for i, line in enumerate(lines):
            if line.strip().startswith("A01:"):
                data_start = i
                break
        else:
            raise ValueError("Measurement data not found in file")
        for l in lines[data_start:]:
            l = l.strip()
            if not l:
                continue
            well, value = l.split(":")
            well_luminosity[well.strip()] = float(value.strip())
        return well_luminosity



class LuminescenceResults(BasePherastarResult):
    luminescence: Dict[str, float]
    sfg_mw: float = 26800  # Da
    dilution_factor: float = ((20+10)/20)*20*20*20*2

    @classmethod
    def from_pherastar_data(cls, file_path: str) -> "LuminescenceResults":
        """
        Create a LuminescenceResults object from a pherastar data.
        """
        well_luminosity = cls.parse_pherastar_well_luminosity(file_path)
        return cls(luminescence=well_luminosity)

    def titer(self) -> Dict[str, float]:
        """
        Calculate the concentration of the sample.
        """
        assay_concentration = self.normalized_assay_concentration()
        for well, conc in assay_concentration.items():
            if conc is None:
                assay_concentration[well] = None
            else:   
                assay_concentration[well] = (conc * self.dilution_factor) * self.sfg_mw * 10**-9 # g/L
        return assay_concentration

    def normalized_assay_concentration(self) -> Dict[str, float]:
        """
        Normalize luminescence values to the mean of hibit standard curve. Result is in nM.
        """
        std_protein_nM_conc_dict = {
            "std_1": 1.000,
            "std_2": 0.500,
            "std_3": 0.100,
            "std_4": 0.050,
            "std_5": 0.010,
            "std_6": 0.005,
            "std_7": 0.001,
            "std_8": 0.0005
        }
        # Map standard names to their corresponding wells
        std_well_map = {
            "std_1": [f"A{col:02d}" for col in range(1, 5)],
            "std_2": [f"B{col:02d}" for col in range(1, 5)],
            "std_3": [f"C{col:02d}" for col in range(1, 5)],
            "std_4": [f"D{col:02d}" for col in range(1, 5)],
            "std_5": [f"E{col:02d}" for col in range(1, 5)],
            "std_6": [f"F{col:02d}" for col in range(1, 5)],
            "std_7": [f"G{col:02d}" for col in range(1, 5)],
            "std_8": [f"H{col:02d}" for col in range(1, 5)],
        }
        # Wells in columns 1-4 for rows I-P are negative controls
        neg_control_wells = []
        for row in [chr(i) for i in range(ord("I"), ord("P")+1)]:
            for col in range(1, 5):
                neg_control_wells.append(f"{row}{col:02d}")

        std_mean_lum = {}
        for std_name, std_wells in std_well_map.items():
            std_mean_lum[std_name] = np.mean([self.luminescence[well] for well in std_wells])
        # Calculate standard curve of best fit (linear regression in log-log space)
        # Prepare x (concentration) and y (mean luminescence) arrays
        std_concs = []
        std_lums = []
        for std_name in std_well_map:
            conc = std_protein_nM_conc_dict[std_name]
            lum = std_mean_lum[std_name]
            std_concs.append(conc)
            std_lums.append(lum)
        # Use only standards with nonzero concentration and positive luminescence
        std_concs = np.array(std_concs)
        std_lums = np.array(std_lums)
        mask = (std_concs > 0) & (std_lums > 0)
        log_concs = np.log10(std_concs[mask])
        log_lums = np.log10(std_lums[mask])
        # Fit linear regression: log_lums = slope * log_concs + intercept
        slope, intercept = np.polyfit(log_concs, log_lums, 1)
        # Store the fit for later use if needed
        self._std_curve_fit = {"slope": slope, "intercept": intercept}
        logger.info(f"Standard curve fit: slope={slope}, intercept={intercept}")

        # Map wells in columns 5-24 (actual samples) to concentrations using the standard curve
        # We'll store the results in a dictionary: {well: calculated_concentration}
        sample_concentrations = {}
        for well, lum in self.luminescence.items():
            if well not in neg_control_wells and lum is not None and lum > 0:
                # Use the fitted standard curve to calculate concentration from luminescence
                # log10(lum) = slope * log10(conc) + intercept
                # => log10(conc) = (log10(lum) - intercept) / slope
                log_lum = np.log10(lum)
                log_conc = (log_lum - intercept) / slope
                conc = 10 ** log_conc
                sample_concentrations[well] = conc
            else:
                sample_concentrations[well] = None
        return sample_concentrations


class FluorescenceResults(BasePherastarResult):
    fluorescence: Dict[str, float]
    dilution_factor: float = 1
    cvs: Dict[str, float] = {}

    @classmethod
    def from_pherastar_data(cls, file_path: str) -> "FluorescenceResults":
        """
        Create a FluorescenceResults object from a pherastar data.
        """
        well_fluorescence = cls.parse_pherastar_well_luminosity(file_path)
        return cls(fluorescence=well_fluorescence)

    def titer(self) -> Dict[str, float]:
        raise NotImplementedError("Titer calculation for fluorescence is not implemented yet")


class ConcentrationResults(BaseModel):
    concentration_g_L: Dict[str, float]


class ReagentFlags(BaseModel):
    flags: Dict[str, list[str]]


class ProteinTarget(str, Enum):
    EGFP_HIBIT_M8975491 = "egfp-hibit-m8975491"
    SFGFP_HIBIT_M8975481 = "sfgfp-hibit-m8975481"
    SFGFP_JEWETT_M9134733 = "sfgfp-jewett-m9134733"
    SFGFP_M8918409 = "sfgfp-m8918409"

    def __str__(self) -> str:
        return self.value


class Plate(BaseModel):
    """
    This class represents a plate of samples. Each sample is a 20 microliter reaction in a well. Besides experimental samples we want to run, we need to have:
    * Positive and negative controls - these are wells that consist our known-good mixture and buffer without necessary reagents. This is meant to show that our experimental conditions are working.
    * HiBit standards - these wells are reserved for HiBit standards. That means required concentrations of protein to create a standard curve of HiBit luminescence.
    * Replicates - both controls and samples themselves should be run in replicates. Currently we run 4 replicates per sample.

    After samples are defined, it's good idea to randomize the plate layout. This would ensure that effects stemming from plate geometry are minimized.
    """

    model_config = {"arbitrary_types_allowed": True, "validate_assignment": True}

    plate_id: str  # unique identifier for the plate
    run_id: Optional[str] = None  # optional run identifier for plate runs
    samples: List[Sample]
    n_rows: int = 16
    n_columns: int = 24

    protein_target: Optional[ProteinTarget] = None
    protein_molecular_weight: Optional[float] = None

    reserved_columns: list[int] = [1, 2, 3, 4]  # negative controls and hibit standards

    replicate_factor: int = 4

    # reserved wells
    positive_controls: int = 0
    positive_controls_mixed: int = 2

    # results
    luminescence_results: Optional[LuminescenceResults] = None
    fluorescence_results: Optional[FluorescenceResults] = None
    concentration_results: Optional[ConcentrationResults] = None

    # qc 
    fluorescence_qc: Optional[FluorescenceResults] = None
    reagent_flags: Optional[ReagentFlags] = None
    standards_metrics: Optional[dict[str, float]] = None
    reagent_lots: Optional[dict[str, str]] = None

    # metadata
    metadata: Optional[dict] = None

    skip_validation: bool = False

    @property
    def coordinates(self) -> List[str]:
        """
        Get the coordinates of the plate.
        """
        rows = [chr(i) for i in range(ord("A"), ord("P") + 1)]
        return [
            f"{row}{column:02d}"
            for column in range(1, self.n_columns + 1)
            if column not in self.reserved_columns
            for row in rows[:self.n_rows]
        ]

    @property
    def n_wells(self) -> int:
        return (self.n_columns - len(self.reserved_columns)) * self.n_rows

    @model_validator(mode="after")
    def validate_sample_count(self) -> bool:
        """
        Validate the number of samples on the plate.
        """
        if self.skip_validation:
            return self
        assert len(self.samples) == self.n_wells, (
            f"Number of samples: {len(self.samples)}, expected: {self.n_wells}"
        )
        return self

    @model_validator(mode="after")
    def validate_replicates(self) -> bool:
        """
        Ensure sufficient replicates are present.
        """
        if self.skip_validation:
            return self
        # Check control replicates
        assert len([s for s in self.samples if s.sample_type == SampleType.POSITIVE_CONTROL_MIXED]) >= (self.positive_controls_mixed) * self.replicate_factor, (
            f"Positive control mixed samples: {len([s for s in self.samples if s.sample_type == SampleType.POSITIVE_CONTROL_MIXED])}, expected: {(self.positive_controls_mixed) * self.replicate_factor}"
        )
        assert len([s for s in self.samples if s.sample_type == SampleType.POSITIVE_CONTROL]) >= (self.positive_controls) * self.replicate_factor, (
            f"Positive control samples: {len([s for s in self.samples if s.sample_type == SampleType.POSITIVE_CONTROL])}, expected: {(self.positive_controls) * self.replicate_factor}"
        )

        # Check experimental samples
        experimental_samples = [
            s for s in self.samples if s.sample_type == SampleType.EXPERIMENTAL
        ]

        # Group samples by their reagent list
        sample_groups = {}
        for sample in experimental_samples:
            # Convert reagent list to tuple for hashing
            reagent_tuple = tuple(sample.reagent_list.model_dump().items())
            if reagent_tuple not in sample_groups:
                sample_groups[reagent_tuple] = []
            sample_groups[reagent_tuple].append(sample)

        # Check if each group has exactly 'replicates' number of samples
        for group in sample_groups.values():
            assert len(group) >= self.replicate_factor, (
                f"Not enough replicates for sample {group[0].sample_id}, expected: {self.replicate_factor}, current: {len(group)}"
            )

        return self

    def get_plate_layout(self) -> list[list[str]]:
        """
        Get the plate layout. Physical plate has 16 rows and 24 columns.
        """
        return dict(zip(self.coordinates, self.samples))

    def randomize_sample_layout(self) -> list[list[str]]:
        """
        Unless model wants to cotnrol plate layout itself, it's good idea to randomize it.
        """
        random.shuffle(self.samples)

    @property
    def available_independent_exp(self) -> int:
        """
        Returns the number of independent experimental samples that can be run on the plate,
        accounting for required controls and replication. All plate sample types, including experimental,
        should be replicated according to replicate_factor.
        """
        return (
            self.n_wells
            - self.positive_controls * self.replicate_factor
            - len(self.reserved_columns) * self.n_rows
        ) / self.replicate_factor

    def get_plate_dataframe(self, include_qc=False) -> pd.DataFrame:
        """
        Get the plate as a dataframe.
        """
        # Exclude non-reagent fields (e.g., control flags) from reagent columns
        reagent_columns = [k for k, f in ReagentList.model_fields.items() if not f.exclude]
        rows = []
        sample_layout = self.get_plate_layout()
        if self.luminescence_results:
            titer = self.luminescence_results.titer()

        for sample_coord, sample in sample_layout.items():
            row = {
                "sample_id": sample.sample_id,
                "sample_type": sample.sample_type,
                "well_coordinate": sample_coord,
            }
            if sample.sample_type == SampleType.EXPERIMENTAL:
                for k, v in sample.reagent_list.model_dump().items():
                    row[k] = v
            elif sample.sample_type == SampleType.POSITIVE_CONTROL_MIXED and sample.reagent_list:
                for k, v in sample.reagent_list.model_dump().items():
                    row[k] = v
            else:
                for k in reagent_columns:
                    row[k] = 0.0
            if sample.sample_type == SampleType.STANDARD:
                if self.standards_metrics:
                    for metric, value in self.standards_metrics.items():
                        row[f"standards {metric}"] = value
            if self.luminescence_results:
                row["luminescence"] = self.luminescence_results.luminescence[sample_coord]
                row["titer"] = titer[sample_coord]
            if self.fluorescence_results:
                row["fluorescence"] = self.fluorescence_results.fluorescence[sample_coord]
                if self.fluorescence_results.cvs:
                    row["fluorescence_cv"] = self.fluorescence_results.cvs.get(sample_coord, None)
            if include_qc:
                if self.fluorescence_qc:
                    row["fluorescence_qc"] = self.fluorescence_qc.fluorescence[sample_coord]
            if self.reagent_flags:
                row["reagent_flags"] = "|".join(self.reagent_flags.flags.get(sample_coord, []))
            if self.concentration_results:
                row["concentration (g/L)"] = self.concentration_results.concentration_g_L.get(sample_coord, None)
            row["sample_cost"] = sample.sample_cost
            if sample.metadata:
                for k, v in sample.metadata.items():
                    row[k] = v
            rows.append(row)
        return pd.DataFrame(rows)