kishorkukreja/awesome-supply-chain/quality-management

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
class ControlCharts:
    """
    Statistical Process Control charts
    X-bar & R charts, p-charts, c-charts, etc.
    """
    def __init__(self, data, subgroup_size=5):
        """
        Parameters:
        - data: array of measurements or DataFrame
        - subgroup_size: sample size per subgroup
        """
        self.data = np.array(data)
        self.subgroup_size = subgroup_size
        self.n_subgroups = len(data) // subgroup_size
    def xbar_r_chart(self):
        """
        X-bar and R (Range) control chart for variables data
        Returns control limits and data for plotting
        """
        # Reshape into subgroups
        subgroups = self.data[:self.n_subgroups * self.subgroup_size].reshape(
            self.n_subgroups, self.subgroup_size
        )
        # Calculate subgroup means and ranges
        xbar = subgroups.mean(axis=1)
        R = subgroups.max(axis=1) - subgroups.min(axis=1)
        # Overall mean and average range
        xbar_mean = xbar.mean()
        R_mean = R.mean()
        # Control chart constants (for n=5, A2=0.577, D3=0, D4=2.114)
        # For other sample sizes, use full table
        constants = {
            2: {'A2': 1.880, 'D3': 0, 'D4': 3.267},
            3: {'A2': 1.023, 'D3': 0, 'D4': 2.574},
            4: {'A2': 0.729, 'D3': 0, 'D4': 2.282},
            5: {'A2': 0.577, 'D3': 0, 'D4': 2.114},
            6: {'A2': 0.483, 'D3': 0, 'D4': 2.004},
            7: {'A2': 0.419, 'D3': 0.076, 'D4': 1.924},
            8: {'A2': 0.373, 'D3': 0.136, 'D4': 1.864},
            9: {'A2': 0.337, 'D3': 0.184, 'D4': 1.816},
            10: {'A2': 0.308, 'D3': 0.223, 'D4': 1.777}
        }
        n = self.subgroup_size
        A2 = constants.get(n, constants[5])['A2']
        D3 = constants.get(n, constants[5])['D3']
        D4 = constants.get(n, constants[5])['D4']
        # X-bar chart limits
        xbar_ucl = xbar_mean + A2 * R_mean
        xbar_lcl = xbar_mean - A2 * R_mean
        # R chart limits
        r_ucl = D4 * R_mean
        r_lcl = D3 * R_mean
        # Detect out-of-control points
        xbar_ooc = (xbar > xbar_ucl) | (xbar < xbar_lcl)
        r_ooc = (R > r_ucl) | (R < r_lcl)
        return {
            'xbar': xbar,
            'xbar_mean': xbar_mean,
            'xbar_ucl': xbar_ucl,
            'xbar_lcl': xbar_lcl,
            'xbar_out_of_control': xbar_ooc,
            'R': R,
            'R_mean': R_mean,
            'R_ucl': r_ucl,
            'R_lcl': r_lcl,
            'R_out_of_control': r_ooc,
            'in_control': not (xbar_ooc.any() or r_ooc.any())
        }
    def p_chart(self, defects, sample_sizes):
        """
        p-chart for proportion defective (attribute data)
        Parameters:
        - defects: array of number of defectives per sample
        - sample_sizes: array of sample sizes (can be variable)
        """
        defects = np.array(defects)
        sample_sizes = np.array(sample_sizes)
        # Proportion defective per sample
        p = defects / sample_sizes
        # Average proportion defective
        p_bar = defects.sum() / sample_sizes.sum()
        # Control limits (3-sigma)
        # For variable sample size, calculate limits for each point
        if len(set(sample_sizes)) == 1:
            # Constant sample size
            n = sample_sizes[0]
            ucl = p_bar + 3 * np.sqrt(p_bar * (1 - p_bar) / n)
            lcl = max(0, p_bar - 3 * np.sqrt(p_bar * (1 - p_bar) / n))
            ucl = np.full_like(p, ucl)
            lcl = np.full_like(p, lcl)
        else:
            # Variable sample size
            ucl = p_bar + 3 * np.sqrt(p_bar * (1 - p_bar) / sample_sizes)
            lcl = np.maximum(0, p_bar - 3 * np.sqrt(p_bar * (1 - p_bar) / sample_sizes))
        # Out of control points
        ooc = (p > ucl) | (p < lcl)
        return {
            'p': p,
            'p_bar': p_bar,
            'ucl': ucl,
            'lcl': lcl,
            'out_of_control': ooc,
            'in_control': not ooc.any()
        }
    def c_chart(self, defects_per_unit):
        """
        c-chart for count of defects per unit
        Parameters:
        - defects_per_unit: array of defect counts per unit
        """
        c = np.array(defects_per_unit)
        # Average defects per unit
        c_bar = c.mean()
        # Control limits (based on Poisson distribution)
        ucl = c_bar + 3 * np.sqrt(c_bar)
        lcl = max(0, c_bar - 3 * np.sqrt(c_bar))
        # Out of control
        ooc = (c > ucl) | (c < lcl)
        return {
            'c': c,
            'c_bar': c_bar,
            'ucl': ucl,
            'lcl': lcl,
            'out_of_control': ooc,
            'in_control': not ooc.any()
        }
    def plot_xbar_r_chart(self, results):
        """Plot X-bar and R control charts"""
        fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
        # X-bar chart
        x = range(1, len(results['xbar']) + 1)
        ax1.plot(x, results['xbar'], 'bo-', label='Subgroup Mean')
        ax1.axhline(results['xbar_mean'], color='green', linestyle='-', linewidth=2, label='Center Line')
        ax1.axhline(results['xbar_ucl'], color='red', linestyle='--', linewidth=2, label='UCL')
        ax1.axhline(results['xbar_lcl'], color='red', linestyle='--', linewidth=2, label='LCL')
        # Mark out-of-control points
        ooc_indices = np.where(results['xbar_out_of_control'])[0]
        if len(ooc_indices) > 0:
            ax1.plot(ooc_indices + 1, results['xbar'][ooc_indices], 'rx', markersize=12, markeredgewidth=3)
        ax1.set_xlabel('Subgroup Number')
        ax1.set_ylabel('Subgroup Mean')
        ax1.set_title('X-bar Control Chart')
        ax1.legend()
        ax1.grid(True, alpha=0.3)
        # R chart
        ax2.plot(x, results['R'], 'bo-', label='Subgroup Range')
        ax2.axhline(results['R_mean'], color='green', linestyle='-', linewidth=2, label='Center Line')
        ax2.axhline(results['R_ucl'], color='red', linestyle='--', linewidth=2, label='UCL')
        ax2.axhline(results['R_lcl'], color='red', linestyle='--', linewidth=2, label='LCL')
        # Mark out-of-control points
        ooc_indices = np.where(results['R_out_of_control'])[0]
        if len(ooc_indices) > 0:
            ax2.plot(ooc_indices + 1, results['R'][ooc_indices], 'rx', markersize=12, markeredgewidth=3)
        ax2.set_xlabel('Subgroup Number')
        ax2.set_ylabel('Subgroup Range')
        ax2.set_title('R Control Chart')
        ax2.legend()
        ax2.grid(True, alpha=0.3)
        plt.tight_layout()
        return fig
