import math

# 1. Finding a problem that affects only part of a lot (e.g. mould in some
#    bags). With n independent units drawn at random, the chance that at
#    least one is affected is 1 − (1 − p)^n.
def units_needed(p, confidence=0.95):
    return math.ceil(math.log(1 - confidence) / math.log(1 - p))

for p in (0.20, 0.05, 0.01):
    print(f"{100 * p:4.0f} % of units affected: sample {units_needed(p):3d} "
          f"units for a 95 % chance of catching at least one")

# 2. Estimating an average (e.g. polysaccharide %) to within ± d with 95 %
#    confidence when units vary with standard deviation s: n = (1.96 s / d)².
for s, d in ((3.0, 1.0), (3.0, 2.0), (6.0, 2.0)):
    n = math.ceil((1.96 * s / d) ** 2)
    print(f"sd {s:.0f} points, margin ± {d:.0f}: {n:3d} units to sample")
#   20 % of units affected: sample  14 units for a 95 % chance of catching at least one
#    5 % of units affected: sample  59 units for a 95 % chance of catching at least one
#    1 % of units affected: sample 299 units for a 95 % chance of catching at least one
# sd 3 points, margin ± 1:  35 units to sample
# sd 3 points, margin ± 2:   9 units to sample
# sd 6 points, margin ± 2:  35 units to sample

import math

# HPLC fingerprints as peak areas at the same retention times (invented).
reference = [12.0, 30.5, 8.2, 44.0, 5.1, 19.6, 3.3]
batches = {
    "batch A": [11.4, 31.9, 7.8, 42.6, 5.5, 20.3, 3.0],
    "batch B": [14.8, 24.1, 9.9, 39.0, 7.2, 25.5, 4.6],
    "batch C": [12.5, 29.8, 0.4, 47.3, 5.0, 18.9, 3.4],   # one peak missing
}

def cosine(a, b):
    dot = sum(x * y for x, y in zip(a, b))
    return dot / math.sqrt(sum(x * x for x in a) * sum(y * y for y in b))

def correlation(a, b):
    ma, mb = sum(a) / len(a), sum(b) / len(b)
    return cosine([x - ma for x in a], [y - mb for y in b])

for name, peaks in batches.items():
    worst = max(range(len(peaks)), key=lambda i: abs(peaks[i] / reference[i] - 1))
    print(f"{name}: cosine {cosine(reference, peaks):.3f}, "
          f"correlation {correlation(reference, peaks):.3f}, "
          f"largest change at peak {worst + 1} "
          f"({100 * (peaks[worst] / reference[worst] - 1):+.0f} %)")
# batch A: cosine 0.999, correlation 0.998, largest change at peak 7 (-9 %)
# batch B: cosine 0.984, correlation 0.967, largest change at peak 5 (+41 %)
# batch C: cosine 0.990, correlation 0.984, largest change at peak 3 (-95 %)

# Minimum water activity for growth (typical values, see the text).
AW_MIN = {
    "most Gram-negative bacteria": 0.97,
    "most Gram-positive bacteria": 0.90,
    "Staphylococcus aureus":       0.86,
    "most yeasts":                 0.88,
    "most moulds":                 0.80,
    "xerophilic moulds":           0.61,
}

def gamma_aw(aw, aw_min):
    """Relative growth rate (0 to 1) in the 'gamma' approach:
    linear from aw_min (no growth) to aw = 1 (no limitation)."""
    return max(0.0, (aw - aw_min) / (1 - aw_min))

levels = (0.95, 0.85, 0.70, 0.55, 0.30)
print("organism group               " + "  ".join(f"aw {a:.2f}" for a in levels))
for name, m in AW_MIN.items():
    print(f"{name:28s} " + "  ".join(f"{gamma_aw(a, m):7.2f}" for a in levels))
print("growth stops below about 0.6; 0.3 is chosen for texture and slower chemistry")
# organism group               aw 0.95  aw 0.85  aw 0.70  aw 0.55  aw 0.30
# most Gram-negative bacteria     0.00     0.00     0.00     0.00     0.00
# most Gram-positive bacteria     0.50     0.00     0.00     0.00     0.00
# Staphylococcus aureus           0.64     0.00     0.00     0.00     0.00
# most yeasts                     0.58     0.00     0.00     0.00     0.00
# most moulds                     0.75     0.25     0.00     0.00     0.00
# xerophilic moulds               0.87     0.62     0.23     0.00     0.00
# growth stops below about 0.6; 0.3 is chosen for texture and slower chemistry