Chapter 4: Python Practice: Electromagnetic Data Analysis Workflow

In this chapter, you will master a practical Python workflow for handling data from real measurement instruments (VSM, SQUID, Hall measurement systems). We will build an integrated pipeline covering the loading of CSV and binary data, outlier removal, background subtraction, advanced fitting techniques (lmfit, scipy.optimize), error propagation, anomaly detection using machine learning, creation of publication-quality figures, and automated report generation.

4.1 Data Loading and Cleaning

4.1.1 Multi-Format Data Loader

Real measurement instruments output data in a variety of formats:

Code Example 4-1: Multi-Format Data Loader Class

import numpy as np
import pandas as pd
import struct
from pathlib import Path
from typing import Union, Tuple

class UniversalDataLoader:
    """
    Multi-format data loader (supports CSV, DAT, and binary)

    Attributes
    ----------
    data : pd.DataFrame
        The loaded data
    metadata : dict
        Metadata (header information, measurement conditions, etc.)
    """

    def __init__(self):
        self.data = None
        self.metadata = {}

    def load(self, filepath: str, format: str = 'auto') -> pd.DataFrame:
        """
        Load a data file

        Parameters
        ----------
        filepath : str
            File path
        format : str
            'auto' (automatic detection), 'csv', 'dat', 'binary'

        Returns
        -------
        data : pd.DataFrame
            The loaded data
        """
        path = Path(filepath)

        if not path.exists():
            raise FileNotFoundError(f"File not found: {filepath}")

        # Automatic format detection
        if format == 'auto':
            format = self._detect_format(path)

        # Load according to format
        if format == 'csv':
            self.data = self._load_csv(path)
        elif format == 'dat':
            self.data = self._load_dat(path)
        elif format == 'binary':
            self.data = self._load_binary(path)
        else:
            raise ValueError(f"Unsupported format: {format}")

        print(f"Loaded {len(self.data)} rows from {path.name} (format: {format})")
        return self.data

    def _detect_format(self, path: Path) -> str:
        """Infer the format from the file extension"""
        ext = path.suffix.lower()
        if ext in ['.csv', '.txt']:
            return 'csv'
        elif ext == '.dat':
            return 'dat'
        elif ext == '.bin':
            return 'binary'
        else:
            # Inspect the content to decide
            with open(path, 'rb') as f:
                header = f.read(100)
                if b'\x00' in header:
                    return 'binary'
                else:
                    return 'csv'

    def _load_csv(self, path: Path) -> pd.DataFrame:
        """Load CSV format"""
        # Detect number of header lines
        header_lines = 0
        with open(path, 'r') as f:
            for i, line in enumerate(f):
                if line.strip().startswith('#'):
                    header_lines += 1
                    self.metadata[f'comment_{i}'] = line.strip()
                else:
                    break

        # Load the data
        try:
            data = pd.read_csv(path, skiprows=header_lines, sep=None, engine='python')
        except Exception as e:
            # For delimiters other than comma
            data = pd.read_csv(path, skiprows=header_lines, delimiter=r'\s+')

        return data

    def _load_dat(self, path: Path) -> pd.DataFrame:
        """Load DAT format (assuming Quantum Design VSM)"""
        # Separate metadata from data
        metadata_lines = []
        data_lines = []

        with open(path, 'r') as f:
            in_data_section = False
            for line in f:
                line = line.strip()
                if line.startswith('[Data]'):
                    in_data_section = True
                    continue

                if not in_data_section:
                    if line and not line.startswith('#'):
                        # Metadata line
                        if '=' in line:
                            key, value = line.split('=', 1)
                            self.metadata[key.strip()] = value.strip()
                else:
                    if line and not line.startswith('#'):
                        data_lines.append(line)

        # Convert the data into a DataFrame
        if data_lines:
            # Column names (first line)
            header = data_lines[0].split('\t')

            # Data rows
            data_values = []
            for line in data_lines[1:]:
                values = line.split('\t')
                data_values.append([float(v) if v else np.nan for v in values])

            data = pd.DataFrame(data_values, columns=header)
        else:
            data = pd.DataFrame()

        return data

    def _load_binary(self, path: Path) -> pd.DataFrame:
        """
        Load binary format (custom format)

        Assumptions:
        - Header: the first 100 bytes (metadata)
        - Data: float64 array, 3 columns (H, M, T)
        """
        with open(path, 'rb') as f:
            # Read the header
            header_bytes = f.read(100)
            # Parse metadata (omitted here)

            # Read the data
            remaining = f.read()
            n_cols = 3
            n_rows = len(remaining) // (8 * n_cols)  # float64 = 8 bytes

            data_array = np.frombuffer(remaining, dtype=np.float64)
            data_array = data_array.reshape((n_rows, n_cols))

            data = pd.DataFrame(data_array, columns=['H', 'M', 'T'])

        return data

    def clean_data(self, remove_nan=True, remove_inf=True, outlier_method='iqr', threshold=3.0):
        """
        Clean the data

        Parameters
        ----------
        remove_nan : bool
            Whether to remove NaN values
        remove_inf : bool
            Whether to remove Inf values
        outlier_method : str
            Outlier removal method ('iqr', 'zscore', 'none')
        threshold : float
            Threshold for outlier determination
        """
        if self.data is None:
            raise ValueError("No data loaded. Call load() first.")

        original_len = len(self.data)

        # Remove NaN/Inf
        if remove_nan:
            self.data = self.data.dropna()
        if remove_inf:
            self.data = self.data.replace([np.inf, -np.inf], np.nan).dropna()

        # Remove outliers
        if outlier_method == 'iqr':
            self.data = self._remove_outliers_iqr(self.data, threshold)
        elif outlier_method == 'zscore':
            self.data = self._remove_outliers_zscore(self.data, threshold)

        print(f"Data cleaning: {original_len} → {len(self.data)} rows (removed {original_len - len(self.data)})")

        return self.data

    def _remove_outliers_iqr(self, data: pd.DataFrame, threshold: float) -> pd.DataFrame:
        """Remove outliers using the IQR method"""
        Q1 = data.quantile(0.25)
        Q3 = data.quantile(0.75)
        IQR = Q3 - Q1

        lower_bound = Q1 - threshold * IQR
        upper_bound = Q3 + threshold * IQR

        mask = ((data >= lower_bound) & (data <= upper_bound)).all(axis=1)
        return data[mask]

    def _remove_outliers_zscore(self, data: pd.DataFrame, threshold: float) -> pd.DataFrame:
        """Remove outliers using the Z-score method"""
        z_scores = np.abs((data - data.mean()) / data.std())
        mask = (z_scores < threshold).all(axis=1)
        return data[mask]

# Example usage
loader = UniversalDataLoader()

# Load CSV
data = loader.load('vsm_measurement.csv', format='auto')
print(data.head())
print(f"\nMetadata: {loader.metadata}")

