import numpy as np

import matplotlib.pyplot as plt

from scipy.stats import spearmanr
import pandas as pd
from scipy import stats
import scipy.optimize as opt
from statsmodels.distributions.empirical_distribution import ECDF
from scipy.optimize import bisect

def empirical_cdf_inverse(data, u):

"""
Compute the inverse empirical CDF (quantile function).
"""
sorted_data = np.sort(data)
quantile_index = np.clip(int(u * (len(sorted_data)-1)), 0, len(sorted_data) -
1)
return sorted_data[quantile_index]

def gaussian_copula_cdf_inverse(u, rho):

"""
Inverse of the partial derivative of the Gaussian copula CDF.
Uses numerical root-finding to approximate the inverse.
"""

def func(x):
return stats.norm.cdf((x - rho * stats.norm.ppf(u)) / np.sqrt(1 - rho **
2)) - u

return bisect(func, -5, 5) # Root finding in reasonable range

def simulate_xt(omega_xt, p01, p11, rho, empirical_data):

"""
Simulates x_t given x_t-1 using a Gaussian copula with empirical marginals.
"""
# Compute the transformed uniform variable
u = ((p01 + p11) * omega_xt - p01) / p11

# Inverse of copula partial derivative

u_cond = gaussian_copula_cdf_inverse(u, rho)

# Transform back using inverse empirical CDF

x_t = empirical_cdf_inverse(empirical_data, u_cond)

return x_t

def gumbel_copula_cdf_inverse(u, theta):

"""
Computes the inverse of the Gumbel copula's conditional CDF using root-finding.
"""

def func(v):
term1 = (-np.log(u)) ** theta
term2 = (-np.log(v)) ** theta
C_uv = np.exp(- (term1 + term2) ** (1 / theta))
print(u,v)
print('***')
 return C_uv - u # Find root where this equals zero

return bisect(func, 1e-6, 1 - 1e-6) # Numerical root-finding in (0,1) range

def simulate_xt_gumbel(x_t_minus_1, omega_xt, p01, p11, theta, empirical_data):

"""
Simulates x_t given x_t-1 using a Gumbel copula with empirical marginals.
"""
# Compute the transformed uniform variable
u = ((p01 + p11) * omega_xt - p01) / p11

# Ensure u is within valid probability range

u = np.clip(u, 1e-6, 1 - 1e-6)

# Inverse of Gumbel copula partial derivative

u_cond = gumbel_copula_cdf_inverse(u, theta)

# Transform back using inverse empirical CDF

x_t = empirical_cdf_inverse(empirical_data, u_cond)

return x_t

def compute_xt_empirical(p01, p11, omega_xt, theta, empirical_data):

"""
Computes X_t from X_{t-1} using the Gumbel copula and empirical marginals.
"""
U= ((p01 + p11) * omega_xt - p01) / p11
print(U,theta)
U_t = gumbel_copula_cdf_inverse(U, theta)
print(U_t)
X_t = empirical_cdf_inverse(empirical_data, U_t)

return X_t

# Gumbel copula CDF

def gumbel_copula_cdf(u, v, theta):
if theta == 1:
return u * v # independent case
inner_term = ((-np.log(u))**theta + (-np.log(v))**theta)**(1/theta)
return np.exp(-inner_term)

# Partial derivative of Gumbel copula with respect to u

def gumbel_copula_partial_u(u, v, theta):
if theta == 1:
return v # independent case
inner_term = (-np.log(u))**theta + (-np.log(v))**theta
c_uv = gumbel_copula_cdf(u, v, theta)
return c_uv * ((-np.log(u))**(theta-1)) / (u * inner_term**(1 - 1/theta))

