Arko007/walnut-rancidity-predictor
0
1"""2Synthetic Walnut Storage Dataset Generator3Simulates Indian storage conditions using Arrhenius-based food chemistry kinetics.4"""5 6import numpy as np7import pandas as pd8from pathlib import Path9import random10 11np.random.seed(42)12random.seed(42)13 14# Physical constants15R = 8.314 # J/(mol·K) — gas constant16Ea = 80000.0 # J/mol — activation energy for lipid oxidation17A = 1.5e12 # pre-exponential factor (1/day)18 19# Rancidity threshold20PV_THRESHOLD = 5.0 # meq/kg21 22# Storage scenario profiles (Indian context)23SCENARIOS = {24 "cold_storage": {25 "temp_range": (2.0, 8.0),26 "hum_range": (40.0, 60.0),27 "weight": 0.20,28 },29 "hill_region": {30 "temp_range": (5.0, 20.0),31 "hum_range": (35.0, 65.0),32 "weight": 0.25,33 },34 "ambient_warehouse": {35 "temp_range": (18.0, 32.0),36 "hum_range": (50.0, 75.0),37 "weight": 0.30,38 },39 "hot_transport": {40 "temp_range": (28.0, 40.0),41 "hum_range": (55.0, 85.0),42 "weight": 0.25,43 },44}45 46SCENARIO_NAMES = list(SCENARIOS.keys())47SCENARIO_WEIGHTS = [SCENARIOS[s]["weight"] for s in SCENARIO_NAMES]48 49 50def sigmoid(x):51 return 1.0 / (1.0 + np.exp(-x))52 53 54def arrhenius_rate(T_celsius: float, humidity: float, moisture: float) -> float:55 """Compute daily oxidation rate constant via Arrhenius kinetics."""56 T_kelvin = T_celsius + 273.1557 k_base = A * np.exp(-Ea / (R * T_kelvin))58 59 # Humidity and moisture accelerate oxidation60 h_factor = 1.0 + 0.015 * (humidity - 40.0) # relative to 40 % RH61 m_factor = 1.0 + 0.20 * (moisture - 4.0) # relative to 4 % moisture62 h_factor = max(h_factor, 0.5)63 m_factor = max(m_factor, 0.5)64 65 return k_base * h_factor * m_factor66 67 68def generate_sequence(69 seq_len: int,70 scenario: str,71 start_day: int = 0,72) -> list[dict]:73 """74 Generate a single walnut storage sequence.75 Returns a list of daily records.76 """77 scen = SCENARIOS[scenario]78 temp_lo, temp_hi = scen["temp_range"]79 hum_lo, hum_hi = scen["hum_range"]80 81 # Initial chemical state82 PV0 = np.random.uniform(0.3, 1.2) # initial peroxide value83 FFA0 = np.random.uniform(0.05, 0.15) # initial free fatty acids84 base_oxy = np.random.uniform(0.18, 0.23)85 86 # Day when PV will exceed threshold (for shelf-life label)87 # We estimate it analytically for the "average" conditions of this scenario88 T_avg = (temp_lo + temp_hi) / 2.089 H_avg = (hum_lo + hum_hi) / 2.090 M_avg = np.random.uniform(3.0, 8.0)91 k_avg = arrhenius_rate(T_avg, H_avg, M_avg)92 if k_avg > 0 and PV0 < PV_THRESHOLD:93 days_to_rancid = np.log(PV_THRESHOLD / PV0) / k_avg94 else:95 days_to_rancid = 0.096 97 records = []98 PV = PV099 FFA = FFA0100 101 for i in range(seq_len):102 day = start_day + i103 104 # Daily environmental fluctuation105 temp = np.clip(np.random.normal((temp_lo + temp_hi) / 2,106 (temp_hi - temp_lo) / 6),107 temp_lo, temp_hi)108 humidity = np.clip(np.random.normal((hum_lo + hum_hi) / 2,109 (hum_hi - hum_lo) / 6),110 hum_lo, hum_hi)111 moisture = np.clip(np.random.normal(5.0, 1.2), 3.0, 8.0)112 oxygen = np.clip(np.random.normal(base_oxy, 0.005), 0.18, 0.23)113 114 k = arrhenius_rate(temp, humidity, moisture)115 116 # Peroxide value (Arrhenius model, Euler step)117 PV = PV * np.exp(k) # PV(t+1) = PV(t) * exp(k)118 PV = max(PV, 0.01)119 120 # Derived oxidation indicators121 FFA = FFA + 0.002 * k * PV + np.random.normal(0, 0.002)122 hexanal = 0.15 * PV + 0.05 * max(PV - 2.0, 0) ** 1.3 \123 + np.random.normal(0, 0.05)124 oxidation_ix = 0.3 * PV + 0.4 * max(FFA, 0) + np.random.normal(0, 0.05)125 126 FFA = max(FFA, 0.01)127 hexanal = max(hexanal, 0.0)128 oxidation_ix = max(oxidation_ix, 0.0)129 130 # Targets131 # rancidity_probability: sigmoid centred at PV = 5132 rand_prob = float(sigmoid(PV - PV_THRESHOLD))133 134 # shelf_life_remaining_days: days until PV > threshold from NOW135 shelf_life = float(max(days_to_rancid - day, 0.0))136 137 # decay_curve_value: normalised PV on [0,1] scale for regression138 decay_curve = float(min(PV / 10.0, 1.0))139 140 records.append({141 "day": day,142 "temperature": round(float(temp), 3),143 "humidity": round(float(humidity), 3),144 "moisture": round(float(moisture), 3),145 "oxygen": round(float(oxygen), 5),146 "peroxide_value": round(float(PV), 4),147 "free_fatty_acids": round(float(FFA), 4),148 "hexanal_level": round(float(hexanal), 4),149 "oxidation_index": round(float(oxidation_ix), 4),150 "rancidity_probability": round(rand_prob, 6),151 "shelf_life_remaining_days": round(shelf_life, 2),152 "decay_curve_value": round(decay_curve, 6),153 })154 155 return records156 157 158def generate_dataset(target_sequences: int = 90000) -> pd.DataFrame:159 """Generate the full dataset as a flat CSV."""160 all_records = []161 seq_count = 0162 163 print(f"Generating {target_sequences} sequences …")164 log_every = target_sequences // 10165 166 while seq_count < target_sequences:167 scenario = random.choices(SCENARIO_NAMES, weights=SCENARIO_WEIGHTS, k=1)[0]168 seq_len = random.randint(30, 90)169 start_day = random.randint(0, 90)170 171 records = generate_sequence(seq_len, scenario, start_day)172 173 # Tag each record with a sequence_id174 seq_id = seq_count175 for rec in records:176 rec["sequence_id"] = seq_id177 all_records.extend(records)178 179 seq_count += 1180 if seq_count % log_every == 0:181 print(f" {seq_count}/{target_sequences} sequences "182 f"({len(all_records):,} rows)")183 184 df = pd.DataFrame(all_records)185 # Reorder columns186 cols = [187 "sequence_id", "day",188 "temperature", "humidity", "moisture", "oxygen",189 "peroxide_value", "free_fatty_acids", "hexanal_level", "oxidation_index",190 "rancidity_probability", "shelf_life_remaining_days", "decay_curve_value",191 ]192 df = df[cols]193 return df194 195 196if __name__ == "__main__":197 out_dir = Path("data")198 out_dir.mkdir(exist_ok=True)199 200 df = generate_dataset(target_sequences=90000)201 202 out_path = out_dir / "walnut_storage_timeseries.csv"203 df.to_csv(out_path, index=False)204 print(f"\nDataset saved → {out_path}")205 print(f"Total rows : {len(df):,}")206 print(f"Total seqs : {df['sequence_id'].nunique():,}")207 print(df.describe().to_string())208