Cherry Blossom Forecasting

def prepare_features(year, temp_data, weather_data, bloom_data):
    """
    Prepare comprehensive feature set for hybrid model
    
    Returns dictionary of 20 features capturing:
    - Temperature patterns (hibernation and spring)
    - Photoperiod effects
    - Atmospheric conditions
    - Temporal trends
    """
    features = {}
    
    # === 1. HIBERNATION PHASE FEATURES ===
    nov_1 = datetime(year - 1, 11, 1)
    feb_1 = datetime(year, 2, 1)
    
    hibernation_temps = []
    for d in range((feb_1 - nov_1).days):
        date = nov_1 + timedelta(days=d)
        temp = get_avg_temp(temp_data, date)
        hibernation_temps.append(temp)
    
    features['hib_mean_temp'] = np.mean(hibernation_temps)
    features['hib_std_temp'] = np.std(hibernation_temps)
    features['hib_min_temp'] = np.min(hibernation_temps)
    
    # DTS-style accumulation
    features['hib_accumulated'] = sum([
        max(0, t) for t in hibernation_temps
    ])
    
    # Yoshimi-style coldest period
    jan_15 = datetime(year, 1, 15)
    features['coldest_period'] = find_coldest_period(
        temp_data, jan_15, 60
    )
    
    # === 2. EARLY SPRING FEATURES ===
    mar_1 = datetime(year, 3, 1)
    spring_temps = []
    for d in range(30):  # March data
        date = mar_1 + timedelta(days=d)
        temp = get_avg_temp(temp_data, date)
        spring_temps.append(temp)
    
    features['spring_temp_mean'] = np.mean(spring_temps)
    features['spring_temp_max'] = np.max(spring_temps)
    
    # Temperature trend (warming rate)
    x = np.arange(len(spring_temps))
    slope = np.polyfit(x, spring_temps, 1)[0]
    features['spring_temp_trend'] = slope
    
    # === 3. PHOTOPERIOD FEATURES ===
    features['feb_1_day_length'] = calculate_day_length(feb_1)
    features['mar_1_day_length'] = calculate_day_length(mar_1)
    features['day_length_change'] = (
        features['mar_1_day_length'] - 
        features['feb_1_day_length']
    )
    
    # === 4. GROWTH PERIOD FEATURES ===
    # Days above growth threshold
    growth_days = sum([1 for t in spring_temps if t > 10])
    features['days_above_10C'] = growth_days
    
    # Heat accumulation rate
    heat_sum = sum([max(0, t-5) for t in spring_temps])
    features['heat_accumulation_rate'] = heat_sum / 30
    
    # === 5. ATMOSPHERIC FEATURES ===
    # Humidity
    humidity_data = get_humidity(weather_data, year, feb_1, mar_1)
    features['avg_humidity'] = np.mean(humidity_data)
    features['std_humidity'] = np.std(humidity_data)
    
    # Precipitation
    precip_data = get_precipitation(weather_data, year, feb_1, mar_1)
    features['total_precipitation'] = np.sum(precip_data)
    features['rainy_days'] = np.sum(precip_data > 1.0)
    
    # Sea-level pressure
    pressure_data = get_pressure(weather_data, year, feb_1, mar_1)
    features['avg_pressure'] = np.mean(pressure_data)
    features['pressure_variability'] = np.std(pressure_data)
    
    # === 6. TEMPORAL CONTEXT ===
    # Recent climate trend
    recent_years = range(year - 5, year)
    recent_temps = []
    for y in recent_years:
        year_data = temp_data[temp_data['Year'] == y]
        winter_spring = year_data[
            (year_data['Month'] >= 11) | (year_data['Month'] <= 4)
        ]
        recent_temps.append(winter_spring['Avg_Temp'].mean())
    
    trend = np.polyfit(range(5), recent_temps, 1)[0]
    features['climate_trend'] = trend
    
    # Previous year's bloom (if available)
    prev_bloom = get_bloom_day_of_year(bloom_data, year - 1)
    features['prev_year_bloom'] = prev_bloom if prev_bloom else 85
    
    return features

from sklearn.ensemble import RandomForestRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import cross_val_score