# Conditional CDF using Gumbel copula for a given x (if x > 0) and rainfall value y
def O_Y_given_X_x_and_X_positive(y, x):
u, v = H_x(x), H_y(y)
c_uv_partial = gumbel_copula_partial_u(u, v, theta)
h_x = np.count_nonzero(rainfall_data[:-1] == x) / len(rainfall_data[:-1]) if
len(rainfall_data[:-1]) > 0 else 0
f_y = np.count_nonzero(rainfall_data[1:] == y) / len(rainfall_data[1:]) if
len(rainfall_data[1:]) > 0 else 0
denom = P_10 * h_x + P_11 * f_y
return (P_10 * h_x + P_11 * f_y * c_uv_partial) / denom if denom > 0 else
H_y(y)

# Conditional CDF for X = 0 (when previous day is dry)

def O_Y_given_X_0(y):
denom = P_00 + P_01
return (P_00 + P_01 * H_y(y)) / denom if denom > 0 else H_y(y)

# Define the Gumbel Copula function

def gumbel_copula(u, v, theta):
term = (-np.log(u)) ** theta + (-np.log(v)) ** theta
return np.exp(-term ** (1 / theta))

# Partial derivative with respect to u

def gumbel_partial_derivative_u(u, v, theta):
C_uv = gumbel_copula(u, v, theta)
term1 = (-np.log(u)) ** (theta - 1) / u
term2 = (-np.log(u)) ** theta + (-np.log(v)) ** theta
return C_uv * (term1 / term2 ** (1 - 1 / theta))

# Inverse function: Solve for u given v and the partial derivative value
def inverse_gumbel_partial(du, v, theta):
"""Compute inverse of Gumbel copula's partial derivative."""
if du <= 0 or du >= 1 or v <= 0: # Check for invalid values
return np.nan # Return NaN so we can detect failure

term = (-np.log(du)) ** theta + (-np.log(v)) ** theta

if term <= 0: # Ensure valid computation
return np.nan
return np.exp(-term ** (1 / theta))

# Simulation function
def simulate_markov_chain(T, p00, p01, p11,p10):

xt =np.zeros(T)# Store results

o1= np.zeros(T)
pt1= np.zeros(T)
# # st=[]
# st = np.zeros(T)
for t in range(1, T - 1):
while True: # Retry mechanism
O = np.random.uniform() # Generate new random value
# O = np.random.uniform(0,1)
if xt[t - 1] == 0 and O > (p00 / (p00 + p01)):
p1 = ((p00 + p01) * O - p00) / p01
xt[t] = empirical_cdf_inverse(rainfall_data, p1)
o1[t] = O
pt1[t] = p1
st = "*****"
break # Exit while loop
 elif xt[t - 1] > 0 and O > (p01 / (p01 + p11)):
du = np.clip(((p01 + p11) * O - p01) / p11, 1e-5, 1 - 1e-5)
v = xt[t - 1]
u_inverse = inverse_gumbel_partial(du, v, theta)

if np.isnan(u_inverse): # Check for invalid copula output

continue # Retry loop with a new O

xt[t] = empirical_cdf_inverse(rainfall_data, u_inverse)

o1[t] = O
pt1[t] = u_inverse
st = "**1**"
break # Exit while loop

elif xt[t - 1] > 0 and O > (p10 / (p10 + p00)):

p1 = 1*((p00 + p10) * O - p10) / p00
# xt[t] = empirical_cdf_inverse(rainfall_data, p1)
xt[t]=0
o1[t] = O
pt1[t] =p1
st = "**2**"
break # Exit while loop

else:
xt[t] = 0
o1[t] = O
pt1[t] = 0
st = "**3**"
break # Exit while loop

print(t,"{:.3f}".format(rainfall_data[t]),"{:.3f}".format(xt[t]),"{:.3f}".format(xt
[t-1]),"{:.3f}".format(pt1[t]),"{:.3f}".format(o1[t]),st)

return xt

df1 = pd.read_csv('D:\\pythonProject1\\venv\\Rain_imd\\IMD_DataSET\\annual.csv')
rainfall_data = df1.iloc[3:19, 7].values # Example selection of data

print(rainfall_data)

n = len(rainfall_data) - 1

