Skip to content

Commit

Permalink
Finish basic closure objective function generator. It is 40x faster t…
Browse files Browse the repository at this point in the history 
…han the previsou version
  • Loading branch information
@KilinW
KilinW committed Nov 16, 2023
1 parent 8f036c2 commit d8a03b5
Show file tree
Hide file tree
Showing 4 changed files with 112 additions and 135 deletions.
25 changes: 13 additions & 12 deletions examples/configuration.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,8 +74,8 @@ def month_start_end():
stats_month[month, scale, :] = [
model.mean(),
model.cvar(),
model.acf(),
model.skewness(),
model.acf(),
model.pDry(0),
]

Expand All @@ -89,8 +89,8 @@ def month_start_end():
stats_month_seperate[month, year, scale, :] = [
model.mean(),
model.cvar(),
model.acf(),
model.skewness(),
model.acf(),
model.pDry(0),
]

Expand All @@ -104,14 +104,14 @@ def month_start_end():
target = BLRPRxConfig.default_target([1, 3, 6, 24])
target[Stat_Props.MEAN] = target_np[:, 0]
target[Stat_Props.CVAR] = target_np[:, 1]
target[Stat_Props.AR1] = target_np[:, 2]
target[Stat_Props.SKEWNESS] = target_np[:, 3]
target[Stat_Props.SKEWNESS] = target_np[:, 2]
target[Stat_Props.AR1] = target_np[:, 3]

weight = BLRPRxConfig.default_weight([1, 3, 6, 24])
weight[Stat_Props.MEAN] = weight_np[:, 0]
weight[Stat_Props.CVAR] = weight_np[:, 1]
weight[Stat_Props.AR1] = weight_np[:, 2]
weight[Stat_Props.SKEWNESS] = weight_np[:, 3]
weight[Stat_Props.SKEWNESS] = weight_np[:, 2]
weight[Stat_Props.AR1] = weight_np[:, 3]

mask = BLRPRxConfig.default_mask([1, 3, 6, 24])
mask[Stat_Props.MEAN] = [1, 0, 0, 0]
Expand All @@ -120,14 +120,15 @@ def month_start_end():
mask[Stat_Props.SKEWNESS] = [1, 1, 1, 1]

config = BLRPRxConfig(target=target, weight=weight, mask=mask)
a = config.get_evaluation_func()
arr = np.array([0.016679733103341976, 0.08270236178820184, 0.34970877070925505, 9.017352714561754, 0.9931496975448589, 1.01, 0.971862948182735])
obj_func = config.get_evaluation_func()
arr = np.array([0.016679733103341976, 0.08270236178820184, 0.34970877070925505, 9.017352714561754, 0.9931496975448589, 1.0, 0.971862948182735])
params = BLRPRx_params(*arr)

old_config = BLRPRxFitter()
model = BLRPRx()

print(a(arr))
print(old_config._evaluate(arr, target_np, weight_np, model))
print(obj_func(arr))
print(old_config.evaluate(params, target_np, weight_np, model))

#print(timeit.timeit(lambda: a(arr), number=389253))
#print(timeit.timeit(lambda: old_config._evaluate(arr, target_np, weight_np, model), number=389253))
print(timeit.timeit(lambda: obj_func(arr), number=100000))
print(timeit.timeit(lambda: old_config._evaluate(arr, target_np, weight_np, model), number=100000))
107 changes: 54 additions & 53 deletions legacy/compare_obj.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,61 +4,51 @@
import pandas as pd

from pyBL.fitting.fitter import BLRPRxFitter
from pyBL.models import Stat_Propsps, BLRPRx, BLRPRx_params
from pyBL.models import Stat_Props, BLRPRx, BLRPRx_params
from pyBL.timeseries import IntensityMRLE

timescale = [1, 3600, 3 * 3600, 6 * 3600, 24 * 3600]
props = [Stat_PropsStat_PropsStat_PropsStat_PropsStat_Props
Stat_Props.MEAN,
Stat_Props.CVAR,
Stat_Props.AR1,
Stat_Props.SKEWNESS,
Stat_Props.pDRY,
]
from Library.Fitting_lib.objectiveFunction import Exponential_func
from Library.BLRPRmodel.BLRPRx import BLRPRx as BLRPRx_legacy

# Set timezone to UTC
os.environ["TZ"] = "UTC"
timescale = [1, 3600, 3*3600, 6*3600, 24*3600]
props = [Stat_Props.MEAN, Stat_Props.CVAR, Stat_Props.AR1, Stat_Props.SKEWNESS, Stat_Props.pDRY]