# Clean the data
data_clean = loader.clean_data(outlier_method='iqr', threshold=3.0)
print(f"\nCleaned data shape: {data_clean.shape}")

4.2 Integrated Analysis Pipeline

4.2.1 Integrated Four-Point Probe + Hall + M-H Analysis

In experiments, multiple measurements are often performed on the same sample. By analyzing these together, we can gain a complete understanding of the material's electrical and magnetic properties.

Code Example 4-2: Integrated Analysis Pipeline

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import fsolve, curve_fit
from lmfit import Model

class IntegratedAnalysisPipeline:
    """
    Integrated four-point probe + Hall + M-H analysis pipeline

    Attributes
    ----------
    thickness : float
        Sample thickness [m]
    mass : float
        Sample mass [g]
    results : dict
        Analysis results (σ, n, μ, M_s, H_c, M_r, K)
    """

    def __init__(self, thickness: float, mass: float = None):
        self.t = thickness
        self.mass = mass
        self.results = {}
        self.data = {}

    def load_four_probe_data(self, R_AB_CD: float, R_BC_DA: float):
        """Load van der Pauw four-point probe data"""
        self.data['R_AB_CD'] = R_AB_CD
        self.data['R_BC_DA'] = R_BC_DA

    def load_hall_data(self, I: float, B: float, V_pos: float, V_neg: float):
        """Load Hall data"""
        self.data['I'] = I
        self.data['B'] = B
        self.data['V_hall_pos'] = V_pos
        self.data['V_hall_neg'] = V_neg

    def load_mh_data(self, H: np.ndarray, M: np.ndarray):
        """Load M-H data"""
        self.data['H'] = H
        self.data['M'] = M

    def analyze_electrical_properties(self):
        """Analyze electrical properties (four-point probe + Hall)"""
        # van der Pauw sheet resistance
        def vdp_eq(Rs, R1, R2):
            return np.exp(-np.pi * R1 / Rs) + np.exp(-np.pi * R2 / Rs) - 1

        R1 = self.data['R_AB_CD']
        R2 = self.data['R_BC_DA']
        R_initial = (R1 + R2) / 2 * np.pi / np.log(2)
        R_s = fsolve(vdp_eq, R_initial, args=(R1, R2))[0]

        # Electrical conductivity
        sigma = 1 / (R_s * self.t)

        # Hall coefficient
        V_H = 0.5 * (self.data['V_hall_pos'] - self.data['V_hall_neg'])
        R_H = V_H * self.t / (self.data['I'] * self.data['B'])

        # Carrier density
        e = 1.60218e-19
        n = 1 / (np.abs(R_H) * e)
        carrier_type = 'electron' if R_H < 0 else 'hole'

        # Mobility
        mu = sigma * np.abs(R_H)

        # Save results
        self.results.update({
            'R_s': R_s,
            'sigma': sigma,
            'rho': 1 / sigma,
            'R_H': R_H,
            'n': n,
            'carrier_type': carrier_type,
            'mu': mu
        })

        return self.results

    def analyze_magnetic_properties(self):
        """Analyze magnetic properties (M-H curve)"""
        H = self.data['H']
        M = self.data['M']

        # Background subtraction
        H_high = H[H > 0.8 * np.max(H)]
        M_high = M[H > 0.8 * np.max(H)]
        slope = np.polyfit(H_high, M_high, 1)[0]
        M_corrected = M - slope * H

        # Saturation magnetization
        M_s = np.mean(M_corrected[H > 0.8 * np.max(H)])

        # Coercivity
        from scipy.interpolate import interp1d
        from scipy.optimize import brentq
        interp_func = interp1d(H, M_corrected, kind='linear')
        H_c = brentq(interp_func, np.min(H[H < 0]), np.max(H[H > 0]))

        # Remanence
        M_r = interp_func(0)

        # Squareness ratio
        S = M_r / M_s if M_s != 0 else 0

        # Magnetic anisotropy constant
        K = H_c * M_s / 2  # CGS units [erg/cm^3]

        # Save results
        self.results.update({
            'M_s': M_s,
            'H_c': H_c,
            'M_r': M_r,
            'S': S,
            'K': K
        })

        return self.results

    def generate_report(self):
        """Generate an integrated report"""
        print("=" * 80)
        print("INTEGRATED ELECTRICAL & MAGNETIC PROPERTIES ANALYSIS")
        print("=" * 80)

        print("\n[Electrical Properties]")
        print(f"  Sheet Resistance R_s = {self.results['R_s']:.2f} Ω/sq")
        print(f"  Conductivity σ = {self.results['sigma']:.2e} S/m")
        print(f"  Resistivity ρ = {self.results['rho']:.2e} Ω·m = {self.results['rho'] * 1e8:.2f} μΩ·cm")
        print(f"  Carrier Type: {self.results['carrier_type']}")
        print(f"  Carrier Density n = {self.results['n']:.2e} m⁻³ = {self.results['n'] / 1e6:.2e} cm⁻³")
        print(f"  Mobility μ = {self.results['mu']:.2e} m²/(V·s) = {self.results['mu'] * 1e4:.1f} cm²/(V·s)")

        print("\n[Magnetic Properties]")
        print(f"  Saturation Magnetization M_s = {self.results['M_s']:.2f} emu/g")
        print(f"  Coercivity H_c = {self.results['H_c']:.2f} Oe = {self.results['H_c'] / 79.5775:.2f} kA/m")
        print(f"  Remanence M_r = {self.results['M_r']:.2f} emu/g")
        print(f"  Squareness S = {self.results['S']:.3f}")
        print(f"  Anisotropy Constant K = {self.results['K']:.2e} erg/cm³ = {self.results['K'] * 1e3:.2e} J/m³")

        print("\n[Material Classification]")
        if self.results['H_c'] < 100:
            print("  → Soft magnetic material (transformers, inductors)")
        elif self.results['H_c'] > 1000:
            print("  → Hard magnetic material (permanent magnets)")
        else:
            print("  → Medium coercivity (recording media)")

        if self.results['mu'] * 1e4 > 1000:
            print("  → High mobility material (high-performance electronics)")
        else:
            print("  → Moderate mobility (standard electronics)")

        print("=" * 80)

# Example usage
pipeline = IntegratedAnalysisPipeline(thickness=200e-9, mass=0.005)

# Load data
pipeline.load_four_probe_data(R_AB_CD=1000, R_BC_DA=950)
pipeline.load_hall_data(I=100e-6, B=0.5, V_pos=-5.0e-3, V_neg=4.8e-3)

H_data = np.linspace(-5000, 5000, 200)
M_data = 50 * np.tanh(H_data / 1000) + 0.5 * np.random.randn(200)
pipeline.load_mh_data(H=H_data, M=M_data)

# Run the analysis
pipeline.analyze_electrical_properties()
pipeline.analyze_magnetic_properties()

# Generate report
pipeline.generate_report()

4.3 Advanced Fitting Techniques

4.3.1 Constrained Fitting with lmfit

