# DATA CLEANING AND DATA PREPROCESSING

import pandas as pd
import numpy as np
import requests

# Data Cleaning

# # Load and inspect data
typhoon_data = pd.read_csv('typhoon_data.csv')
earthquake_data = pd.read_csv('earthquake_data.csv')
socioeconomic_data = pd.read_csv('socioeconomic_data.csv')

# Data Preprocessing

# 1. Data Cleaning
# Handling missing values
typhoon_data = typhoon_data.dropna(subset=['Date', 'Intensity'])  # Drop rows with critical missing values
typhoon_data['Intensity'] = typhoon_data['Intensity'].fillna(typhoon_data['Intensity'].mean())  # Impute average value

# Handling outliers in numerical data (example: Intensity)
q1 = typhoon_data['Intensity'].quantile(0.25)
q3 = typhoon_data['Intensity'].quantile(0.75)
iqr = q3 - q1
upper_bound = q3 + 1.5 * iqr
lower_bound = q1 - 1.5 * iqr
typhoon_data = typhoon_data[(typhoon_data['Intensity'] >= lower_bound) & (typhoon_data['Intensity'] <= upper_bound)]

# 2. Data Transformation
# Standardize numerical columns (example: Population Density)
socioeconomic_data['Population_Density'] = (socioeconomic_data['Population_Density'] - socioeconomic_data['Population_Density'].mean()) / socioeconomic_data['Population_Density'].std()

# Encode categorical columns (if necessary)
typhoon_data['Region'] = typhoon_data['Region'].astype('category').cat.codes

# 3. Data Integration and Merging
# Merge datasets on common key (e.g., Region)
combined_data = pd.merge(typhoon_data, earthquake_data, on='Region', how='outer')
combined_data = pd.merge(combined_data, socioeconomic_data, on='Region', how='left')

# 4. Feature Engineering
# Creating a lag feature for disaster frequency by region
combined_data['Lag_Disaster_Count'] = combined_data.groupby('Region')['Date_x'].transform(lambda x: x.diff().dt.days)

# 5. Data Reduction (Optional: Reduce dimensionality)
# Example: Dropping columns that are not needed for analysis
combined_data = combined_data.drop(columns=['Unnecessary_Column'])

# 6. Data Validation
# Check for data consistency and types
print("\nCombined Data Information:")
print(combined_data.info())

# Check value distributions
print("\nSummary Statistics:")
print(combined_data.describe())

# 7. Documentation and Logging
# Save processed data to CSV for traceability
combined_data.to_csv('processed_combined_data.csv', index=False)
print("\nData preprocessing completed and saved.")


# EXPLORATORY DATA ANALYSIS (EDA)

1. Basic Data Exploration
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# Load the dataset
combined_data = pd.read_csv('processed_combined_data.csv')

# 1. Summary Statistics
print("Summary Statistics:\n", combined_data.describe())

# 2. Data Distribution Analysis
numerical_columns = ['Intensity', 'Population_Density', 'Economic_Impact_Index', 'Historical_Frequency']
for col in numerical_columns:
    plt.figure(figsize=(8, 4))
    sns.histplot(combined_data[col], kde=True)
    plt.title(f'Distribution of {col}')
    plt.show()

# 3. Missing Values Analysis
missing_values = combined_data.isnull().sum()
print("\nMissing Values:\n", missing_values)
sns.heatmap(combined_data.isnull(), cbar=False, cmap='viridis')
plt.title('Heatmap of Missing Values')
plt.show()

# 4. Outlier Detection
for col in numerical_columns:
    plt.figure(figsize=(8, 4))
    sns.boxplot(x=combined_data[col])
    plt.title(f'Box Plot for {col}')
    plt.show()

2. Temporal Analysis
# 5. Time Series Analysis
combined_data['Date'] = pd.to_datetime(combined_data['Date'])
plt.figure(figsize=(14, 6))
sns.lineplot(data=combined_data, x='Date', y='Intensity', hue='Disaster_Type')
plt.title('Disaster Intensity Over Time')
plt.xlabel('Date')
plt.ylabel('Intensity')
plt.show()

# 6. Seasonal Decomposition
import statsmodels.api as sm

