Chapter 10: Advanced OOP Patterns - Architecting Scientific Software

# Can you identify the OOP patterns in this NumPy/Matplotlib code?
import numpy as np
import matplotlib.pyplot as plt

# Pattern 1: Factory
arr = np.array([1, 2, 3])  # What pattern?

# Pattern 2: Context Manager  
with plt.style.context('dark_background'):  # What pattern?
    fig, ax = plt.subplots()  # Pattern 3: ?

# Pattern 4: Method chaining
result = arr.reshape(3, 1).mean(axis=0)  # What pattern?

from abc import ABC, abstractmethod
import numpy as np

# First, see the problem ABCs solve in astronomy
class BadCCDArray:
    def read_counts(self): 
        return np.array([100, 200, 300])  # photon counts

class BadIRArray:
    def get_flux(self):  # Different method name!
        return np.array([50, 60, 70])  # flux in erg/s/cm²

# Without ABCs, incompatible interfaces cause runtime errors
# You've seen this problem when different instruments use different formats

from abc import ABC, abstractmethod

class AstronomicalDetector(ABC):
    """Abstract base defining detector interface in CGS units."""
    
    @abstractmethod
    def read_frame(self) -> np.ndarray:
        """Read one frame of data in erg/s/cm²."""
        pass
    
    @abstractmethod
    def get_noise(self) -> float:
        """Return noise in erg/s/cm²."""
        pass
    
    @abstractmethod
    def get_area_cm2(self) -> float:
        """Return detector area in cm²."""
        pass
    
    # Concrete methods provide shared functionality
    def photons_to_flux(self, counts, wavelength_cm):
        """Convert photon counts to flux in erg/s/cm².
        
        E = hν = hc/λ where:
        h = 6.626e-27 erg·s (Planck in CGS)
        c = 2.998e10 cm/s (speed of light in CGS)
        """
        h = 6.626e-27  # erg·s
        c = 2.998e10   # cm/s
        energy_per_photon = h * c / wavelength_cm
        return counts * energy_per_photon

# Now implementations MUST follow the interface
class CCDDetector(AstronomicalDetector):
    """CCD implementation with CGS units."""
    
    def __init__(self, area_cm2=4.0, gain=1.0):
        self.area = area_cm2  # typically 2x2 cm for astronomy CCDs
        self.gain = gain
        self.read_noise = 5.0  # electrons
    
    def read_frame(self) -> np.ndarray:
        """Read CCD frame, return flux in erg/s/cm²."""
        # Simulate Poisson photon noise
        photons = np.random.poisson(1000, size=(1024, 1024))
        # Convert to flux at 550nm (5.5e-5 cm)
        return self.photons_to_flux(photons, 5.5e-5)
    
    def get_noise(self) -> float:
        """RMS noise in erg/s/cm²."""
        # Convert read noise to flux
        return self.photons_to_flux(self.read_noise, 5.5e-5)
    
    def get_area_cm2(self) -> float:
        """CCD area in cm²."""
        return self.area

class InfraredArray(AstronomicalDetector):
    """IR array implementation with CGS units."""
    
    def __init__(self, area_cm2=1.0):
        self.area = area_cm2  # smaller pixels for IR
        self.dark_current = 0.1  # e-/s at 77K
    
    def read_frame(self) -> np.ndarray:
        """Read IR frame at 2.2 microns."""
        # IR arrays have different noise characteristics
        signal = np.random.randn(256, 256) * 10 + 100
        # Convert to flux at 2.2 microns (2.2e-4 cm)
        return self.photons_to_flux(signal, 2.2e-4)
    
    def get_noise(self) -> float:
        """Dark current noise in erg/s/cm²."""
        return self.photons_to_flux(self.dark_current, 2.2e-4)
    
    def get_area_cm2(self) -> float:
        return self.area

# Process any detector uniformly
def calculate_sensitivity(detector: AstronomicalDetector) -> float:
    """Calculate sensitivity for any detector in erg/s/cm²."""
    signal = detector.read_frame()
    noise = detector.get_noise()
    area = detector.get_area_cm2()
    
    # Signal-to-noise ratio calculation
    snr = np.mean(signal) / noise if noise > 0 else 0
    sensitivity = noise * 3 / area  # 3-sigma detection limit
    
    return sensitivity

# Works with ANY detector implementing the interface!
ccd = CCDDetector(area_cm2=4.0)
ir = InfraredArray(area_cm2=1.0)

print(f"CCD sensitivity: {calculate_sensitivity(ccd):.2e} erg/s/cm³")
print(f"IR sensitivity: {calculate_sensitivity(ir):.2e} erg/s/cm³")

from abc import ABC, abstractmethod
from typing import Tuple
import numpy as np

class SpectralAnalyzer(ABC):
    """Advanced ABC using NumPy arrays and CGS units."""
    
    @property
    @abstractmethod
    def wavelength_range_cm(self) -> Tuple[float, float]:
        """Wavelength coverage in cm."""
        pass
    
    @abstractmethod
    def find_lines(self, flux: np.ndarray, 
                   threshold_sigma: float = 3.0) -> np.ndarray:
        """Find emission lines in flux [erg/s/cm²/cm].
        
        Returns wavelengths in cm of detected lines.
        """
        pass
    
    def flux_to_photons(self, flux_cgs, wavelength_cm):
        """Convert flux to photon rate using CGS units.
        
        flux: erg/s/cm²/cm
        returns: photons/s/cm²/cm
        """
        h = 6.626e-27  # erg·s
        c = 2.998e10   # cm/s
        return flux_cgs * wavelength_cm / (h * c)