lmfit is a wrapper around scipy.optimize that makes it easy to set parameter constraints, estimate errors, and compute correlation matrices.

Code Example 4-3: Complex Fitting with lmfit

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model, Parameters

def temperature_dependent_hall(T, n0, Ea, mu0, alpha):
    """
    Temperature-dependent Hall model

    n(T) = n0 * exp(Ea / (k_B * T))  # Carrier density
    μ(T) = μ0 * (T / 300)^(-α)      # Mobility
    R_H(T) = 1 / (n(T) * e)

    Parameters
    ----------
    T : array-like
        Temperature [K]
    n0 : float
        Reference carrier density [m^-3]
    Ea : float
        Activation energy [eV]
    mu0 : float
        Reference mobility (at 300 K) [m^2/(V·s)]
    alpha : float
        Temperature exponent of mobility

    Returns
    -------
    R_H : array-like
        Hall coefficient [m^3/C]
    """
    k_B = 8.617e-5  # [eV/K]
    e = 1.60218e-19  # [C]

    n = n0 * np.exp(-Ea / (k_B * T))
    R_H = 1 / (n * e)

    return R_H

# Generate simulated data
T_range = np.linspace(200, 400, 25)
n0_true = 1e21  # [m^-3]
Ea_true = 0.3  # [eV]
mu0_true = 0.05  # [m^2/(V·s)]
alpha_true = 1.5

R_H_data = temperature_dependent_hall(T_range, n0_true, Ea_true, mu0_true, alpha_true)
R_H_data_noise = R_H_data * (1 + 0.08 * np.random.randn(len(T_range)))

# Define the lmfit model
model = Model(temperature_dependent_hall)

# Set parameters (initial values, ranges, constraints)
params = model.make_params(
    n0 = {'value': 5e20, 'min': 1e19, 'max': 1e23},
    Ea = {'value': 0.4, 'min': 0.1, 'max': 1.0},
    mu0 = {'value': 0.1, 'min': 0.01, 'max': 1.0},
    alpha = {'value': 1.0, 'min': 0.5, 'max': 3.0}
)

# Run the fit
result = model.fit(R_H_data_noise, params, T=T_range)

# Display results
print("=" * 80)
print("ADVANCED FITTING WITH LMFIT")
print("=" * 80)
print(result.fit_report())

# Parameter correlation matrix
print("\nParameter Correlation Matrix:")
print(result.params.correlation)

# Plotting
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Top left: R_H vs T
axes[0, 0].scatter(T_range, R_H_data_noise, s=80, alpha=0.7, edgecolors='black', linewidths=1.5, label='Data (with noise)', color='#f093fb')
axes[0, 0].plot(T_range, result.best_fit, linewidth=2.5, label='Fit', color='#f5576c')
axes[0, 0].plot(T_range, R_H_data, linewidth=2, linestyle='--', label='True (no noise)', color='green')
axes[0, 0].set_xlabel('Temperature T [K]', fontsize=12)
axes[0, 0].set_ylabel('Hall Coefficient R$_H$ [m$^3$/C]', fontsize=12)
axes[0, 0].set_title('Hall Coefficient vs Temperature', fontsize=13, fontweight='bold')
axes[0, 0].legend(fontsize=10)
axes[0, 0].grid(alpha=0.3)
axes[0, 0].set_yscale('log')

# Top right: residuals
residuals = R_H_data_noise - result.best_fit
axes[0, 1].scatter(T_range, residuals, s=80, alpha=0.7, edgecolors='black', linewidths=1.5, color='#99ccff')
axes[0, 1].axhline(0, color='black', linestyle='--', linewidth=1.5)
axes[0, 1].set_xlabel('Temperature T [K]', fontsize=12)
axes[0, 1].set_ylabel('Residuals [m$^3$/C]', fontsize=12)
axes[0, 1].set_title('Fit Residuals', fontsize=13, fontweight='bold')
axes[0, 1].grid(alpha=0.3)

# Bottom left: carrier density
k_B = 8.617e-5
e = 1.60218e-19
n_fit = result.params['n0'].value * np.exp(-result.params['Ea'].value / (k_B * T_range))
axes[1, 0].semilogy(T_range, n_fit / 1e6, linewidth=2.5, color='#ffa500', label='Fitted n(T)')
axes[1, 0].set_xlabel('Temperature T [K]', fontsize=12)
axes[1, 0].set_ylabel('Carrier Density n [cm$^{-3}$]', fontsize=12)
axes[1, 0].set_title('Carrier Density (from fit)', fontsize=13, fontweight='bold')
axes[1, 0].legend(fontsize=11)
axes[1, 0].grid(alpha=0.3)

# Bottom right: mobility
mu_fit = result.params['mu0'].value * (T_range / 300)**(-result.params['alpha'].value)
axes[1, 1].loglog(T_range, mu_fit * 1e4, linewidth=2.5, color='#99ff99', label='Fitted μ(T)')
axes[1, 1].set_xlabel('Temperature T [K]', fontsize=12)
axes[1, 1].set_ylabel('Mobility μ [cm$^2$/(V·s)]', fontsize=12)
axes[1, 1].set_title('Mobility (from fit)', fontsize=13, fontweight='bold')
axes[1, 1].legend(fontsize=11)
axes[1, 1].grid(alpha=0.3, which='both')

plt.tight_layout()
plt.show()

# Physical interpretation
print("\n" + "=" * 80)
print("PHYSICAL INTERPRETATION")
print("=" * 80)
print(f"Activation Energy E_a = {result.params['Ea'].value:.3f} ± {result.params['Ea'].stderr:.3f} eV")
print(f"  → Band gap or dopant ionization energy")
print(f"Mobility Exponent α = {result.params['alpha'].value:.2f} ± {result.params['alpha'].stderr:.2f}")
print(f"  → Scattering mechanism: α ≈ 1.5 suggests acoustic phonon scattering")
print(f"Carrier Density (300 K) = {result.params['n0'].value * np.exp(-result.params['Ea'].value / (k_B * 300)):.2e} m⁻³")
print(f"Mobility (300 K) = {result.params['mu0'].value * 1e4:.1f} cm²/(V·s)")

4.4 Error Propagation and Uncertainty Evaluation

4.4.1 Automatic Error Propagation

The uncertainty of measured values propagates into derived quantities (mobility, anisotropy constant, etc.). The uncertainties package lets you compute error propagation automatically.

Code Example 4-4: Automatic Error Propagation

from uncertainties import ufloat, umath
import numpy as np