def train_hybrid_model(training_years, temp_data, 
                      weather_data, bloom_data):
    """
    Train the hybrid Random Forest model
    
    Uses 5-fold cross-validation for parameter tuning
    Returns trained model and scaler
    """
    X_train = []
    y_train = []
    
    # Prepare training data
    print("Preparing features for training...")
    for year in training_years:
        features = prepare_features(
            year, temp_data, weather_data, bloom_data
        )
        feature_vector = list(features.values())
        X_train.append(feature_vector)
        
        # Target: day of year when bloom occurs
        actual_bloom = get_bloom_day_of_year(bloom_data, year)
        y_train.append(actual_bloom)
    
    X_train = np.array(X_train)
    y_train = np.array(y_train)
    
    print(f"Training set: {len(X_train)} samples, "
          f"{X_train.shape[1]} features")
    
    # Feature scaling (helpful even for trees)
    scaler = StandardScaler()
    X_train_scaled = scaler.fit_transform(X_train)
    
    # Hyperparameter tuning through cross-validation
    best_score = -np.inf
    best_params = None
    
    for n_est in [50, 100, 150]:
        for max_d in [8, 10, 12]:
            for min_split in [3, 5, 7]:
                model = RandomForestRegressor(
                    n_estimators=n_est,
                    max_depth=max_d,
                    min_samples_split=min_split,
                    random_state=42,
                    n_jobs=-1
                )
                
                # 5-fold cross-validation
                scores = cross_val_score(
                    model, X_train_scaled, y_train,
                    cv=5, scoring='neg_mean_absolute_error'
                )
                mean_score = scores.mean()
                
                if mean_score > best_score:
                    best_score = mean_score
                    best_params = {
                        'n_estimators': n_est,
                        'max_depth': max_d,
                        'min_samples_split': min_split
                    }
    
    print(f"Best params: {best_params}")
    print(f"CV MAE: {-best_score:.2f} days")
    
    # Train final model with best parameters
    final_model = RandomForestRegressor(
        **best_params,
        random_state=42,
        n_jobs=-1
    )
    final_model.fit(X_train_scaled, y_train)
    
    return final_model, scaler

# Train the model
hybrid_model, scaler = train_hybrid_model(
    training_years, temp_data, weather_data, bloom_data
)

# Extract feature importances
feature_names = [
    'hib_mean_temp', 'hib_std_temp', 'hib_min_temp',
    'hib_accumulated', 'coldest_period',
    'spring_temp_mean', 'spring_temp_max', 'spring_temp_trend',
    'feb_1_day_length', 'mar_1_day_length', 'day_length_change',
    'days_above_10C', 'heat_accumulation_rate',
    'avg_humidity', 'std_humidity',
    'total_precipitation', 'rainy_days',
    'avg_pressure', 'pressure_variability',
    'climate_trend', 'prev_year_bloom'
]

importances = hybrid_model.feature_importances_
importance_df = pd.DataFrame({
    'Feature': feature_names,
    'Importance': importances
}).sort_values('Importance', ascending=False)

print("Top 10 Most Important Features:")
print("=" * 50)
for idx, row in importance_df.head(10).iterrows():
    print(f"{row['Feature']:30s} {row['Importance']:.4f}")

def predict_hybrid_model(model, scaler, year, 
                        temp_data, weather_data, bloom_data):
    """
    Make bloom date prediction using hybrid model
    """
    # Extract features for target year
    features = prepare_features(
        year, temp_data, weather_data, bloom_data
    )
    feature_vector = np.array(list(features.values())).reshape(1, -1)
    
    # Scale features
    feature_vector_scaled = scaler.transform(feature_vector)
    
    # Predict day of year
    predicted_day = model.predict(feature_vector_scaled)[0]
    
    # Convert to actual date
    predicted_bloom = day_of_year_to_date(year, int(predicted_day))
    
    return predicted_bloom, predicted_day

# Evaluate on test years
print("\nHybrid Model Test Results")
print("=" * 70)
hybrid_predictions = []

for year in test_years:
    pred_date, pred_day = predict_hybrid_model(
        hybrid_model, scaler, year, 
        temp_data, weather_data, bloom_data
    )
    actual_date = get_actual_bloom(bloom_data, year)
    actual_day = actual_date.timetuple().tm_yday
    
    error = pred_day - actual_day
    
    hybrid_predictions.append({
        'Year': year,
        'Predicted': pred_date.strftime('%b %d'),
        'Actual': actual_date.strftime('%b %d'),
        'Error': int(error)
    })
    
    print(f"{year}: Pred={pred_date.strftime('%b %d')}, "
          f"Actual={actual_date.strftime('%b %d')}, "
          f"Error={error:+.0f} days")

# Calculate metrics
errors = [abs(p['Error']) for p in hybrid_predictions]
mae = np.mean(errors)
rmse = np.sqrt(np.mean([e**2 for e in errors]))

print(f"\nMAE: {mae:.2f} days")
print(f"RMSE: {rmse:.2f} days")