function Y = Select_Best_Ca_Approximant(lambda2, nu2, epsilon, gamma)
    %
    % Uses a pre-trained shallow net with 13 hidden nodes to identify the 
    % most accurate Ca2+ 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.9247;  -1.672;  -1.933;  0.311;  -0.06191;  -0.8043;  2.186;  -1.886;  2.72];

    W1 = [ 0.1058 -1.12 0.3318 0.2996;  0.4088 -0.486 1.684 -0.2093; ...
           0.104 0.6083 1.476 -0.2241;  -0.08185 -0.1843 -0.08538 -2.121; ...
           0.04953 -1.097 0.2677 -1.012;  -0.1052 0.9986 -0.4409 -0.1182; ...
          -0.1082 -0.1143 -0.5137 -2.203;  0.00209 1.142 0.3669 1.564; ...
          -0.1677 -0.3878 -1.074 -3.092];

    % Layer 2 biases and weights
    b2 = [ -600.1;  183.9;  253.9;  160.1];
    
    W2 = [ 321.5 -43.95 -19.91 -35.14 56.8 365 582.2 73.38 49.38;  ...
          -337.8 35.92 -131.8 43.13 -52.01 -402 -476.3 71.82 224.8; ...
          -39.95 -10.93 91.5 0.2355 -11.75 -41.87 -152.8 -98.19 -117.9; ...
           56.58 17.95 61.06 -9.99 8.303 76.36 48.51 -47.53 -158.9];

        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 (argmax = best method)
end
