#ifndef TMVA_SOFIE_ROPERATOR_LSTM_I
#define TMVA_SOFIE_ROPERATOR_LSTM_I

namespace TMVA {
namespace Experimental {
namespace SOFIE {

template<typename T>
auto ROperator_LSTM<T>::TypeInference(std::vector<ETensorType> input)
-> std::vector<ETensorType> {
	ETensorType out = input[0];
	return {out, out};
}

template<typename T>
auto ROperator_LSTM<T>::ShapeInference(std::vector<std::vector<size_t>> input)
-> std::vector<std::vector<size_t>> {
	size_t num_directions = input[1][0];
	size_t hidden_size = input[1][1] / 4;
	if (fAttrLayout == 0) {
		size_t seq_length = input[0][0];
		size_t batch_size = input[0][1];
		std::vector<std::vector<size_t>> ret(
			{{seq_length, num_directions, batch_size, hidden_size},
			 {num_directions, batch_size, hidden_size},
			 {num_directions, batch_size, hidden_size}});
		return ret;
	} else {
		size_t batch_size = input[0][0];
		size_t seq_length = input[0][1];
		std::vector<std::vector<size_t>> ret(
			{{batch_size, seq_length, num_directions, hidden_size},
			 {batch_size, num_directions, hidden_size},
			 {batch_size, num_directions, hidden_size}});
		return ret;
	}
}

template<typename T>
auto ROperator_LSTM<T>::Initialize(RModel &model)
-> void {
   fUseSession = model.UseSession();
   // Check the input and output tensors
   if (!model.CheckIfTensorAlreadyExist(fNX)) {
		throw std::runtime_error("TMVA SOFIE LSTM Op input tensor " + fNX + "  is not found in model.");
	}
	fShapeX = model.GetTensorShape(fNX);
	if (fShapeX.size() != 3) {
		throw std::runtime_error("TMVA SOFIE LSTM Op input tensor " + fNX + " is not of 3 dimensions.");
	}
	if (!model.CheckIfTensorAlreadyExist(fNW)) {
		throw std::runtime_error("TMVA SOFIE LSTM Op input tensor " + fNW + "  is not found in model.");
	}
	fShapeW = model.GetTensorShape(fNW);
	if (fShapeW.size() != 3) {
		throw std::runtime_error("TMVA SOFIE LSTM Op input tensor " + fNW + " is not of 3 dimensions.");
	}
	if (!model.CheckIfTensorAlreadyExist(fNR)) {
		throw std::runtime_error("TMVA SOFIE LSTM Op input tensor " + fNR + "  is not found in model.");
	}
	fShapeR = model.GetTensorShape(fNR);
	if (fShapeR.size() != 3) {
		throw std::runtime_error("TMVA SOFIE LSTM Op input tensor " + fNR + " is not of 3 dimensions.");
	}
	if (!fNB.empty()) {
		if (!model.CheckIfTensorAlreadyExist(fNB)) {
			throw std::runtime_error("TMVA SOFIE LSTM op input tensor " + fNB + " is not  found in model.");
		}
		fShapeB = model.GetTensorShape(fNB);
		if (fShapeB.size() != 2 && fShapeB.size() != 5) {
			throw std::runtime_error("TMVA SOFIE LSTM op input tensor " + fNB + " is not of 2 or 5 dimensions.");
		}
		if (fShapeB.size() == 2) {
			// Broadcasting the bias
			auto original_data = model.GetInitializedTensorData(fNB);
			size_t num_directions = fShapeW[0];
			size_t seq_length = (fAttrLayout == 0)? fShapeX[0] : fShapeX[1];
			size_t batch_size = (fAttrLayout == 0)? fShapeX[1] : fShapeX[0];
			if (fType == "float") {
				float *original_bias = static_cast<float*>(original_data.get());
				float *new_bias = new float[4 * num_directions * seq_length * batch_size * fAttrHiddenSize];
				for (size_t gate = 0; gate < 4; gate++) {
					float sum[fAttrHiddenSize];
					for (size_t direction = 0; direction < num_directions; direction++) {
						size_t offset = direction * 8 * fAttrHiddenSize + gate * fAttrHiddenSize;
						for (size_t h = 0; h < fAttrHiddenSize; h++) {
							sum[h] = original_bias[offset + h] + original_bias[offset + h + 4 * fAttrHiddenSize];
						}
						for (size_t seq = 0; seq < seq_length; seq++) {
							for (size_t batch = 0; batch < batch_size; batch++) {
								size_t bias_offset = gate * num_directions * seq_length * batch_size * fAttrHiddenSize
									+ direction * seq_length * batch_size * fAttrHiddenSize
									+ seq * batch_size * fAttrHiddenSize + batch * fAttrHiddenSize;
								std::copy(sum, sum + fAttrHiddenSize, new_bias + bias_offset);
							}
						}
					}
				}
				std::vector<size_t> new_bias_shape = {4, num_directions, seq_length, batch_size, fAttrHiddenSize};
				std::shared_ptr<void> new_bias_ptr(new_bias, std::default_delete<float[]>());
				model.UpdateInitializedTensor(fNB, model.GetTensorType(fNB), new_bias_shape, new_bias_ptr);
				fShapeB = model.GetTensorShape(fNB);
			}
		}
	}
	if (!fNSequence_lens.empty()) {
		if (!model.CheckIfTensorAlreadyExist(fNSequence_lens)) {
			throw std::runtime_error("TMVA SOFIE LSTM Op input tensor " +
                                  fNSequence_lens +
                                  "is not found in model.");
      }
      fShapeSequence_lens = model.GetTensorShape(fNSequence_lens);
      if (fShapeSequence_lens.size() != 1) {
         throw std::runtime_error("TMVA SOFIE LSTM Op input tensor " +
                                  fNSequence_lens +
                                  " is not of 1 dimension.");
      }
   }
   if (!fNInitial_h.empty()) {
      if (!model.CheckIfTensorAlreadyExist(fNInitial_h)) {
        throw std::runtime_error("TMVA SOFIE LSTM Op input tensor " +
                                 fNInitial_h + " is not found in model.");
      }
      fShapeInitial_h = model.GetTensorShape(fNInitial_h);
      if (fShapeInitial_h.size() != 3) {
        throw std::runtime_error("TMVA SOFIE LSTM Op input tensor " +
                                 fNInitial_h + " is not of 3 dimensions.");
      }
   }
   if (!fNInitial_c.empty()) {
      if (!model.CheckIfTensorAlreadyExist(fNInitial_c)) {
         throw std::runtime_error("TMVA SOFIE LSTM Op input tensor " +
                                  fNInitial_c + " is not found in model.");
      }
      fShapeInitial_c = model.GetTensorShape(fNInitial_c);
      if (fShapeInitial_c.size() != 3) {
         throw std::runtime_error("TMVA SOFIE LSTM Op input tensor " +
                                  fNInitial_c + " is not of 3 dimensions.");
      }
   }
   if (!fNP.empty()) {
      if (!model.CheckIfTensorAlreadyExist(fNP)) {
         throw std::runtime_error("TMVA SOFIE LSTM op input tensor " + fNP + " is not  found in model.");
      }
      fShapeP = model.GetTensorShape(fNP);
      if (fShapeP.size() != 2 && fShapeP.size() != 4) {
         throw std::runtime_error("TMVA SOFIE LSTM op input tensor " + fNP + " is not of 2 or 4 dimensions.");
      }
      if (fShapeP.size() == 2) {
         // Broadcasting the weight for peepholes
         auto original_data = model.GetInitializedTensorData(fNP);
         size_t num_directions = fShapeW[0];
         size_t batch_size = (fAttrLayout == 0)? fShapeX[1] : fShapeX[0];
         if (fType == "float") {
            float *original_p = static_cast<float*>(original_data.get());
            float *new_p = new float[num_directions * 3 * batch_size * fAttrHiddenSize];
            for (size_t direction = 0; direction < num_directions; direction++) {
               for (size_t gate = 0; gate < 3; gate++) {
                  size_t p_offset = direction * 3 * fAttrHiddenSize + gate * fAttrHiddenSize;
                  for (size_t batch = 0; batch < batch_size; batch++) {
                     size_t offset = direction * 3 * batch_size * fAttrHiddenSize
                        + gate * batch_size * fAttrHiddenSize + batch * fAttrHiddenSize;
                     std::copy(original_p + p_offset, original_p + p_offset + fAttrHiddenSize,
                        new_p + offset);
                  }
               }
            }
            std::vector<size_t> new_p_shape = {num_directions, 3, batch_size, fAttrHiddenSize};
            std::shared_ptr<void> new_p_ptr(new_p, std::default_delete<float[]>());
            model.UpdateInitializedTensor(fNP, model.GetTensorType(fNP), new_p_shape, new_p_ptr);
            fShapeP = model.GetTensorShape(fNP);
         }
      }
   }
	if (!fNY.empty()) {
      fShapeY = ShapeInference({fShapeX, fShapeW})[0];
      if (!model.CheckIfTensorAlreadyExist(fNY)) {
         model.AddIntermediateTensor(fNY, model.GetTensorType(fNX), fShapeY);
      }
	}
	if (!fNY_h.empty()) {
      fShapeY_h = ShapeInference({fShapeX, fShapeW})[1];
      if (!model.CheckIfTensorAlreadyExist(fNY_h)) {
         model.AddIntermediateTensor(fNY_h, model.GetTensorType(fNX), fShapeY_h);
      }
	}
	if (!fNY_c.empty()) {
      fShapeY_c = ShapeInference({fShapeX, fShapeW})[2];
      if (!model.CheckIfTensorAlreadyExist(fNY_c)) {
         model.AddIntermediateTensor(fNY_c, model.GetTensorType(fNX), fShapeY_c);
      }
	}
	// Check the attributes
	for (auto &activation : fAttrActivations) {
     if (activation != "Relu" && activation != "Tanh" &&
          activation != "Sigmoid" && activation != "Affine" &&
          activation != "LeakyRelu" && activation != "ThresholdRelu" &&
          activation != "ScaledTanh" && activation != "HardSigmoid" &&
          activation != "Elu" && activation != "Softsign" &&
          activation != "Softplus") {
         throw std::runtime_error("TMVA SOFIE - Activation function " +
                                 activation + " not implemented");
      }
	}
   if (fAttrDirection != "forward" && fAttrDirection != "backward" &&
       fAttrDirection != "bidirectional") {
      throw std::runtime_error(
          "TMVA SOFIE - Invalid LSTM direction fAttrDirection = " +
          fAttrDirection);
   }
   if (4 * fAttrHiddenSize != fShapeW[1]) {
      throw std::runtime_error(
          "TMVA SOFIE - fAttrHiddenSize must be equal to " +
          std::to_string(fShapeW[1] / 4));
   }
   if (fAttrInputForget > 1) {
      throw std::runtime_error(
         "TMVA SOFIE - fAttrInputForget = " + std::to_string(fAttrInputForget)
         + " must be 0 or 1.");
   }
   if (fAttrLayout > 1) {
      throw std::runtime_error("TMVA SOFIE - Layout fAttrLayout = " +
                               std::to_string(fAttrLayout) +
                               " must be 0 (timewise) or 1 (batchwise)");
   }
	if (fAttrActivations.empty()) {
      if (fAttrDirection == "bidirectional") {
         fAttrActivations = {"Sigmoid", "Tanh", "Tanh", "Sigmoid", "Tanh", "Tanh"};
      } else {
         fAttrActivations = {"Sigmoid", "Tanh", "Tanh"};
      }
	}
}

// generate code for Session data members (e.g. internal vectors)
template <typename T>
std::string ROperator_LSTM<T>::GenerateSessionMembersCode(std::string opName)
{
   opName = "op_" + opName;
   std::stringstream out;

   size_t num_directions = fShapeW[0];
   size_t seq_length = (fAttrLayout == 0) ? fShapeX[0] : fShapeX[1];
   size_t batch_size = (fAttrLayout == 0) ? fShapeX[1] : fShapeX[0];
   size_t input_size = fShapeX[2];

   if (fAttrLayout != 0) {
      out << "std::vector<" << fType << "> fVec_" << opName << "_input = std::vector<" << fType << ">("
          << seq_length * batch_size * input_size << ");\n";
      out << "std::vector<" << fType << "> fVec_" << opName << "_initial_hidden_state = std::vector<" << fType << ">("
          << num_directions * batch_size * fAttrHiddenSize << ");\n";
      out << "std::vector<" << fType << "> fVec_" << opName << "_initial_cell_state = std::vector<" << fType << ">("
          << num_directions * batch_size * fAttrHiddenSize << ");\n";
   }
   // Set the feedforward
   size_t ff_size = seq_length * batch_size * fAttrHiddenSize;
   out << "std::vector<" << fType << "> fVec_" << opName << "_ff_input_gate = std::vector<" << fType << ">(" << ff_size << ");\n";
   out << "std::vector<" << fType << "> fVec_" << opName << "_ff_output_gate = std::vector<" << fType << ">(" << ff_size << ");\n";
   out << "std::vector<" << fType << "> fVec_" << opName << "_ff_cell_gate = std::vector<" << fType << ">(" << ff_size << ");\n";
   if (fAttrInputForget == 0)
      out << "std::vector<" << fType << "> fVec_" << opName << "_ff_forget_gate = std::vector<" << fType << ">(" << ff_size << ");\n";
   // gate results
   size_t hs_size = seq_length * num_directions * batch_size * fAttrHiddenSize;
   out << "std::vector<" << fType << "> fVec_" << opName << "_input_gate = std::vector<" << fType << ">(" << hs_size << ");\n";
   out << "std::vector<" << fType << "> fVec_" << opName << "_output_gate = std::vector<" << fType << ">(" << hs_size << ");\n";
   out << "std::vector<" << fType << "> fVec_" << opName << "_cell_gate = std::vector<" << fType << ">(" << hs_size << ");\n";
   if (fAttrInputForget == 0)
      out << "std::vector<" << fType << "> fVec_" << opName << "_forget_gate = std::vector<" << fType << ">(" << hs_size << ");\n";
   // cell state     
   out << "std::vector<" << fType << "> fVec_" << opName << "_cell_state = std::vector<" << fType << ">(" << hs_size << ");\n";
   out << "std::vector<" << fType << "> fVec_" << opName << "_new_cell_state = std::vector<" << fType << ">(" << hs_size << ");\n";
   // hiddden state
   if (fAttrLayout != 0 || fNY.empty()) {
      out << "std::vector<" << fType << "> fVec_" << opName << "_hidden_state = std::vector<" << fType << ">(" << hs_size << ");\n";
   }

   out << "\n";

   return out.str();
}

template<typename T>
auto ROperator_LSTM<T>::Generate(std::string OpName)
-> std::string {
	OpName = "op_" + OpName;
	std::stringstream out;

	size_t seq_length = (fAttrLayout == 0) ? fShapeX[0] : fShapeX[1];
	size_t batch_size = (fAttrLayout == 0) ? fShapeX[1] : fShapeX[0];
	size_t input_size = fShapeX[2];
	size_t num_directions = fShapeW[0];

   // set the input
   if (fAttrLayout == 0) {
      out << SP << fType << " *" << OpName << "_input = tensor_" << fNX << ";\n";
   } else {
      if (fUseSession) 
         out << SP << fType << " * " << OpName << "_input = fVec_" << OpName << "_input.data();\n";
      else
         out << SP << fType << "  " << OpName << "_input[" << seq_length * batch_size * input_size << "] = {0};\n";

      out << SP << "for(size_t seq = 0; seq < " << seq_length << "; seq++) {\n";
      out << SP << SP << "for(size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
      out << SP << SP << SP << "for(size_t i = 0; i < " << input_size << "; i++) {\n";
      out << SP << SP << SP << SP << OpName << "_input[seq * " << batch_size * input_size
         << " + batch * " << input_size << " + i] = " << "tensor_" << fNX << "[batch * "
         << seq_length * input_size << " + seq * " << input_size << " + i];\n";
      out << SP << SP << SP << "}\n";
      out << SP << SP << "}\n";
      out << SP << "}\n";
   }

   // Set the initial hidden state
   if (!fNInitial_h.empty()) {
      if (fAttrLayout == 0) {
         out << SP << fType << " *" << OpName << "_initial_hidden_state = " << " tensor_"
                << fNInitial_h << ";\n";
      } else {
         if (fUseSession)
            out << SP << fType << " * " << OpName << "_initial_hidden_state = fVec_" << OpName
                   << "_initial_hidden_state.data();\n";
         else
            out << SP << fType << "  " << OpName << "_initial_hidden_state[" << num_directions * batch_size *
                fAttrHiddenSize << "] = {0};\n";

         for (size_t direction = 0; direction < num_directions; direction++) {
            out << SP << "for(size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
            out << SP << SP << "for(size_t h = 0; h < " << fAttrHiddenSize << "; h++) {\n";
            out << SP << SP << SP << OpName << "_initial_hidden_state["
                << direction * batch_size * fAttrHiddenSize << " + batch * " << fAttrHiddenSize
                << " + h] = tensor_" << fNInitial_h << "[batch * " << num_directions * fAttrHiddenSize
                << " + " << direction * fAttrHiddenSize << " + h];\n";
            out << SP << SP << "}\n";
            out << SP << "}\n";
         }
      }
   }

   // Set the initial cell state
   if (!fNInitial_c.empty()) {
      if (fAttrLayout == 0) {
         out << SP << fType << " *" << OpName << "_initial_cell_state = " << " tensor_"
                << fNInitial_c << ";\n";
      } else {
         if (fUseSession)
            out << SP << fType << " * " << OpName << "_initial_cell_state = fVec_" << OpName
                << "_initial_cell_state.data();\n";
         else
            out << SP << fType << "  " << OpName << "_initial_cell_state[" << num_directions * batch_size *
                fAttrHiddenSize << "] = {0};\n";

         for (size_t direction = 0; direction < num_directions; direction++) {
            out << SP << "for(size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
            out << SP << SP << "for(size_t h = 0; h < " << fAttrHiddenSize << "; h++) {\n";
            out << SP << SP << SP << OpName << "_initial_cell_state["
                << direction * batch_size * fAttrHiddenSize << " + batch * " << fAttrHiddenSize
                << " + h] = tensor_" << fNInitial_c << "[batch * " << num_directions * fAttrHiddenSize
                << " + " << direction * fAttrHiddenSize << " + h];\n";
            out << SP << SP << "}\n";
            out << SP << "}\n";
         }
      }
   }

   // Set the feedforward
   size_t ff_size = seq_length * batch_size * fAttrHiddenSize;
   if (fUseSession) {
      out << SP << fType << " * " << OpName << "_ff_input_gate = fVec_" << OpName << "_ff_input_gate.data();\n";
      out << SP << fType << " * " << OpName << "_ff_output_gate = fVec_" << OpName << "_ff_output_gate.data();\n";
      out << SP << fType << " * " << OpName << "_ff_cell_gate = fVec_" << OpName << "_ff_cell_gate.data();\n";
      if (fAttrInputForget == 0) {
         out << SP << fType << " * " << OpName << "_ff_forget_gate = fVec_" << OpName << "_ff_forget_gate.data();\n";
      }
   } else {
      out << SP << fType << "  " << OpName << "_ff_input_gate[" << ff_size << "] = {0};\n";
      out << SP << fType << "  " << OpName << "_ff_output_gate[" << ff_size << "] = {0};\n";
      out << SP << fType << "  " << OpName << "_ff_cell_gate[" << ff_size << "] = {0};\n";
      if (fAttrInputForget == 0) {
         out << SP << fType << "  " << OpName << "_ff_forget_gate[" << ff_size << "] = {0};\n";
      }
   }
   // Set the gates
   size_t hidden_state_size = seq_length * num_directions * batch_size * fAttrHiddenSize;
   if (fUseSession) {
      out << SP << fType << " * " << OpName << "_input_gate = fVec_" << OpName << "_input_gate.data();\n";
      out << SP << fType << " * " << OpName << "_output_gate = fVec_" << OpName << "_output_gate.data();\n";
      out << SP << fType << " * " << OpName << "_cell_gate = fVec_" << OpName << "_cell_gate.data();\n";
      if (fAttrInputForget == 0) {
         out << SP << fType << " * " << OpName << "_forget_gate = fVec_" << OpName << "_forget_gate.data();\n";
      }
   } else {
      out << SP << fType << "  " << OpName << "_input_gate[" << hidden_state_size << "] = {0};\n";
      out << SP << fType << "  " << OpName << "_output_gate[" << hidden_state_size << "] = {0};\n";
      out << SP << fType << "  " << OpName << "_cell_gate[" << hidden_state_size << "] = {0};\n";
      if (fAttrInputForget == 0) {
         out << SP << fType << "  " << OpName << "_forget_gate[" << hidden_state_size << "] = {0};\n";
      }
   }
   // Set the cell state and the new cell state = h(cell state)
   if (fUseSession) {
      out << SP << fType << " * " << OpName << "_cell_state = fVec_" << OpName << "_cell_state.data();\n";
      out << SP << fType << " * " << OpName << "_new_cell_state = fVec_" << OpName << "_new_cell_state.data();\n";
   } else {
      out << SP << fType << "  " << OpName << "_cell_state[" << hidden_state_size << "] = {0};\n";
      out << SP << fType << "  " << OpName << "_new_cell_state[" << hidden_state_size << "] = {0};\n";
   }

   // Set the hidden state
   if (fAttrLayout == 0 && !fNY.empty()) {
      out << SP << fType << " *" << OpName << "_hidden_state = tensor_" << fNY << ";\n";
   } else {
      if (fUseSession) {
         out << SP << fType << " * " << OpName << "_hidden_state = fVec_" << OpName << "_hidden_state.data();\n";
      } else {
         out << SP << fType << "  " << OpName << "_hidden_state[" << hidden_state_size << "] = {0};\n";
      }
   }

   out << SP << "char " << OpName << "_transA = 'N';\n";
   out << SP << "char " << OpName << "_transB = 'T';\n";
   out << SP << "int " << OpName << "_m = " << seq_length * batch_size << ";\n";
   out << SP << "int " << OpName << "_n = " << fAttrHiddenSize << ";\n";
   out << SP << "int " << OpName << "_k = " << input_size << ";\n";
   if (fType == "float") {
      out << SP << fType << "  " << OpName << "_alpha = 1.;\n";
      out << SP << fType << "  " << OpName << "_beta = 0.;\n";
   }
   if (!fNB.empty()) {
      out << SP << "int " << OpName << "_bias_size = " << seq_length * batch_size * fAttrHiddenSize << ";\n";
      out << SP << "int " << OpName << "_incx = 1;\n";
      out << SP << "int " << OpName << "_incy = 1;\n";
   }

   for (size_t direction = 0; direction < num_directions; direction++) {
      if (direction == 0) {
         if (fType == "float") {
            // input_gate = input * weight_i^T
            out << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                << OpName <<"_n, &" << OpName << "_m, &" << OpName << "_k, &" << OpName << "_alpha, tensor_"
                << fNW << ", &" << OpName << "_k, " << OpName << "_input, &" << OpName << "_k, &"
               << OpName << "_beta, " << OpName << "_ff_input_gate, &" << OpName << "_n);\n";
            // output_gate = input * weight_o^T
            size_t wo_offset = fAttrHiddenSize * input_size;
            out << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                << OpName <<"_n, &" << OpName << "_m, &" << OpName << "_k, &" << OpName << "_alpha, tensor_"
               << fNW << " + " << wo_offset << ", &" << OpName << "_k, " << OpName << "_input, &"
               << OpName << "_k, &" << OpName << "_beta, " << OpName << "_ff_output_gate, &" << OpName << "_n);\n";
            // cell_gate = input * weight_c^T
            size_t wc_offset = 3 * fAttrHiddenSize * input_size;
            out << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                << OpName <<"_n, &" << OpName << "_m, &" << OpName << "_k, &" << OpName << "_alpha, tensor_"
               << fNW << " + " << wc_offset << ", &" << OpName << "_k, " << OpName << "_input, &"
               << OpName << "_k, &" << OpName << "_beta, " << OpName << "_ff_cell_gate, &" << OpName << "_n);\n";
         }
      } else {
         if (fType == "float") {
            // input_gate = input * weight_i^T
            size_t wi_offset = 4 * fAttrHiddenSize * input_size;
            out << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                << OpName <<"_n, &" << OpName << "_m, &" << OpName << "_k, &" << OpName << "_alpha, tensor_"
               << fNW << " + " << wi_offset << ", &" << OpName << "_k, " << OpName << "_input, &"
               << OpName << "_k, &" << OpName << "_beta, " << OpName << "_ff_input_gate, &" << OpName << "_n);\n";
            // output_gate = input * weight_o^T
            size_t wo_offset = 4 * fAttrHiddenSize * input_size + 1 * fAttrHiddenSize * input_size;
            out << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                << OpName <<"_n, &" << OpName << "_m, &" << OpName << "_k, &" << OpName << "_alpha, tensor_"
               << fNW << " + " << wo_offset << ", &" << OpName << "_k, " << OpName << "_input, &"
               << OpName << "_k, &" << OpName << "_beta, " << OpName << "_ff_output_gate, &" << OpName << "_n);\n";
            // cell_gate = input * weight_c^T
            size_t wc_offset = 4 * fAttrHiddenSize * input_size + 3 * fAttrHiddenSize * input_size;
            out << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                << OpName <<"_n, &" << OpName << "_m, &" << OpName << "_k, &" << OpName << "_alpha, tensor_"
               << fNW << " + " << wc_offset << ", &" << OpName << "_k, " << OpName << "_input, &"
               << OpName << "_k, &" << OpName << "_beta, " << OpName << "_ff_cell_gate, &" << OpName << "_n);\n";
         }
      }
      if (fAttrInputForget == 0) {
         // forget_gate = input * weight_f^T
         if (direction == 0) {
            if (fType == "float") {
               size_t wf_offset = 2 * fAttrHiddenSize * input_size;
               out << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                   << OpName <<"_n, &" << OpName << "_m, &" << OpName << "_k, &" << OpName << "_alpha, tensor_"
                   << fNW << " + " << wf_offset << ", &" << OpName << "_k, " << OpName << "_input, &"
                   << OpName << "_k, &" << OpName << "_beta, " << OpName << "_ff_forget_gate, &" << OpName << "_n);\n";
            }
         } else {
            if (fType == "float") {
               size_t wf_offset = 4 * fAttrHiddenSize * input_size + 2 * fAttrHiddenSize * input_size;
               out << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                   << OpName <<"_n, &" << OpName << "_m, &" << OpName << "_k, &" << OpName << "_alpha, tensor_"
                   << fNW << " + " << wf_offset << ", &" << OpName << "_k, " << OpName << "_input, &"
                   << OpName << "_k, &" << OpName << "_beta, " << OpName << "_ff_forget_gate, &" << OpName << "_n);\n";
            }
         }
      }

      // Add the bias
      if (!fNB.empty()) {
         if (direction == 0) {
            if (fType == "float") {
               // ff_input_gate += bias_i
               out << SP << "BLAS::saxpy_(&" << OpName << "_bias_size, &" << OpName << "_alpha, tensor_"
                   << fNB << ", &" << OpName << "_incx, " << OpName << "_ff_input_gate, &" << OpName << "_incy);\n";
               // ff_output_gate += bias_o
               size_t bo_offset =  seq_length * batch_size * fAttrHiddenSize;
               out << SP << "BLAS::saxpy_(&" << OpName << "_bias_size, &" << OpName << "_alpha, tensor_"
                   << fNB << " + " << bo_offset << ", &" << OpName << "_incx, " << OpName << "_ff_output_gate, &"
                   << OpName << "_incy);\n";
               // ff_cell_gate += bias_c
               size_t bc_offset = 3 * seq_length * batch_size * fAttrHiddenSize;
               out << SP << "BLAS::saxpy_(&" << OpName << "_bias_size, &" << OpName << "_alpha, tensor_"
                   << fNB << " + " << bc_offset << ", &" << OpName << "_incx, " << OpName << "_ff_cell_gate, &"
                   << OpName << "_incy);\n";
            }
         } else {
            if (fType == "float") {
               // ff_input_gate += bias_i
               size_t bi_offset = 4 * seq_length * batch_size * fAttrHiddenSize;
               out << SP << "BLAS::saxpy_(&" << OpName << "_bias_size, &" << OpName << "_alpha, tensor_"
                   << fNB << " + " << bi_offset << ", &" << OpName << "_incx, " << OpName << "_ff_input_gate, &"
                   << OpName << "_incy);\n";
               // ff_output_gate += bias_o
               size_t bo_offset = 4 * seq_length * batch_size * fAttrHiddenSize
                  + seq_length * batch_size * fAttrHiddenSize;
               out << SP << "BLAS::saxpy_(&" << OpName << "_bias_size, &" << OpName << "_alpha, tensor_"
                   << fNB << " + " << bo_offset << ", &" << OpName << "_incx, " << OpName << "_ff_output_gate, &"
                   << OpName << "_incy);\n";
               // ff_cell_gate += bias_c
               size_t bc_offset = 4 * num_directions * seq_length * batch_size * fAttrHiddenSize
                  + 3 * seq_length * batch_size * fAttrHiddenSize;
               out << SP << "BLAS::saxpy_(&" << OpName << "_bias_size, &" << OpName << "_alpha, tensor_"
                   << fNB << " + " << bc_offset << ", &" << OpName << "_incx, " << OpName << "_ff_cell_gate, &"
                   << OpName << "_incy);\n";
            }
         }
         if (fAttrInputForget == 0) {
            // ff_forget_gate += bias_f
            if (direction == 0) {
               if (fType == "float") {
                  size_t bo_offset = 2 * seq_length * batch_size * fAttrHiddenSize;
                  out << SP << "BLAS::saxpy_(&" << OpName << "_bias_size, &" << OpName << "_alpha, tensor_"
                      << fNB << " + " << bo_offset << ", &" << OpName << "_incx, " << OpName << "_ff_forget_gate, &"
                      << OpName << "_incy);\n";
               }
            } else {
               if (fType == "float") {
                  size_t bo_offset = 4 * seq_length * batch_size * fAttrHiddenSize
                     + 2 * seq_length * batch_size * fAttrHiddenSize;
                  out << SP << "BLAS::saxpy_(&" << OpName << "_bias_size, &" << OpName << "_alpha, tensor_"
                      << fNB << " + " << bo_offset << ", &" << OpName << "_incx, " << OpName << "_ff_forget_gate, &"
                      << OpName << "_incy);\n";
               }
            }
         }
      }

      // Copy ff_input_gate, ff_output_gate, ff_cell_gate and ff_forget_gate into input_gate, output_gate,
      //   cell_gate and forget_gate
      out << SP << "for (size_t seq = 0; seq < " << seq_length << "; seq++) {\n";
      out << SP << SP << "size_t ff_offset = seq * " << batch_size * fAttrHiddenSize << ";\n";
      if (direction == 0) {
         out << SP << SP << "size_t gate_offset = seq * " << num_directions * batch_size * fAttrHiddenSize
            << ";\n";
      } else {
         out << SP << SP << "size_t gate_offset = seq * " << num_directions * batch_size * fAttrHiddenSize
             << " + " << batch_size * fAttrHiddenSize << ";\n";
      }
      size_t ff_seq_size = batch_size * fAttrHiddenSize;
      out << SP << SP << "std::copy(" << OpName << "_ff_input_gate + ff_offset, " << OpName
          << "_ff_input_gate + ff_offset + " << ff_seq_size << ", " << OpName << "_input_gate + gate_offset);\n";
      out << SP << SP << "std::copy(" << OpName << "_ff_output_gate + ff_offset, " << OpName
          << "_ff_output_gate + ff_offset + " << ff_seq_size << ", " << OpName << "_output_gate + gate_offset);\n";
      out << SP << SP << "std::copy(" << OpName << "_ff_cell_gate + ff_offset, " << OpName
          << "_ff_cell_gate + ff_offset + " << ff_seq_size << ", " << OpName << "_cell_gate + gate_offset);\n";
      if (fAttrInputForget == 0) {
         out << SP << SP << "std::copy(" << OpName << "_ff_forget_gate + ff_offset, " << OpName
             << "_ff_forget_gate + ff_offset + " << ff_seq_size << ", " << OpName << "_forget_gate + gate_offset);\n";
      }
      out << SP << "}\n";

      out << SP << "for (size_t seq = 0; seq < " << seq_length << "; seq++) {\n";
      if (fAttrDirection == "backward" || direction == 1) {
         out << SP << SP << "size_t index = " << seq_length - 1 << " - seq;\n";
      } else {
         out << SP << SP << "size_t index = seq;\n";
      }
      out << SP << SP << "int m2 = " << batch_size << ";\n";
      if (direction == 0) {
         out << SP << SP << "size_t offset = index * " << num_directions * batch_size * fAttrHiddenSize
              << ";\n";
      } else {
         out << SP << SP << "size_t offset = index * " << num_directions * batch_size * fAttrHiddenSize
             << " + " << batch_size * fAttrHiddenSize << ";\n";
      }
      size_t size = batch_size * fAttrHiddenSize;
      // gate = gate + initial_hidden_state * Recurrence^T
      out << SP << SP << "if (seq == 0) {\n";
      if (!fNInitial_h.empty()) {
         if (direction == 0) {
            if (fType == "float") {
               out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                   << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << ", &"
                   << OpName << "_n, " << OpName << "_initial_hidden_state, &" << OpName << "_n, &" << OpName
                   << "_alpha, " << OpName << "_input_gate + offset, &" << OpName << "_n);\n";
               size_t ro_offset = fAttrHiddenSize * fAttrHiddenSize;
               out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                   << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << " + "
                   << ro_offset << ", &" << OpName << "_n, " << OpName << "_initial_hidden_state, &" << OpName
                   << "_n, &" << OpName << "_alpha, " << OpName << "_output_gate + offset, &" << OpName << "_n);\n";
               size_t rc_offset = 3 * fAttrHiddenSize * fAttrHiddenSize;
               out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                   << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << " + "
                   << rc_offset << ", &" << OpName << "_n, " << OpName << "_initial_hidden_state, &" << OpName
                   << "_n, &" << OpName << "_alpha, " << OpName << "_cell_gate + offset, &" << OpName << "_n);\n";
               if (fAttrInputForget == 0) {
                  size_t rf_offset = 2 * fAttrHiddenSize * fAttrHiddenSize;
                  out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                      << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << " + "
                      << rf_offset << ", &" << OpName << "_n, " << OpName << "_initial_hidden_state, &" << OpName
                      << "_n, &" << OpName << "_alpha, " << OpName << "_forget_gate + offset, &" << OpName << "_n);\n";
               }
            }
         } else { // direction=1
            if (fType == "float") {
               size_t ri_offset = 4 * fAttrHiddenSize * fAttrHiddenSize;
               out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                   << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << " + "
                   << ri_offset << ", &" << OpName << "_n, " << OpName << "_initial_hidden_state, &" << OpName
                   << "_n, &" << OpName << "_alpha, " << OpName << "_input_gate + offset, &" << OpName << "_n);\n";
               size_t ro_offset = 4 * fAttrHiddenSize * fAttrHiddenSize + 1 * fAttrHiddenSize * fAttrHiddenSize;
               out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                   << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << " + "
                   << ro_offset << ", &" << OpName << "_n, " << OpName << "_initial_hidden_state, &" << OpName
                   << "_n, &" << OpName << "_alpha, " << OpName << "_output_gate + offset, &" << OpName << "_n);\n";
               size_t rc_offset = 4 * fAttrHiddenSize * fAttrHiddenSize + 3 * fAttrHiddenSize * fAttrHiddenSize;
               out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                   << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << " + "
                   << rc_offset << ", &" << OpName << "_n, " << OpName << "_initial_hidden_state, &" << OpName
                   << "_n, &" << OpName << "_alpha, " << OpName << "_cell_gate + offset, &" << OpName << "_n);\n";
               if (fAttrInputForget == 0) {
                  size_t rf_offset = 4 * fAttrHiddenSize * fAttrHiddenSize + 2 * fAttrHiddenSize * fAttrHiddenSize;
                  out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                      << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << " + "
                      << rf_offset << ", &" << OpName << "_n, " << OpName << "_initial_hidden_state, &" << OpName
                      << "_n, &" << OpName << "_alpha, " << OpName << "_forget_gate + offset, &" << OpName << "_n);\n";
               }
            }
         }
      }
      out << SP << SP << "} else {\n";
      // gate = gate + previous_hidden_state * Recurrence^T
      if (direction == 0) {
         if (fAttrDirection == "backward") {
            out << SP << SP << SP << "size_t previous_offset = (index + 1) * "
                << num_directions * batch_size * fAttrHiddenSize << ";\n";
         } else {
            out << SP << SP << SP << "size_t previous_offset = (seq - 1) * "
                << num_directions * batch_size * fAttrHiddenSize << ";\n";
         }
         if (fType == "float") {
            out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
             << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << ", &"
             << OpName << "_n, " << OpName << "_hidden_state + previous_offset, &" << OpName << "_n, &"
             << OpName << "_alpha, " << OpName << "_input_gate + offset, &" << OpName << "_n);\n";
            size_t ro_offset = 1 * fAttrHiddenSize * fAttrHiddenSize;
            out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
             << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << " + "
             << ro_offset << ", &" << OpName << "_n, " << OpName << "_hidden_state + previous_offset, &"
             << OpName << "_n, &" << OpName << "_alpha, " << OpName << "_output_gate + offset, &"
             << OpName << "_n);\n";
            size_t rc_offset = 3 * fAttrHiddenSize * fAttrHiddenSize;
            out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
             << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << " + "
             << rc_offset << ", &" << OpName << "_n, " << OpName << "_hidden_state + previous_offset, &"
             << OpName << "_n, &" << OpName << "_alpha, " << OpName << "_cell_gate + offset, &"
             << OpName << "_n);\n";
            if (fAttrInputForget == 0) {
               size_t rf_offset = 2 * fAttrHiddenSize * fAttrHiddenSize;
               out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                   << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << " + "
                   << rf_offset << ", &" << OpName << "_n, " << OpName << "_hidden_state + previous_offset, &"
                   << OpName << "_n, &" << OpName << "_alpha, " << OpName << "_forget_gate + offset, &"
                   << OpName << "_n);\n";
            }
         }
      } else {
         out << SP << SP << SP << "size_t previous_offset = (index + 1) * "
             << num_directions * batch_size * fAttrHiddenSize << " + " << batch_size * fAttrHiddenSize << ";\n";
         if (fType == "float") {
            size_t ri_offset = 4 * fAttrHiddenSize * fAttrHiddenSize;
            out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
             << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << " + "
             << ri_offset << ", &" << OpName << "_n, " << OpName << "_hidden_state + previous_offset, &"
             << OpName << "_n, &" << OpName << "_alpha, " << OpName << "_input_gate + offset, &"
             << OpName << "_n);\n";
            size_t ro_offset = 4 * fAttrHiddenSize * fAttrHiddenSize + fAttrHiddenSize * fAttrHiddenSize;
            out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
             << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << " + "
             << ro_offset << ", &" << OpName << "_n, " << OpName << "_hidden_state + previous_offset, &"
             << OpName << "_n, &" << OpName << "_alpha, " << OpName << "_output_gate + offset, &"
             << OpName << "_n);\n";
            size_t rc_offset = 4 * fAttrHiddenSize * fAttrHiddenSize + 3 * fAttrHiddenSize * fAttrHiddenSize;
            out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
             << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << " + "
             << rc_offset << ", &" << OpName << "_n, " << OpName << "_hidden_state + previous_offset, &"
             << OpName << "_n, &" << OpName << "_alpha, " << OpName << "_cell_gate + offset, &"
             << OpName << "_n);\n";
            if (fAttrInputForget == 0) {
               size_t rf_offset = 4 * fAttrHiddenSize * fAttrHiddenSize + 2 * fAttrHiddenSize * fAttrHiddenSize;
               out << SP << SP << SP << "BLAS::sgemm_(&" << OpName << "_transB, &" << OpName << "_transA, &"
                   << OpName << "_n, &m2, &" << OpName << "_n, &" << OpName << "_alpha, tensor_" << fNR << " + "
                   << rf_offset << ", &" << OpName << "_n, " << OpName << "_hidden_state + previous_offset, &"
                   << OpName << "_n, &" << OpName << "_alpha, " << OpName << "_forget_gate + offset, &"
                   << OpName << "_n);\n";
            }
         }
      }
      out << SP << SP << "}\n";

      // Clip the elements of the cell gate into the range [-fAttrClip, fAttrClip]
      if (fAttrClip > .0) {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float x = (" << OpName << "_cell_gate[i] > " << -fAttrClip << ") ? "
                << OpName << "_cell_gate[i] : " << -fAttrClip << ";\n";
         }
         out << SP << SP << SP << OpName << "_cell_gate[i] = (x < " << fAttrClip << ") ? x : "
             << fAttrClip << ";\n";
         out << SP << SP << "}\n";
      }
      // Apply the activation function to the cell gate, cell_gate = g(cell_gate)
      if (fAttrActivations[direction * 3 + 1] == "Relu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_cell_gate[i] < 0.)\n";
         out << SP << SP << SP << SP << OpName << "_cell_gate[i] = 0.;\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 1] == "Tanh") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float ex = exp(-2 * " << OpName << "_cell_gate[i]);\n";
         }
         out << SP << SP << SP << SP << OpName << "_cell_gate[i] = (1. - ex) / (1. + ex);\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 1] == "Sigmoid") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_cell_gate[i] = 1. / (1. + exp(-" << OpName
             << "_cell_gate[i]));\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 1] == "Affine") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_cell_gate[i] = "
             << fAttrActivationAlpha[direction * 3 + 1] << " * " << OpName << "_cell_gate[i] + "
             << fAttrActivationBeta[direction * 3 + 1] << ";\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 1] == "ScaledTanh") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float ex = exp(-2 * " << fAttrActivationBeta[direction * 3 + 1]
                << " * "<< OpName << "_cell_gate[i]);\n";
            }
            out << SP << SP << SP << SP << OpName << "_cell_gate[i] = "
                << fAttrActivationAlpha[direction * 3 + 1] << " * (1. - ex) / (1. + ex);\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 1] == "HardSigmoid") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float a = " << fAttrActivationAlpha[direction * 3 + 1] << " * "
                << OpName << "_cell_gate[i] + " << fAttrActivationBeta[direction * 3 + 1] << ";\n";
            out << SP << SP << SP << "float b = (a > 0.) ? a : 0.;\n";
         }
         out << SP << SP << SP << SP << OpName << "_cell_gate[i] = (b < 1.) ? b : 1.;\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 1] == "LeakyRelu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_cell_gate[i] < 0.)\n";
         out << SP << SP << SP << SP << OpName << "_cell_gate[i] = "
             << fAttrActivationAlpha[direction * 3 + 1] << " * " << OpName << "_cell_gate[i];\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 1] == "ThresholdRelu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_cell_gate[i] < "
             << fAttrActivationAlpha[direction * 3 + 1] << ")\n";
         out << SP << SP << SP << SP << OpName << "_cell_gate[i] = 0.;\n";
         out << SP << SP << "}";
      } else if (fAttrActivations[direction * 3 + 1] == "Elu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_cell_gate[i] < 0.)\n";
         out << SP << SP << SP << SP << OpName << "_cell_gate[i] = "
             << fAttrActivationAlpha[direction * 3 + 1] << " * exp(" << OpName << "_cell_gate[i] - 1.);\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 1] == "Softsign") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_cell_gate[i] = " << OpName
             << "_cell_gate[i] / (1. + abs(" << OpName << "_cell_gate[i]));\n";
         out << SP << SP << "}\n";
      } else { // fAttrActivations[direction * 3 + 1] = Softplus
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_cell_gate[i] = log(1. + exp("
             << OpName << "_cell_gate[i]));\n";
         out << SP << SP << "}\n";
      }

      // Peephole connections for the input gate and the forget gate
      if (!fNP.empty()) {
         // gate = 1.0 * gate + previous_cell_state * P^T
         out << SP << SP << "if (seq == 0) {\n";
         if (!fNInitial_c.empty()) {
            if (direction == 0) {
               out << SP << SP << SP << "for (size_t i = 0; i < " << size << "; i++) {\n";
               out << SP << SP << SP << SP << OpName << "_input_gate[i + offset] += tensor_" << fNP
                   << "[i] * " << OpName << "_initial_cell_state[i];\n";
               out << SP << SP << SP << "}\n";
               if (fAttrInputForget == 0) {
                  size_t pf_offset = batch_size * fAttrHiddenSize;
                  out << SP << SP << SP << "for (size_t i = 0; i < " << size << "; i++) {\n";
                  out << SP << SP << SP << SP << OpName << "_forget_gate[i + offset] += tensor_" << fNP
                      << "[i + " << pf_offset << "] * " << OpName << "_initial_cell_state[i];\n";
                  out << SP << SP << SP << "}\n";
               }
            } else {
               size_t pi_offset = 3 * batch_size * fAttrHiddenSize;
               size_t initial_c_offset = batch_size * fAttrHiddenSize;
               out << SP << SP << SP << "for (size_t i = 0; i < " << size << "; i++) {\n";
               out << SP << SP << SP << SP << OpName << "_input_gate[i + offset] += tensor_" << fNP
                   << "[i + " << pi_offset << "] * " << OpName << "_initial_cell_state[i + " << initial_c_offset
                   << "];\n";
               out << SP << SP << SP << "}\n";
               if (fAttrInputForget == 0) {
                  size_t pf_offset = 3 * batch_size * fAttrHiddenSize + batch_size * fAttrHiddenSize;
                  out << SP << SP << SP << "for (size_t i = 0; i < " << size << "; i++) {\n";
                  out << SP << SP << SP << SP << OpName << "_forget_gate[i + offset] += tensor_" << fNP
                      << "[i + " << pf_offset << "] * " << OpName << "_initial_cell_state[i + " << initial_c_offset
                      << "];\n";
                  out << SP << SP << SP << "}\n";
               }
            }
         }
         out << SP << SP << "} else {\n";
         if (direction == 0) {
            if (fAttrDirection == "backward") {
               out << SP << SP << SP << "size_t c_offset = (index + 1) * "
                   << num_directions * batch_size * fAttrHiddenSize << ";\n";
            } else {
               out << SP << SP << SP << "size_t c_offset = (seq - 1) * "
                   << num_directions * batch_size * fAttrHiddenSize << ";\n";
            }
            out << SP << SP << SP << "for (size_t i = 0; i < " << size << "; i++) {\n";
            out << SP << SP << SP << SP << OpName << "_input_gate[i + offset] += tensor_" << fNP
                << "[i] * " << OpName << "_cell_state[i + c_offset];\n";
            out << SP << SP << SP << "}\n";
            if (fAttrInputForget == 0) {
               size_t pf_offset = batch_size * fAttrHiddenSize;
               out << SP << SP << SP << "for (size_t i = 0; i < " << size << "; i++) {\n";
               out << SP << SP << SP << SP << OpName << "_forget_gate[i + offset] += tensor_" << fNP
                   << "[i + " << pf_offset << "] * " << OpName << "_cell_state[i + c_offset];\n";
               out << SP << SP << SP << "}\n";
            }
         } else { // direction=1
            size_t pi_offset = 3 * batch_size * fAttrHiddenSize;
            out << SP << SP << SP << "size_t c_offset = (index + 1) * "
                << num_directions * batch_size * fAttrHiddenSize << " + " << batch_size * fAttrHiddenSize << ";\n";
            out << SP << SP << SP << "for (size_t i = 0; i < " << size << "; i++) {\n";
            out << SP << SP << SP << SP << OpName << "_input_gate[i + offset] += tensor_" << fNP
                << "[i + " << pi_offset << "] * " << OpName << "_cell_state[i + c_offset];\n";
            out << SP << SP << SP << "}\n";
            if (fAttrInputForget == 0) {
               size_t pf_offset = 3 * batch_size * fAttrHiddenSize + batch_size * fAttrHiddenSize;
               out << SP << SP << SP << "for (size_t i = 0; i < " << size << "; i++) {\n";
               out << SP << SP << SP << SP << OpName << "_forget_gate[i + offset] += tensor_" << fNP
                   << "[i + " << pf_offset << "] * " << OpName << "_cell_state[i + c_offset];\n";
               out << SP << SP << SP << "}\n";
            }
         }
         out << SP << SP << "}\n";
      }

      // Clip the elements of the input gate into the range [-fAttrClip, fAttrClip]
      if (fAttrClip > .0) {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float x = (" << OpName << "_input_gate[i] > " << -fAttrClip << ") ? "
                << OpName << "_input_gate[i] : " << -fAttrClip << ";\n";
         }
         out << SP << SP << SP << OpName << "_input_gate[i] = (x < " << fAttrClip << ") ? x : "
             << fAttrClip << ";\n";
         out << SP << SP << "}\n";
      }
      // Apply the activation function to the input gate
      if (fAttrActivations[direction * 3] == "Relu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_input_gate[i] < 0.)\n";
         out << SP << SP << SP << SP << OpName << "_input_gate[i] = 0.;\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "Tanh") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float ex = exp(-2 * " << OpName << "_input_gate[i]);\n";
         }
         out << SP << SP << SP << SP << OpName << "_input_gate[i] = (1. - ex) / (1. + ex);\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "Sigmoid") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_input_gate[i] = 1. / (1. + exp(-" << OpName
             << "_input_gate[i]));\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "Affine") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_input_gate[i] = "
             << fAttrActivationAlpha[direction * 3] << " * " << OpName << "_input_gate[i] + "
             << fAttrActivationBeta[direction * 3] << ";\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "ScaledTanh") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float ex = exp(-2 * " << fAttrActivationBeta[direction * 3]
                << " * "<< OpName << "_input_gate[i]);\n";
            }
            out << SP << SP << SP << SP << OpName << "_input_gate[i] = "
                << fAttrActivationAlpha[direction * 3] << " * (1. - ex) / (1. + ex);\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "HardSigmoid") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float a = " << fAttrActivationAlpha[direction * 3] << " * "
                << OpName << "_input_gate[i] + " << fAttrActivationBeta[direction * 3] << ";\n";
            out << SP << SP << SP << "float b = (a > 0.) ? a : 0.;\n";
         }
         out << SP << SP << SP << SP << OpName << "_input_gate[i] = (b < 1.) ? b : 1.;\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "LeakyRelu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_input_gate[i] < 0.)\n";
         out << SP << SP << SP << SP << OpName << "_input_gate[i] = "
             << fAttrActivationAlpha[direction * 3] << " * " << OpName << "_input_gate[i];\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "ThresholdRelu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_input_gate[i] < "
             << fAttrActivationAlpha[direction * 3] << ")\n";
         out << SP << SP << SP << SP << OpName << "_input_gate[i] = 0.;\n";
         out << SP << SP << "}";
      } else if (fAttrActivations[direction * 3] == "Elu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_input_gate[i] < 0.)\n";
         out << SP << SP << SP << SP << OpName << "_input_gate[i] = "
             << fAttrActivationAlpha[direction * 3] << " * exp(" << OpName << "_input_gate[i] - 1.);\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "Softsign") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_input_gate[i] = " << OpName
             << "_input_gate[i] / (1. + abs(" << OpName << "_input_gate[i]));\n";
         out << SP << SP << "}\n";
      } else { // fAttrActivations[direction * 3] = Softplus
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_input_gate[i] = log(1. + exp("
             << OpName << "_input_gate[i]));\n";
         out << SP << SP << "}\n";
      }

      if (fAttrInputForget == 0) {
         // Clip the elements of the forget gate into the range [-fAttrClip, fAttrClip]
         if (fAttrClip > .0) {
            out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
            if (fType == "float") {
               out << SP << SP << SP << "float x = (" << OpName << "_forget_gate[i] > "
                   << -fAttrClip << ") ? " << OpName << "_forget_gate[i] : " << -fAttrClip << ";\n";
            }
            out << SP << SP << SP << OpName << "_forget_gate[i] = (x < " << fAttrClip
                << ") ? x : " << fAttrClip << ";\n";
            out << SP << SP << "}\n";
         }
         // Apply the activation function to the forget gate, cell_gate = g(cell_gate)
         if (fAttrActivations[direction * 3] == "Relu") {
            out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
            out << SP << SP << SP << "if (" << OpName << "_forget_gate[i] < 0.)\n";
            out << SP << SP << SP << SP << OpName << "_forget_gate[i] = 0.;\n";
            out << SP << SP << "}\n";
         } else if (fAttrActivations[direction * 3] == "Tanh") {
            out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
            if (fType == "float") {
               out << SP << SP << SP << "float ex = exp(-2 * " << OpName << "_forget_gate[i]);\n";
            }
            out << SP << SP << SP << SP << OpName << "_forget_gate[i] = (1. - ex) / (1. + ex);\n";
            out << SP << SP << "}\n";
         } else if (fAttrActivations[direction * 3] == "Sigmoid") {
            out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
            out << SP << SP << SP << SP << OpName << "_forget_gate[i] = 1. / (1. + exp(-"
                << OpName << "_forget_gate[i]));\n";
            out << SP << SP << "}\n";
         } else if (fAttrActivations[direction * 3] == "Affine") {
            out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
            out << SP << SP << SP << SP << OpName << "_forget_gate[i] = "
                << fAttrActivationAlpha[direction * 3] << " * " << OpName << "_forget_gate[i] + "
               << fAttrActivationBeta[direction * 3] << ";\n";
            out << SP << SP << "}\n";
         } else if (fAttrActivations[direction * 3] == "ScaledTanh") {
            out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
            if (fType == "float") {
               out << SP << SP << SP << "float ex = exp(-2 * " << fAttrActivationBeta[direction * 3]
                   << " * "<< OpName << "_forget_gate[i]);\n";
            }
            out << SP << SP << SP << SP << OpName << "_forget_gate[i] = "
                << fAttrActivationAlpha[direction * 3] << " * (1. - ex) / (1. + ex);\n";
            out << SP << SP << "}\n";
         } else if (fAttrActivations[direction * 3] == "HardSigmoid") {
            out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
            if (fType == "float") {
               out << SP << SP << SP << "float a = " << fAttrActivationAlpha[direction * 3] << " * "
                   << OpName << "_forget_gate[i] + " << fAttrActivationBeta[direction * 3] << ";\n";
               out << SP << SP << SP << "float b = (a > 0.) ? a : 0.;\n";
            }
            out << SP << SP << SP << SP << OpName << "_forget_gate[i] = (b < 1.) ? b : 1.;\n";
            out << SP << SP << "}\n";
         } else if (fAttrActivations[direction * 3] == "LeakyRelu") {
            out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
            out << SP << SP << SP << "if (" << OpName << "_forget_gate[i] < 0.)\n";
            out << SP << SP << SP << SP << OpName << "_forget_gate[i] = "
                << fAttrActivationAlpha[direction * 3] << " * " << OpName << "_forget_gate[i];\n";
            out << SP << SP << "}\n";
         } else if (fAttrActivations[direction * 3] == "ThresholdRelu") {
            out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
            out << SP << SP << SP << "if (" << OpName << "_forget_gate[i] < "
                << fAttrActivationAlpha[direction * 3] << ")\n";
            out << SP << SP << SP << SP << OpName << "_forget_gate[i] = 0.;\n";
            out << SP << SP << "}";
         } else if (fAttrActivations[direction * 3] == "Elu") {
            out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
            out << SP << SP << SP << "if (" << OpName << "_forget_gate[i] < 0.)\n";
            out << SP << SP << SP << SP << OpName << "_forget_gate[i] = "
                << fAttrActivationAlpha[direction * 3] << " * exp(" << OpName << "_forget_gate[i] - 1.);\n";
            out << SP << SP << "}\n";
         } else if (fAttrActivations[direction * 3] == "Softsign") {
            out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
            out << SP << SP << SP << SP << OpName << "_forget_gate[i] = " << OpName
                << "_forget_gate[i] / (1. + abs(" << OpName << "_forget_gate[i]));\n";
            out << SP << SP << "}\n";
         } else { // fAttrActivations[direction * 3] = Softplus
            out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
            out << SP << SP << SP << SP << OpName << "_forget_gate[i] = log(1. + exp("
                << OpName << "_forget_gate[i]));\n";
            out << SP << SP << "}\n";
         }
      }

      // cell_state = input_gate o cell_gate
      out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
      out << SP << SP << SP << OpName << "_cell_state[i] = " << OpName << "_input_gate[i] * "
          << OpName << "_cell_gate[i];\n";
      out << SP << SP << "}\n";

      if (fAttrInputForget == 0) {
         out << SP << SP << "if (seq == 0) {\n";
         if (!fNInitial_c.empty()) {
            // cell_state += forget_gate o initial_cell_state
            out << SP << SP << SP << "for (size_t i = 0; i < " << size << "; i++) {\n";
            out << SP << SP << SP << SP << OpName << "_cell_state[i + offset] += "
                << OpName << "_forget_gate[i + offset] * " << OpName << "_initial_cell_state[i];\n";
            out << SP << SP << SP << "}\n";
         }
         out << SP << SP << "} else {\n";
         // cell_state += forget_gate o previous_cell_state
         if (direction == 0) {
            if (fAttrDirection == "backward") {
               out << SP << SP << SP << "size_t previous_offset = (index + 1) * "
                   << num_directions * batch_size * fAttrHiddenSize << ";\n";
            } else {
               out << SP << SP << SP << "size_t previous_offset = (seq - 1) * "
                   << num_directions * batch_size * fAttrHiddenSize << ";\n";
            }
         } else { // direction=1
            out << SP << SP << SP << "size_t previous_offset = (index + 1) * "
                << num_directions * batch_size * fAttrHiddenSize << " + " << batch_size * fAttrHiddenSize << ";\n";
         }
         out << SP << SP << SP << "for (size_t i = 0; i < " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_cell_state[i + offset] += "
             << OpName << "_forget_gate[i + offset] * " << OpName << "_cell_state[i + previous_offset];\n";
         out << SP << SP << SP << "}\n";
         out << SP << SP << "}\n";
      }

      if (!fNP.empty()) {
         // Peephole connection for the output gate
         if (direction == 0) {
            size_t p_offset = 2 * batch_size * fAttrHiddenSize;
            out << SP << SP << SP << "for (size_t i = 0; i < " << size << "; i++) {\n";
            out << SP << SP << SP << SP << OpName << "_output_gate[i + offset] += tensor_"
                << fNP << "[i + " << p_offset << "] * " << OpName << "_cell_state[i + offset];\n";
            out << SP << SP << SP << "}\n";
         } else { // direction=1
            size_t p_offset = 3 * batch_size * fAttrHiddenSize + 2 * batch_size * fAttrHiddenSize;
            out << SP << SP << SP << "for (size_t i = 0; i < " << size << "; i++) {\n";
            out << SP << SP << SP << SP << OpName << "_output_gate[i + offset] += tensor_"
                << fNP << "[i + " << p_offset << "] * " << OpName << "_cell_state[i + offset];\n";
            out << SP << SP << SP << "}\n";
         }
      }

      // Clip the elements of the output gate into the range [-fAttrClip, fAttrClip]
      if (fAttrClip > .0) {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float x = (" << OpName << "_output_gate[i] > " << -fAttrClip
                << ") ? " << OpName << "_output_gate[i] : " << -fAttrClip << ";\n";
         }
         out << SP << SP << SP << OpName << "_output_gate[i] = (x < " << fAttrClip << ") ? x : "
             << fAttrClip << ";\n";
         out << SP << SP << "}\n";
      }
      // Apply the activation function to the output gate
      if (fAttrActivations[direction * 3] == "Relu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_output_gate[i] < 0.)\n";
         out << SP << SP << SP << SP << OpName << "_output_gate[i] = 0.;\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "Tanh") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float ex = exp(-2 * " << OpName << "_output_gate[i]);\n";
         }
         out << SP << SP << SP << SP << OpName << "_output_gate[i] = (1. - ex) / (1. + ex);\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "Sigmoid") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_output_gate[i] = 1. / (1. + exp(-" << OpName
             << "_output_gate[i]));\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "Affine") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_output_gate[i] = "
             << fAttrActivationAlpha[direction * 3] << " * " << OpName << "_output_gate[i] + "
             << fAttrActivationBeta[direction * 3] << ";\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "ScaledTanh") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float ex = exp(-2 * " << fAttrActivationBeta[direction * 3]
                << " * "<< OpName << "_output_gate[i]);\n";
            }
            out << SP << SP << SP << SP << OpName << "_output_gate[i] = "
                << fAttrActivationAlpha[direction * 3] << " * (1. - ex) / (1. + ex);\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "HardSigmoid") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float a = " << fAttrActivationAlpha[direction * 3] << " * "
                << OpName << "_output_gate[i] + " << fAttrActivationBeta[direction * 3] << ";\n";
            out << SP << SP << SP << "float b = (a > 0.) ? a : 0.;\n";
         }
         out << SP << SP << SP << SP << OpName << "_output_gate[i] = (b < 1.) ? b : 1.;\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "LeakyRelu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_output_gate[i] < 0.)\n";
         out << SP << SP << SP << SP << OpName << "_output_gate[i] = "
             << fAttrActivationAlpha[direction * 3] << " * " << OpName << "_output_gate[i];\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "ThresholdRelu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_output_gate[i] < "
             << fAttrActivationAlpha[direction * 3] << ")\n";
         out << SP << SP << SP << SP << OpName << "_output_gate[i] = 0.;\n";
         out << SP << SP << "}";
      } else if (fAttrActivations[direction * 3] == "Elu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_output_gate[i] < 0.)\n";
         out << SP << SP << SP << SP << OpName << "_output_gate[i] = "
             << fAttrActivationAlpha[direction * 3] << " * exp(" << OpName << "_output_gate[i] - 1.);\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3] == "Softsign") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_output_gate[i] = " << OpName
             << "_output_gate[i] / (1. + abs(" << OpName << "_output_gate[i]));\n";
         out << SP << SP << "}\n";
      } else { // fAttrActivations[direction * 3] = Softplus
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_output_gate[i] = log(1. + exp("
             << OpName << "_output_gate[i]));\n";
         out << SP << SP << "}\n";
      }

      // copy cell_state into new_cell_state
      out << SP << SP << "std::copy(" << OpName << "_cell_state + offset, " << OpName
          << "_cell_state + offset + " << size << ", "<< OpName << "_new_cell_state + offset);\n";
      // Clip the elements of the new_cell_state into the range [-fAttrClip, fAttrClip]
      if (fAttrClip > .0) {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float x = (" << OpName << "_new_cell_state[i] > " << -fAttrClip
                << ") ? " << OpName << "_new_cell_state[i] : " << -fAttrClip << ";\n";
         }
         out << SP << SP << SP << OpName << "_new_cell_state[i] = (x < " << fAttrClip << ") ? x : "
             << fAttrClip << ";\n";
         out << SP << SP << "}\n";
      }
      // Apply the activation function to the new cell state
      if (fAttrActivations[direction * 3 + 2] == "Relu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_new_cell_state[i] < 0.)\n";
         out << SP << SP << SP << SP << OpName << "_new_cell_state[i] = 0.;\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 2] == "Tanh") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float ex = exp(-2 * " << OpName << "_new_cell_state[i]);\n";
         }
         out << SP << SP << SP << SP << OpName << "_new_cell_state[i] = (1. - ex) / (1. + ex);\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 2] == "Sigmoid") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_new_cell_state[i] = 1. / (1. + exp(-" << OpName
             << "_new_cell_state[i]));\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 2] == "Affine") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_new_cell_state[i] = "
             << fAttrActivationAlpha[direction * 3 + 2] << " * " << OpName << "_new_cell_state[i] + "
             << fAttrActivationBeta[direction * 3 + 2] << ";\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 2] == "ScaledTanh") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float ex = exp(-2 * " << fAttrActivationBeta[direction * 3 + 2]
                << " * "<< OpName << "_new_cell_state[i]);\n";
            }
            out << SP << SP << SP << SP << OpName << "_new_cell_state[i] = "
                << fAttrActivationAlpha[direction * 3 + 2] << " * (1. - ex) / (1. + ex);\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 2] == "HardSigmoid") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         if (fType == "float") {
            out << SP << SP << SP << "float a = " << fAttrActivationAlpha[direction * 3 + 2] << " * "
                << OpName << "_new_cell_state[i] + " << fAttrActivationBeta[direction * 3 + 2] << ";\n";
            out << SP << SP << SP << "float b = (a > 0.) ? a : 0.;\n";
         }
         out << SP << SP << SP << SP << OpName << "_new_cell_state[i] = (b < 1.) ? b : 1.;\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 2] == "LeakyRelu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_new_cell_state[i] < 0.)\n";
         out << SP << SP << SP << SP << OpName << "_new_cell_state[i] = "
             << fAttrActivationAlpha[direction * 3 + 2] << " * " << OpName << "_new_cell_state[i];\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 2] == "ThresholdRelu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_new_cell_state[i] < "
             << fAttrActivationAlpha[direction * 3 + 2] << ")\n";
         out << SP << SP << SP << SP << OpName << "_new_cell_state[i] = 0.;\n";
         out << SP << SP << "}";
      } else if (fAttrActivations[direction * 3 + 2] == "Elu") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << "if (" << OpName << "_new_cell_state[i] < 0.)\n";
         out << SP << SP << SP << SP << OpName << "_new_cell_state[i] = "
             << fAttrActivationAlpha[direction * 3 + 2] << " * exp(" << OpName << "_new_cell_state[i] - 1.);\n";
         out << SP << SP << "}\n";
      } else if (fAttrActivations[direction * 3 + 2] == "Softsign") {
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_new_cell_state[i] = " << OpName
             << "_new_cell_state[i] / (1. + abs(" << OpName << "_new_cell_state[i]));\n";
         out << SP << SP << "}\n";
      } else { // fAttrActivations[direction * 3 + 2] = Softplus
         out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
         out << SP << SP << SP << SP << OpName << "_new_cell_state[i] = log(1. + exp("
             << OpName << "_new_cell_state[i]));\n";
         out << SP << SP << "}\n";
      }

      // hidden_state = output_gate o new_cell_state
      out << SP << SP << "for (size_t i = offset; i < offset + " << size << "; i++) {\n";
      out << SP << SP << SP << OpName << "_hidden_state[i] = " << OpName << "_output_gate[i] * "
          << OpName << "_new_cell_state[i];\n";
      out << SP << SP << "}\n";
      out << SP << "}\n";
   }

   // Padding the hidden state for LSTM with different sequence lengths
   if (!fNSequence_lens.empty()) {
      out << SP << "for (size_t seq = 0; seq < " << seq_length << "; seq++) {\n";
      out << SP << SP << "for (size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
      out << SP << SP << SP << "if (seq >= tensor_" << fNSequence_lens << "[batch]) {\n";
      for (size_t direction = 0; direction < num_directions; direction++) {
         out << SP << SP << SP << SP << SP << "for (size_t h = 0; h < " << fAttrHiddenSize << "; h++) {\n";
         out << SP << SP << SP << SP << SP << SP << "size_t idx = seq * "
             << num_directions * batch_size * fAttrHiddenSize + direction * batch_size * fAttrHiddenSize
             << " + batch * " << fAttrHiddenSize << " + h;\n";
         out << SP << SP << SP << SP << SP << SP << OpName << "_cell_state[idx] = 0.;\n";
         out << SP << SP << SP << SP << SP << SP << OpName << "_hidden_state[idx] = 0.;\n";
         out << SP << SP << SP << SP << SP << "}\n";
      }
      out << SP << SP << SP << "}\n";
      out << SP << SP << "}\n";
      out << SP << "}\n";
   }

   // Copy the hidden state into y and y_h and copy cell_state into y_c
   if (fAttrLayout == 0) {
      if (!fNY_h.empty()) {
         // Copy hidden_state into Y_h
         if (fNSequence_lens.empty()) {
            size_t y_h_size = batch_size * fAttrHiddenSize;
            if (fAttrDirection == "backward") {
               out << SP << "std::copy(" << OpName << "_hidden_state, " << OpName << "_hidden_state + "
                   << y_h_size << ", tensor_" << fNY_h << ");\n";
            } else {
               size_t offset = (seq_length - 1) * num_directions * batch_size * fAttrHiddenSize;
               out << SP << "std::copy(" << OpName << "_hidden_state + " << offset << ", " << OpName
                   << "_hidden_state + " << offset << " + " << y_h_size << ", tensor_" << fNY_h << ");\n";
            }
            if (num_directions == 2) {
               out << SP << "std::copy(" << OpName << "_hidden_state + " << y_h_size << ", " << OpName
                   << "_hidden_state + " << 2 * y_h_size << ", tensor_" << fNY_h << " + " << y_h_size << ");\n";
            }
         } else { // LSTM with different sequence lengths
            if (fAttrDirection == "backward") {
               out << SP << "for (size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
               out << SP << SP << "size_t offset = batch * " << fAttrHiddenSize << ";\n";
               out << SP << SP << "std::copy(" << OpName << "_hidden_state + offset, " << OpName
                   << "_hidden_state + offset + " << fAttrHiddenSize << ", tensor_" << fNY_h << " + offset);\n";
               out << SP << "}\n";
            } else {
               out << SP << "for (size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
               out << SP << SP << "size_t seq = " << "tensor_" << fNSequence_lens << "[batch] - 1;\n";
               out << SP << SP << "size_t offset = seq * " << num_directions * batch_size * fAttrHiddenSize
                   << " + batch * " << fAttrHiddenSize << ";\n";
               out << SP << SP << "size_t y_h_offset = batch * " << fAttrHiddenSize << ";\n";
               out << SP << SP << "std::copy(" << OpName << "_hidden_state + offset, " << OpName
                   << "_hidden_state + offset + " << fAttrHiddenSize << ", tensor_" << fNY_h << " + y_h_offset);\n";
               out << SP << "}\n";
            }
            if (num_directions == 2) {
               out << SP << "for (size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
               out << SP << SP << "size_t offset = " << batch_size * fAttrHiddenSize
                   << " + batch * " << fAttrHiddenSize << ";\n";
               out << SP << SP << "size_t y_h_offset = " << batch_size * fAttrHiddenSize
                   << " + batch * " << fAttrHiddenSize << ";\n";
               out << SP << SP << "std::copy(" << OpName << "_hidden_state + offset, " << OpName
                   << "_hidden_state + offset + " << fAttrHiddenSize << ", tensor_" << fNY_h << " + y_h_offset);\n";
               out << SP << "}\n";
            }
         }
      }
      if (!fNY_c.empty()) {
         // Copy cell_state into Y_c
         if (fNSequence_lens.empty()) {
            size_t y_h_size = batch_size * fAttrHiddenSize;
            if (fAttrDirection == "backward") {
               out << SP << "std::copy(" << OpName << "_cell_state, " << OpName << "_hidden_state + "
                   << y_h_size << ", tensor_" << fNY_c << ");\n";
            } else {
               size_t offset = (seq_length - 1) * num_directions * batch_size * fAttrHiddenSize;
               out << SP << "std::copy(" << OpName << "_cell_state + " << offset << ", " << OpName
                   << "_cell_state + " << offset << " + " << y_h_size << ", tensor_" << fNY_c << ");\n";
            }
            if (num_directions == 2) {
               out << SP << "std::copy(" << OpName << "_cell_state + " << y_h_size << ", " << OpName
                   << "_cell_state + " << 2 * y_h_size << ", tensor_" << fNY_c << " + " << y_h_size << ");\n";
            }
         } else { // LSTM with different sequence lengths
            if (fAttrDirection == "backward") {
               out << SP << "for (size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
               out << SP << SP << "size_t offset = batch * " << fAttrHiddenSize << ";\n";
               out << SP << SP << "std::copy(" << OpName << "_cell_state + offset, " << OpName
                   << "_cell_state + offset + " << fAttrHiddenSize << ", tensor_" << fNY_c << " + offset);\n";
               out << SP << "}\n";
            } else {
               out << SP << "for (size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
               out << SP << SP << "size_t seq = " << "tensor_" << fNSequence_lens << "[batch] - 1;\n";
               out << SP << SP << "size_t offset = seq * " << num_directions * batch_size * fAttrHiddenSize
                   << " + batch * " << fAttrHiddenSize << ";\n";
               out << SP << SP << "size_t y_h_offset = batch * " << fAttrHiddenSize << ";\n";
               out << SP << SP << "std::copy(" << OpName << "_cell_state + offset, " << OpName
                   << "_cell_state + offset + " << fAttrHiddenSize << ", tensor_" << fNY_c << " + y_h_offset);\n";
               out << SP << "}\n";
            }
            if (num_directions == 2) {
               out << SP << "for (size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
               out << SP << SP << "size_t offset = " << batch_size * fAttrHiddenSize
                   << " + batch * " << fAttrHiddenSize << ";\n";
               out << SP << SP << "size_t y_h_offset = " << batch_size * fAttrHiddenSize
                   << " + batch * " << fAttrHiddenSize << ";\n";
               out << SP << SP << "std::copy(" << OpName << "_cell_state + offset, " << OpName
                   << "_cell_state + offset + " << fAttrHiddenSize << ", tensor_" << fNY_c << " + y_h_offset);\n";
               out << SP << "}\n";
            }
         }
      }
   } else { // fAttrLayout=1
      if (!fNY.empty()) {
         // Copy hidden_state into Y
         for (size_t direction = 0; direction < num_directions; direction++) {
            out << SP << "for (size_t seq = 0; seq < " << seq_length << "; seq++) {\n";
            out << SP << SP << "for (size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
            out << SP << SP << SP << "size_t offset = seq * " << num_directions * batch_size * fAttrHiddenSize
                << " + " << direction * batch_size * fAttrHiddenSize << " + batch * " << fAttrHiddenSize << ";\n";
            out << SP << SP << SP << "size_t y_offset = batch * " << seq_length * num_directions * fAttrHiddenSize
                << " + seq * " << num_directions * fAttrHiddenSize << " + " << direction * fAttrHiddenSize << ";\n";
            out << SP << SP << SP << "std::copy(" << OpName << "_hidden_state + offset, " << OpName
                << "_hidden_state + offset + " << fAttrHiddenSize << ", tensor_" << fNY << " + y_offset);\n";
            out << SP << SP << "}\n";
            out << SP << "}\n";
         }
      }
      if (!fNY_h.empty()) {
         // Copy the hidden_state into Y_h
         if (fAttrDirection == "backward") {
            out << SP << "for (size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
            out << SP << SP << "size_t offset = batch * " << fAttrHiddenSize << ";\n";
            out << SP << SP << "size_t y_h_offset = batch * " << num_directions * fAttrHiddenSize << ";\n";
            out << SP << SP << "std::copy(" << OpName << "_hidden_state + offset, " << OpName
                << "_hidden_state + offset + " << fAttrHiddenSize << ", tensor_" << fNY_h << " + y_h_offset);\n";
            out << SP << "}\n";
         } else {
            out << SP << "for (size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
            if (fNSequence_lens.empty()) {
               out << SP << SP << "size_t seq = " << seq_length - 1 << ";\n";
            } else {
               out << SP << SP << "size_t seq = " << "tensor_" << fNSequence_lens << "[batch] - 1;\n";
            }
            out << SP << SP << "size_t offset = seq * " << num_directions * batch_size * fAttrHiddenSize
                << " + batch * " << fAttrHiddenSize << ";\n";
            out << SP << SP << "size_t y_h_offset = batch * " << num_directions * fAttrHiddenSize << ";\n";
            out << SP << SP << "std::copy(" << OpName << "_hidden_state + offset, " << OpName
                << "_hidden_state + offset + " << fAttrHiddenSize << ", tensor_" << fNY_h << " + y_h_offset);\n";
            out << SP << "}\n";
         }
         if (num_directions == 2) {
            out << SP << "for (size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
            out << SP << SP << "size_t offset = " << batch_size * fAttrHiddenSize << " + batch * "
                << fAttrHiddenSize << ";\n";
            out << SP << SP << "size_t y_h_offset = batch * " << num_directions * fAttrHiddenSize << " + "
                << fAttrHiddenSize << ";\n";
            out << SP << SP << "std::copy(" << OpName << "_hidden_state + offset, " << OpName
                << "_hidden_state + offset + " << fAttrHiddenSize << ", tensor_" << fNY_h << " + y_h_offset);\n";
            out << SP << "}\n";
         }
      }

      if (!fNY_c.empty()) {
         // copy the cell_state into Y_c
         if (fAttrDirection == "backward") {
            out << SP << "for (size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
            out << SP << SP << "size_t offset = batch * " << fAttrHiddenSize << ";\n";
            out << SP << SP << "size_t y_h_offset = batch * " << num_directions * fAttrHiddenSize << ";\n";
            out << SP << SP << "std::copy(" << OpName << "_cell_state + offset, " << OpName
                << "_cell_state + offset + " << fAttrHiddenSize << ", tensor_" << fNY_c << " + y_h_offset);\n";
            out << SP << "}\n";
         } else {
            out << SP << "for (size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
            if (fNSequence_lens.empty()) {
               out << SP << SP << "size_t seq = " << seq_length - 1 << ";\n";
            } else {
               out << SP << SP << "size_t seq = " << "tensor_" << fNSequence_lens << "[batch] - 1;\n";
            }
            out << SP << SP << "size_t offset = seq * " << num_directions * batch_size * fAttrHiddenSize
                << " + batch * " << fAttrHiddenSize << ";\n";
            out << SP << SP << "size_t y_h_offset = batch * " << num_directions * fAttrHiddenSize << ";\n";
            out << SP << SP << "std::copy(" << OpName << "_cell_state + offset, " << OpName
                << "_cell_state + offset + " << fAttrHiddenSize << ", tensor_" << fNY_c << " + y_h_offset);\n";
            out << SP << "}\n";
         }
         if (num_directions == 2) {
            out << SP << "for (size_t batch = 0; batch < " << batch_size << "; batch++) {\n";
            out << SP << SP << "size_t offset = " << batch_size * fAttrHiddenSize << " + batch * "
                << fAttrHiddenSize << ";\n";
            out << SP << SP << "size_t y_h_offset = batch * " << num_directions * fAttrHiddenSize << " + "
                << fAttrHiddenSize << ";\n";
            out << SP << SP << "std::copy(" << OpName << "_cell_state + offset, " << OpName
                << "_cell_state + offset + " << fAttrHiddenSize << ", tensor_" << fNY_c << " + y_h_offset);\n";
            out << SP << "}\n";
         }
      }
   }

   return out.str();
}

} // namespace SOFIE
} // namespace Experimental
} // namespace TMVA

#endif
