/*

 * Copyright (c) 2023, Alliance for Open Media. All rights reserved.

 *

 * This source code is subject to the terms of the BSD 2 Clause License and

 * the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License

 * was not distributed with this source code in the LICENSE file, you can

 * obtain it at www.aomedia.org/license/software. If the Alliance for Open

 * Media Patent License 1.0 was not distributed with this source code in the

 * PATENTS file, you can obtain it at www.aomedia.org/license/patent.

 */

#include <stdbool.h>

#include <assert.h>

#include <immintrin.h>

#include "config/av1_rtcd.h"

#include "av1/encoder/ml.h"

#include "av1/encoder/x86/ml_sse3.h"

#define CALC_OUTPUT_FOR_2ROWS                                               \

  const int index = weight_idx + (2 * i * tot_num_inputs);                  \

  const __m256 weight0 = _mm256_loadu_ps(&weights[index]);                  \

  const __m256 weight1 = _mm256_loadu_ps(&weights[index + tot_num_inputs]); \

  const __m256 mul0 = _mm256_mul_ps(inputs256, weight0);                    \

  const __m256 mul1 = _mm256_mul_ps(inputs256, weight1);                    \

  hadd[i] = _mm256_hadd_ps(mul0, mul1);

static inline void nn_propagate_8to1(

    const float *const inputs, const float *const weights,

    const float *const bias, int num_inputs_to_process, int tot_num_inputs,

    int num_outputs, float *const output_nodes, int is_clip_required) {

  // Process one output row at a time.

  for (int out = 0; out < num_outputs; out++) {

    __m256 in_result = _mm256_setzero_ps();

    float bias_val = bias[out];

    for (int in = 0; in < num_inputs_to_process; in += 8) {

      const __m256 inputs256 = _mm256_loadu_ps(&inputs[in]);

      const int weight_idx = in + (out * tot_num_inputs);

      const __m256 weight0 = _mm256_loadu_ps(&weights[weight_idx]);

      const __m256 mul0 = _mm256_mul_ps(inputs256, weight0);

      in_result = _mm256_add_ps(in_result, mul0);

    }

    const __m128 low_128 = _mm256_castps256_ps128(in_result);

    const __m128 high_128 = _mm256_extractf128_ps(in_result, 1);

    const __m128 sum_par_0 = _mm_add_ps(low_128, high_128);

    const __m128 sum_par_1 = _mm_hadd_ps(sum_par_0, sum_par_0);

    const __m128 sum_tot =

        _mm_add_ps(_mm_shuffle_ps(sum_par_1, sum_par_1, 0x99), sum_par_1);

    bias_val += _mm_cvtss_f32(sum_tot);

    if (is_clip_required) bias_val = AOMMAX(bias_val, 0);

    output_nodes[out] = bias_val;

  }

}

static inline void nn_propagate_8to4(

    const float *const inputs, const float *const weights,

    const float *const bias, int num_inputs_to_process, int tot_num_inputs,

    int num_outputs, float *const output_nodes, int is_clip_required) {

  __m256 hadd[2];

  for (int out = 0; out < num_outputs; out += 4) {

    __m128 bias_reg = _mm_loadu_ps(&bias[out]);

    __m128 in_result = _mm_setzero_ps();

    for (int in = 0; in < num_inputs_to_process; in += 8) {

      const __m256 inputs256 = _mm256_loadu_ps(&inputs[in]);

      const int weight_idx = in + (out * tot_num_inputs);

      // Process two output row at a time.

      for (int i = 0; i < 2; i++) {

        CALC_OUTPUT_FOR_2ROWS

      }

      const __m256 sum_par = _mm256_hadd_ps(hadd[0], hadd[1]);

      const __m128 low_128 = _mm256_castps256_ps128(sum_par);

      const __m128 high_128 = _mm256_extractf128_ps(sum_par, 1);

      const __m128 result = _mm_add_ps(low_128, high_128);

      in_result = _mm_add_ps(in_result, result);

    }

    in_result = _mm_add_ps(in_result, bias_reg);

    if (is_clip_required) in_result = _mm_max_ps(in_result, _mm_setzero_ps());

    _mm_storeu_ps(&output_nodes[out], in_result);

  }

}