disaster_counts = combined_data['Date'].dt.to_period('M').value_counts().sort_index()
decomposition = sm.tsa.seasonal_decompose(disaster_counts, model='additive')
decomposition.plot()
plt.show()

# 7. Year-on-Year Comparison
combined_data['Year'] = combined_data['Date'].dt.year
sns.countplot(data=combined_data, x='Year', hue='Disaster_Type')
plt.title('Year-on-Year Disaster Count Comparison')
plt.show()

# 8. Monthly/Seasonal Patterns
combined_data['Month'] = combined_data['Date'].dt.month
sns.countplot(data=combined_data, x='Month', hue='Disaster_Type')
plt.title('Monthly Disaster Count')
plt.show()

3. Geospatial Analysis
import folium

# 9. Mapping Disaster Incidents
map_center = [13.41, 122.56]
disaster_map = folium.Map(location=map_center, zoom_start=6)

for _, row in combined_data.iterrows():
    folium.CircleMarker(
        location=[row['Latitude'], row['Longitude']],
        radius=3,
        color='red' if row['Disaster_Type'] == 'Typhoon' else 'blue',
        fill=True,
        fill_opacity=0.6,
        popup=f"{row['Disaster_Type']} - {row['Date'].year}"
    ).add_to(disaster_map)

disaster_map.save("disaster_map.html")

# 10. Heatmaps for disaster frequency can be done using libraries like geopandas if more complex visualizations are needed.

4. Multivariate Analysis
# 11. Correlation Matrix
plt.figure(figsize=(10, 8))
sns.heatmap(combined_data.corr(), annot=True, cmap='coolwarm')
plt.title('Correlation Matrix')
plt.show()

# 12. Pair Plot Analysis
sns.pairplot(combined_data[numerical_columns + ['Disaster_Type']], hue='Disaster_Type', diag_kind='kde')
plt.show()

# 13. Scatter Plot with Regression Line
sns.lmplot(data=combined_data, x='Intensity', y='Economic_Impact_Index', hue='Disaster_Type')
plt.title('Intensity vs Economic Impact')
plt.show()

5. Statistical Analysis
from scipy import stats

# 14. Hypothesis Testing (t-test between typhoon and earthquake intensities)
typhoon_intensity = combined_data[combined_data['Disaster_Type'] == 'Typhoon']['Intensity']
earthquake_intensity = combined_data[combined_data['Disaster_Type'] == 'Earthquake']['Intensity']
t_stat, p_value = stats.ttest_ind(typhoon_intensity, earthquake_intensity)
print(f"T-test result: t-stat={t_stat}, p-value={p_value}")

# 15. Chi-Square Test (Region vs Disaster Type)
contingency_table = pd.crosstab(combined_data['Region'], combined_data['Disaster_Type'])
chi2, p, dof, ex = stats.chi2_contingency(contingency_table)
print(f"Chi-Square Test: chi2={chi2}, p-value={p}")

# 16. Feature Importance Analysis using tree-based models can be implemented in predictive modeling.

6. Categorical Data Analysis
# 18. Bar Plots for counts
sns.countplot(data=combined_data, x='Disaster_Type', order=combined_data['Disaster_Type'].value_counts().index)
plt.title('Disaster Count by Type')
plt.show()

# 19. Pie Chart for proportion of disaster types
combined_data['Disaster_Type'].value_counts().plot.pie(autopct='%1.1f%%', figsize=(8, 8))
plt.title('Proportion of Disaster Types')
plt.show()

# 20. Count Plots for regional analysis
sns.countplot(data=combined_data, y='Region', order=combined_data['Region'].value_counts().index, hue='Disaster_Type')
plt.title('Disaster Count by Region')
plt.show()

7. Advanced Visualizations
# 21. Time-Based Heatmaps
monthly_yearly_counts = combined_data.groupby([combined_data['Date'].dt.year, combined_data['Date'].dt.month]).size().unstack()
plt.figure(figsize=(12, 8))
sns.heatmap(monthly_yearly_counts, cmap='Blues', annot=True, fmt='d')
plt.title('Disaster Occurrences Heatmap')
plt.xlabel('Month')
plt.ylabel('Year')
plt.show()

# 22. Animated Plots can be implemented using Plotly or matplotlib animation for more dynamic visualizations.

# 23. Sunburst Chart (using Plotly for hierarchical visualization)
import plotly.express as px