# Example usage - Variable data (measurements)
np.random.seed(42)
# Simulate process data (mean=100, std=2)
data = np.random.normal(100, 2, 100)
spc = ControlCharts(data, subgroup_size=5)
results = spc.xbar_r_chart()
print("X-bar Chart:")
print(f"  Center Line (X-bar-bar): {results['xbar_mean']:.2f}")
print(f"  UCL: {results['xbar_ucl']:.2f}")
print(f"  LCL: {results['xbar_lcl']:.2f}")
print(f"  Process in control: {results['in_control']}")
print("\nR Chart:")
print(f"  Center Line (R-bar): {results['R_mean']:.2f}")
print(f"  UCL: {results['R_ucl']:.2f}")
print(f"  LCL: {results['R_lcl']:.2f}")
# Plot
fig = spc.plot_xbar_r_chart(results)
plt.show()
# Example - Attribute data (p-chart)
defects = np.array([5, 8, 3, 6, 12, 4, 7, 5, 9, 6])
sample_sizes = np.array([100] * 10)
p_chart_results = spc.p_chart(defects, sample_sizes)
print(f"\np-chart:")
print(f"  Average proportion defective: {p_chart_results['p_bar']:.4f} ({p_chart_results['p_bar']*100:.2f}%)")
print(f"  UCL: {p_chart_results['ucl'][0]:.4f}")
print(f"  LCL: {p_chart_results['lcl'][0]:.4f}")
print(f"  Process in control: {p_chart_results['in_control']}")

class ProcessCapability:
    """
    Process capability analysis (Cp, Cpk, Pp, Ppk)
    Measures how well process meets specifications
    """
    def __init__(self, data, lsl, usl, target=None):
        """
        Parameters:
        - data: array of process measurements
        - lsl: lower specification limit
        - usl: upper specification limit
        - target: target value (optional, default is midpoint)
        """
        self.data = np.array(data)
        self.lsl = lsl
        self.usl = usl
        self.target = target if target is not None else (lsl + usl) / 2
    def calculate_capability(self):
        """Calculate process capability indices"""
        # Process statistics
        mean = self.data.mean()
        std = self.data.std(ddof=1)  # Sample standard deviation
        # Within-subgroup variation (short-term capability)
        # For simplicity, using overall std; ideally use R-bar/d2 method
        sigma_within = std
        # Specification width
        spec_width = self.usl - self.lsl
        # Cp: Potential capability (assumes process centered)
        # Cp = (USL - LSL) / (6 * sigma)
        cp = spec_width / (6 * sigma_within)
        # Cpk: Actual capability (accounts for centering)
        # Cpk = min[(USL - mean)/(3*sigma), (mean - LSL)/(3*sigma)]
        cpu = (self.usl - mean) / (3 * sigma_within)
        cpl = (mean - self.lsl) / (3 * sigma_within)
        cpk = min(cpu, cpl)