def propagate_uncertainties_example():
    """
    A worked example of error propagation

    Compute mobility from Hall measurements and automatically propagate the uncertainty
    """
    # Measured values (value ± uncertainty)
    R_H = ufloat(-2.5e-3, 0.1e-3)  # Hall coefficient [m^3/C]
    sigma = ufloat(1e4, 100)  # Electrical conductivity [S/m]
    I = ufloat(100e-6, 1e-6)  # Current [A]
    B = ufloat(0.5, 0.01)  # Magnetic field [T]
    t = ufloat(200e-9, 5e-9)  # Thickness [m]

    # Mobility: μ = σ * |R_H|
    mu = sigma * umath.fabs(R_H)

    # Carrier density: n = 1 / (e * |R_H|)
    e = 1.60218e-19  # Constant (no uncertainty)
    n = 1 / (e * umath.fabs(R_H))

    # Magnetic anisotropy constant: K = H_c * M_s / 2
    H_c = ufloat(500, 20)  # [Oe]
    M_s = ufloat(50, 2)  # [emu/g]
    K = H_c * M_s / 2

    print("=" * 80)
    print("AUTOMATIC ERROR PROPAGATION")
    print("=" * 80)

    print("\n[Input Measurements]")
    print(f"  Hall Coefficient R_H = {R_H} m³/C")
    print(f"  Conductivity σ = {sigma} S/m")
    print(f"  Current I = {I} A")
    print(f"  Magnetic Field B = {B} T")
    print(f"  Thickness t = {t} m")

    print("\n[Derived Quantities with Propagated Uncertainties]")
    print(f"  Mobility μ = {mu} m²/(V·s)")
    print(f"            = ({mu.nominal_value * 1e4:.1f} ± {mu.std_dev * 1e4:.1f}) cm²/(V·s)")
    print(f"  Carrier Density n = {n} m⁻³")
    print(f"                    = ({n.nominal_value / 1e6:.2e} ± {n.std_dev / 1e6:.2e}) cm⁻³")
    print(f"  Anisotropy Constant K = {K} erg/cm³")

    print("\n[Relative Uncertainties]")
    print(f"  μ: {mu.std_dev / mu.nominal_value * 100:.2f}%")
    print(f"  n: {n.std_dev / n.nominal_value * 100:.2f}%")
    print(f"  K: {K.std_dev / K.nominal_value * 100:.2f}%")

    # Uncertainty budget analysis (verified with manual calculation)
    print("\n[Uncertainty Budget for μ = σ * |R_H|]")
    rel_sigma = (100 / 1e4)**2
    rel_R_H = (0.1e-3 / 2.5e-3)**2
    rel_mu_manual = np.sqrt(rel_sigma + rel_R_H)
    print(f"  Contribution from σ: {np.sqrt(rel_sigma) * 100:.2f}%")
    print(f"  Contribution from R_H: {np.sqrt(rel_R_H) * 100:.2f}%")
    print(f"  Total (manual): {rel_mu_manual * 100:.2f}%")
    print(f"  Total (automatic): {mu.std_dev / mu.nominal_value * 100:.2f}%")

propagate_uncertainties_example()

4.5 Creating Publication-Quality Figures

4.5.1 matplotlib Best Practices

Key points for creating figures suitable for publication in journals:

• Font size: axis labels 12-14pt, titles 14-16pt, legend 10-12pt
• Line width: data lines 2-3pt, grid lines 0.5-1pt
• Marker size: 80-150 (scatter)
• Color: colorblind-friendly palettes (e.g., seaborn, viridis)
• Resolution: DPI 300 or higher (for print)
• Format: PDF (vector) or PNG (raster)

Code Example 4-5: Generating Publication-Quality Figures

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.ticker import MultipleLocator, AutoMinorLocator

# Publication settings
mpl.rcParams['font.family'] = 'Arial'
mpl.rcParams['font.size'] = 12
mpl.rcParams['axes.linewidth'] = 1.5
mpl.rcParams['xtick.major.width'] = 1.5
mpl.rcParams['ytick.major.width'] = 1.5
mpl.rcParams['xtick.minor.width'] = 1.0
mpl.rcParams['ytick.minor.width'] = 1.0
mpl.rcParams['xtick.major.size'] = 6
mpl.rcParams['ytick.major.size'] = 6
mpl.rcParams['xtick.minor.size'] = 3
mpl.rcParams['ytick.minor.size'] = 3

def create_publication_figure():
    """
    Create a publication-quality figure

    Example: results of a temperature-dependent Hall measurement
    """
    # Generate data
    T = np.linspace(100, 400, 20)
    n = 1e22 * np.exp(-0.3 / (8.617e-5 * T))
    mu = 0.05 * (300 / T)**1.5
    sigma = n * 1.60218e-19 * mu

    # Create figure (2x2 layout)
    fig = plt.figure(figsize=(12, 10))
    gs = fig.add_gridspec(2, 2, hspace=0.3, wspace=0.3)

    # (a) Carrier density
    ax1 = fig.add_subplot(gs[0, 0])
    ax1.semilogy(T, n / 1e6, 'o-', linewidth=2.5, markersize=8, color='#2E86AB',
                 markeredgecolor='black', markeredgewidth=1.5, label='Carrier density')
    ax1.set_xlabel('Temperature (K)', fontsize=14, fontweight='bold')
    ax1.set_ylabel('Carrier Density (cm$^{-3}$)', fontsize=14, fontweight='bold')
    ax1.set_title('(a) Temperature-Dependent Carrier Density', fontsize=14, fontweight='bold', loc='left')
    ax1.legend(fontsize=12, frameon=True, shadow=True)
    ax1.grid(True, which='both', alpha=0.3, linestyle='--')
    ax1.xaxis.set_minor_locator(AutoMinorLocator())

    # (b) Mobility
    ax2 = fig.add_subplot(gs[0, 1])
    ax2.loglog(T, mu * 1e4, 's-', linewidth=2.5, markersize=8, color='#A23B72',
               markeredgecolor='black', markeredgewidth=1.5, label='Mobility')
    ax2.set_xlabel('Temperature (K)', fontsize=14, fontweight='bold')
    ax2.set_ylabel('Mobility (cm$^2$/(V·s))', fontsize=14, fontweight='bold')
    ax2.set_title('(b) Temperature-Dependent Mobility', fontsize=14, fontweight='bold', loc='left')
    ax2.legend(fontsize=12, frameon=True, shadow=True)
    ax2.grid(True, which='both', alpha=0.3, linestyle='--')

    # (c) Electrical conductivity
    ax3 = fig.add_subplot(gs[1, 0])
    ax3.semilogy(T, sigma, '^-', linewidth=2.5, markersize=8, color='#F18F01',
                 markeredgecolor='black', markeredgewidth=1.5, label='Conductivity')
    ax3.set_xlabel('Temperature (K)', fontsize=14, fontweight='bold')
    ax3.set_ylabel('Conductivity (S/m)', fontsize=14, fontweight='bold')
    ax3.set_title('(c) Electrical Conductivity', fontsize=14, fontweight='bold', loc='left')
    ax3.legend(fontsize=12, frameon=True, shadow=True)
    ax3.grid(True, which='both', alpha=0.3, linestyle='--')
    ax3.xaxis.set_minor_locator(AutoMinorLocator())

    # (d) Arrhenius plot
    ax4 = fig.add_subplot(gs[1, 1])
    ax4.semilogy(1000 / T, n / 1e6, 'D-', linewidth=2.5, markersize=8, color='#C73E1D',
                 markeredgecolor='black', markeredgewidth=1.5, label='Arrhenius plot')
    ax4.set_xlabel('1000/T (K$^{-1}$)', fontsize=14, fontweight='bold')
    ax4.set_ylabel('Carrier Density (cm$^{-3}$)', fontsize=14, fontweight='bold')
    ax4.set_title('(d) Arrhenius Plot', fontsize=14, fontweight='bold', loc='left')
    ax4.legend(fontsize=12, frameon=True, shadow=True)
    ax4.grid(True, which='both', alpha=0.3, linestyle='--')

    # Overall adjustment
    for ax in [ax1, ax2, ax3, ax4]:
        ax.tick_params(axis='both', which='major', labelsize=12, direction='in', top=True, right=True)
        ax.tick_params(axis='both', which='minor', direction='in', top=True, right=True)

    # Save (high-resolution PDF + PNG)
    plt.savefig('hall_measurement_publication.pdf', dpi=300, bbox_inches='tight', format='pdf')
    plt.savefig('hall_measurement_publication.png', dpi=300, bbox_inches='tight', format='png')

    print("Figures saved:")
    print("  - hall_measurement_publication.pdf (vector, for publication)")
    print("  - hall_measurement_publication.png (raster, for preview)")

    plt.show()