class OpticalSpectrograph(SpectralAnalyzer):
    """Optical spectrograph with CGS units."""
    
    @property
    def wavelength_range_cm(self):
        # 380-750 nm in cm
        return (3.8e-5, 7.5e-5)
    
    def find_lines(self, flux, threshold_sigma=3.0):
        """Find emission lines using median filtering."""
        # Simple line finding with NumPy
        median = np.median(flux)
        std = np.std(flux)
        peaks = np.where(flux > median + threshold_sigma * std)[0]
        
        # Convert indices to wavelengths in cm
        wave_start, wave_end = self.wavelength_range_cm
        wavelengths = np.linspace(wave_start, wave_end, len(flux))
        
        return wavelengths[peaks]

# Using the abstract interface with NumPy
spectrograph = OpticalSpectrograph()
mock_flux = np.random.randn(1000) * 1e-12 + 5e-12  # erg/s/cm²/cm
lines = spectrograph.find_lines(mock_flux)
print(f"Found {len(lines)} emission lines")
print(f"Wavelength range: {spectrograph.wavelength_range_cm} cm")

from dataclasses import dataclass, field
from typing import List, Optional
import numpy as np

# Before dataclasses - lots of boilerplate
class OldStyleObservation:
    def __init__(self, target: str, ra_deg: float, dec_deg: float,
                 exposure_s: float, airmass: float, 
                 flux_cgs: Optional[float] = None):
        self.target = target
        self.ra_deg = ra_deg
        self.dec_deg = dec_deg
        self.exposure_s = exposure_s
        self.airmass = airmass
        self.flux_cgs = flux_cgs
    
    def __repr__(self):
        return (f"OldStyleObservation(target={self.target}, "
                f"ra_deg={self.ra_deg}, ...)")
    
    def __eq__(self, other):
        return (self.target == other.target and 
                self.ra_deg == other.ra_deg and ...)

# With dataclasses - automatic and clean!
@dataclass
class Observation:
    """Astronomical observation with automatic methods."""
    target: str
    ra_deg: float
    dec_deg: float
    exposure_s: float
    airmass: float
    flux_cgs: Optional[float] = None  # erg/s/cm²
    
    def __post_init__(self):
        """Validate after auto-generated __init__."""
        if not 0 <= self.ra_deg <= 360:
            raise ValueError(f"RA {self.ra_deg} outside [0,360]")
        if not -90 <= self.dec_deg <= 90:
            raise ValueError(f"Dec {self.dec_deg} outside [-90,90]")
        if self.airmass < 1.0:
            raise ValueError(f"Airmass {self.airmass} < 1.0")
    
    def extinction_magnitudes(self) -> float:
        """Atmospheric extinction in magnitudes.
        
        Uses Rayleigh scattering at sea level.
        """
        k_extinction = 0.15  # mag/airmass at 550nm
        return k_extinction * (self.airmass - 1.0)

# Automatic __init__, __repr__, __eq__ and more!
obs = Observation(
    target="M31",
    ra_deg=10.68,
    dec_deg=41.27,
    exposure_s=300,
    airmass=1.2,
    flux_cgs=1.5e-12
)

print(obs)  # Nice automatic __repr__
print(f"Extinction: {obs.extinction_magnitudes():.3f} mag")

from __future__ import annotations  # Enables forward references
from dataclasses import dataclass, field, asdict
import numpy as np
from typing import ClassVar, Optional

@dataclass
class SpectralLine:
    """Spectral line with CGS units and validation."""
    
    # Class variable (shared by all instances)
    c_cgs: ClassVar[float] = 2.998e10  # cm/s
    
    # Instance variables with types
    wavelength_cm: float
    flux_cgs: float  # erg/s/cm²
    width_cm: float
    name: str = "Unknown"
    
    # Computed field with default_factory
    measurements: List[float] = field(default_factory=list)
    
    # Field with metadata
    snr: float = field(default=0.0, metadata={"unit": "dimensionless"})
    
    def __post_init__(self):
        """Calculate SNR after initialization."""
        if self.flux_cgs > 0 and self.width_cm > 0:
            # Simple SNR estimate
            self.snr = self.flux_cgs / (self.width_cm * 1e-13)
    
    @property
    def velocity_width_km_s(self) -> float:
        """Line width in velocity space (km/s)."""
        return (self.width_cm / self.wavelength_cm) * self.c_cgs / 1e5
    
    @property
    def frequency_hz(self) -> float:
        """Frequency in Hz."""
        return self.c_cgs / self.wavelength_cm

# Create H-alpha line
h_alpha = SpectralLine(
    wavelength_cm=6.563e-5,
    flux_cgs=1e-13,
    width_cm=1e-7,
    name="H-alpha"
)

