clc;
clear;
close all;

rand("seed",42);

% ============================================
% 1. Generate training samples using LHS
% ============================================

Ntrain = 500;
d = 2;

Utrain = zeros(Ntrain,d);

for j = 1:d

    p = randperm(Ntrain);

    Utrain(:,j) = ...
        (p'-rand(Ntrain,1))/Ntrain;

end

Xtrain = -2+4*Utrain;

x = Xtrain(:,1);
y = Xtrain(:,2);

% Static function
Ytrain = sin(x+0.5).*cos(y-0.5).* ...
    exp(-((x-0.5).^2+(y+0.5).^2));

% ============================================
% 2. Generate independent test samples
% ============================================

Ntest = 2000;

Utest = zeros(Ntest,d);

for j = 1:d

    p = randperm(Ntest);

    Utest(:,j) = ...
        (p'-rand(Ntest,1))/Ntest;

end

Xtest = -2+4*Utest;

xt = Xtest(:,1);
yt = Xtest(:,2);

Ytest = sin(xt+0.5).*cos(yt-0.5).* ...
    exp(-((xt-0.5).^2+(yt+0.5).^2));

% ============================================
% 3. Prepare training and test matrices
% ============================================

P = Xtrain';
T = Ytrain';

Ptest = Xtest';
Ttest = Ytest';

% ============================================
% 4. Train the model
% ============================================

spread = 1;
goal = 0;
maxNeurons = 20;

[w,c,s,sr,er,gError] = ...
    newllm(P,T,spread,goal,maxNeurons);

% ============================================
% 5. Predict test samples
% ============================================

Ypred = llmsim(w,c,s,Ptest);

% ============================================
% 6. Evaluation metrics
% ============================================

error = Ttest-Ypred;

MSE = mean(error.^2);
RMSE = sqrt(MSE);
MAE = mean(abs(error));

R2 = 1-sum(error.^2)/ ...
    sum((Ttest-mean(Ttest)).^2);

fprintf('\n========== TEST RESULTS ==========\n');

fprintf('Training samples = %d\n',Ntrain);
fprintf('Test samples = %d\n',Ntest);
fprintf('Number of neurons = %d\n',size(c,2));

fprintf('Training SSE = %.8f\n',gError);
fprintf('Test MSE = %.8f\n',MSE);
fprintf('Test RMSE = %.8f\n',RMSE);
fprintf('Test MAE = %.8f\n',MAE);
fprintf('Test R2 = %.8f\n',R2);

% ============================================
% 7. Save datasets and results
% ============================================

csvwrite('TrainData.csv',[Xtrain,Ytrain]);
csvwrite('TestData.csv',[Xtest,Ytest]);

csvwrite('TestPrediction.csv', ...
    [Ttest',Ypred']);

% ============================================
% 8. Plot actual versus predicted values
% ============================================

figure;

plot(Ttest,Ypred,'o');
hold on;

plot([min(Ttest),max(Ttest)], ...
     [min(Ttest),max(Ttest)],'r-');

xlabel('Actual output');
ylabel('Predicted output');

title('LOLIMOT: Actual vs Predicted');

grid on;