% problem 1, part c
clc, clear

% define domain for x1, x2, x3, x4
domain = [0, 2];

% define different resolutions for sampling
h_vals = [1, 0.8, 0.75, 0.5, 0.25, 0.1];%, 0.05];

% create empty vector to store K values
K_vals = zeros(size(h_vals));

% loop over h values
for hi = 1:length(h_vals)
    h = h_vals(hi);
    max_K = 0;

    % generate sample points within domain using h
    x_range = domain(1):h:domain(2);

    fprintf('For h = %.4f, we have %d^4 = %d total samples.\n', h, length(x_range), length(x_range)^4);
    
    % loop over 4D grid of values (for x1, x2, x3, x4)
    for x1 = x_range
        for x2 = x_range
            for x3 = x_range
                for x4 = x_range
                    % evaluate Jacobian at this sampled point
                    J = [-2-sin(x1), 10.*x2, 0, -2.*x4;
                        0, -1, 1, 0;
                        0, cos(x2), -8.*x3, 0;
                        0, 0, 0, -1];

                    % compute K = 2-norm of Jacobian
                    K = norm(J, 2);

                    % update maximum Lipschitz constant
                    if K > max_K
                        max_K = K;
                    end
                end
            end
        end
    end
    % store result for this h
    K_vals(hi) = max_K;
end

% plot results of K for various h
f1 = figure();
plot(h_vals, K_vals, "-o", "Linewidth", 3, "Color", "r", "MarkerFaceColor", "r");
set(gca, "FontSize", 12, "XDir", "reverse");
grid;
title("\textbf{Approximation of Lipschitz Constant}", "\textbf{vs Sampling Resolution}", ...
    "FontSize", 16, "Interpreter", "latex");
xlabel("Step Size $(h)$", "FontSize", 16, "Interpreter", "latex");
ylabel("Numerical Lipschitz Constant $(K)$", "FontSize", 16, "Interpreter", "latex");

% section for randomly sampling points
% set seed
rng(66);
% try increasing sample sizes
num_samples = [10, 100, 1000, 5000, 10000, 50000]; 

K_random = zeros(size(num_samples));

for k = 1:length(num_samples)
    N = num_samples(k);
    max_K = 0;

    fprintf('Generating %d random samples.\n', N);
    
    for i = 1:N
        % Uniformly sample each x_i in [0, 2]
        x1 = 2 * rand();
        x2 = 2 * rand();
        x3 = 2 * rand();
        x4 = 2 * rand();
        
        % evaluate Jacobian at this randomly sampled point
        J = [
            -2 - sin(x1), 10 * x2, 0, -2 * x4;
            0, -1, 1, 0;
            0, cos(x2), -8 * x3, 0;
            0, 0, 0, -1
        ];
        
        % compute K = 2-norm of Jacobian
        K = norm(J, 2);
        
        if K > max_K
            max_K = K;
        end
    end
    % store result
    K_random(k) = max_K;
end

% Plot Lipschitz estimate vs number of random samples
f2 = figure;
plot(num_samples, K_random, "-o", "Linewidth", 3, "Color", "b", "MarkerFaceColor", "b");
set(gca, "FontSize", 12);
grid;
title("\textbf{Approximation of Lipschitz Constant}", "\textbf{via Random Sampling}", ...
    "FontSize", 16, "Interpreter", "latex");
xlabel("Number of Random Samples", "FontSize", 16, "Interpreter", "latex");
ylabel("Numerical Lipschitz Constant $(K)$", "FontSize", 16, "Interpreter", "latex");