print(f"H-alpha: {h_alpha.frequency_hz:.2e} Hz")
print(f"Velocity width: {h_alpha.velocity_width_km_s:.1f} km/s")
print(f"SNR: {h_alpha.snr:.1f}")

# Convert to dictionary (useful for saving)
line_dict = asdict(h_alpha)
print(f"As dict: {line_dict}")

# Frozen (immutable) dataclass for constants
@dataclass(frozen=True)
class PhysicalConstants:
    """Immutable physical constants in CGS."""
    c: float = 2.998e10      # cm/s
    h: float = 6.626e-27     # erg·s
    k: float = 1.381e-16     # erg/K
    m_e: float = 9.109e-28   # g
    m_p: float = 1.673e-24   # g
    
    def photon_energy(self, wavelength_cm: float) -> float:
        """Photon energy in ergs."""
        return self.h * self.c / wavelength_cm

# Constants are immutable
constants = PhysicalConstants()
# constants.c = 3e10  # This would raise FrozenInstanceError
print(f"Photon at 550nm: {constants.photon_energy(5.5e-5):.2e} erg")

import asyncio
import time
import numpy as np

# Synchronous version - everything waits
def sync_observe(target: str) -> dict:
    """Traditional blocking observation."""
    print(f"Starting observation of {target}")
    
    # Each step blocks
    print("Slewing telescope...")
    time.sleep(2)  # Simulate 2 second slew
    
    print("Exposing CCD...")
    time.sleep(3)  # Simulate 3 second exposure
    
    print("Reading CCD...")
    time.sleep(1)  # Simulate 1 second readout
    
    return {"target": target, "counts": 1000}

# Time the synchronous version
start = time.perf_counter()
result = sync_observe("M31")
sync_time = time.perf_counter() - start
print(f"Synchronous time: {sync_time:.1f}s\n")

# Asynchronous version - operations can overlap
async def slew_telescope(ra_deg: float, dec_deg: float):
    """Async telescope slewing."""
    print(f"Slewing to RA={ra_deg:.1f}, Dec={dec_deg:.1f}")
    await asyncio.sleep(2)  # Non-blocking sleep
    print("Slew complete")
    return True

async def expose_ccd(exposure_s: float):
    """Async CCD exposure."""
    print(f"Starting {exposure_s}s exposure")
    await asyncio.sleep(3)  # Simulated exposure
    
    # Generate mock data in CGS units
    photons = np.random.poisson(1000, size=(100, 100))
    flux_cgs = photons * 3.6e-12  # erg/s/cm² at 550nm
    
    print("Exposure complete")
    return flux_cgs

async def configure_spectrograph(resolution: int):
    """Async spectrograph configuration."""
    print(f"Configuring spectrograph R={resolution}")
    await asyncio.sleep(1)
    print("Spectrograph ready")
    return resolution

async def async_observe(target: str, ra: float, dec: float):
    """Parallel observation with multiple instruments."""
    print(f"Starting async observation of {target}")
    
    # Launch operations concurrently!
    tasks = [
        slew_telescope(ra, dec),
        expose_ccd(30.0),
        configure_spectrograph(5000)
    ]
    
    # Wait for all to complete
    results = await asyncio.gather(*tasks)
    
    return {
        "target": target,
        "pointed": results[0],
        "flux": results[1].mean(),
        "resolution": results[2]
    }

# Run async version
async def main():
    start = time.perf_counter()
    result = await async_observe("M31", 10.68, 41.27)
    async_time = time.perf_counter() - start
    
    print(f"\nAsync time: {async_time:.1f}s")
    print(f"Speedup: {sync_time/async_time:.1f}x")
    print(f"Mean flux: {result['flux']:.2e} erg/s/cm²")

# Note: In Jupyter, use await directly
# In scripts, use asyncio.run(main())
await main()

from dataclasses import dataclass
from typing import List, Optional
import asyncio

# Note: In Jupyter notebooks, use 'await' directly in cells
# In Python scripts, use asyncio.run(main())
# Some Jupyter versions may require nest_asyncio:
# import nest_asyncio
# nest_asyncio.apply()

@dataclass
class TelescopeState:
    """Current telescope state."""
    ra_deg: float = 0.0
    dec_deg: float = 0.0
    filter: str = "V"
    focus_mm: float = 0.0
    dome_azimuth_deg: float = 0.0

