function [netb,accuracy,imdsValidation,indexes] = ...
    IAS_CNN_digits_best_kernel(net,H,kernels)
% IAS_CNN_digits_best_kernel - find best matching kernels for weights in net
% On input:
%     net (Matlab neural net struct): base neural net with learned kernels
%     H (int): filter size (to get HxH)
%     kernels (kernels struct): set of standard image processing kernels
% On output:
%     netb (neural net structure): net trained using best matched kernels
%     accuracy (float): classification accuracy
%     imdsValidation (Validation stucture): validation set
%     indexes (1x4 vector): best match kernel indexes
% Call:
%     [netb,ab,valb,ind] = IAS_CNN_digits_best_kernel(net,5,kernels);
% Author:
%     T. Henderson
%     UU
%     Spring 2025
%

num_kernels = length(kernels);
W = IAS_extract_weights(net,2);
W1 = IAS_W2K(W,1);
W2 = IAS_W2K(W,2);
W3 = IAS_W2K(W,3);
W4 = IAS_W2K(W,4);
err1 = Inf;
err2 = Inf;
err3 = Inf;
err4 = Inf;
ind1 = 0;
ind2 = 0;
ind3 = 0;
ind4 = 0;
for k = 1:num_kernels
    % kernel
    K = kernels(k).kernel;
    K = IAS_scale(W1,K);
    err = mean(mean(abs(K-W1)));
    if err<err1
        err1 = err;
        ind1 = k;
    end
    err = mean(mean(abs(K-W2)));
    if err<err2
        err2 = err;
        ind2 = k;
    end
    err = mean(mean(abs(K-W3)));
    if err<err3
        err3 = err;
        ind3 = k;
    end
    err = mean(mean(abs(K-W4)));
    if err<err4
        err4 = err;
        ind4 = k;
    end
end
indexes = [ind1,ind2,ind3,ind4];
K1 = kernels(ind1).kernel;
K1 = IAS_scale(W1,K1);
K2 = kernels(ind2).kernel;
K2 = IAS_scale(W1,K2);
K3 = kernels(ind3).kernel;
K3 = IAS_scale(W1,K3);
K4 = kernels(ind4).kernel;
K4 = IAS_scale(W1,K4);
W = zeros(5,5,1,4);
W = IAS_K2W(K1,W,1);
W = IAS_K2W(K2,W,2);
W = IAS_K2W(K3,W,3);
W = IAS_K2W(K4,W,4);

% create image data store
digitDatasetPath = fullfile(matlabroot,"toolbox","nnet","nndemos", ...
    "nndatasets","DigitDataset");
imds = imageDatastore(digitDatasetPath, ...
    IncludeSubfolders=true,LabelSource="foldernames");

% divide into training and validation sets
numTrainFiles = 750; [imdsTrain,imdsValidation] = splitEachLabel(imds,...
    numTrainFiles,"randomize");

% create network layers
layers = [
    imageInputLayer([28 28 1])
	
    convolution2dLayer(H,4,'Weights',W,'WeightLearnRateFactor',0,...
    Padding="same")
    batchNormalizationLayer
    reluLayer
	
    maxPooling2dLayer(2,Stride=2)
	
    convolution2dLayer(H,4,Padding="same")
    batchNormalizationLayer
    reluLayer
	
    maxPooling2dLayer(2,Stride=2)
		
    fullyConnectedLayer(10)
    softmaxLayer
    classificationLayer
    ];

% train the network
options = trainingOptions("sgdm", ...
    InitialLearnRate=0.01, ...
    MaxEpochs=8, ...
    Plots="training-progress");    

netb = trainNetwork(imdsTrain,layers,options);

% Check accuracy
YPred = classify(netb,imdsValidation);
YValidation = imdsValidation.Labels;
accuracy = sum(YPred == YValidation)/numel(YValidation);

tch = 0;