# Set timezone to UTC
os.environ['TZ'] = 'UTC'

def rain_timeseries():
data_path = os.path.join(os.path.dirname(__file__), "01 Input_data", "elmdon.csv")
data_path = os.path.join(os.path.dirname(__file__), '01 Input_data', 'elmdon.csv')
data = pd.read_csv(data_path, parse_dates=["datatime"])
data["datatime"] = data["datatime"].astype("int64") // 10**9
time = data["datatime"].to_numpy()
intensity = data["Elmdon"].to_numpy()
data['datatime'] = data['datatime'].astype("int64") // 10 ** 9
time = data['datatime'].to_numpy()
intensity = data['Elmdon'].to_numpy()
return time, intensity


def month_start_end():
# Generate first day of each month from 1980 to 2010
day = pd.date_range(start="1980-01-01", end="2000-01-01", freq="MS")
day = pd.date_range(start='1980-01-01', end='2000-01-01', freq='MS')
# Convert to unix time
month_srt = day.astype("int64") // 10**9
month_srt = day.astype("int64") // 10 ** 9
month_end = month_srt
# Stack them together
month_interval = np.stack((month_srt[:-1], month_end[1:]), axis=1)
# Group the month_interval by month
month_interval = np.reshape(month_interval, (-1, 12, 2))
return month_interval


time, intensity = rain_timeseries()
mrle = IntensityMRLE(time, intensity / 3600)
mrle = IntensityMRLE(time, intensity/3600)
month_interval_each_year = month_start_end()

# Segment the mrle timeseries into months from 1900 to 2100
mrle_month_each = np.empty(
(12, len(month_interval_each_year), 5), dtype=IntensityMRLE
) # (month, year, scale)
mrle_month_each = np.empty((12, len(month_interval_each_year), 5), dtype=IntensityMRLE) # (month, year, scale)
for i, year in enumerate(month_interval_each_year):
for j, month in enumerate(year):
for k, scale in enumerate(timescale):
mrle_month_each[j, i, k] = mrle[month[0] : month[1]].rescale(
scale
) # 1s scale
mrle_month_each[j, i, k] = mrle[month[0]:month[1]].rescale(scale) # 1s scale

# MRLE that stores the total of each month
mrle_month_total = np.empty((12, 5), dtype=IntensityMRLE) # (month, scale)
mrle_month_total = np.empty((12, 5), dtype=IntensityMRLE) # (month, scale)
for i in range(12):
for j in range(len(mrle_month_each[0])):
for k, scale in enumerate(timescale):
Expand All @@ -70,42 +60,53 @@ def month_start_end():
for month in range(12):
for scale in range(5):
model = mrle_month_total[month, scale]
stats_month[month, scale, :] = [
model.mean(),
model.cvar(),
model.acf(),
model.skewness(),
model.pDry(0),
]

stats_month_seperate = np.zeros(
(12, len(month_interval_each_year), 5, 5)
) # (month, year, scale, stats)
stats_month[month, scale, :] = [model.mean(), model.cvar(), model.acf(), model.skewness(), model.pDry(0)]

stats_month_seperate = np.zeros((12, len(month_interval_each_year), 5, 5)) # (month, year, scale, stats)
for month in range(12):
for year in range(len(month_interval_each_year)):
for scale in range(5):
model = mrle_month_each[month, year, scale]
stats_month_seperate[month, year, scale, :] = [
model.mean(),
model.cvar(),
model.acf(),
model.skewness(),
model.pDry(0),
]

stats_weight = 1 / np.nanvar(
stats_month_seperate, axis=1
) # (month, scale, stats) (12, 5, 5)
stats_month_seperate[month, year, scale, :] = [model.mean(), model.cvar(), model.acf(), model.skewness(), model.pDry(0)]

stats_weight = 1/np.nanvar(stats_month_seperate, axis=1) # (month, scale, stats) (12, 5, 5)

target = stats_month[:, 1:, :]
weight = stats_weight[:, 1:, :]


# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

params = BLRPRx_params(
lambda_=0.01, phi=0.1, kappa=0.1, alpha=2.1, nu=0.3, sigmax_mux=1, iota=0.1
)
params = BLRPRx_params(lambda_=0.016679733103341976,
phi=0.08270236178820184,
kappa=0.34970877070925505,
alpha=9.017352714561754,
nu=0.9931496975448589,
sigmax_mux=1.0,
iota=0.971862948182735)
legacy_params = [params.lambda_, params.iota, params.sigmax_mux, 0, params.alpha, params.alpha/params.nu, params.kappa, params.phi, 1.0]
fitter = BLRPRxFitter()
model = BLRPRx(params=params)