class AsyncTelescope:
    """Async telescope controller with realistic operations."""
    
    def __init__(self, name: str):
        self.name = name
        self.state = TelescopeState()
        self.is_tracking = False
    
    async def slew_to(self, ra: float, dec: float) -> float:
        """Slew to coordinates, return time taken."""
        # Calculate slew time based on distance
        distance = np.sqrt((ra - self.state.ra_deg)**2 + 
                          (dec - self.state.dec_deg)**2)
        slew_time = min(distance / 2.0, 30.0)  # 2 deg/s, max 30s
        
        print(f"{self.name}: Slewing {distance:.1f}° ({slew_time:.1f}s)")
        await asyncio.sleep(slew_time)
        
        self.state.ra_deg = ra
        self.state.dec_deg = dec
        self.is_tracking = True
        
        return slew_time
    
    async def change_filter(self, filter_name: str):
        """Change filter wheel."""
        if filter_name != self.state.filter:
            print(f"{self.name}: Changing to {filter_name} filter")
            await asyncio.sleep(5.0)  # Filter wheel rotation
            self.state.filter = filter_name
    
    async def auto_focus(self) -> float:
        """Autofocus routine, return FWHM in arcsec."""
        print(f"{self.name}: Starting autofocus")
        
        best_focus = 0.0
        best_fwhm = float('inf')
        
        # Simulate focus sweep
        for focus in np.linspace(-2, 2, 5):
            await asyncio.sleep(1.0)  # Move focus
            
            # Simulate FWHM measurement
            fwhm = abs(focus) + np.random.random() + 0.8
            if fwhm < best_fwhm:
                best_fwhm = fwhm
                best_focus = focus
        
        self.state.focus_mm = best_focus
        print(f"{self.name}: Focus complete, FWHM={best_fwhm:.2f}\"")
        return best_fwhm

# Async observation sequence
async def observe_targets(telescope: AsyncTelescope, 
                         targets: List[tuple]):
    """Observe multiple targets efficiently."""
    
    for name, ra, dec, filter_name in targets:
        print(f"\n--- Observing {name} ---")
        
        # Parallel operations where possible
        tasks = [
            telescope.slew_to(ra, dec),
            telescope.change_filter(filter_name)
        ]
        
        # Slew and filter change happen in parallel
        await asyncio.gather(*tasks)
        
        # Focus if needed (every 5 targets)
        if np.random.random() < 0.2:
            await telescope.auto_focus()
        
        # Simulate exposure
        print(f"Exposing {name} in {filter_name}")
        await asyncio.sleep(10.0)
        
        print(f"Completed {name}")

# Run observation sequence
async def night_observations():
    telescope = AsyncTelescope("10m Keck")
    
    targets = [
        ("M31", 10.68, 41.27, "V"),
        ("M42", 83.82, -5.39, "R"),
        ("M51", 202.47, 47.20, "B"),
    ]
    
    start = time.perf_counter()
    await observe_targets(telescope, targets)
    total_time = time.perf_counter() - start
    
    print(f"\nTotal observation time: {total_time:.1f}s")

await night_observations()

# Descriptors for unit-safe astronomy
class UnitProperty:
    """Descriptor ensuring CGS units are maintained."""
    
    def __init__(self, unit_name: str, cgs_conversion: float = 1.0,
                 min_val: Optional[float] = None):
        self.unit_name = unit_name
        self.cgs_conversion = cgs_conversion
        self.min_val = min_val
    
    def __set_name__(self, owner, name):
        """Called when descriptor is assigned to class."""
        self.name = f"_{name}"
        self.public_name = name
    
    def __get__(self, obj, objtype=None):
        """Get value in CGS units."""
        if obj is None:
            return self
        return getattr(obj, self.name, None)
    
    def __set__(self, obj, value):
        """Set value with validation and conversion."""
        # Handle tuple of (value, unit)
        if isinstance(value, tuple):
            val, unit = value
            if unit != self.unit_name:
                # Convert to CGS
                if unit == "km" and self.unit_name == "cm":
                    val *= 1e5
                elif unit == "AU" and self.unit_name == "cm":
                    val *= 1.496e13
                else:
                    raise ValueError(f"Cannot convert {unit} to {self.unit_name}")
        else:
            val = value
        
        # Validate
        if self.min_val is not None and val < self.min_val:
            raise ValueError(f"{self.public_name} = {val} below minimum {self.min_val}")
        
        setattr(obj, self.name, val)

class CelestialObject:
    """Object with unit-safe properties."""
    
    # Descriptors enforce CGS units
    distance = UnitProperty("cm", min_val=0)
    mass = UnitProperty("g", min_val=0)
    radius = UnitProperty("cm", min_val=0)
    temperature = UnitProperty("K", min_val=0)
    
    def __init__(self, name: str):
        self.name = name
    
    @property
    def escape_velocity_cm_s(self) -> float:
        """Escape velocity in cm/s."""
        if self.mass and self.radius:
            G = 6.674e-8  # cm³/g/s²
            return np.sqrt(2 * G * self.mass / self.radius)
        return 0.0

# Unit-safe usage
earth = CelestialObject("Earth")
earth.distance = (1, "AU")  # Converts to cm automatically
earth.mass = 5.972e27  # grams
earth.radius = (6371, "km")  # Converts to cm

print(f"Earth distance: {earth.distance:.2e} cm")
print(f"Earth radius: {earth.radius:.2e} cm")
print(f"Escape velocity: {earth.escape_velocity_cm_s/1e5:.1f} km/s")

# Metaclass for instrument registry (like Astropy's registry)
class InstrumentMeta(type):
    """Metaclass that auto-registers instruments."""
    
    _registry = {}
    
    def __new__(mcs, name, bases, namespace):
        # Create the class
        cls = super().__new__(mcs, name, bases, namespace)
        
        # Register if it has instrument_code
        if 'instrument_code' in namespace:
            code = namespace['instrument_code']
            mcs._registry[code] = cls
            print(f"Registered {name} as '{code}'")
        
        return cls
    
    @classmethod
    def create(mcs, code: str, **kwargs):
        """Factory method using registry."""
        if code not in mcs._registry:
            available = list(mcs._registry.keys())
            raise ValueError(f"Unknown instrument '{code}'. "
                           f"Available: {available}")
        return mcs._registry[code](**kwargs)
    
    @classmethod
    def list_instruments(mcs):
        """List all registered instruments."""
        return list(mcs._registry.keys())

