Artificial Intelligence
— Basics & Essentials

% Breadth-First Search on an adjacency list
% Graph represented as cell array of neighbour lists

function path = bfs(adj, start, goal)
  n = length(adj);
  visited = false(1, n);
  parent  = zeros(1, n);   % 0 = no parent
  queue   = start;
  visited(start) = true;

  while ~isempty(queue)
    node  = queue(1);
    queue = queue(2:end);   % dequeue

    if node == goal
      path = reconstruct(parent, start, goal);
      return;
    end

    for nb = adj{node}
      if ~visited(nb)
        visited(nb) = true;
        parent(nb)  = node;
        queue       = [queue, nb];   % enqueue
      end
    end
  end
  path = [];   % no path found
end

function path = reconstruct(parent, start, goal)
  path = goal;
  node = goal;
  while node ~= start
    node = parent(node);
    path = [node, path];
  end
end

% --- Example usage ---
adj = {[2,3], [1,4,5], [1,6], [2], [2], [3]};
p   = bfs(adj, 1, 6);
disp(p)   % → 1  3  6

% Generate full truth table for n propositional variables

function T = truth_table(n)
  rows = 2^n;
  T = false(rows, n);
  for col = 1:n
    period = 2^(n - col);
    for row = 1:rows
      T(row, col) = mod(floor((row-1) / period), 2) == 0;
    end
  end
end

% Truth table for: (P ∧ Q) → R
T = truth_table(3);   % cols: P, Q, R
P = T(:,1); Q = T(:,2); R = T(:,3);
formula = (~(P & Q)) | R;   % (P∧Q)→R ≡ ¬(P∧Q) ∨ R

disp('P  Q  R  | (P∧Q)→R');
for i = 1:rows(T)
  fprintf('%d  %d  %d  |  %d\n', P(i), Q(i), R(i), formula(i));
end

% Naive Bayes classifier (Gaussian likelihood)
% Train on (X, y) where X is n×d feature matrix, y is n×1 class labels

function model = nb_train(X, y)
  classes = unique(y);
  model.classes = classes;
  for k = 1:length(classes)
    idx = (y == classes(k));
    model.prior(k)  = mean(idx);            % P(Cₖ)
    model.mu(k,:)   = mean(X(idx,:), 1);    % μₖ per feature
    model.sigma(k,:)= var(X(idx,:), 0, 1);  % σ²ₖ per feature
  end
end

function ypred = nb_predict(model, X)
  n = rows(X); K = length(model.classes);
  logpost = zeros(n, K);
  for k = 1:K
    % log P(Cₖ) + Σ log N(xⱼ; μₖⱼ, σ²ₖⱼ)
    logpost(:,k) = log(model.prior(k)) + ...
      sum(-0.5*log(2*pi*model.sigma(k,:)) - ...
          (bsxfun(@minus, X, model.mu(k,:))).^2 ./ ...
          (2*model.sigma(k,:)), 2);
  end
  [~, idx] = max(logpost, [], 2);
  ypred = model.classes(idx);
end

% --- Synthetic data example ---
randn('state', 42);
X = [randn(50,2); randn(50,2) + 3];
y = [ones(50,1); 2*ones(50,1)];
model   = nb_train(X, y);
ypred   = nb_predict(model, X);
accuracy = mean(ypred == y);
fprintf('Accuracy: %.1f%%\n', accuracy*100);

% Principal Component Analysis — implemented from scratch