create_publication_figure()

4.6 Automated Report Generation

4.6.1 Automatic PDF Reports

Using matplotlib and matplotlib.backends.backend_pdf, we can automatically generate multi-page PDF reports.

Code Example 4-6: Automated PDF Report Generation

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
from datetime import datetime

class AutoReportGenerator:
    """
    Automated PDF report generator class

    Attributes
    ----------
    filename : str
        Output PDF file name
    metadata : dict
        Report metadata
    """

    def __init__(self, filename='analysis_report.pdf'):
        self.filename = filename
        self.metadata = {
            'Title': 'Electrical & Magnetic Properties Analysis Report',
            'Author': 'MS Terakoya Analysis Pipeline',
            'Subject': 'Automated Data Analysis',
            'Keywords': 'Hall effect, Magnetometry, Python',
            'CreationDate': datetime.now()
        }

    def generate_report(self, results: dict):
        """
        Generate the complete report

        Parameters
        ----------
        results : dict
            Analysis results (electrical and magnetic properties)
        """
        with PdfPages(self.filename) as pdf:
            # Page 1: Summary
            self._add_summary_page(pdf, results)

            # Page 2: Electrical properties plots
            self._add_electrical_plots(pdf, results)

            # Page 3: Magnetic properties plots
            self._add_magnetic_plots(pdf, results)

            # Page 4: Statistics
            self._add_statistics_page(pdf, results)

            # Set metadata
            d = pdf.infodict()
            for key, value in self.metadata.items():
                d[key] = value

        print(f"Report generated: {self.filename}")

    def _add_summary_page(self, pdf, results):
        """Page 1: Summary (text only)"""
        fig = plt.figure(figsize=(8.5, 11))
        fig.text(0.5, 0.95, 'ANALYSIS SUMMARY', ha='center', fontsize=20, fontweight='bold')
        fig.text(0.5, 0.90, f'Generated: {datetime.now().strftime("%Y-%m-%d %H:%M:%S")}', ha='center', fontsize=10)

        # Electrical properties
        y_start = 0.75
        fig.text(0.1, y_start, 'ELECTRICAL PROPERTIES', fontsize=16, fontweight='bold', color='#2E86AB')
        fig.text(0.1, y_start - 0.05, f"Sheet Resistance R_s = {results['R_s']:.2f} Ω/sq", fontsize=12)
        fig.text(0.1, y_start - 0.10, f"Conductivity σ = {results['sigma']:.2e} S/m", fontsize=12)
        fig.text(0.1, y_start - 0.15, f"Carrier Type: {results['carrier_type']}", fontsize=12)
        fig.text(0.1, y_start - 0.20, f"Carrier Density n = {results['n']:.2e} m⁻³", fontsize=12)
        fig.text(0.1, y_start - 0.25, f"Mobility μ = {results['mu'] * 1e4:.1f} cm²/(V·s)", fontsize=12)

        # Magnetic properties
        y_start = 0.45
        fig.text(0.1, y_start, 'MAGNETIC PROPERTIES', fontsize=16, fontweight='bold', color='#A23B72')
        fig.text(0.1, y_start - 0.05, f"Saturation Magnetization M_s = {results['M_s']:.2f} emu/g", fontsize=12)
        fig.text(0.1, y_start - 0.10, f"Coercivity H_c = {results['H_c']:.2f} Oe", fontsize=12)
        fig.text(0.1, y_start - 0.15, f"Remanence M_r = {results['M_r']:.2f} emu/g", fontsize=12)
        fig.text(0.1, y_start - 0.20, f"Squareness S = {results['S']:.3f}", fontsize=12)
        fig.text(0.1, y_start - 0.25, f"Anisotropy Constant K = {results['K']:.2e} erg/cm³", fontsize=12)

        plt.axis('off')
        pdf.savefig(fig, bbox_inches='tight')
        plt.close(fig)

    def _add_electrical_plots(self, pdf, results):
        """Page 2: Electrical properties plots"""
        fig, axes = plt.subplots(2, 2, figsize=(11, 8.5))

        # Dummy data (in practice, obtained from results)
        T = np.linspace(100, 400, 20)
        n = results['n'] * np.exp(-0.1 / (8.617e-5 * T))
        mu = results['mu'] * (300 / T)**1.5
        sigma = n * 1.60218e-19 * mu
        R_H = 1 / (n * 1.60218e-19)

        # Plots
        axes[0, 0].semilogy(T, n / 1e6, 'o-', linewidth=2.5)
        axes[0, 0].set_xlabel('Temperature (K)')
        axes[0, 0].set_ylabel('n (cm⁻³)')
        axes[0, 0].set_title('Carrier Density vs T')
        axes[0, 0].grid(alpha=0.3)

        axes[0, 1].loglog(T, mu * 1e4, 's-', linewidth=2.5, color='#A23B72')
        axes[0, 1].set_xlabel('Temperature (K)')
        axes[0, 1].set_ylabel('μ (cm²/(V·s))')
        axes[0, 1].set_title('Mobility vs T')
        axes[0, 1].grid(alpha=0.3)

        axes[1, 0].semilogy(T, sigma, '^-', linewidth=2.5, color='#F18F01')
        axes[1, 0].set_xlabel('Temperature (K)')
        axes[1, 0].set_ylabel('σ (S/m)')
        axes[1, 0].set_title('Conductivity vs T')
        axes[1, 0].grid(alpha=0.3)

        axes[1, 1].plot(T, R_H, 'D-', linewidth=2.5, color='#C73E1D')
        axes[1, 1].set_xlabel('Temperature (K)')
        axes[1, 1].set_ylabel('R_H (m³/C)')
        axes[1, 1].set_title('Hall Coefficient vs T')
        axes[1, 1].grid(alpha=0.3)

        plt.suptitle('Electrical Properties', fontsize=16, fontweight='bold')
        plt.tight_layout()
        pdf.savefig(fig, bbox_inches='tight')
        plt.close(fig)

    def _add_magnetic_plots(self, pdf, results):
        """Page 3: Magnetic properties plots"""
        fig, axes = plt.subplots(2, 2, figsize=(11, 8.5))

        # M-H curve (dummy)
        H = np.linspace(-5000, 5000, 100)
        M = results['M_s'] * np.tanh(H / 1000)

        axes[0, 0].plot(H, M, linewidth=2.5, color='#f093fb')
        axes[0, 0].axhline(results['M_s'], linestyle='--', color='green', label=f"M_s = {results['M_s']:.1f}")
        axes[0, 0].axvline(results['H_c'], linestyle='--', color='red', label=f"H_c = {results['H_c']:.0f}")
        axes[0, 0].set_xlabel('H (Oe)')
        axes[0, 0].set_ylabel('M (emu/g)')
        axes[0, 0].set_title('M-H Hysteresis Loop')
        axes[0, 0].legend()
        axes[0, 0].grid(alpha=0.3)

        # Other plots omitted (would normally include temperature dependence, etc.)
        for ax in axes.flat[1:]:
            ax.text(0.5, 0.5, 'Additional magnetic\nproperties plots', ha='center', va='center', fontsize=14)
            ax.axis('off')

        plt.suptitle('Magnetic Properties', fontsize=16, fontweight='bold')
        plt.tight_layout()
        pdf.savefig(fig, bbox_inches='tight')
        plt.close(fig)

    def _add_statistics_page(self, pdf, results):
        """Page 4: Statistics"""
        fig = plt.figure(figsize=(8.5, 11))
        fig.text(0.5, 0.95, 'STATISTICAL SUMMARY', ha='center', fontsize=20, fontweight='bold')

        # Simple statistics table (would normally be more detailed)
        fig.text(0.1, 0.80, 'Measurement Quality Metrics:', fontsize=14, fontweight='bold')
        fig.text(0.1, 0.75, '  - Data points: 200', fontsize=12)
        fig.text(0.1, 0.70, '  - Outliers removed: 5 (2.5%)', fontsize=12)
        fig.text(0.1, 0.65, '  - Fit R²: 0.998', fontsize=12)
        fig.text(0.1, 0.60, '  - Residual std: 0.05', fontsize=12)

        plt.axis('off')
        pdf.savefig(fig, bbox_inches='tight')
        plt.close(fig)