sunburst_data = combined_data.groupby(['Disaster_Type', 'Region']).size().reset_index(name='Count')
fig = px.sunburst(sunburst_data, path=['Disaster_Type', 'Region'], values='Count', title='Disaster Distribution by Type and Region')
fig.show()

8. Additional Exploratory Analyses
# 24. Lag Analysis
combined_data['Lag_Disaster_Intensity'] = combined_data.groupby('Region')['Intensity'].shift(1)
sns.scatterplot(data=combined_data, x='Lag_Disaster_Intensity', y='Intensity', hue='Disaster_Type')
plt.title('Lag Analysis: Previous vs Current Intensity')
plt.show()

# 25. Moving Averages
combined_data['Rolling_Intensity'] = combined_data['Intensity'].rolling(window=3).mean()
plt.figure(figsize=(14, 6))
sns.lineplot(data=combined_data, x='Date', y='Rolling_Intensity')
plt.title('3-Month Moving Average of Intensity')
plt.show()

# 26. Clustering Analysis (e.g., K-Means)
from sklearn.cluster import KMeans

kmeans = KMeans(n_clusters=3)
combined_data['Cluster'] = kmeans.fit_predict(combined_data[numerical_columns])
sns.scatterplot(data=combined_data, x='Intensity', y='Economic_Impact_Index', hue='Cluster', palette='viridis')
plt.title('Clustering Analysis')
plt.show()

# 27. PCA (Principal Component Analysis)
from sklearn.decomposition import PCA

pca = PCA(n_components=2)
pca_result = pca.fit_transform(combined_data[numerical_columns])
combined_data['PCA1'] = pca_result[:, 0]
combined_data['PCA2'] = pca_result[:, 1]
sns.scatterplot(x='PCA1', y='PCA2', hue='Disaster_Type', data=combined_data)
plt.title('PCA Result')
plt.show()

9. Insights and Recommendation
# 28. Trend Identification
plt.figure(figsize=(14, 6))
sns.lineplot(data=combined_data, x='Date', y='Intensity', hue='Region')
plt.title('Trends in Disaster Intensity by Region')
plt.show()

# 29. Regional Vulnerability Analysis
vulnerability = combined_data.groupby('Region')[['Intensity', 'Economic_Impact_Index']].mean().sort_values(by='Intensity', ascending=False)
print("Regional Vulnerability Analysis:\n", vulnerability)

# 30. Anomaly Detection
combined_data['Anomaly'] = (combined_data['Intensity'] > combined_data['Intensity'].quantile(0.95)).astype(int)
sns.scatterplot(data=combined_data, x='Date', y='Intensity', hue='Anomaly')
plt.title('Anomaly Detection in Disaster Intensity')
plt.show()


# ADVANCED PREDICTIVE MODELING

# -*- coding: utf-8 -*-
"""
Advanced Predictive Modeling Script for Disaster Preparedness and Response
Author: Jahaziel Titular
Description: This script performs advanced predictive modeling using Python, covering
the entire process from data preparation to post-modeling analysis.
"""

# Section A: Data Preparation for Modeling
# Import necessary libraries
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split, GridSearchCV, TimeSeriesSplit
from sklearn.preprocessing import StandardScaler, LabelEncoder
from sklearn.impute import SimpleImputer
import matplotlib.pyplot as plt
import seaborn as sns

# Load dataset
print("Loading data...")
data = pd.read_csv('processed_combined_data.csv', parse_dates=['Date'])
print("Data loaded successfully.\n")

# 1. Feature Selection and Engineering
print("Performing feature engineering...")
data['Year'] = data['Date'].dt.year
data['Month'] = data['Date'].dt.month
data['Lag_Intensity'] = data.groupby('Region')['Intensity'].shift(1)

# Drop rows with NaN in lag features
data.dropna(subset=['Lag_Intensity'], inplace=True)

# 2. Data Splitting
print("Splitting data into training and test sets...")
X = data[['Intensity', 'Population_Density', 'Historical_Frequency', 'Lag_Intensity', 'Month']]
y_classification = data['Disaster_Type_Label']
y_regression = data['Economic_Impact_Index']