static inline void nn_propagate_8to8(

    const float *const inputs, const float *const weights,

    const float *const bias, int num_inputs_to_process, int tot_num_inputs,

    int num_outputs, float *const output_nodes, int is_clip_required) {

  __m256 hadd[4];

  for (int out = 0; out < num_outputs; out += 8) {

    __m256 bias_reg = _mm256_loadu_ps(&bias[out]);

    __m256 in_result = _mm256_setzero_ps();

    for (int in = 0; in < num_inputs_to_process; in += 8) {

      const __m256 inputs256 = _mm256_loadu_ps(&inputs[in]);

      const int weight_idx = in + (out * tot_num_inputs);

      // Process two output rows at a time.

      for (int i = 0; i < 4; i++) {

        CALC_OUTPUT_FOR_2ROWS

      }

      const __m256 hh0 = _mm256_hadd_ps(hadd[0], hadd[1]);

      const __m256 hh1 = _mm256_hadd_ps(hadd[2], hadd[3]);

      __m256 ht_0 = _mm256_permute2f128_ps(hh0, hh1, 0x20);

      __m256 ht_1 = _mm256_permute2f128_ps(hh0, hh1, 0x31);

      __m256 result = _mm256_add_ps(ht_0, ht_1);

      in_result = _mm256_add_ps(in_result, result);

    }

    in_result = _mm256_add_ps(in_result, bias_reg);

    if (is_clip_required)

      in_result = _mm256_max_ps(in_result, _mm256_setzero_ps());

    _mm256_storeu_ps(&output_nodes[out], in_result);

  }

}

static inline void nn_propagate_input_multiple_of_8(

    const float *const inputs, const float *const weights,

    const float *const bias, int num_inputs_to_process, int tot_num_inputs,

    bool is_output_layer, int num_outputs, float *const output_nodes) {

  // The saturation of output is considered for hidden layer which is not equal

  // to final hidden layer.

  const int is_clip_required =

      !is_output_layer && num_inputs_to_process == tot_num_inputs;

  if (num_outputs % 8 == 0) {

    nn_propagate_8to8(inputs, weights, bias, num_inputs_to_process,

                      tot_num_inputs, num_outputs, output_nodes,

                      is_clip_required);

  } else if (num_outputs % 4 == 0) {

    nn_propagate_8to4(inputs, weights, bias, num_inputs_to_process,

                      tot_num_inputs, num_outputs, output_nodes,

                      is_clip_required);

  } else {

    nn_propagate_8to1(inputs, weights, bias, num_inputs_to_process,

                      tot_num_inputs, num_outputs, output_nodes,

                      is_clip_required);

  }

}