class Instrument(metaclass=InstrumentMeta):
    """Base class with metaclass."""
    
    def __init__(self, **config):
        self.config = config

# Classes auto-register on definition!
class WFC3(Instrument):
    """Hubble Wide Field Camera 3."""
    instrument_code = "WFC3"
    
    def __init__(self, channel="UVIS", **config):
        super().__init__(**config)
        self.channel = channel  # UVIS or IR

class MIRI(Instrument):
    """JWST Mid-Infrared Instrument."""
    instrument_code = "MIRI"
    
    def __init__(self, mode="imaging", **config):
        super().__init__(**config)
        self.mode = mode

class HARPS(Instrument):
    """High Accuracy Radial velocity Planet Searcher."""
    instrument_code = "HARPS"
    
    def __init__(self, resolution=115000, **config):
        super().__init__(**config)
        self.resolution = resolution

# Use the registry
print(f"Available: {InstrumentMeta.list_instruments()}")

# Create instruments by code
miri = InstrumentMeta.create("MIRI", mode="spectroscopy")
harps = InstrumentMeta.create("HARPS")

print(f"Created {miri.__class__.__name__} in {miri.mode} mode")
print(f"HARPS resolution: R={harps.resolution}")

import sys
import time
from functools import lru_cache

# Memory optimization with __slots__
class RegularStar:
    """Normal class - flexible but memory-hungry."""
    
    def __init__(self, ra, dec, mag, parallax):
        self.ra = ra  # degrees
        self.dec = dec  # degrees  
        self.mag = mag  # magnitude
        self.parallax = parallax  # milliarcsec

class OptimizedStar:
    """Using __slots__ for memory efficiency."""
    
    __slots__ = ['ra', 'dec', 'mag', 'parallax', '_distance_cache']
    
    def __init__(self, ra, dec, mag, parallax):
        self.ra = ra
        self.dec = dec
        self.mag = mag
        self.parallax = parallax
        self._distance_cache = None
    
    @property
    def distance_pc(self):
        """Distance in parsecs (cached)."""
        if self._distance_cache is None:
            if self.parallax > 0:
                self._distance_cache = 1000.0 / self.parallax
            else:
                self._distance_cache = float('inf')
        return self._distance_cache

# Compare memory usage
regular = RegularStar(10.68, 41.27, 3.44, 24.36)
optimized = OptimizedStar(10.68, 41.27, 3.44, 24.36)

print(f"Regular: {sys.getsizeof(regular.__dict__)} bytes")
print(f"Optimized: {sys.getsizeof(optimized)} bytes")

# For Gaia catalog (2 billion stars), the difference matters!
n_stars = 100000  # Subset for demo
regular_mem = n_stars * sys.getsizeof(regular.__dict__)
optimized_mem = n_stars * 56  # Approximate slotted size

print(f"\n{n_stars:,} stars:")
print(f"Regular: {regular_mem/1e6:.1f} MB")
print(f"Optimized: {optimized_mem/1e6:.1f} MB")
print(f"Savings: {(regular_mem-optimized_mem)/1e6:.1f} MB")

# LRU cache for expensive calculations
class GalaxySpectrum:
    """Galaxy spectrum with expensive computations."""
    
    def __init__(self, wavelengths_cm, fluxes_cgs):
        self.wavelengths = wavelengths_cm  # cm
        self.fluxes = fluxes_cgs  # erg/s/cm²/cm
    
    @lru_cache(maxsize=128)
    def find_redshift(self, template_lines_cm):
        """Find redshift by template matching (expensive)."""
        print("Computing redshift...")  # Shows when cache misses
        
        # Simplified cross-correlation
        best_z = 0.0
        best_corr = -float('inf')
        
        for z in np.linspace(0, 0.1, 100):
            # Shift template lines
            shifted = [line * (1 + z) for line in template_lines_cm]
            
            # Simple correlation metric
            corr = sum(1 for line in shifted 
                      if any(abs(line - w) < 1e-7 
                            for w in self.wavelengths))
            
            if corr > best_corr:
                best_corr = corr
                best_z = z
        
        return best_z
    
    def velocity_km_s(self, template_lines_cm):
        """Recession velocity from cached redshift."""
        z = self.find_redshift(template_lines_cm)
        c_km_s = 2.998e5
        return z * c_km_s

# Test caching
wavelengths = np.linspace(4e-5, 7e-5, 1000)  # 400-700 nm in cm
fluxes = np.random.randn(1000) * 1e-13

galaxy = GalaxySpectrum(wavelengths, fluxes)

# H-alpha and H-beta rest wavelengths in cm
h_lines = tuple([4.861e-5, 6.563e-5])