X_train, X_test, y_train_class, y_test_class = train_test_split(X, y_classification, test_size=0.3, random_state=42)
X_train_reg, X_test_reg, y_train_reg, y_test_reg = train_test_split(X, y_regression, test_size=0.3, random_state=42)
print("Data splitting completed.\n")

# 3. Data Imputation
print("Imputing missing values...")
imputer = SimpleImputer(strategy='mean')
X_train = pd.DataFrame(imputer.fit_transform(X_train), columns=X.columns)
X_test = pd.DataFrame(imputer.transform(X_test), columns=X.columns)
print("Imputation completed.\n")

# Section B: Model Selection and Training
from sklearn.ensemble import RandomForestClassifier, GradientBoostingRegressor
from xgboost import XGBClassifier, XGBRegressor

# 4. Model Initialization
print("Initializing models...")
rf_clf = RandomForestClassifier(random_state=42)
xgb_clf = XGBClassifier(random_state=42)
gbr_reg = GradientBoostingRegressor(random_state=42)
xgb_reg = XGBRegressor(random_state=42)
print("Models initialized.\n")

# 5. Hyperparameter Tuning
print("Performing hyperparameter tuning...")
param_grid_clf = {
    'n_estimators': [50, 100],
    'max_depth': [3, 5, 7]
}
grid_clf = GridSearchCV(xgb_clf, param_grid_clf, cv=3, scoring='accuracy')
grid_clf.fit(X_train, y_train_class)
best_clf = grid_clf.best_estimator_
print("Best classifier parameters found:", grid_clf.best_params_, "\n")

param_grid_reg = {
    'n_estimators': [100, 150],
    'learning_rate': [0.01, 0.1],
    'max_depth': [3, 5]
}
grid_reg = GridSearchCV(gbr_reg, param_grid_reg, cv=3, scoring='neg_mean_squared_error')
grid_reg.fit(X_train_reg, y_train_reg)
best_reg = grid_reg.best_estimator_
print("Best regressor parameters found:", grid_reg.best_params_, "\n")

# Section C: Model Evaluation
from sklearn.metrics import accuracy_score, mean_squared_error, mean_absolute_error, classification_report

# 8. Performance Metrics - Classification
print("Evaluating classification model...")
y_pred_class = best_clf.predict(X_test)
print("Classification Report:\n", classification_report(y_test_class, y_pred_class))
print("Accuracy:", accuracy_score(y_test_class, y_pred_class), "\n")

# 8. Performance Metrics - Regression
print("Evaluating regression model...")
y_pred_reg = best_reg.predict(X_test_reg)
print("Mean Absolute Error:", mean_absolute_error(y_test_reg, y_pred_reg))
print("Mean Squared Error:", mean_squared_error(y_test_reg, y_pred_reg), "\n")

# Section D: Model Optimization and Tuning
from sklearn.ensemble import StackingClassifier, VotingClassifier, StackingRegressor

# 11. Ensemble Techniques
print("Creating ensemble models...")
stacking_clf = StackingClassifier(
    estimators=[('rf', rf_clf), ('xgb', xgb_clf)],
    final_estimator=RandomForestClassifier(random_state=42)
)
stacking_clf.fit(X_train, y_train_class)

stacking_reg = StackingRegressor(
    estimators=[('gbr', gbr_reg), ('xgb', xgb_reg)],
    final_estimator=GradientBoostingRegressor(random_state=42)
)
stacking_reg.fit(X_train_reg, y_train_reg)
print("Ensemble models created.\n")

# 13. Refinement and Retraining
print("Refining models with cross-validation...")
cv = TimeSeriesSplit(n_splits=5)
for train_idx, val_idx in cv.split(X_train):
    X_cv_train, X_cv_val = X_train.iloc[train_idx], X_train.iloc[val_idx]
    y_cv_train, y_cv_val = y_train_class.iloc[train_idx], y_train_class.iloc[val_idx]
    stacking_clf.fit(X_cv_train, y_cv_train)
print("Cross-validation completed.\n")

# Section E: Time Series Forecasting Specifics
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense, Dropout
from tensorflow.keras.preprocessing.sequence import TimeseriesGenerator

# 14. Building the LSTM Model
print("Building LSTM model...")
lstm_model = Sequential([
    LSTM(50, activation='relu', input_shape=(12, 1)),
    Dropout(0.2),
    Dense(1)
])
lstm_model.compile(optimizer='adam', loss='mean_squared_error')