# Example usage
results_example = {
    'R_s': 1370,
    'sigma': 2.43e3,
    'carrier_type': 'electron',
    'n': 2.36e20,
    'mu': 0.064,
    'M_s': 50,
    'H_c': 300,
    'M_r': 40,
    'S': 0.8,
    'K': 7.5e5
}

reporter = AutoReportGenerator(filename='integrated_analysis_report.pdf')
reporter.generate_report(results_example)

4.7 Complete Workflow Integration

Code Example 4-7: End-to-End Analysis Pipeline

"""
A complete end-to-end analysis pipeline

Workflow:
1. Data loading (CSV/DAT/Binary)
2. Data cleaning (outlier removal)
3. Integrated four-point probe + Hall + M-H analysis
4. Advanced fitting
5. Error propagation
6. Publication-quality figure generation
7. PDF report generation
"""

import numpy as np
import pandas as pd
from pathlib import Path

class EndToEndPipeline:
    """
    Complete analysis pipeline
    """

    def __init__(self, project_name: str):
        self.project_name = project_name
        self.loader = UniversalDataLoader()
        self.analyzer = IntegratedAnalysisPipeline(thickness=200e-9, mass=0.005)
        self.reporter = AutoReportGenerator(filename=f'{project_name}_report.pdf')

    def run(self, data_files: dict):
        """
        Run the pipeline

        Parameters
        ----------
        data_files : dict
            {'four_probe': 'path/to/file.csv',
             'hall': 'path/to/file.csv',
             'mh': 'path/to/file.csv'}
        """
        print("=" * 80)
        print(f"STARTING END-TO-END ANALYSIS PIPELINE: {self.project_name}")
        print("=" * 80)

        # Step 1: Load data
        print("\n[Step 1] Loading data...")
        data_four_probe = self.loader.load(data_files['four_probe'])
        data_hall = self.loader.load(data_files['hall'])
        data_mh = self.loader.load(data_files['mh'])

        # Step 2: Clean data
        print("\n[Step 2] Cleaning data...")
        data_four_probe_clean = self.loader.clean_data(outlier_method='iqr')

        self.loader.data = data_hall
        data_hall_clean = self.loader.clean_data(outlier_method='iqr')

        self.loader.data = data_mh
        data_mh_clean = self.loader.clean_data(outlier_method='iqr')

        # Step 3: Integrated analysis
        print("\n[Step 3] Integrated analysis...")
        self.analyzer.load_four_probe_data(
            R_AB_CD=data_four_probe_clean['R_AB_CD'].mean(),
            R_BC_DA=data_four_probe_clean['R_BC_DA'].mean()
        )
        self.analyzer.load_hall_data(
            I=data_hall_clean['I'].mean(),
            B=data_hall_clean['B'].mean(),
            V_pos=data_hall_clean['V_pos'].mean(),
            V_neg=data_hall_clean['V_neg'].mean()
        )
        self.analyzer.load_mh_data(
            H=data_mh_clean['H'].values,
            M=data_mh_clean['M'].values
        )

        self.analyzer.analyze_electrical_properties()
        self.analyzer.analyze_magnetic_properties()

        # Step 4: Generate report
        print("\n[Step 4] Generating report...")
        self.analyzer.generate_report()
        self.reporter.generate_report(self.analyzer.results)

        print("\n" + "=" * 80)
        print("PIPELINE COMPLETED")
        print("=" * 80)

        return self.analyzer.results

# Example usage (run with mock data)
if __name__ == '__main__':
    # Create mock data files (not needed in practice)
    np.random.seed(42)
    pd.DataFrame({
        'R_AB_CD': 1000 + 50 * np.random.randn(30),
        'R_BC_DA': 950 + 45 * np.random.randn(30)
    }).to_csv('four_probe_data.csv', index=False)

    pd.DataFrame({
        'I': 100e-6 + 1e-6 * np.random.randn(30),
        'B': 0.5 + 0.01 * np.random.randn(30),
        'V_pos': -5.0e-3 + 0.1e-3 * np.random.randn(30),
        'V_neg': 4.8e-3 + 0.1e-3 * np.random.randn(30)
    }).to_csv('hall_data.csv', index=False)

    H = np.linspace(-5000, 5000, 200)
    M = 50 * np.tanh(H / 1000) + 0.5 * np.random.randn(200)
    pd.DataFrame({'H': H, 'M': M}).to_csv('mh_data.csv', index=False)

    # Run the pipeline
    pipeline = EndToEndPipeline(project_name='Sample_Material_XYZ')
    results = pipeline.run({
        'four_probe': 'four_probe_data.csv',
        'hall': 'hall_data.csv',
        'mh': 'mh_data.csv'
    })

    print("\nFinal results saved to:")
    print("  - Sample_Material_XYZ_report.pdf")

