function [optimized_params, model, weights] = exponential_column_fit(data, x)
    % EXPONENTIAL_COLUMN_FIT Fits data using exponential scaling and column weights
    %
    % Inputs:
    %   data - 32000×8 matrix of measurements
    %   x    - 32000×1 vector of independent variable values (optional)
    %
    % Outputs:
    %   optimized_params - Final parameters [a, b, w1_raw, ..., w8_raw]
    %   model            - Best fit model matrix
    %   weights          - Optimized column weights [0-1]
    
    % Default x values if not provided
    if nargin < 2
        x = (1:size(data,1))';
    end
    
    % Normalize x to [0 1] for numerical stability
    x_norm = (x - min(x)) / (max(x) - min(x));
    
    % Initial parameters [a, b, w1_raw, ..., w8_raw]
    init_params = [1, 0, zeros(1, 8)]; % Raw weights will be sigmoid-transformed
    
    % Optimization options
    options = optimset('Display', 'final', 'MaxIter', 1000, 'MaxFunEvals', 3000);
    
    % Perform optimization
    optimized_raw = fminsearch(@(p) objective_fn(p, data, x_norm), init_params, options);
    
    % Extract and transform parameters
    a = optimized_raw(1);
    b = optimized_raw(2);
    weights = 1./(1 + exp(-optimized_raw(3:end))); % Sigmoid to [0,1]
    
    % Create final model
    exp_vector = a * exp(b * x_norm);
    model = exp_vector * weights;
    
    % Combine parameters for output
    optimized_params = [a, b, weights];
end

function total_error = objective_fn(params, data, x_norm)
    % Weight regularization strength
    lambda = 0.1; % Adjust to control weight sparsity
    
    % Extract parameters
    a = params(1);
    b = params(2);
    raw_weights = params(3:end);
    
    % Transform weights and ensure they stay in (0,1)
    weights = 1./(1 + exp(-raw_weights));
    
    % Calculate exponential component
    exp_component = a * exp(b * x_norm);
    
    % Create weighted model
    model = exp_component * weights;
    
    % Calculate main error term
    fit_error = sum((data - model).^2, 'all');
    
    % Regularization to prevent extreme weights
    reg_error = lambda * sum(weights.*(1 - weights));
    
    % Total error to minimize
    total_error = fit_error + reg_error;
end