# Prepare data for LSTM
sequence_data = data['Intensity'].values.reshape(-1, 1)
generator = TimeseriesGenerator(sequence_data, sequence_data, length=12, batch_size=1)
lstm_model.fit(generator, epochs=20)
print("LSTM model trained.\n")

# 16. Model Forecasting
print("Forecasting with LSTM model...")
lstm_predictions = lstm_model.predict(generator)
plt.figure(figsize=(12, 6))
plt.plot(sequence_data, label='Actual')
plt.plot(np.arange(12, len(lstm_predictions) + 12), lstm_predictions, label='Predicted', color='red')
plt.title('LSTM Forecasting')
plt.legend()
plt.show()

# Section F: Post-Modeling Analysis
# 17. Residual Analysis
residuals = y_test_reg - y_pred_reg
sns.histplot(residuals, kde=True)
plt.title('Residual Analysis')
plt.show()

# 18. Comparison of Model Performance
print("Comparison of model performance completed.\n")

# Visualizations for Predictions
print("Generating visualizations...")

# 1. Actual vs. Predicted for regression
plt.figure(figsize=(10, 6))
plt.scatter(y_test_reg, y_pred_reg, alpha=0.6, edgecolor='k')
plt.plot([y_test_reg.min(), y_test_reg.max()], [y_test_reg.min(), y_test_reg.max()], 'r--')
plt.xlabel('Actual Economic Impact Index')
plt.ylabel('Predicted Economic Impact Index')
plt.title('Actual vs Predicted: Regression')
plt.grid(True)
plt.show()

# 2. Residual Plot for regression
residuals = y_test_reg - y_pred_reg
plt.figure(figsize=(10, 6))
sns.histplot(residuals, kde=True)
plt.title('Residual Analysis: Regression Model')
plt.xlabel('Residuals')
plt.ylabel('Frequency')
plt.show()

# 3. Feature Importance Plot for the best regressor (e.g., XGBoost)
xgb_feature_importance = pd.Series(best_reg.feature_importances_, index=X.columns)
xgb_feature_importance.sort_values().plot(kind='barh', figsize=(8, 6))
plt.title('Feature Importance: XGBoost Regressor')
plt.xlabel('Importance Score')
plt.show()

# 4. Classification Model Visualization
plt.figure(figsize=(10, 6))
plt.scatter(range(len(y_test_class)), y_test_class, label='Actual', alpha=0.6)
plt.scatter(range(len(y_pred_class)), y_pred_class, label='Predicted', alpha=0.6, color='red')
plt.title('Classification Model Predictions')
plt.xlabel('Sample Index')
plt.ylabel('Disaster Type Label')
plt.legend()
plt.grid(True)
plt.show()

# 5. Time Series Forecasting Visualizations
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense, Dropout
from tensorflow.keras.preprocessing.sequence import TimeseriesGenerator

print("Building and training LSTM model...")
lstm_model = Sequential([
    LSTM(50, activation='relu', input_shape=(12, 1)),
    Dropout(0.2),
    Dense(1)
])
lstm_model.compile(optimizer='adam', loss='mean_squared_error')

sequence_data = data['Intensity'].values.reshape(-1, 1)
generator = TimeseriesGenerator(sequence_data, sequence_data, length=12, batch_size=1)
lstm_model.fit(generator, epochs=20)

# 6. Forecasting Visualization
print("Forecasting with LSTM model...")
lstm_predictions = lstm_model.predict(generator)
plt.figure(figsize=(12, 6))
plt.plot(sequence_data, label='Actual')
plt.plot(np.arange(12, len(lstm_predictions) + 12), lstm_predictions, label='Predicted', color='red')
plt.fill_between(np.arange(12, len(lstm_predictions) + 12), lstm_predictions.flatten() - 0.1 * lstm_predictions.flatten(),
                 lstm_predictions.flatten() + 0.1 * lstm_predictions.flatten(), color='red', alpha=0.2)
plt.title('LSTM Time Series Forecast with Confidence Interval')
plt.legend()
plt.grid(True)
plt.show()

# Documentation and Reporting
print("Modeling process completed. Visualizations generated for insights.")