from AlgorithmImports import *
import numbers
from math import log, ceil
import pandas as pd
import numpy as np
from sklearn.covariance import OAS

# ==============================================================================
# MLFINLAB: CRITICAL LINE ALGORITHM (CLA) CLASSES
# ==============================================================================

class CLA:
    def __init__(self, weight_bounds=(0, 1), calculate_returns="mean"):
        self.weight_bounds = weight_bounds
        self.calculate_returns = calculate_returns
        self.weights = list()
        self.lambdas = list()
        self.gammas = list()
        self.free_weights = list()
        self.expected_returns = None
        self.cov_matrix = None
        self.lower_bounds = None
        self.upper_bounds = None
        self.max_sharpe = None
        self.min_var = None
        self.efficient_frontier_means = None
        self.efficient_frontier_sigma = None

    @staticmethod
    def _infnone(number):
        return float("-inf") if number is None else number

    def _init_algo(self):
        structured_array = np.zeros((self.expected_returns.shape[0]), dtype=[("id", int), ("mu", float)])
        expected_returns = [self.expected_returns[i][0] for i in range(self.expected_returns.shape[0])]
        structured_array[:] = list(zip(list(range(self.expected_returns.shape[0])), expected_returns))
        expected_returns = np.sort(structured_array, order="mu")
        
        index, weights = expected_returns.shape[0], np.copy(self.lower_bounds)
        while np.sum(weights) < 1:
            index -= 1
            weights[expected_returns[index][0]] = self.upper_bounds[expected_returns[index][0]]
        weights[expected_returns[index][0]] += 1 - np.sum(weights)
        return [expected_returns[index][0]], weights

    @staticmethod
    def _compute_bi(c_final, asset_bounds_i):
        if c_final > 0:
            return asset_bounds_i[1][0]
        return asset_bounds_i[0][0]

    def _compute_w(self, covar_f_inv, covar_fb, mean_f, w_b):
        ones_f = np.ones(mean_f.shape)
        g_1 = np.dot(np.dot(ones_f.T, covar_f_inv), mean_f)
        g_2 = np.dot(np.dot(ones_f.T, covar_f_inv), ones_f)
        if w_b is None:
            g_final, w_1 = float(-self.lambdas[-1] * g_1 / g_2 + 1 / g_2), 0
        else:
            ones_b = np.ones(w_b.shape)
            g_3 = np.dot(ones_b.T, w_b)
            g_4 = np.dot(covar_f_inv, covar_fb)
            w_1 = np.dot(g_4, w_b)
            g_4 = np.dot(ones_f.T, w_1)
            g_final = float(-self.lambdas[-1] * g_1 / g_2 + (1 - g_3 + g_4) / g_2)

        w_2 = np.dot(covar_f_inv, ones_f)
        w_3 = np.dot(covar_f_inv, mean_f)
        free_asset_weights = -1*w_1 + g_final * w_2 + self.lambdas[-1] * w_3
        return free_asset_weights, g_final

    def _compute_lambda(self, covar_f_inv, covar_fb, mean_f, w_b, asset_index, b_i):
        ones_f = np.ones(mean_f.shape)
        c_1 = np.dot(np.dot(ones_f.T, covar_f_inv), ones_f)
        c_2 = np.dot(covar_f_inv, mean_f)
        c_3 = np.dot(np.dot(ones_f.T, covar_f_inv), mean_f)
        c_4 = np.dot(covar_f_inv, ones_f)
        c_final = -1*c_1 * c_2[asset_index] + c_3 * c_4[asset_index]
        if c_final == 0:
            return None, None

        if isinstance(b_i, list):
            b_i = self._compute_bi(c_final, b_i)

        if w_b is None:
            return float((c_4[asset_index] - c_1 * b_i) / c_final), b_i

        ones_b = np.ones(w_b.shape)
        l_1 = np.dot(ones_b.T, w_b)
        l_2 = np.dot(covar_f_inv, covar_fb)
        l_3 = np.dot(l_2, w_b)
        l_2 = np.dot(ones_f.T, l_3)
        lambda_value = float(((1 - l_1 + l_2) * c_4[asset_index] - c_1 * (b_i + l_3[asset_index])) / c_final)
        return lambda_value, b_i

    def _get_matrices(self, free_weights):
        covar_f = self._reduce_matrix(self.cov_matrix, free_weights, free_weights)
        mean_f = self._reduce_matrix(self.expected_returns, free_weights, [0])
        bounded_weights = self._get_bounded_weights(free_weights)
        covar_fb = self._reduce_matrix(self.cov_matrix, free_weights, bounded_weights)
        w_b = self._reduce_matrix(self.weights[-1], bounded_weights, [0])
        return covar_f, covar_fb, mean_f, w_b

    def _get_bounded_weights(self, free_weights):
        return self._diff_lists(list(range(self.expected_returns.shape[0])), free_weights)

    @staticmethod
    def _diff_lists(list_1, list_2):
        return list(set(list_1) - set(list_2))

    @staticmethod
    def _reduce_matrix(matrix, row_indices, col_indices):
        return matrix[np.ix_(row_indices, col_indices)]

    def _purge_num_err(self, tol):
        index_1 = 0
        while True:
            flag = False
            if index_1 == len(self.weights):
                break
            if abs(sum(self.weights[index_1]) - 1) > tol:
                flag = True
            else:
                for index_2 in range(self.weights[index_1].shape[0]):
                    if (self.weights[index_1][index_2] - self.lower_bounds[index_2] < -tol
                            or self.weights[index_1][index_2] - self.upper_bounds[index_2] > tol):
                        flag = True
                        break
            if flag is True:
                del self.weights[index_1]
                del self.lambdas[index_1]
                del self.gammas[index_1]
                del self.free_weights[index_1]
            else:
                index_1 += 1

    def _purge_excess(self):
        index_1, repeat = 0, False
        while True:
            if repeat is False:
                index_1 += 1
            if index_1 >= len(self.weights) - 1:
                break
            weights = self.weights[index_1]
            mean = np.dot(weights.T, self.expected_returns)[0, 0]
            index_2, repeat = index_1 + 1, False
            while True:
                if index_2 == len(self.weights):
                    break
                weights = self.weights[index_2]
                mean_ = np.dot(weights.T, self.expected_returns)[0, 0]
                if mean < mean_:
                    del self.weights[index_1]
                    del self.lambdas[index_1]
                    del self.gammas[index_1]
                    del self.free_weights[index_1]
                    repeat = True
                    break
                index_2 += 1

    @staticmethod
    def _golden_section(obj, left, right, **kwargs):
        tol, sign, args = 1.0e-9, -1, None
        args = kwargs.get("args", None)
        num_iterations = int(ceil(-2.078087 * log(tol / abs(right - left))))
        gs_ratio = 0.618033989
        complementary_gs_ratio = 1.0 - gs_ratio

        x_1 = gs_ratio * left + complementary_gs_ratio * right
        x_2 = complementary_gs_ratio * left + gs_ratio * right
        f_1 = sign * obj(x_1, *args)
        f_2 = sign * obj(x_2, *args)

        for _ in range(num_iterations):
            if f_1 > f_2:
                left = x_1
                x_1 = x_2
                f_1 = f_2
                x_2 = complementary_gs_ratio * left + gs_ratio * right
                f_2 = sign * obj(x_2, *args)
            else:
                right = x_2
                x_2 = x_1
                f_2 = f_1
                x_1 = gs_ratio * left + complementary_gs_ratio * right
                f_1 = sign * obj(x_1, *args)

        if f_1 < f_2:
            return x_1, sign * f_1
        return x_2, sign * f_2

    def _eval_sr(self, alpha, w_0, w_1):
        weights = alpha * w_0 + (1 - alpha) * w_1
        returns = np.dot(weights.T, self.expected_returns)[0, 0]
        volatility = np.dot(np.dot(weights.T, self.cov_matrix), weights)[0, 0] ** 0.5
        return returns / volatility

    def _bound_free_weight(self, free_weights):
        lambda_in = None
        i_in = None
        bi_in = None
        if len(free_weights) > 1:
            covar_f, covar_fb, mean_f, w_b = self._get_matrices(free_weights)
            covar_f_inv = np.linalg.inv(covar_f)
            j = 0
            for i in free_weights:
                lambda_i, b_i = self._compute_lambda(
                    covar_f_inv, covar_fb, mean_f, w_b, j, [self.lower_bounds[i], self.upper_bounds[i]]
                )
                if self._infnone(lambda_i) > self._infnone(lambda_in):
                    lambda_in, i_in, bi_in = lambda_i, i, b_i
                j += 1
        return lambda_in, i_in, bi_in

    def _free_bound_weight(self, free_weights):
        lambda_out = None
        i_out = None
        if len(free_weights) < self.expected_returns.shape[0]:
            bounded_weight_indices = self._get_bounded_weights(free_weights)
            for i in bounded_weight_indices:
                covar_f, covar_fb, mean_f, w_b = self._get_matrices(free_weights + [i])
                covar_f_inv = np.linalg.inv(covar_f)
                lambda_i, _ = self._compute_lambda(
                    covar_f_inv, covar_fb, mean_f, w_b, mean_f.shape[0] - 1, self.weights[-1][i]
                )
                if (self.lambdas[-1] is None or lambda_i < self.lambdas[-1]) and lambda_i > self._infnone(lambda_out):
                    lambda_out, i_out = lambda_i, i
        return lambda_out, i_out

    def _initialise(self, asset_prices, resample_by):
        pass

    @staticmethod
    def _calculate_mean_historical_returns(asset_prices, frequency=252):
        returns = asset_prices.pct_change().dropna(how="all")
        returns = returns.mean() * frequency
        return returns

    @staticmethod
    def _calculate_exponential_historical_returns(asset_prices, frequency=252, span=500):
        returns = asset_prices.pct_change().dropna(how="all")
        returns = returns.ewm(span=span).mean().iloc[-1] * frequency
        return returns

    def allocate(self, asset_prices, solution="cla_turning_points", resample_by="B"):
        pass

    def _compute_solution(self, assets, solution):
        if solution == "max_sharpe":
            self.max_sharpe, self.weights = self._max_sharpe()
            self.weights = pd.DataFrame(self.weights)
            self.weights.index = assets
            self.weights = self.weights.T
        elif solution == "min_volatility":
            self.min_var, self.weights = self._min_volatility()
            self.weights = pd.DataFrame(self.weights)
            self.weights.index = assets
            self.weights = self.weights.T
        elif solution == "efficient_frontier":
            self.efficient_frontier_means, self.efficient_frontier_sigma, self.weights = self._efficient_frontier()
            weights_copy = self.weights.copy()
            for i, turning_point in enumerate(weights_copy):
                self.weights[i] = turning_point.reshape(1, -1)[0]
            self.weights = pd.DataFrame(self.weights, columns=assets)
        elif solution == "cla_turning_points":
            weights_copy = self.weights.copy()
            for i, turning_point in enumerate(weights_copy):
                self.weights[i] = turning_point.reshape(1, -1)[0]
            self.weights = pd.DataFrame(self.weights, columns=assets)
        else:
            raise ValueError("Unknown solution string specified.")

    def _max_sharpe(self):
        w_sr, sharpe_ratios = [], []
        for i in range(len(self.weights) - 1):
            w_0 = np.copy(self.weights[i])
            w_1 = np.copy(self.weights[i + 1])
            kwargs = {"minimum": False, "args": (w_0, w_1)}
            alpha, sharpe_ratio = self._golden_section(self._eval_sr, 0, 1, **kwargs)
            w_sr.append(alpha * w_0 + (1 - alpha) * w_1)
            sharpe_ratios.append(sharpe_ratio)
        maximum_sharp_ratio = max(sharpe_ratios)
        weights_with_max_sharpe_ratio = w_sr[sharpe_ratios.index(maximum_sharp_ratio)]
        return maximum_sharp_ratio, weights_with_max_sharpe_ratio

    def _min_volatility(self):
        var = []
        for weights in self.weights:
            volatility = np.dot(np.dot(weights.T, self.cov_matrix), weights)
            var.append(volatility)
        min_var = min(var)
        return min_var ** .5, self.weights[var.index(min_var)]

    def _efficient_frontier(self, points=100):
        means, sigma, weights = [], [], []
        partitions = np.linspace(0, 1, points // len(self.weights))[:-1]
        b = list(range(len(self.weights) - 1))
        for i in b:
            w_0, w_1 = self.weights[i], self.weights[i + 1]
            if i == b[-1]:
                partitions = np.linspace(0, 1, points // len(self.weights))
            for j in partitions:
                w = w_1 * j + (1 - j) * w_0
                weights.append(np.copy(w))
                means.append(np.dot(w.T, self.expected_returns)[0, 0])
                sigma.append(np.dot(np.dot(w.T, self.cov_matrix), w)[0, 0] ** 0.5)
        return means, sigma, weights


class CLAModified(CLA):
    def _initialise(self, asset_prices, covariance, resample_by):
        if not isinstance(asset_prices, pd.DataFrame):
            raise ValueError("Asset prices matrix must be a dataframe")
        if not isinstance(asset_prices.index, pd.DatetimeIndex):
            raise ValueError("Asset prices dataframe must be indexed by date.")

        asset_prices = asset_prices.resample(resample_by).last()

        if self.calculate_returns == "mean":
            self.expected_returns = self._calculate_mean_historical_returns(asset_prices=asset_prices)
        elif self.calculate_returns == "exponential":
            self.expected_returns = self._calculate_exponential_historical_returns(asset_prices=asset_prices)
        else:
            raise ValueError("Unknown returns specified. Supported returns - mean, exponential")
            
        self.expected_returns = np.array(self.expected_returns).reshape((len(self.expected_returns), 1))
        if (self.expected_returns == np.ones(self.expected_returns.shape) * self.expected_returns.mean()).all():
            self.expected_returns[-1, 0] += 1e-5

        self.cov_matrix = np.asarray(covariance)

        if isinstance(self.weight_bounds[0], numbers.Real):
            self.lower_bounds = np.ones(self.expected_returns.shape) * self.weight_bounds[0]
        else:
            self.lower_bounds = np.array(self.weight_bounds[0]).reshape(self.expected_returns.shape)

        if isinstance(self.weight_bounds[0], numbers.Real):
            self.upper_bounds = np.ones(self.expected_returns.shape) * self.weight_bounds[1]
        else:
            self.upper_bounds = np.array(self.weight_bounds[1]).reshape(self.expected_returns.shape)

        self.weights = []
        self.lambdas = []
        self.gammas = []
        self.free_weights = []

    def allocate(self, asset_prices, covariance, solution="cla_turning_points", resample_by="B"):
        self._initialise(asset_prices=asset_prices, covariance=covariance, resample_by=resample_by)
        assets = asset_prices.columns

        free_weights, weights = self._init_algo()
        self.weights.append(np.copy(weights))  
        self.lambdas.append(None)
        self.gammas.append(None)
        self.free_weights.append(free_weights[:])
        
        while True:
            lambda_in, i_in, bi_in = self._bound_free_weight(free_weights)
            lambda_out, i_out = self._free_bound_weight(free_weights)

            if (lambda_in is None or lambda_in < 0) and (lambda_out is None or lambda_out < 0):
                self.lambdas.append(0)
                covar_f, covar_fb, mean_f, w_b = self._get_matrices(free_weights)
                covar_f_inv = np.linalg.inv(covar_f)
                mean_f = np.zeros(mean_f.shape)
            else:
                if self._infnone(lambda_in) > self._infnone(lambda_out):
                    self.lambdas.append(lambda_in)
                    free_weights.remove(i_in)
                    weights[i_in] = bi_in 
                else:
                    self.lambdas.append(lambda_out)
                    free_weights.append(i_out)
                covar_f, covar_fb, mean_f, w_b = self._get_matrices(free_weights)
                covar_f_inv = np.linalg.inv(covar_f)

            w_f, gamma = self._compute_w(covar_f_inv, covar_fb, mean_f, w_b)
            for i in range(len(free_weights)):
                weights[free_weights[i]] = w_f[i]
            self.weights.append(np.copy(weights)) 
            self.gammas.append(gamma)
            self.free_weights.append(free_weights[:])
            if self.lambdas[-1] == 0:
                break

        self._purge_num_err(10e-10)
        self._purge_excess()
        self._compute_solution(assets=assets, solution=solution)

# ==============================================================================
# CUSTOM SECURITY INITIALIZER (For Slippage & Brokerage Models)
# ==============================================================================

class CustomSecurityInitializer(BrokerageModelSecurityInitializer):
    def __init__(self, brokerage_model, security_seeder):
        super().__init__(brokerage_model, security_seeder)

    def Initialize(self, security):
        super().Initialize(security)
        security.SetSlippageModel(ConstantSlippageModel(0.001))

# ==============================================================================
# QUANTCONNECT ALGORITHM
# ==============================================================================

class UltimateFactorCLA(QCAlgorithm):
    def Initialize(self):
        self.SetStartDate(2021, 1, 1)
        self.SetCash(100000)
        
        self.UniverseSettings.Resolution = Resolution.Minute
        self.AddUniverse(self.FundamentalSelection)
        self.spy = self.AddEquity("SPY", Resolution.Minute).Symbol
        
        self.max_candidates = 50 
        self.final_count = 25     
        self.candidates = []
        
        self.next_universe_time = self.Time
        self.pending_weights = {}
        self.pending_liquidations = []

        self.SetWarmUp(60)
        
        self.SetBrokerageModel(BrokerageName.ALPACA)
        self.SetSecurityInitializer(CustomSecurityInitializer(self.BrokerageModel, SecuritySeeder.Null))

        self.Schedule.On(self.DateRules.WeekStart(self.spy), 
                         self.TimeRules.AfterMarketOpen(self.spy, 30), 
                         self.QueueTrades)

    def FundamentalSelection(self, fundamental):
        if self.Time < self.next_universe_time:
            return Universe.Unchanged
            
        self.next_universe_time = self.Time + timedelta(days=7)

        # 1. Base Filter: Must have Price/Cap and a valid P/B ratio
        filtered = [f for f in fundamental if f.HasFundamentalData 
                                          and f.Price > 5 
                                          and f.MarketCap > 1e8 
                                          and f.ValuationRatios.PBRatio > 0]
        
        # Ensure we have enough data to bucket safely
        if len(filtered) < 1000: return Universe.Unchanged

        # 2. Sort by Market Cap to define our universe boundaries
        sorted_by_size = sorted(filtered, key=lambda x: x.MarketCap, reverse=True)
        
        # 3. Stratify into Cap Buckets
        large_caps = sorted_by_size[:200]        # Top 200 Largest Companies
        mid_caps = sorted_by_size[200:500]       # Next 300 Companies
        small_caps = sorted_by_size[500:1000]    # Next 500 Companies

        # Helper function: Sorts a given bucket by Value (High Book-to-Market)
        def get_value_stocks(bucket, count):
            # Book-to-Market = 1 / PBRatio
            sorted_bucket = sorted(bucket, key=lambda x: 1 / x.ValuationRatios.PBRatio, reverse=True)
            return [x.Symbol for x in sorted_bucket[:count]]

        # 4. Extract Value Stocks from each strata (Total = 50 Candidates)
        large_value_symbols = get_value_stocks(large_caps, 20)  # 40%
        mid_value_symbols = get_value_stocks(mid_caps, 15)      # 30%
        small_value_symbols = get_value_stocks(small_caps, 15)  # 30%
        
        # 5. Combine into the final candidate list
        self.candidates = large_value_symbols + mid_value_symbols + small_value_symbols
        
        self.Debug(f"Stratified Universe: 20 Large-Value, 15 Mid-Value, 15 Small-Value generated.")
        return self.candidates

    def QueueTrades(self):
        if self.IsWarmingUp or not self.candidates: return

        # Risk-Off Check
        spy_history = self.History(self.spy, 200, Resolution.Daily)
        if not spy_history.empty:
            spy_current = spy_history['close'].iloc[-1]
            spy_sma = spy_history['close'].mean()
            if spy_current < spy_sma:
                self.Liquidate() 
                self.pending_weights.clear() 
                self.Debug("Market Risk-Off: Liquidating to Cash")
                return

        # Fetch candidate history
        history = self.History(self.candidates, 60, Resolution.Daily)
        if history.empty: return
        
        prices = history['close'].unstack(level=0).ffill().dropna(axis=1)
        if prices.empty or len(prices.columns) < self.final_count: return

        # Target Top 25 Momentum Stocks from our Fama-French Universe
        mom_scores = (prices.iloc[-1] / prices.iloc[0]) - 1
        top_symbols = mom_scores.sort_values(ascending=False).head(self.final_count).index.tolist()
        
        target_prices = prices[top_symbols]
        target_prices.index = pd.to_datetime(target_prices.index) 
        
        try:
            # Generate Sklearn Shrinkage Covariance Matrix
            target_returns = target_prices.pct_change().dropna(how='all')
            oas = OAS()
            oas.fit(target_returns)
            target_cov = pd.DataFrame(oas.covariance_, index=target_returns.columns, columns=target_returns.columns)
            
            # Execute Critical Line Algorithm (CLA) Optimization
            cla = CLAModified(weight_bounds=(0, 1))
            cla.allocate(asset_prices=target_prices, covariance=target_cov, resample_by='B', solution="max_sharpe")
            
            # Extract optimal target weights
            cla_weights_series = cla.weights.iloc[0]
            
            # Clean floating point dust to prevent QC target percentage errors
            clean_weights = {}
            for sym, w in cla_weights_series.items():
                if w > 0.0001:  
                    clean_weights[sym] = round(w, 4) 
            
            # Queue Liquidations
            current_holdings = [x.Key for x in self.Portfolio if x.Value.Invested]
            self.pending_liquidations = [sym for sym in current_holdings if sym not in clean_weights]
            
            # Set Target Weights Directly
            self.pending_weights = clean_weights
            
            # Trigger first attempt
            self.ExecutePendingTrades()
                
        except Exception as e:
            self.Debug(f"Trade Execution Error: {e}")

    def OnData(self, data):
        """Fallback loop: Retries stuck trades every 10 minutes"""
        if not self.pending_weights and not self.pending_liquidations: return
            
        if self.Time.minute % 10 == 0:
            self.ExecutePendingTrades()

    def ExecutePendingTrades(self):
        """Safely processes the queues only if valid live price data exists"""
        completed_liquidations = []
        for symbol in self.pending_liquidations:
            if self.Securities.ContainsKey(symbol) and self.Securities[symbol].Price > 0:
                self.SetHoldings(symbol, 0)
                completed_liquidations.append(symbol)
                
        for symbol in completed_liquidations:
            self.pending_liquidations.remove(symbol)

        completed_allocations = []
        for symbol, weight in self.pending_weights.items():
            if self.Securities.ContainsKey(symbol) and self.Securities[symbol].Price > 0:
                self.SetHoldings(symbol, weight)
                completed_allocations.append(symbol)

        for symbol in completed_allocations:
            del self.pending_weights[symbol]

from AlgorithmImports import *
import pandas as pd
import numpy as np
from scipy.cluster.hierarchy import dendrogram, linkage
from scipy.spatial.distance import squareform
from sklearn.covariance import OAS

# ==============================================================================
# MARCOS LÓPEZ DE PRADO HRP CLASSES (with Shrinkage & Modifications)
# ==============================================================================

class HierarchicalRiskParity:
    def __init__(self):
        self.weights = list()
        self.seriated_correlations = None
        self.seriated_distances = None
        self.ordered_indices = None
        self.clusters = None

    @staticmethod
    def _tree_clustering(correlation, method='single'):
        distances = np.sqrt((1 - correlation).round(5) / 2)
        clusters = linkage(squareform(distances.values), method=method)
        return distances, clusters

    def _quasi_diagnalization(self, num_assets, curr_index):
        if curr_index < num_assets:
            return [curr_index]
        left = int(self.clusters[curr_index - num_assets, 0])
        right = int(self.clusters[curr_index - num_assets, 1])
        return (self._quasi_diagnalization(num_assets, left) + self._quasi_diagnalization(num_assets, right))

    def _get_seriated_matrix(self, assets, distances, correlations):
        ordering = assets[self.ordered_indices]
        seriated_distances = distances.loc[ordering, ordering]
        seriated_correlations = correlations.loc[ordering, ordering]
        return seriated_distances, seriated_correlations

    def _recursive_bisection(self, covariances, assets):
        self.weights = pd.Series(1.0, index=self.ordered_indices)
        clustered_alphas = [self.ordered_indices]

        while clustered_alphas:
            clustered_alphas = [cluster[start:end]
                                for cluster in clustered_alphas
                                for start, end in ((0, len(cluster) // 2), (len(cluster) // 2, len(cluster)))
                                if len(cluster) > 1]

            for subcluster in range(0, len(clustered_alphas), 2):
                left_cluster = clustered_alphas[subcluster]
                right_cluster = clustered_alphas[subcluster + 1]

                left_subcovar = covariances.iloc[left_cluster, left_cluster]
                inv_diag = 1 / np.diag(left_subcovar.values)
                parity_w = inv_diag * (1 / np.sum(inv_diag))
                left_cluster_var = np.dot(parity_w, np.dot(left_subcovar, parity_w))

                right_subcovar = covariances.iloc[right_cluster, right_cluster]
                inv_diag = 1 / np.diag(right_subcovar.values)
                parity_w = inv_diag * (1 / np.sum(inv_diag))
                right_cluster_var = np.dot(parity_w, np.dot(right_subcovar, parity_w))

                alloc_factor = 1 - left_cluster_var / (left_cluster_var + right_cluster_var)
                self.weights[left_cluster] *= alloc_factor
                self.weights[right_cluster] *= 1 - alloc_factor

        self.weights.index = assets[self.ordered_indices]
        self.weights = pd.DataFrame(self.weights).T

    @staticmethod
    def _calculate_returns(asset_prices, resample_by):
        if resample_by is not None:
            asset_prices = asset_prices.resample(resample_by).last()
        asset_returns = asset_prices.pct_change()
        asset_returns = asset_returns.dropna(how='all')
        return asset_returns

    @staticmethod
    def _shrink_covariance(covariance):
        oas = OAS()
        oas.fit(covariance)
        shrinked_covariance = oas.covariance_
        return pd.DataFrame(shrinked_covariance, index=covariance.columns, columns=covariance.columns)

    @staticmethod
    def _cov2corr(covariance):
        d_matrix = np.zeros_like(covariance)
        diagnoal_sqrt = np.sqrt(np.diag(covariance))
        np.fill_diagonal(d_matrix, diagnoal_sqrt)
        d_inv = np.linalg.inv(d_matrix)
        corr = np.dot(np.dot(d_inv, covariance), d_inv)
        corr = pd.DataFrame(corr, index=covariance.columns, columns=covariance.columns)
        return corr


class HierarchicalRiskParityModified(HierarchicalRiskParity):
    def allocate(self, asset_prices, covariance, resample_by='B', use_shrinkage=False):
        if not isinstance(asset_prices, pd.DataFrame):
            raise ValueError("Asset prices matrix must be a dataframe")
        if not isinstance(asset_prices.index, pd.DatetimeIndex):
            raise ValueError("Asset prices dataframe must be indexed by date.")

        asset_returns = self._calculate_returns(asset_prices, resample_by=resample_by)
        num_assets = asset_returns.shape[1]
        assets = asset_returns.columns
        
        cov = pd.DataFrame(covariance, columns=assets, index=assets)
        
        if use_shrinkage:
            cov = self._shrink_covariance(covariance=cov)
        corr = self._cov2corr(covariance=cov)

        distances, self.clusters = self._tree_clustering(correlation=corr)
        self.ordered_indices = self._quasi_diagnalization(num_assets, 2 * num_assets - 2)
        self.seriated_distances, self.seriated_correlations = self._get_seriated_matrix(assets=assets, distances=distances, correlations=corr)
        self._recursive_bisection(covariances=cov, assets=assets)


# ==============================================================================
# CUSTOM SECURITY INITIALIZER (For Slippage & Brokerage Models)
# ==============================================================================

class CustomSecurityInitializer(BrokerageModelSecurityInitializer):
    def __init__(self, brokerage_model, security_seeder):
        super().__init__(brokerage_model, security_seeder)

    def Initialize(self, security):
        # 1. Apply the default Alpaca brokerage models (fees, margin, etc.)
        super().Initialize(security)
        # 2. Apply our custom 0.1% slippage penalty on top
        security.SetSlippageModel(ConstantSlippageModel(0.001))


# ==============================================================================
# QUANTCONNECT ALGORITHM
# ==============================================================================

class UltimateFactorHRP(QCAlgorithm):
    def Initialize(self):
        self.SetStartDate(2021, 1, 1)
        self.SetCash(100000)
        
        # 1. Environment & Universe
        self.UniverseSettings.Resolution = Resolution.Minute
        self.AddUniverse(self.FundamentalSelection)
        self.spy = self.AddEquity("SPY", Resolution.Minute).Symbol
        
        # 2. Strategy Variables
        self.max_candidates = 50 
        self.final_count = 15     
        self.candidates = []
        self.weight_buffer = 0.02  # 2% weight delta buffer
        
        # 3. State Variables
        self.next_universe_time = self.Time
        self.pending_weights = {}
        self.pending_liquidations = []

        self.SetWarmUp(60)
        
        # 4. Brokerage & Custom Security Initializer
        self.SetBrokerageModel(BrokerageName.ALPACA)
        self.SetSecurityInitializer(CustomSecurityInitializer(self.BrokerageModel, SecuritySeeder.Null))

        # 5. Scheduling (Daily Rebalance)
        self.Schedule.On(self.DateRules.EveryDay(self.spy), 
                         self.TimeRules.AfterMarketOpen(self.spy, 30), 
                         self.QueueTrades)

    def FundamentalSelection(self, fundamental):
        """Locks universe updates to a weekly cycle"""
        if self.Time < self.next_universe_time:
            return Universe.Unchanged
            
        self.next_universe_time = self.Time + timedelta(days=7)

        filtered = [f for f in fundamental if f.HasFundamentalData and f.Price > 5 and f.MarketCap > 1e8]
        sorted_by_cap = sorted(filtered, key=lambda x: x.MarketCap, reverse=True)
        self.candidates = [x.Symbol for x in sorted_by_cap[:self.max_candidates]]
        
        return self.candidates

    def QueueTrades(self):
        """Calculates daily targets and pushes them into the execution queues"""
        if self.IsWarmingUp or not self.candidates: return

        # Risk-Off Check
        spy_history = self.History(self.spy, 200, Resolution.Daily)
        if not spy_history.empty:
            spy_current = spy_history['close'].iloc[-1]
            spy_sma = spy_history['close'].mean()
            if spy_current < spy_sma:
                self.Liquidate() 
                self.pending_weights.clear() 
                self.Debug("Market Risk-Off: Liquidating to Cash")
                return

        # Fetch candidate history
        history = self.History(self.candidates, 60, Resolution.Daily)
        if history.empty: return
        
        prices = history['close'].unstack(level=0).ffill().dropna(axis=1)
        if prices.empty or len(prices.columns) < self.final_count: return

        # Target Top 15 Momentum Stocks
        mom_scores = (prices.iloc[-1] / prices.iloc[0]) - 1
        top_symbols = mom_scores.sort_values(ascending=False).head(self.final_count).index.tolist()
        
        target_prices = prices[top_symbols]
        target_prices.index = pd.to_datetime(target_prices.index) 
        
        try:
            # HRP Math & Allocation
            target_returns = target_prices.pct_change().dropna(how='all')
            target_cov = target_returns.cov()
            
            hrp = HierarchicalRiskParityModified()
            hrp.allocate(asset_prices=target_prices, covariance=target_cov, resample_by='B', use_shrinkage=True)
            
            hrp_weights_series = hrp.weights.iloc[0]
            
            # Queue Liquidations
            current_holdings = [x.Key for x in self.Portfolio if x.Value.Invested]
            self.pending_liquidations = [sym for sym in current_holdings if sym not in hrp_weights_series.index]
            
            # Apply Weight Delta Filter
            self.pending_weights.clear()
            for symbol, target_weight in hrp_weights_series.items():
                current_weight = 0
                if self.Portfolio.ContainsKey(symbol) and self.Portfolio[symbol].Invested:
                    current_weight = self.Portfolio[symbol].HoldingsValue / self.Portfolio.TotalPortfolioValue
                
                weight_delta = abs(target_weight - current_weight)
                
                if weight_delta >= self.weight_buffer:
                    self.pending_weights[symbol] = target_weight
            
            # Trigger first attempt
            self.ExecutePendingTrades()
                
        except Exception as e:
            self.Debug(f"Trade Execution Error: {e}")

    def OnData(self, data):
        """Fallback loop: Retries stuck trades every 10 minutes"""
        if not self.pending_weights and not self.pending_liquidations: return
            
        if self.Time.minute % 10 == 0:
            self.ExecutePendingTrades()

    def ExecutePendingTrades(self):
        """Safely processes the queues only if valid live price data exists"""
        completed_liquidations = []
        for symbol in self.pending_liquidations:
            if self.Securities.ContainsKey(symbol) and self.Securities[symbol].Price > 0:
                self.SetHoldings(symbol, 0)
                completed_liquidations.append(symbol)
                
        for symbol in completed_liquidations:
            self.pending_liquidations.remove(symbol)

        completed_allocations = []
        for symbol, weight in self.pending_weights.items():
            if self.Securities.ContainsKey(symbol) and self.Securities[symbol].Price > 0:
                self.SetHoldings(symbol, weight)
                completed_allocations.append(symbol)

        for symbol in completed_allocations:
            del self.pending_weights[symbol]