print(fitter.evaluate(x=model.params(), model=model, target=target, weight=weight))
new_score = fitter.evaluate(x=model.get_params(), model=model, target=target[0], weight=weight[0])
legacy_target = [target[0, 0, 0], target[0, 0, 1], target[0, 0, 2], target[0, 0, 2],
target[0, 1, 1], target[0, 1, 3], target[0, 1, 2],
target[0, 2, 1], target[0, 2, 3], target[0, 2, 2],
target[0, 3, 1], target[0, 3, 3], target[0, 3, 2]]
legacy_weight = [weight[0, 0, 0], weight[0, 0, 1], weight[0, 0, 3], weight[0, 0, 2],
weight[0, 1, 1], weight[0, 1, 3], weight[0, 1, 2],
weight[0, 2, 1], weight[0, 2, 3], weight[0, 2, 2],
weight[0, 3, 1], weight[0, 3, 3], weight[0, 3, 2]]

legacy_target_cor = [target[0, 0, 0], target[0, 0, 1], target[0, 0, 2], target[0, 0, 3],
target[0, 1, 1], target[0, 1, 2], target[0, 1, 3],
target[0, 2, 1], target[0, 2, 2], target[0, 2, 3],
target[0, 3, 1], target[0, 3, 2], target[0, 3, 3]]
legacy_weight_cor = [weight[0, 0, 0], weight[0, 0, 1], weight[0, 0, 2], weight[0, 0, 3],
weight[0, 1, 1], weight[0, 1, 2], weight[0, 1, 3],
weight[0, 2, 1], weight[0, 2, 2], weight[0, 2, 3],
weight[0, 3, 1], weight[0, 3, 2], weight[0, 3, 3]]
old_score = Exponential_func(legacy_params, 0, [1, 3, 6, 24], legacy_target, legacy_weight, model=BLRPRx_legacy(mode='exponential'))
cor_score = Exponential_func(legacy_params, 0, [1, 3, 6, 24], legacy_target_cor, legacy_weight_cor, model=BLRPRx_legacy(mode='exponential'))
print(new_score, old_score, cor_score)
print('>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>')
102 changes: 40 additions & 62 deletions src/pyBL/fitting/fitter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -271,7 +271,6 @@ def _evaluate(
for prop in enabled_props:
if self._mask[i, prop.value]:
result[i, prop.value] = model.get_prop(prop, ts)
print(result)
diff = np.sum(np.power(result - target, 2) * weight)
return diff

Expand Down Expand Up @@ -330,74 +329,53 @@ def get_evaluation_func(self) -> Callable[[npt.NDArray], np.float64]:

f1_func = self._rci_model.get_f1
f2_func = self._rci_model.get_f2

@nb.njit
def evaluation_func(x: npt.NDArray[np.float64]) -> np.float64:
predict = np.zeros((len(scales), stats_len), dtype=np.float64)
lambda_, phi, kappa, alpha, nu, sigmax_mux, iota = x
f1 = f1_func(sigmax_mux)
f2 = f2_func(sigmax_mux)
for t, scale in enumerate(scales):
mean = None
variance = None
moment_3rd = None
for stat in range(stats_len):
if not mask_np[t, stat]:
continue
if stat == 0: # Stat_Props.MEAN.value
if mean is None:
mean = _blrprx_mean(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota)
predict[t, stat] = mean
elif stat == 1: # Stat_Props.CVAR.value
if mean is None:
mean = _blrprx_mean(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota)
if variance is None:
variance = _blrprx_variance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1)
predict[t, stat] = np.sqrt(variance) / mean
elif stat == 2: # Stat_Props.SKEWNESS.value
if moment_3rd is None:
moment_3rd = _blrprx_moment_3rd(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1, f2)
if variance is None:
variance = _blrprx_variance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1)
predict[t, stat] = moment_3rd / np.power(variance, 1.5)
elif stat == 3: # Stat_Props.AR1.value
if variance is None:
variance = _blrprx_variance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1)
predict[t, stat] = _blrprx_covariance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1, 1.0) / variance
elif stat == 4: # Stat_Props.AR2.value
if variance is None:
variance = _blrprx_variance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1)
predict[t, stat] = _blrprx_covariance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1, 2.0) / variance
elif stat == 5: # Stat_Props.AR3.value
if variance is None:
variance = _blrprx_variance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1)
predict[t, stat] = _blrprx_covariance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1, 3.0) / variance
elif stat == 6: # Stat_Props.pDRY.value
predict[t, stat] = 0
elif stat == 7: # Stat_Props.MSIT.value
predict[t, stat] = 0
elif stat == 8: # Stat_Props.MSD.value
predict[t, stat] = 0
elif stat == 9: # Stat_Props.MCIT.value
predict[t, stat] = 0
elif stat == 10: # Stat_Props.MCD.value
predict[t, stat] = 0
elif stat == 11: # Stat_Props.MCS.value
predict[t, stat] = 0
elif stat == 12: # Stat_Props.MPC.value
predict[t, stat] = 0
elif stat == 13: # Stat_Props.VAR.value
if variance is None:
variance = _blrprx_variance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1)
predict[t, stat] = variance
elif stat == 14: # Stat_Props.AC1.value
predict[t, stat] = _blrprx_covariance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1, 1.0)
elif stat == 15: # Stat_Props.AC2.value
predict[t, stat] = _blrprx_covariance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1, 2.0)
elif stat == 16: # Stat_Props.AC3.value
predict[t, stat] = _blrprx_covariance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1, 3.0)

