function Y = Select_Best_Buffer_Approximant(lambda2, nu2, epsilon, gamma)
    %
    % Uses a pre-trained shallow net with 13 hidden nodes to identify the 
    % most accurate buffer approximant for the given values of model parameters
    %
    % Arguments are the 4 model parameters: lambda2, nu2, epsilon, gamma
    %
    % Output Y is an integer specifying the best approximant:
    % Y=1: RBA; Y=2: ExpPadeA; Y=3: ExpPadeB; Y=4: ExpExp.
    %
    %          Victor Matveev and Yinbo Chen (NJIT) March 2021
    %  (GPL) Free to use, share, modify, with appropriate attribution
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    X = log( [gamma ; epsilon; lambda2; nu2] ) / log(10);
    
    % Input offest and gain
    Offset = [ 1.301;  1.301;  3;   3];
    Gain   = [ 1.537; 1.537; 0.5; 1/3];

    % Layer 1 biases and weights
    b1 = [ 0.522; 1.657; 0.1008; 0.3276; 0.4301; -1.382; 0.006693; 1.88; -0.2144; 0.1921; -0.006686; 4.175; 3.225];

    W1 = [ 0.001029 0.79 -0.01037 0.9306; -0.07892 -0.1066 -1.133 -1.959; ...
           0.01768 -0.0984 0.3685 0.485; -0.06547 -0.1331 0.2821 -0.2874; ...
           -0.02164 0.002288 -0.5498 0.87; -0.04598 1.489 0.134 0.73;     ...
            0.06377 -0.1128 0.1973 0.6311; -0.1308 -0.09162 -1.062 -0.9823; ...
            0.1002 -0.2121 0.8472 -0.1529; -0.06318 0.915 -0.5073 0.7691; ...
            -0.061 1.029 -1.171 1.282; 0.1637 -1.409 2.978 1.638; ...
            0.05212 -0.4442 -1.535 -2.208];
        
    % Layer 2 biases and weights
    b2 = [ 760.5; -1729; 236.1; 732.4];
    
    W2 = [ -245.1 245.5 380.1 -247.3 116.3 -16.88 -425.2 -607.2 40.26 139 20.41 -62.52 -303.8; ...
           -129.8 -338 -527.9 281 236.5 34.02 310.8 802.9 284.6 321.2 -141.6 -17.33 1271; ...
            327.9 6.26 -907.6 612.3 -291.4 33.42 1098 195.9 -218.4 -219 44.29 51 -605.9; ...
            46.67 85.69 1056 -645.5 -61.77 -50.9 -983.6 -394.5 -106.6 -242.5 76.61 28.41 -363];

        X  = Gain .* ( X + Offset ) - 1; % Input -> Layer 1 affine transform
        Z  = tanh( b1 + W1 * X );        % Layer 1 Transform
    [~, Y] =  max( b2 + W2 * Z );        % Layer 2 Transform (index of max = best method)
end