function [Z, V, explained] = pca_scratch(X, r)
  % X: n×d data matrix (samples × features)
  % r: number of principal components to keep

  % 1. Centre the data
  mu  = mean(X, 1);
  Xc  = bsxfun(@minus, X, mu);

  % 2. Covariance matrix
  Cov = (Xc' * Xc) / (rows(X) - 1);

  % 3. Eigendecomposition
  [V, D] = eig(Cov);

  % 4. Sort descending by eigenvalue
  eigvals = diag(D);
  [eigvals, idx] = sort(eigvals, 'descend');
  V = V(:, idx);

  % 5. Variance explained
  explained = eigvals / sum(eigvals);

  % 6. Project onto top-r components
  Z = Xc * V(:, 1:r);
end

% --- Example: 4D iris-like data → 2D ---
randn('state', 0);
X = randn(150, 4);
X(:,1:2) = X(:,1:2) * 2;   % inflate first two dims

[Z, V, exp_var] = pca_scratch(X, 2);
fprintf('Variance explained by PC1: %.1f%%\n', exp_var(1)*100);
fprintf('Variance explained by PC2: %.1f%%\n', exp_var(2)*100);
scatter(Z(:,1), Z(:,2), 20, 'filled');
xlabel('PC1'); ylabel('PC2'); title('PCA Projection');

% MLP Forward Pass — 2 hidden layers, sigmoid activation

function sigma = sigmoid(z)
  sigma = 1 ./ (1 + exp(-z));
end

function [a_out, cache] = forward_pass(X, W1, b1, W2, b2, W3, b3)
  % Layer 1
  z1 = X  * W1' + repmat(b1, rows(X), 1);
  a1 = sigmoid(z1);
  % Layer 2
  z2 = a1 * W2' + repmat(b2, rows(X), 1);
  a2 = sigmoid(z2);
  % Output layer
  z3 = a2 * W3' + repmat(b3, rows(X), 1);
  a_out = sigmoid(z3);   % binary classification
  cache = {z1,a1,z2,a2,z3};
end

% --- Initialise a 2→4→4→1 network randomly ---
randn('state', 1);
W1 = randn(4,2) * 0.1;  b1 = zeros(1,4);
W2 = randn(4,4) * 0.1;  b2 = zeros(1,4);
W3 = randn(1,4) * 0.1;  b3 = zeros(1,1);

X = [0 0; 0 1; 1 0; 1 1];   % XOR inputs
[output, ~] = forward_pass(X, W1,b1, W2,b2, W3,b3);
disp(output)

% Linear Regression: Normal Equations vs Gradient Descent

%% Generate synthetic data
rand('state', 42);
n = 100;
x = 2*rand(n,1);
y = 3*x + 1 + 0.4*randn(n,1);   % y = 3x + 1 + noise
X = [ones(n,1), x];                 % design matrix [1, x]

%% 1. Normal Equations (closed-form)
theta_ne = (X'*X) \ (X'*y);
fprintf('Normal Equations: θ₀=%.4f, θ₁=%.4f\n', theta_ne(1), theta_ne(2));

%% 2. Gradient Descent
theta = zeros(2,1);
alpha = 0.1;
iters = 500;
J_hist = zeros(iters,1);

for t = 1:iters
  err  = X*theta - y;
  grad = (X'*err) / n;       % ∇J = (1/n)Xᵀ(Xθ−y)
  theta = theta - alpha*grad;
  J_hist(t) = mean(err.^2) / 2;
end
fprintf('Gradient Descent:  θ₀=%.4f, θ₁=%.4f\n', theta(1), theta(2));

%% Plot loss curve
plot(1:iters, J_hist, 'b-', 'LineWidth', 2);
xlabel('Iteration'); ylabel('MSE Loss');
title('Gradient Descent Convergence');

% Solve XOR with 2-layer network trained via backprop
% Architecture: 2 → 4 → 1  (sigmoid throughout)

randn('state', 7);
X = [0 0; 0 1; 1 0; 1 1];
y = [0; 1; 1; 0];          % XOR truth table

% Init weights (Xavier-like)
W1 = randn(4,2) / sqrt(2);  b1 = zeros(4,1);
W2 = randn(1,4) / sqrt(4);  b2 = zeros(1,1);

alpha = 1.0;   % learning rate
losses = zeros(5000,1);

for epoch = 1:5000
  %% Forward
  Z1 = W1*X' + b1;          % 4×4
  A1 = 1./(1+exp(-Z1));     % sigmoid
  Z2 = W2*A1 + b2;          % 1×4
  A2 = 1./(1+exp(-Z2));     % output

  %% Loss (binary cross-entropy)
  L = -mean(y'.*log(A2+1e-9) + (1-y').*log(1-A2+1e-9));
  losses(epoch) = L;

  %% Backward
  dA2 = -(y'./(A2+1e-9) - (1-y')./(1-A2+1e-9)) / 4;
  dZ2 = dA2 .* A2 .* (1 - A2);        % σ'(z) = σ(1−σ)
  dW2 = dZ2 * A1';                       % 1×4
  db2 = sum(dZ2, 2);

  dA1 = W2' * dZ2;                       % 4×4
  dZ1 = dA1 .* A1 .* (1 - A1);
  dW1 = dZ1 * X;                         % 4×2
  db1 = sum(dZ1, 2);

  %% Update
  W1 = W1 - alpha * dW1;  b1 = b1 - alpha * db1;
  W2 = W2 - alpha * dW2;  b2 = b2 - alpha * db2;
