function [best_samples, optimized_params, best_fit] = optimize_exponential_fit(data_matrix)
    % OPTIMIZE_EXPONENTIAL_FIT Performs exponential fit and selects best samples
    % 
    % Inputs:
    %   data_matrix - 5x10 matrix where rows are samples, columns are variables
    % 
    % Outputs:
    %   best_samples - Logical vector indicating selected samples (1×5)
    %   optimized_params - Optimized exponential parameters [a, b, c]
    %   best_fit - Best fit values at each variable position
    
    % Validate input
    [n_samples, n_vars] = size(data_matrix);
    if n_samples ~= 5 || n_vars ~= 10
        error('Input must be a 5×10 matrix');
    end
    
    % Prepare data vectors
    x_data = repmat(1:n_vars, n_samples, 1)';  % Variable indices
    y_data = data_matrix';                      % Transpose for correct ordering
    x_data = x_data(:);                         % Column vector of variable indices
    y_data = y_data(:);                         % Column vector of all values
    
    % Define exponential model
    exp_model = @(p, x) p(1) * exp(p(2) * x) + p(3);
    
    % Initialize parameters [a, b, c, w1, w2, w3, w4, w5]
    init_params = [2, -0.1, 0, zeros(1, n_samples)];  % Start with all samples active
    
    % Optimization options
    options = optimset('Display', 'final', 'MaxIter', 1000, 'MaxFunEvals', 3000, 'TolX', 1e-6);
    
    % Perform optimization
    optimized_params = fminsearch(@(p) objective_fn(p, x_data, y_data, exp_model, n_vars),...
                                 init_params, options);
    
    % Determine best samples (weights > 0 after sigmoid transform)
    sample_weights = 1./(1 + exp(-optimized_params(4:end)));
    best_samples = sample_weights > 0.5;
    
    % Calculate best fit
    best_fit = exp_model(optimized_params(1:3), 1:n_vars);
    
    % Visualization
    figure;
    plot(1:n_vars, data_matrix', 'b.'); hold on;
    plot(1:n_vars, data_matrix(best_samples,:)', 'r.', 'MarkerSize', 12);
    plot(1:n_vars, best_fit, 'k-', 'LineWidth', 2);
    xlabel('Variable Index'); ylabel('Value');
    legend('All Data', 'Selected Samples', 'Exponential Fit');
    title('Exponential Fit with Sample Selection');
    grid on;
end

function error = objective_fn(params, x_data, y_data, exp_model, n_vars)
    % Extract parameters and transform weights
    exp_params = params(1:3);
    sample_weights = 1./(1 + exp(-params(4:end)));  % Sigmoid transform
    
    % Create weight vector (repeat weights for each variable)
    weights = repelem(sample_weights, n_vars)';
    
    % Calculate model predictions
    y_pred = exp_model(exp_params, x_data);
    
    % Calculate weighted error
    residuals = y_data - y_pred;
    weighted_error = sum(weights .* residuals.^2);
    
    % Regularization term (encourage sample selection)
    regularization = 0.5 * sum(1 - sample_weights);
    
    % Total error
    error = weighted_error + regularization;
end