void av1_nn_predict_avx2(const float *input_nodes,

                         const NN_CONFIG *const nn_config, int reduce_prec,

                         float *const output) {

  float buf[2][NN_MAX_NODES_PER_LAYER];

  int buf_index = 0;

  int num_inputs = nn_config->num_inputs;

  assert(num_inputs > 0 && num_inputs <= NN_MAX_NODES_PER_LAYER);

  for (int layer = 0; layer <= nn_config->num_hidden_layers; layer++) {

    const float *layer_weights = nn_config->weights[layer];

    const float *layer_bias = nn_config->bias[layer];

    bool is_output_layer = layer == nn_config->num_hidden_layers;

    float *const output_nodes = is_output_layer ? output : &buf[buf_index][0];

    const int num_outputs = is_output_layer

                                ? nn_config->num_outputs

                                : nn_config->num_hidden_nodes[layer];

    assert(num_outputs > 0 && num_outputs <= NN_MAX_NODES_PER_LAYER);

    // Process input multiple of 8 using AVX2 intrinsic.

    if (num_inputs % 8 == 0) {

      nn_propagate_input_multiple_of_8(input_nodes, layer_weights, layer_bias,

                                       num_inputs, num_inputs, is_output_layer,

                                       num_outputs, output_nodes);

    } else {

      // When number of inputs is not multiple of 8, use hybrid approach of AVX2

      // and SSE3 based on the need.

      const int in_mul_8 = num_inputs / 8;

      const int num_inputs_to_process = in_mul_8 * 8;

      int bias_is_considered = 0;

      if (in_mul_8) {

        nn_propagate_input_multiple_of_8(

            input_nodes, layer_weights, layer_bias, num_inputs_to_process,

            num_inputs, is_output_layer, num_outputs, output_nodes);

        bias_is_considered = 1;

      }

      const float *out_temp = bias_is_considered ? output_nodes : layer_bias;

      const int input_remaining = num_inputs % 8;

      if (input_remaining % 4 == 0 && num_outputs % 8 == 0) {

        for (int out = 0; out < num_outputs; out += 8) {

          __m128 out_h = _mm_loadu_ps(&out_temp[out + 4]);

          __m128 out_l = _mm_loadu_ps(&out_temp[out]);

          for (int in = in_mul_8 * 8; in < num_inputs; in += 4) {

            av1_nn_propagate_4to8_sse3(&input_nodes[in],

                                       &layer_weights[out * num_inputs + in],

                                       &out_h, &out_l, num_inputs);

          }

          if (!is_output_layer) {

            const __m128 zero = _mm_setzero_ps();

            out_h = _mm_max_ps(out_h, zero);

            out_l = _mm_max_ps(out_l, zero);

          }

          _mm_storeu_ps(&output_nodes[out + 4], out_h);

          _mm_storeu_ps(&output_nodes[out], out_l);

        }

      } else if (input_remaining % 4 == 0 && num_outputs % 4 == 0) {

        for (int out = 0; out < num_outputs; out += 4) {

          __m128 outputs = _mm_loadu_ps(&out_temp[out]);

          for (int in = in_mul_8 * 8; in < num_inputs; in += 4) {

            av1_nn_propagate_4to4_sse3(&input_nodes[in],

                                       &layer_weights[out * num_inputs + in],

                                       &outputs, num_inputs);

          }

          if (!is_output_layer) outputs = _mm_max_ps(outputs, _mm_setzero_ps());

          _mm_storeu_ps(&output_nodes[out], outputs);

        }

      } else if (input_remaining % 4 == 0) {

        for (int out = 0; out < num_outputs; out++) {

          __m128 outputs = _mm_load1_ps(&out_temp[out]);

          for (int in = in_mul_8 * 8; in < num_inputs; in += 4) {

            av1_nn_propagate_4to1_sse3(&input_nodes[in],

                                       &layer_weights[out * num_inputs + in],

                                       &outputs);

          }

          if (!is_output_layer) outputs = _mm_max_ps(outputs, _mm_setzero_ps());

          output_nodes[out] = _mm_cvtss_f32(outputs);

        }

      } else {

        // Use SSE instructions for scalar operations to avoid the latency

        // of swapping between SIMD and FPU modes.

        for (int out = 0; out < num_outputs; out++) {

          __m128 outputs = _mm_load1_ps(&out_temp[out]);

          for (int in_node = in_mul_8 * 8; in_node < num_inputs; in_node++) {

            __m128 input = _mm_load1_ps(&input_nodes[in_node]);

            __m128 weight =

                _mm_load1_ps(&layer_weights[num_inputs * out + in_node]);

            outputs = _mm_add_ps(outputs, _mm_mul_ps(input, weight));

          }

          if (!is_output_layer) outputs = _mm_max_ps(outputs, _mm_setzero_ps());

          output_nodes[out] = _mm_cvtss_f32(outputs);

        }

      }

    }

    // Before processing the next layer, treat the output of current layer as

    // input to next layer.

    input_nodes = output_nodes;

    num_inputs = num_outputs;

    buf_index = 1 - buf_index;

  }

  if (reduce_prec) av1_nn_output_prec_reduce(output, nn_config->num_outputs);

}