function [X, H, C, Ws] = internal_lstm(X, h0, c0, W, R, b, stateActivation, gateActivation, isTraining)
% This function is undocumented and will change in a future release.

% The help for this method is similar to the CPU method internal_lstm.
%
% This method also outputs Ws, the cuDNN workspace that is required for the
% computation of the backward.

%   Copyright 2019-2023 The MathWorks, Inc.

if ~isreal(X) || ~isreal(h0) || ~isreal(c0) || ~isreal(W) || ~isreal(R) || ~isreal(b)
    error(message('deep:dlarray:ComplexNotSupported'))
end

if nargin < 9
    isTraining = false;
end

if iHasDefaultActivations(stateActivation, gateActivation)
    % Ensure all inputs have the same data type
    [X, h0, c0, W, R, b] = nnet.internal.cnngpu.castToSuperiorFloat(X, h0, c0, W, R, b);

    % Convert to cuDNN weights format
    Wc = nnet.internal.cnngpu.util.LSTMWeightsConverter.toCudnn(W, R, b);

    % Call cuDNN method
    if isTraining
        [X, H, C, Ws] = nnet.internal.cnngpu.lstmForwardTrain(X, h0, c0, Wc);
    else
        [X, H, C] = nnet.internal.cnngpu.lstmForwardPredict(X, h0, c0, Wc);
        Ws = [];        
    end
else
    learnable.W = W;
    learnable.R = R;
    learnable.b = b;
    state.c0 = c0;
    state.y0 = h0;
    options.ReturnLast = false;
    options.Activation = iGetActivationFcn(stateActivation);
    options.RecurrentActivation = iGetRecurrentActivationFcn(gateActivation);
    [X, H, C] = nnet.internal.cnngpu.lstmForwardGeneral(X, learnable, state, options);
    Ws = [];
end
end

function tf = iHasDefaultActivations(stateActivationFunction, gateActivationFunction)
tf = (stateActivationFunction == "tanh") && (gateActivationFunction == "sigmoid");
end

function fcn = iGetActivationFcn(fcnStr)
if fcnStr == "tanh"
    fcn = @tanh;
elseif fcnStr == "relu"
    fcn = @iRelu;
else
    fcn = @iSoftSign;
end
end

function fcn = iGetRecurrentActivationFcn(fcnStr)
if fcnStr == "sigmoid"
    fcn = @iSigmoid;
else
    fcn = @nnet.internal.cnnhost.hardSigmoidForward;
end
end

function y = iSigmoid(x)
y = 1./(1 + exp(-x));
end

function y = iRelu(x)
y = max(0,x);
end

function y = iSoftSign(x)
y = x./(1 + abs(x));
end