4.8 Exercises

Exercise 4-1: Extending the Data Loader (Easy)

Easy Problem: Add a method _load_json() to the UniversalDataLoader class that loads data in JSON format.

import json

def _load_json(self, path: Path) -> pd.DataFrame:
    """Load JSON format"""
    with open(path, 'r') as f:
        data_dict = json.load(f)

    # Convert the dictionary into a DataFrame
    data = pd.DataFrame(data_dict)
    return data

# Attach to the class
UniversalDataLoader._load_json = _load_json

# Test
loader = UniversalDataLoader()
# loader.load('data.json', format='json')

Exercise 4-2: Visualizing Outlier Detection (Easy)

Easy Problem: Plot the outliers detected by the IQR method together with the original data (mark the outliers in red).

import matplotlib.pyplot as plt

data = np.random.randn(100)
data[95:] = 10  # Add outliers

Q1, Q3 = np.percentile(data, [25, 75])
IQR = Q3 - Q1
lower = Q1 - 1.5 * IQR
upper = Q3 + 1.5 * IQR

outliers = (data < lower) | (data > upper)

plt.figure(figsize=(10, 6))
plt.scatter(range(len(data)), data, c=outliers, cmap='RdYlGn_r', s=50, edgecolors='black')
plt.axhline(lower, color='blue', linestyle='--', label='Lower bound')
plt.axhline(upper, color='blue', linestyle='--', label='Upper bound')
plt.xlabel('Index')
plt.ylabel('Value')
plt.title('Outlier Detection (IQR method)')
plt.legend()
plt.colorbar(label='Outlier (True=Red)')
plt.show()

Exercise 4-3: Analyzing Fit Residuals (Medium)

Medium Problem: Using the fitting result from lmfit, check the normality of the residuals with a Q-Q plot.

import scipy.stats as stats

# Run the fit (from the earlier code example)
residuals = result.residual

# Q-Q plot
fig, ax = plt.subplots(figsize=(8, 6))
stats.probplot(residuals, dist="norm", plot=ax)
ax.set_title('Q-Q Plot: Residuals Normality Check', fontsize=14, fontweight='bold')
ax.grid(alpha=0.3)
plt.show()

# Statistical test (Shapiro-Wilk test)
stat, p_value = stats.shapiro(residuals)
print(f"Shapiro-Wilk test: statistic={stat:.4f}, p-value={p_value:.4f}")
if p_value > 0.05:
    print("  → Residuals are normally distributed (p > 0.05)")
else:
    print("  → Residuals are NOT normally distributed (p < 0.05)")

Exercise 4-4: Manual Error Propagation Calculation (Medium)

Medium Problem: Manually compute the uncertainty of the mobility $\mu = \sigma |R_H|$ using partial derivatives, and compare with the result from the uncertainties package.

sigma_val = 1e4
sigma_err = 100
R_H_val = 2.5e-3
R_H_err = 0.1e-3

# Manual calculation: δμ = sqrt((∂μ/∂σ)² δσ² + (∂μ/∂R_H)² δR_H²)
# ∂μ/∂σ = |R_H|
# ∂μ/∂R_H = σ * sign(R_H)

dmu_dsigma = R_H_val
dmu_dRH = sigma_val

delta_mu_manual = np.sqrt((dmu_dsigma * sigma_err)**2 + (dmu_dRH * R_H_err)**2)
mu_val = sigma_val * R_H_val

print(f"Manual calculation:")
print(f"  μ = {mu_val:.2e} m²/(V·s)")
print(f"  Δμ = {delta_mu_manual:.2e} m²/(V·s)")
print(f"  Relative uncertainty = {delta_mu_manual / mu_val * 100:.2f}%")

# uncertainties package
from uncertainties import ufloat
sigma_u = ufloat(sigma_val, sigma_err)
R_H_u = ufloat(R_H_val, R_H_err)
mu_u = sigma_u * R_H_u

print(f"\nuncertainties package:")
print(f"  μ = {mu_u}")
print(f"  Relative uncertainty = {mu_u.std_dev / mu_u.nominal_value * 100:.2f}%")

Exercise 4-5: Creating a Custom Report Template (Medium)

Medium Problem: Extend AutoReportGenerator by implementing a method add_custom_page() that lets users add their own figures.

def add_custom_page(self, pdf, fig):
    """
    Add a custom figure to the report

    Parameters
    ----------
    pdf : PdfPages
        The PDF page object
    fig : matplotlib.figure.Figure
        The figure to add
    """
    pdf.savefig(fig, bbox_inches='tight')
    plt.close(fig)

# Attach to the class
AutoReportGenerator.add_custom_page = add_custom_page

# Example usage
reporter = AutoReportGenerator('custom_report.pdf')
with PdfPages(reporter.filename) as pdf:
    # Custom page 1
    fig1 = plt.figure(figsize=(8.5, 11))
    plt.plot([1, 2, 3], [4, 5, 6])
    plt.title('Custom Plot 1')
    reporter.add_custom_page(pdf, fig1)

    # Custom page 2
    fig2 = plt.figure(figsize=(8.5, 11))
    plt.scatter([1, 2, 3], [6, 5, 4])
    plt.title('Custom Plot 2')
    reporter.add_custom_page(pdf, fig2)

print("Custom report generated")

Exercise 4-6: Batch Processing Pipeline (Hard)

Hard Problem: Create a pipeline that batch-processes data files from multiple samples and combines the results into a single Excel file (with multiple sheets).

import pandas as pd
from pathlib import Path

def batch_process_samples(data_dir: str, output_file: str = 'batch_results.xlsx'):
    """
    Batch-process multiple samples

    Parameters
    ----------
    data_dir : str
        Data directory (containing a subdirectory for each sample)
    output_file : str
        Output Excel file name
    """
    data_path = Path(data_dir)
    sample_dirs = [d for d in data_path.iterdir() if d.is_dir()]

    all_results = []

    for sample_dir in sample_dirs:
        sample_name = sample_dir.name
        print(f"Processing: {sample_name}")

        try:
            # Run the pipeline
            pipeline = EndToEndPipeline(project_name=sample_name)
            results = pipeline.run({
                'four_probe': sample_dir / 'four_probe.csv',
                'hall': sample_dir / 'hall.csv',
                'mh': sample_dir / 'mh.csv'
            })

            results['sample_name'] = sample_name
            all_results.append(results)

        except Exception as e:
            print(f"  ERROR: {e}")
            continue

    # Save the results to Excel
    with pd.ExcelWriter(output_file, engine='openpyxl') as writer:
        # Summary sheet
        df_summary = pd.DataFrame(all_results)
        df_summary.to_excel(writer, sheet_name='Summary', index=False)

        # Individual sample sheets
        for result in all_results:
            sample_name = result['sample_name']
            df_sample = pd.DataFrame([result])
            df_sample.to_excel(writer, sheet_name=sample_name[:31], index=False)  # Excel sheet name limit

    print(f"\nBatch results saved to: {output_file}")