# First call - computes
z1 = galaxy.find_redshift(h_lines)
print(f"Redshift: {z1:.4f}")

# Second call - cached!
z2 = galaxy.find_redshift(h_lines)
print(f"Cached: {z2:.4f}")

# Cache info
print(f"Cache stats: {galaxy.find_redshift.cache_info()}")

"""
Part A: Abstract base and mixins (10 minutes)
Design flexible instrument system with CGS units
"""

from abc import ABC, abstractmethod
import asyncio
import numpy as np
from dataclasses import dataclass
from typing import Optional

# Abstract base for all instruments
class ObservatoryInstrument(ABC):
    """Base class for observatory instruments."""
    
    @abstractmethod
    async def initialize(self) -> bool:
        """Initialize instrument."""
        pass
    
    @abstractmethod
    async def acquire(self, integration_s: float) -> np.ndarray:
        """Acquire data for given time."""
        pass
    
    @abstractmethod
    def get_sensitivity_cgs(self) -> float:
        """Return sensitivity in erg/s/cm²."""
        pass

# Mixins for capabilities
class TemperatureControlMixin:
    """Temperature control for cooled instruments."""
    
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.temperature_k = 293  # Room temp
        self.target_temp_k = 150  # Target for CCDs
    
    async def cool_down(self):
        """Cool to operating temperature."""
        print(f"Cooling from {self.temperature_k}K to {self.target_temp_k}K")
        
        while self.temperature_k > self.target_temp_k:
            await asyncio.sleep(0.1)  # Simulate cooling
            self.temperature_k -= 10
            
        print(f"Reached {self.temperature_k}K")
        return True

class CalibrationMixin:
    """Calibration capability."""
    
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.dark_frame = None
        self.flat_frame = None
    
    async def take_dark(self, integration_s: float):
        """Take dark frame."""
        print(f"Taking {integration_s}s dark")
        await asyncio.sleep(integration_s)
        self.dark_frame = np.random.randn(100, 100) * 10
    
    async def take_flat(self):
        """Take flat field."""
        print("Taking flat field")
        await asyncio.sleep(1)
        self.flat_frame = np.ones((100, 100)) + np.random.randn(100, 100) * 0.01

# Concrete implementations
class CCDCamera(ObservatoryInstrument, TemperatureControlMixin, 
                CalibrationMixin):
    """CCD camera with cooling and calibration."""
    
    def __init__(self, name: str, pixel_size_cm: float = 1.5e-3):
        super().__init__()
        self.name = name
        self.pixel_size = pixel_size_cm  # 15 microns typical
        self.initialized = False
    
    async def initialize(self) -> bool:
        """Initialize CCD with cooling."""
        print(f"Initializing {self.name}")
        
        # Cool down first
        await self.cool_down()
        
        # Take calibration frames
        await asyncio.gather(
            self.take_dark(1.0),
            self.take_flat()
        )
        
        self.initialized = True
        return True
    
    async def acquire(self, integration_s: float) -> np.ndarray:
        """Acquire CCD frame."""
        if not self.initialized:
            raise RuntimeError("CCD not initialized")
        
        print(f"Exposing for {integration_s}s")
        await asyncio.sleep(integration_s)
        
        # Simulate data with Poisson noise
        signal = np.random.poisson(1000 * integration_s, (100, 100))
        
        # Apply calibration if available
        if self.dark_frame is not None:
            signal = signal - self.dark_frame
        if self.flat_frame is not None:
            signal = signal / self.flat_frame
        
        # Convert to flux in erg/s/cm²
        # Assume 550nm, QE=0.9
        h = 6.626e-27  # erg·s
        c = 2.998e10   # cm/s
        wavelength = 5.5e-5  # cm
        
        photon_energy = h * c / wavelength
        area_per_pixel = self.pixel_size**2
        
        flux_cgs = signal * photon_energy / (integration_s * area_per_pixel)
        
        return flux_cgs
    
    def get_sensitivity_cgs(self) -> float:
        """3-sigma sensitivity in erg/s/cm²."""
        # Depends on temperature
        read_noise = 5 if self.temperature_k < 200 else 20
        
        # Convert to flux
        h = 6.626e-27
        c = 2.998e10
        photon_energy = h * c / 5.5e-5
        
        return 3 * read_noise * photon_energy / self.pixel_size**2

# Test the system
async def test_observatory():
    ccd = CCDCamera("Main CCD", pixel_size_cm=1.5e-3)
    
    # Initialize (cooling + calibration)
    await ccd.initialize()
    
    # Take science frame
    data = await ccd.acquire(30.0)
    
    print(f"Mean flux: {data.mean():.2e} erg/s/cm²")
    print(f"Sensitivity: {ccd.get_sensitivity_cgs():.2e} erg/s/cm²")

await test_observatory()

"""
Part B: Complete observatory with async control (15 minutes)
Multiple instruments operating in parallel
"""

@dataclass
class ObservationRequest:
    """Request for observation with units."""
    target: str
    ra_deg: float
    dec_deg: float
    exposures_s: list
    filters: list
    priority: int = 5