end

fprintf('\nFinal predictions (should be ~[0,1,1,0]):\n');
disp(round(A2', 3));
fprintf('Final loss: %.6f\n', losses(end));

% Adam optimiser applied to a simple 2D loss surface
% f(x,y) = x² + 10y²  (Rosenbrock-like elongated bowl)

function [g] = grad_f(params)
  g = [2*params(1); 20*params(2)];
end

% Hyperparameters
alpha = 0.01;
beta1 = 0.9;  beta2 = 0.999;  eps = 1e-8;

theta = [-3; 2];   % start far from optimum (0,0)
m = zeros(2,1);  v = zeros(2,1);
path = theta';

for t = 1:200
  g  = grad_f(theta);
  m  = beta1*m + (1-beta1)*g;
  v  = beta2*v + (1-beta2)*g.^2;
  mh = m / (1 - beta1^t);
  vh = v / (1 - beta2^t);
  theta = theta - alpha * mh ./ (sqrt(vh) + eps);
  path  = [path; theta'];
end

fprintf('Final θ: [%.6f, %.6f]\n', theta(1), theta(2));

% Plot convergence path
[gx,gy] = meshgrid(linspace(-4,4,100), linspace(-3,3,100));
fval = gx.^2 + 10*gy.^2;
contourf(gx, gy, fval, 30, 'LineColor', 'none'); hold on
plot(path(:,1), path(:,2), 'w-o', 'MarkerSize', 3);
title('Adam Optimisation Path');

% Compute Shannon entropy and Information Gain

function H = entropy(labels)
  if isempty(labels), H = 0; return; end
  classes = unique(labels);
  H = 0;
  n = length(labels);
  for c = classes'
    p = sum(labels == c) / n;
    if p > 0
      H = H - p * log2(p);
    end
  end
end

function IG = info_gain(y, feature, threshold)
  left  = y(feature <= threshold);
  right = y(feature >  threshold);
  n = length(y);
  IG = entropy(y) ...
     - (length(left)/n)  * entropy(left) ...
     - (length(right)/n) * entropy(right);
end

% --- Example: find best split threshold ---
feature = [1;2;3;4;5;6;7;8];
labels  = [0;0;0;1;1;1;0;1];

fprintf('Parent entropy: %.4f bits\n', entropy(labels));
thresholds = 1:7;
for t = thresholds
  ig = info_gain(labels, feature, t);
  fprintf('  Split ≤ %d : IG = %.4f\n', t, ig);
end

% Compute fairness metrics for a binary classifier
% y_true: ground truth, y_pred: predictions, group: protected attribute (0/1)

function metrics = fairness_metrics(y_true, y_pred, group)
  g = unique(group);
  for k = 1:length(g)
    idx = group == g(k);
    yt = y_true(idx); yp = y_pred(idx);

    TP = sum(yt==1 & yp==1);
    FP = sum(yt==0 & yp==1);
    TN = sum(yt==0 & yp==0);
    FN = sum(yt==1 & yp==0);

    metrics(k).group       = g(k);
    metrics(k).accuracy    = (TP+TN)/length(yt);
    metrics(k).TPR         = TP/(TP+FN+1e-9);  % true positive rate
    metrics(k).FPR         = FP/(FP+TN+1e-9);  % false positive rate
    metrics(k).pred_pos    = mean(yp);          % P(ŷ=1) — demographic parity

    fprintf('Group %d: Acc=%.2f, TPR=%.2f, FPR=%.2f, P(ŷ=1)=%.2f\n', ...
      g(k), metrics(k).accuracy, metrics(k).TPR, metrics(k).FPR, metrics(k).pred_pos);
  end
  dp_gap = abs(metrics(1).pred_pos - metrics(2).pred_pos);
  fprintf('Demographic Parity Gap: %.4f\n', dp_gap);
end

% --- Example ---
y_true = [1;0;1;1;0;1;0;0;1;1];
y_pred = [1;0;1;0;0;1;1;0;1;0];
group  = [0;0;0;0;0;1;1;1;1;1];
fairness_metrics(y_true, y_pred, group);