print(predict)
mean = _blrprx_mean(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota)
variance = _blrprx_variance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1)
moment_3rd = _blrprx_moment_3rd(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1, f2)
#if not mask_np[t, stat]:
#continue
if mask_np[t, 0]: # Stat_Props.MEAN.value
predict[t, 0] = mean
if mask_np[t,1]: # Stat_Props.CVAR.value
predict[t, 1] = np.sqrt(variance) / mean
if mask_np[t,2]: # Stat_Props.SKEWNESS.value
predict[t, 2] = moment_3rd / np.power(variance, 1.5)
if mask_np[t,3]: # Stat_Props.AR1.value
predict[t, 3] = _blrprx_covariance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1, 1.0) / variance
if mask_np[t,4]: # Stat_Props.AR2.value
predict[t, 4] = _blrprx_covariance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1, 2.0) / variance
if mask_np[t,5]: # Stat_Props.AR3.value
predict[t, 5] = _blrprx_covariance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1, 3.0) / variance
if mask_np[t,6]: # Stat_Props.pDRY.value
predict[t, 6] = 0
if mask_np[t,7]: # Stat_Props.MSIT.value
predict[t, 7] = 0
if mask_np[t,8]: # Stat_Props.MSD.value
predict[t, 8] = 0
if mask_np[t,9]: # Stat_Props.MCIT.value
predict[t, 9] = 0
if mask_np[t,10]: # Stat_Props.MCD.value
predict[t, 10] = 0
if mask_np[t,11]: # Stat_Props.MCS.value
predict[t, 11] = 0
if mask_np[t,12]: # Stat_Props.MPC.value
predict[t, 12] = 0
if mask_np[t,13]: # Stat_Props.VAR.value
predict[t, 13] = variance
if mask_np[t,14]: # Stat_Props.AC1.value
predict[t, 14] = _blrprx_covariance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1, 1.0)
if mask_np[t,15]: # Stat_Props.AC2.value
predict[t, 15] = _blrprx_covariance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1, 2.0)
if mask_np[t,16]: # Stat_Props.AC3.value
predict[t, 16] = _blrprx_covariance(scale, lambda_, phi, kappa, alpha, nu, sigmax_mux, iota, f1, 3.0)
return np.nansum(np.power(predict - target_np, 2) * weight_np * mask_np)
return evaluation_func

Expand Down
13 changes: 5 additions & 8 deletions src/pyBL/models/blrprx.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -245,6 +245,10 @@ def sample(self, duration_hr: float) -> IntensityMRLE:

@nb.njit
def _blrprx_kernel(k: float, u: float, nu: float, alpha: float) -> float:
'''
k: lag
u: timescale
'''
# Modelling rainfall with a Bartlett–Lewis process: new developments(2020) Formula (5)

## TODO: Check if this is still required.
Expand Down Expand Up @@ -360,14 +364,7 @@ def _blrprx_moment_3rd(
):
mu_c = 1.0 + kappa / phi

phi2 = phi**2
phi3 = phi**3
phi4 = phi**4
phi5 = phi**5
phi6 = phi**6
phi7 = phi**7
phi8 = phi**8
phi9 = phi**9
phi2, phi3, phi4, phi5, phi6, phi7, phi8, phi9 = np.power(phi, np.arange(2, 10))

kappa2 = kappa**2

Expand Down

0 comments on commit d8a03b5

Please sign in to comment.