class Observatory:
    """Complete observatory with multiple instruments."""
    
    def __init__(self, name: str):
        self.name = name
        self.instruments = {}
        self.queue = []
        self.completed = []
    
    def add_instrument(self, name: str, instrument: ObservatoryInstrument):
        """Register instrument."""
        self.instruments[name] = instrument
    
    async def initialize_all(self):
        """Initialize all instruments in parallel."""
        print(f"Initializing {self.name} observatory")
        
        tasks = [
            inst.initialize() 
            for inst in self.instruments.values()
        ]
        
        results = await asyncio.gather(*tasks)
        
        if all(results):
            print("All instruments ready")
        else:
            print("Some instruments failed initialization")
    
    async def observe(self, request: ObservationRequest):
        """Execute observation request."""
        print(f"\nObserving {request.target}")
        
        results = {}
        
        for exp_time, filter_name in zip(request.exposures_s, 
                                         request.filters):
            print(f"  {filter_name}: {exp_time}s")
            
            # Use primary CCD
            ccd = self.instruments.get('ccd')
            if ccd:
                data = await ccd.acquire(exp_time)
                results[filter_name] = {
                    'data': data,
                    'mean_flux': data.mean(),
                    'peak_flux': data.max()
                }
        
        self.completed.append((request.target, results))
        return results
    
    async def observe_queue(self):
        """Process observation queue."""
        while self.queue:
            request = self.queue.pop(0)
            await self.observe(request)
    
    def add_request(self, request: ObservationRequest):
        """Add to queue sorted by priority."""
        self.queue.append(request)
        self.queue.sort(key=lambda r: r.priority, reverse=True)

# Create and run observatory
async def full_observatory_demo():
    # Create observatory
    obs = Observatory("Mauna Kea")
    
    # Add instruments
    obs.add_instrument('ccd', CCDCamera("Primary CCD"))
    
    # Initialize everything
    await obs.initialize_all()
    
    # Queue observations
    obs.add_request(ObservationRequest(
        target="M31",
        ra_deg=10.68,
        dec_deg=41.27,
        exposures_s=[30, 60, 60],
        filters=['B', 'V', 'R'],
        priority=10
    ))
    
    obs.add_request(ObservationRequest(
        target="Calibration Star",
        ra_deg=0.0,
        dec_deg=0.0,
        exposures_s=[5, 5],
        filters=['V', 'R'],
        priority=15  # Higher priority
    ))
    
    # Process queue (high priority first)
    await obs.observe_queue()
    
    # Summary
    print(f"\nCompleted {len(obs.completed)} observations")
    for target, results in obs.completed:
        print(f"  {target}: {list(results.keys())}")

await full_observatory_demo()

"""
Implement Factory, Strategy, and Observer patterns
for a data reduction pipeline
"""

from abc import ABC, abstractmethod
from typing import List, Callable
import numpy as np

# Strategy pattern for reduction algorithms
class ReductionStrategy(ABC):
    """Abstract strategy for data reduction."""
    
    @abstractmethod
    def reduce(self, data: np.ndarray) -> np.ndarray:
        """Reduce data using specific algorithm."""
        pass

class MedianCombine(ReductionStrategy):
    """Median combination strategy."""
    
    def reduce(self, data: np.ndarray) -> np.ndarray:
        """Median combine along first axis."""
        return np.median(data, axis=0)

class SigmaClipping(ReductionStrategy):
    """Sigma-clipped mean strategy."""
    
    def __init__(self, sigma: float = 3.0):
        self.sigma = sigma
    
    def reduce(self, data: np.ndarray) -> np.ndarray:
        """Sigma-clipped mean combination."""
        mean = np.mean(data, axis=0)
        std = np.std(data, axis=0)
        
        # Mask outliers
        mask = np.abs(data - mean) < self.sigma * std
        
        # Recalculate with mask
        masked_data = np.where(mask, data, np.nan)
        return np.nanmean(masked_data, axis=0)

# Observer pattern for pipeline monitoring
class PipelineObserver(ABC):
    """Abstract observer for pipeline events."""
    
    @abstractmethod
    def update(self, event: str, data: dict):
        """Called when pipeline state changes."""
        pass

class LogObserver(PipelineObserver):
    """Logs pipeline events."""
    
    def update(self, event: str, data: dict):
        """Log the event."""
        print(f"[LOG] {event}: {data}")

class QualityObserver(PipelineObserver):
    """Monitors data quality."""
    
    def __init__(self):
        self.quality_metrics = []
    
    def update(self, event: str, data: dict):
        """Check quality metrics."""
        if event == "reduction_complete":
            result = data.get('result')
            if result is not None:
                snr = np.mean(result) / np.std(result)
                self.quality_metrics.append(snr)
                
                if snr < 10:
                    print(f"[QUALITY WARNING] Low SNR: {snr:.1f}")