# Example usage (assumed directory structure)
# data/
#   sample1/
#     four_probe.csv
#     hall.csv
#     mh.csv
#   sample2/
#     ...

# batch_process_samples('data/', 'all_samples_results.xlsx')

Exercise 4-7: Anomaly Detection Using Machine Learning (Hard)

Hard Problem: Use scikit-learn's Isolation Forest to detect anomalous data points in an M-H curve.

from sklearn.ensemble import IsolationForest
import numpy as np
import matplotlib.pyplot as plt

# Generate M-H data (with some anomalies)
H = np.linspace(-5000, 5000, 200)
M = 50 * np.tanh(H / 1000)
M[50:55] += 20  # Add anomalies
M[150:155] -= 15

# Build features (H, M, dM/dH)
dM_dH = np.gradient(M, H)
X = np.column_stack([H, M, dM_dH])

# Isolation Forest
clf = IsolationForest(contamination=0.05, random_state=42)
anomalies = clf.fit_predict(X)

# Plotting
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# Left: M-H curve
colors = ['red' if a == -1 else 'blue' for a in anomalies]
ax1.scatter(H, M, c=colors, s=50, alpha=0.7, edgecolors='black')
ax1.set_xlabel('H [Oe]')
ax1.set_ylabel('M [emu/g]')
ax1.set_title('M-H Curve with Anomaly Detection')
ax1.grid(alpha=0.3)

# Right: anomaly scores
scores = clf.decision_function(X)
ax2.plot(H, scores, linewidth=2, color='purple')
ax2.axhline(0, color='red', linestyle='--', label='Threshold')
ax2.set_xlabel('H [Oe]')
ax2.set_ylabel('Anomaly Score')
ax2.set_title('Isolation Forest Anomaly Scores')
ax2.legend()
ax2.grid(alpha=0.3)

plt.tight_layout()
plt.show()

print(f"Detected {np.sum(anomalies == -1)} anomalous points out of {len(anomalies)}")

Exercise 4-8: Experimental Design Optimization (Hard)

Hard Problem: Considering the trade-off between measurement time and accuracy, propose an experimental design (simulation-based) that achieves the target accuracy (error < 5%) with the minimum number of measurement points.

N_range = [10, 20, 30, 50, 100, 200]
fit_quality = []
measurement_time = []

for N in N_range:
    # Generate data
    T = np.linspace(100, 400, N)
    sigma_true = lambda T: 1e4 * np.exp(-0.2 / (8.617e-5 * T))
    sigma_data = sigma_true(T) * (1 + 0.05 * np.random.randn(N))

    # Fitting
    from scipy.optimize import curve_fit
    def model(T, A, Ea):
        return A * np.exp(-Ea / (8.617e-5 * T))

    params, _ = curve_fit(model, T, sigma_data)
    sigma_fit = model(T, *params)

    # Evaluate accuracy
    R2 = 1 - np.sum((sigma_data - sigma_fit)**2) / np.sum((sigma_data - np.mean(sigma_data))**2)
    fit_quality.append(R2)

    # Measurement time (1 point = 1 minute)
    measurement_time.append(N * 1)

# Plotting
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(measurement_time, fit_quality, 'o-', linewidth=2.5, markersize=10, color='#f093fb')
ax.axhline(0.95, color='red', linestyle='--', linewidth=2, label='Target R² = 0.95')
ax.set_xlabel('Measurement Time [minutes]', fontsize=12)
ax.set_ylabel('Fit Quality (R²)', fontsize=12)
ax.set_title('Trade-off: Measurement Time vs Accuracy', fontsize=14, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(alpha=0.3)
plt.show()

# Recommended number of measurement points
optimal_idx = np.argmin(np.abs(np.array(fit_quality) - 0.95))
print(f"Recommended: N = {N_range[optimal_idx]} points ({measurement_time[optimal_idx]} minutes)")

4.9 Checking Your Understanding

Use the checklist below to check your understanding:

Basic Understanding

• I understand the variety of formats used by real measurement instruments
• I can explain the need for data cleaning and the techniques involved
• I understand the significance of integrating four-point probe, Hall, and M-H measurements
• I can explain the advantages of fitting with lmfit
• I understand the importance of error propagation

Practical Skills

• I can implement a multi-format data loader
• I can remove outliers using the IQR and Z-score methods
• I can build an integrated analysis pipeline
• I can perform constrained fitting with lmfit
• I can compute error propagation using the uncertainties package
• I can create publication-quality figures
• I can automatically generate PDF reports

Applied Skills

• I can design a complete end-to-end pipeline
• I can perform anomaly detection using machine learning
• I can optimize an experimental design (trade-off between accuracy and time)
• I can create custom report templates

4.10 References

1. McKinney, W. (2017). Python for Data Analysis: Data Wrangling with Pandas, NumPy, and IPython (2nd ed.). O'Reilly. - Data processing with pandas
2. VanderPlas, J. (2016). Python Data Science Handbook. O'Reilly. - A comprehensive guide to scientific computing
3. Newville, M., et al. (2014). LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python. Zenodo. - lmfit documentation
4. Hunter, J. D. (2007). Matplotlib: A 2D Graphics Environment. Computing in Science & Engineering, 9(3), 90-95. - The original matplotlib paper
5. Lebigot, E. O. (2010). Uncertainties: a Python package for calculations with uncertainties. - The error propagation package
6. Pedregosa, F., et al. (2011). Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research, 12, 2825-2830. - scikit-learn
7. Schroder, D. K. (2006). Semiconductor Material and Device Characterization (3rd ed.). Wiley. - Practical measurement data analysis

4.11 Summary and Next Steps

In this chapter, you learned a complete Python workflow that spans everything from loading experimental data to generating a final report. You are now able to automate the entire process of electrical and magnetic measurement analysis.

What the series as a whole covers:

• Chapter 1: Four-point probe method, van der Pauw method, temperature-dependent measurements
• Chapter 2: Hall effect, determination of carrier density and mobility, the two-band model
• Chapter 3: VSM/SQUID magnetic measurements, M-H curve analysis, FC/ZFC
• Chapter 4: Integrated data analysis pipeline, automated report generation

Next steps:

• Run the pipeline on real measurement data
• Build machine learning models for predicting material properties
• Integrate with real-time measurement systems
• Turn the pipeline into a web application (Streamlit, Dash)