dry_prev = rainfall_data[:-1] == 0
wet_prev = rainfall_data[:-1] > 0
dry_curr = rainfall_data[1:] == 0
wet_curr = rainfall_data[1:] > 0

P_00 = np.sum(dry_prev & dry_curr) / n

P_01 = np.sum(dry_prev & wet_curr) / n
P_10 = np.sum(wet_prev & dry_curr) / n
P_11 = np.sum(wet_prev & wet_curr) / n
print("\nTransition Probabilities:")
print(f"P_00 (Dry-Dry): {P_00:.3f}")
print(f"P_01 (Dry-Wet): {P_01:.3f}")
print(f"P_10 (Wet-Dry): {P_10:.3f}")
print(f"P_11 (Wet-Wet): {P_11:.3f}")

H_x= np.arange(1, len(rainfall_data[1:]) + 1) / len(rainfall_data[1:])

H_y= np.arange(1, len(rainfall_data[:-1]) + 1) / len(rainfall_data[:-1])

wet_indices = (rainfall_data[:-1] > 0) & (rainfall_data[1:] > 0)

wet_y_data = rainfall_data[:-1][wet_indices]
wet_x_data = rainfall_data[1:][wet_indices]

rho, _ = spearmanr(wet_x_data, wet_y_data)

theta = 1 / (1 - rho) # Gumbel copula parameter

print(n,theta)

T = len(rainfall_data)

xt_values = simulate_markov_chain(T, P_00, P_01, P_11,P_10)

simulated_rainfall=xt_values

n = len(simulated_rainfall) - 1
dry_prev = simulated_rainfall[:-1] == 0
wet_prev = simulated_rainfall[:-1] > 0
dry_curr = simulated_rainfall[1:] == 0
wet_curr = simulated_rainfall[1:] > 0

P_00 = np.sum(dry_prev & dry_curr) / n

P_01 = np.sum(dry_prev & wet_curr) / n
P_10 = np.sum(wet_prev & dry_curr) / n
P_11 = np.sum(wet_prev & wet_curr) / n

print("\nTransition Probabilities for simulated rainfall:")

print(f"P_00 (Dry-Dry): {P_00:.3f}")
print(f"P_01 (Dry-Wet): {P_01:.3f}")
print(f"P_10 (Wet-Dry): {P_10:.3f}")
print(f"P_11 (Wet-Wet): {P_11:.3f}")

# y=0 # y is previous day rain i.e. xt-1

# if y==0:
# O1=(np.random.uniform())
# print(O1)
# f1=P_00/(P_00+P_01)
# if O1>f1:
# p1 = ((P_00+P_01)*O1-P_00)/P_01
# x = inverse_ecdf(rainfall_data[1:], p1) # x is todays rain
# else:
# x =0
# y=x
# elif y>0:
# O2 = (np.random.uniform())
# # O2=0.23
# # print(O2)
# f2 = P_01 / (P_01 + P_11)
# # print(f2)
# if O2 > f2:
# p2 = ((P_01 + P_11) * O2 - P_01) / P_11
# U=(np.random.uniform())
# p3 = gumbel_copula_inverse(U, theta, p2)
# # print(p3)
# x = inverse_ecdf(rainfall_data[1:], p3) # x is todays rain
# else:
# x = 0
# y = x

# # print(simulated_rainfall)
# plt.figure(figsize=(10, 6))
# plt.scatter(rainfall_data, simulated_rainfall)
# print(rainfall_data,simulated_rainfall)
plt.figure(figsize=(10, 6))
plt.plot(rainfall_data, label="Observed Rainfall", marker="o", linestyle="-",
alpha=0.7)
plt.plot(simulated_rainfall, label="Simulated Rainfall", marker="x",
linestyle="--", alpha=0.7)
plt.xlabel("Time (days)")
plt.ylabel("Rainfall (mm)")
plt.title("Observed vs Simulated Rainfall")
plt.legend()
plt.grid(True)
plt.show()