# Factory pattern for pipeline creation
class PipelineFactory:
    """Factory for creating appropriate pipelines."""
    
    @staticmethod
    def create_pipeline(obs_type: str) -> 'ReductionPipeline':
        """Create pipeline based on observation type."""
        
        if obs_type == "imaging":
            pipeline = ReductionPipeline()
            pipeline.set_strategy(MedianCombine())
            pipeline.attach(LogObserver())
            return pipeline
        
        elif obs_type == "spectroscopy":
            pipeline = ReductionPipeline()
            pipeline.set_strategy(SigmaClipping(sigma=2.5))
            pipeline.attach(LogObserver())
            pipeline.attach(QualityObserver())
            return pipeline
        
        elif obs_type == "photometry":
            pipeline = ReductionPipeline()
            pipeline.set_strategy(SigmaClipping(sigma=3.0))
            pipeline.attach(QualityObserver())
            return pipeline
        
        else:
            raise ValueError(f"Unknown observation type: {obs_type}")

# Main pipeline using all patterns
class ReductionPipeline:
    """Data reduction pipeline with patterns."""
    
    def __init__(self):
        self.strategy: Optional[ReductionStrategy] = None
        self.observers: List[PipelineObserver] = []
    
    def set_strategy(self, strategy: ReductionStrategy):
        """Set reduction strategy."""
        self.strategy = strategy
        self.notify("strategy_changed", 
                   {"strategy": strategy.__class__.__name__})
    
    def attach(self, observer: PipelineObserver):
        """Attach observer."""
        self.observers.append(observer)
    
    def notify(self, event: str, data: dict):
        """Notify all observers."""
        for observer in self.observers:
            observer.update(event, data)
    
    def process(self, frames: np.ndarray) -> np.ndarray:
        """Process frames using current strategy."""
        if not self.strategy:
            raise RuntimeError("No reduction strategy set")
        
        self.notify("reduction_started", {"n_frames": len(frames)})
        
        # Apply strategy
        result = self.strategy.reduce(frames)
        
        self.notify("reduction_complete", {
            "result_shape": result.shape,
            "result": result
        })
        
        return result

# Demonstrate all patterns together
print("=== Design Patterns Demo ===\n")

# Create different pipelines using factory
imaging_pipeline = PipelineFactory.create_pipeline("imaging")
spectro_pipeline = PipelineFactory.create_pipeline("spectroscopy")

# Generate mock data
frames = np.random.randn(10, 100, 100) * 100 + 1000

# Process with different pipelines
print("Imaging Pipeline:")
imaging_result = imaging_pipeline.process(frames)

print("\nSpectroscopy Pipeline:")
spectro_result = spectro_pipeline.process(frames)

# Change strategy at runtime
print("\nChanging strategy at runtime:")
imaging_pipeline.set_strategy(SigmaClipping(sigma=2.0))
new_result = imaging_pipeline.process(frames)

print("\n✅ Patterns demonstrated:")
print("  - Factory: Created specialized pipelines")
print("  - Strategy: Swappable reduction algorithms")
print("  - Observer: Monitoring and logging")

"""
Debug This! Can you find and fix the performance issue?
"""

# Inefficient implementation
class SlowGalaxy:
    def __init__(self, name, redshift):
        self.name = name
        self.redshift = redshift
    
    def distance_mpc(self):
        """Calculate distance in Mpc (SLOW!)."""
        # Hubble constant
        H0 = 70  # km/s/Mpc
        c = 3e5  # km/s
        
        # This recalculates every time!
        return c * self.redshift / H0
    
    def luminosity_distance(self):
        """Luminosity distance (calls distance repeatedly)."""
        d = self.distance_mpc()
        
        # Cosmological correction (simplified)
        for i in range(1000):  # Unnecessary loop!
            d = d * (1 + self.redshift * 0.001)
        
        return d

# Test slow version
import time

galaxies = [SlowGalaxy(f"Galaxy{i}", 0.01 * i) 
           for i in range(100)]

start = time.perf_counter()
for g in galaxies:
    for _ in range(10):  # Multiple calls
        d = g.luminosity_distance()
slow_time = time.perf_counter() - start

print(f"Slow version: {slow_time*1000:.1f} ms")

# SOLUTION: Optimized implementation
from functools import lru_cache

class FastGalaxy:
    __slots__ = ['name', 'redshift', '_distance_cache']
    
    def __init__(self, name, redshift):
        self.name = name
        self.redshift = redshift
        self._distance_cache = None
    
    @property
    def distance_mpc(self):
        """Cached distance calculation."""
        if self._distance_cache is None:
            H0 = 70
            c = 3e5
            self._distance_cache = c * self.redshift / H0
        return self._distance_cache
    
    @lru_cache(maxsize=1)
    def luminosity_distance(self):
        """Cached luminosity distance."""
        d = self.distance_mpc
        
        # Vectorized calculation instead of loop
        correction = (1 + self.redshift) ** 1.0
        
        return d * correction

# Test fast version
galaxies_fast = [FastGalaxy(f"Galaxy{i}", 0.01 * i) 
                for i in range(100)]

start = time.perf_counter()
for g in galaxies_fast:
    for _ in range(10):
        d = g.luminosity_distance()
fast_time = time.perf_counter() - start

print(f"Fast version: {fast_time*1000:.1f} ms")
print(f"Speedup: {slow_time/fast_time:.1f}x")

print("\n🐛 Performance issues found and fixed:")
print("  1. Recalculating distance every call")
print("  2. Unnecessary loop in luminosity_distance")
print("  3. No caching of expensive calculations")
print("  4. Using dict instead of __slots__")