add lpnorm

CMakeLists.txt

@@ -203,6 +203,7 @@ set(CC_SOURCES
     ${SOURCE_DIR}/math/inverse_fast_fourier_transform.cc
     ${SOURCE_DIR}/math/inverse_fourier_transform.cc
     ${SOURCE_DIR}/math/levinson_durbin_recursion.cc
+    ${SOURCE_DIR}/math/lp_norm_calculation.cc
     ${SOURCE_DIR}/math/matrix.cc
     ${SOURCE_DIR}/math/matrix2d.cc
     ${SOURCE_DIR}/math/minmax_accumulation.cc
@@ -354,6 +355,7 @@ set(MAIN_SOURCES
     ${SOURCE_DIR}/main/lpc2lsp.cc
     ${SOURCE_DIR}/main/lpc2par.cc
     ${SOURCE_DIR}/main/lpccheck.cc
+    ${SOURCE_DIR}/main/lpnorm.cc
     ${SOURCE_DIR}/main/lsp2lpc.cc
     ${SOURCE_DIR}/main/lspcheck.cc
     ${SOURCE_DIR}/main/lspdf.cc

doc/main/lpnorm.rst

@@ -0,0 +1,13 @@
+.. _lpnorm:
+
+lpnorm
+======
+
+.. doxygenfile:: lpnorm.cc
+
+.. seealso::
+
+    :ref:`vopr` :ref:`vsum`
+
+.. doxygenclass:: sptk::LpNormCalculation
+    :members:

include/SPTK/math/lp_norm_calculation.h

@@ -0,0 +1,78 @@
+// ------------------------------------------------------------------------ //
+// Copyright 2021 SPTK Working Group                                        //
+//                                                                          //
+// Licensed under the Apache License, Version 2.0 (the "License");          //
+// you may not use this file except in compliance with the License.         //
+// You may obtain a copy of the License at                                  //
+//                                                                          //
+//     http://www.apache.org/licenses/LICENSE-2.0                           //
+//                                                                          //
+// Unless required by applicable law or agreed to in writing, software      //
+// distributed under the License is distributed on an "AS IS" BASIS,        //
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. //
+// See the License for the specific language governing permissions and      //
+// limitations under the License.                                           //
+// ------------------------------------------------------------------------ //
+
+#ifndef SPTK_MATH_LP_NORM_CALCULATION_H_
+#define SPTK_MATH_LP_NORM_CALCULATION_H_
+
+#include <vector>  // std::vector
+
+#include "SPTK/utils/sptk_utils.h"
+
+namespace sptk {
+
+/**
+ * Calculate Lp-norm.
+ *
+ * The input is the @f$M@f$-th order vector:
+ * @f[
+ *   \begin{array}{cccc}
+ *     x(0), & x(1), & \ldots, & x(M),
+ *   \end{array}
+ * @f]
+ * The output is the Lp-norm represented as
+ * @f[
+ *   \begin{array}{cccc}
+ *     \|x\|_p = \left(\displaystyle\sum_{m=0}^M |x(m)|^p\right)^{\frac{1}{p}}.
+ *   \end{array}
+ * @f]
+ */
+class LpNormCalculation {
+ public:
+  /**
+   * @param[in] num_order Order of input, @f$M@f$.
+   * @param[in] p Power, @f$p@f$.
+   */
+  LpNormCalculation(int num_order, double p);
+
+  virtual ~LpNormCalculation() {
+  }
+
+  /**
+   * @return True if this object is valid.
+   */
+  bool IsValid() const {
+    return is_valid_;
+  }
+
+  /**
+   * @param[in] data @f$M@f$-th order vector.
+   * @param[out] norm Lp-norm.
+   * @return True on success, false on failure.
+   */
+  bool Run(const std::vector<double>& data, double* norm) const;
+
+ private:
+  const int num_order_;
+  const double p_;
+
+  bool is_valid_;
+
+  DISALLOW_COPY_AND_ASSIGN(LpNormCalculation);
+};
+
+}  // namespace sptk
+
+#endif  // SPTK_MATH_LP_NORM_CALCULATION_H_

src/main/lpnorm.cc

@@ -0,0 +1,202 @@
+// ------------------------------------------------------------------------ //
+// Copyright 2021 SPTK Working Group                                        //
+//                                                                          //
+// Licensed under the Apache License, Version 2.0 (the "License");          //
+// you may not use this file except in compliance with the License.         //
+// You may obtain a copy of the License at                                  //
+//                                                                          //
+//     http://www.apache.org/licenses/LICENSE-2.0                           //
+//                                                                          //
+// Unless required by applicable law or agreed to in writing, software      //
+// distributed under the License is distributed on an "AS IS" BASIS,        //
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. //
+// See the License for the specific language governing permissions and      //
+// limitations under the License.                                           //
+// ------------------------------------------------------------------------ //
+
+#include <algorithm>  // std::transform
+#include <fstream>    // std::ifstream
+#include <iomanip>    // std::setw
+#include <iostream>   // std::cerr, std::cin, std::cout, std::endl, etc.
+#include <sstream>    // std::ostringstream
+#include <vector>     // std::vector
+
+#include "GETOPT/ya_getopt.h"
+#include "SPTK/math/lp_norm_calculation.h"
+#include "SPTK/utils/sptk_utils.h"
+
+namespace {
+
+const int kDefaultNumOrder(25);
+const double kDefaultNormOrder(2.0);
+
+void PrintUsage(std::ostream* stream) {
+  // clang-format off
+  *stream << std::endl;
+  *stream << " lpnorm - normalize by Lp-norm" << std::endl;
+  *stream << std::endl;
+  *stream << "  usage:" << std::endl;
+  *stream << "       lpnorm [ options ] [ infile ] > stdout" << std::endl;
+  *stream << "  options:" << std::endl;
+  *stream << "       -l l  : length of vector   (   int)[" << std::setw(5) << std::right << kDefaultNumOrder + 1 << "][   1 <= l <=   ]" << std::endl;  // NOLINT
+  *stream << "       -m m  : order of vector    (   int)[" << std::setw(5) << std::right << "l-1"                << "][   0 <= m <=   ]" << std::endl;  // NOLINT
+  *stream << "       -p p  : order of norm      (double)[" << std::setw(5) << std::right << kDefaultNormOrder    << "][ 1.0 <= p <=   ]" << std::endl;  // NOLINT
+  *stream << "       -h    : print this message" << std::endl;
+  *stream << "  infile:" << std::endl;
+  *stream << "       vector                     (double)[stdin]" << std::endl;
+  *stream << "  stdout:" << std::endl;
+  *stream << "       normalized vector          (double)" << std::endl;
+  *stream << "  notice:" << std::endl;
+  *stream << "       for L-infinity norm, specify -p inf" << std::endl;
+  *stream << std::endl;
+  *stream << " SPTK: version " << sptk::kVersion << std::endl;
+  *stream << std::endl;
+  // clang-format on
+}
+
+}  // namespace
+
+/**
+ * @a lpnorm [ @e option ] [ @e infile ]
+ *
+ * - @b -l @e int
+ *   - length of vector
+ * - @b -m @e int
+ *   - order of vector
+ * - @b -p @e double or @e str
+ *   - order of norm
+ * - @b infile @e str
+ *   - double-type vector sequence
+ * - @b stdout
+ *   - double-type normalized vector sequence
+ *
+ * The below example calculates L2-normalized vector sequence:
+ *
+ * @code{.sh}
+ *   echo 3 4 | x2x +ad | lpnorm -l 2 -p 2 | x2x +da
+ *   # 0.6 0.8
+ * @endcode
+ *
+ * @param[in] argc Number of arguments.
+ * @param[in] argv Argument vector.
+ * @return 0 on success, 1 on failure.
+ */
+int main(int argc, char* argv[]) {
+  int num_order(kDefaultNumOrder);
+  double norm_order(kDefaultNormOrder);
+
+  for (;;) {
+    const int option_char(getopt_long(argc, argv, "l:m:p:h", NULL, NULL));
+    if (-1 == option_char) break;
+
+    switch (option_char) {
+      case 'l': {
+        if (!sptk::ConvertStringToInteger(optarg, &num_order) ||
+            num_order <= 0) {
+          std::ostringstream error_message;
+          error_message
+              << "The argument for the -l option must be a positive integer";
+          sptk::PrintErrorMessage("lpnorm", error_message);
+          return 1;
+        }
+        --num_order;
+        break;
+      }
+      case 'm': {
+        if (!sptk::ConvertStringToInteger(optarg, &num_order) ||
+            num_order < 0) {
+          std::ostringstream error_message;
+          error_message << "The argument for the -m option must be a "
+                        << "non-negative integer";
+          sptk::PrintErrorMessage("lpnorm", error_message);
+          return 1;
+        }
+        break;
+      }
+      case 'p': {
+        if ((!sptk::ConvertSpecialStringToDouble(optarg, &norm_order) &&
+             !sptk::ConvertStringToDouble(optarg, &norm_order)) ||
+            norm_order < 1.0) {
+          std::ostringstream error_message;
+          error_message << "The argument for the -p option must be equal to or "
+                        << "greater than 1.0";
+          sptk::PrintErrorMessage("lpnorm", error_message);
+          return 1;
+        }
+        break;
+      }
+      case 'h': {
+        PrintUsage(&std::cout);
+        return 0;
+      }
+      default: {
+        PrintUsage(&std::cerr);
+        return 1;
+      }
+    }
+  }
+
+  const int num_input_files(argc - optind);
+  if (1 < num_input_files) {
+    std::ostringstream error_message;
+    error_message << "Too many input files";
+    sptk::PrintErrorMessage("lpnorm", error_message);
+    return 1;
+  }
+  const char* input_file(0 == num_input_files ? NULL : argv[optind]);
+
+  if (!sptk::SetBinaryMode()) {
+    std::ostringstream error_message;
+    error_message << "Cannot set translation mode";
+    sptk::PrintErrorMessage("lpnorm", error_message);
+    return 1;
+  }
+
+  std::ifstream ifs;
+  if (NULL != input_file) {
+    ifs.open(input_file, std::ios::in | std::ios::binary);
+    if (ifs.fail()) {
+      std::ostringstream error_message;
+      error_message << "Cannot open file " << input_file;
+      sptk::PrintErrorMessage("lpnorm", error_message);
+      return 1;
+    }
+  }
+  std::istream& input_stream(ifs.is_open() ? ifs : std::cin);
+
+  sptk::LpNormCalculation lp_norm_calculation(num_order, norm_order);
+  if (!lp_norm_calculation.IsValid()) {
+    std::ostringstream error_message;
+    error_message << "Failed to initialize LpNormCalculation";
+    sptk::PrintErrorMessage("lpnorm", error_message);
+    return 1;
+  }
+
+  const int vector_length(num_order + 1);
+  std::vector<double> vector(vector_length);
+  std::vector<double> normalized_vector(vector_length);
+
+  while (sptk::ReadStream(false, 0, 0, vector_length, &vector, &input_stream,
+                          NULL)) {
+    double norm;
+    if (!lp_norm_calculation.Run(vector, &norm)) {
+      std::ostringstream error_message;
+      error_message << "Failed to calculate Lp-norm";
+      sptk::PrintErrorMessage("lpnorm", error_message);
+      return 1;
+    }
+
+    std::transform(vector.begin(), vector.end(), normalized_vector.begin(),
+                   [norm](double x) { return x / norm; });
+
+    if (!sptk::WriteStream(0, vector_length, normalized_vector, &std::cout,
+                           NULL)) {
+      std::ostringstream error_message;
+      error_message << "Failed to write normalized vector";
+      sptk::PrintErrorMessage("lpnorm", error_message);
+      return 1;
+    }
+  }
+
+  return 0;
+}

src/main/sopr.cc

@@ -116,6 +116,8 @@ void PrintUsage(std::ostream* stream) {
   *stream << "           sqrtX    :      sqrt(X)  [ 0.0 <= X <=   ]" << std::endl;  // NOLINT
   *stream << "           lnX      :        ln(X)  [ 0.0 <  X <=   ]" << std::endl;  // NOLINT
   *stream << "           expX     :       exp(X)  [     <= X <=   ]" << std::endl;  // NOLINT
+  *stream << "           inf      : infinity" << std::endl;
+  *stream << "           nan      : not a number" << std::endl;
   *stream << "" << std::endl;
   *stream << "       they are case-insensitive" << std::endl;
   *stream << "" << std::endl;

src/math/lp_norm_calculation.cc

@@ -0,0 +1,75 @@
+// ------------------------------------------------------------------------ //
+// Copyright 2021 SPTK Working Group                                        //
+//                                                                          //
+// Licensed under the Apache License, Version 2.0 (the "License");          //
+// you may not use this file except in compliance with the License.         //
+// You may obtain a copy of the License at                                  //
+//                                                                          //
+//     http://www.apache.org/licenses/LICENSE-2.0                           //
+//                                                                          //
+// Unless required by applicable law or agreed to in writing, software      //
+// distributed under the License is distributed on an "AS IS" BASIS,        //
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. //
+// See the License for the specific language governing permissions and      //
+// limitations under the License.                                           //
+// ------------------------------------------------------------------------ //
+
+#include "SPTK/math/lp_norm_calculation.h"
+
+#include <algorithm>  // std::max_element
+#include <cmath>      // std::abs, std::pow, std::sqrt
+#include <cstddef>    // std::size_t
+#include <limits>     // std::numeric_limits
+#include <numeric>    // std::inner_product
+#include <vector>     // std::vector
+
+namespace sptk {
+
+LpNormCalculation::LpNormCalculation(int num_order, double p)
+    : num_order_(num_order), p_(p), is_valid_(true) {
+  if (num_order_ < 0 || p_ < 1.0) {
+    is_valid_ = false;
+    return;
+  }
+}
+
+bool LpNormCalculation::Run(const std::vector<double>& data,
+                            double* norm) const {
+  if (!is_valid_ || data.size() != static_cast<std::size_t>(num_order_ + 1) ||
+      NULL == norm) {
+    return false;
+  }
+
+  // Specialize for p = 2 as it is the most common case and can be computed
+  // faster than the general case.
+  if (p_ == 2.0) {
+    *norm = std::sqrt(
+        std::inner_product(data.begin(), data.end(), data.begin(), 0.0));
+    return true;
+  }
+
+  const double max_value(std::abs(*std::max_element(
+      data.begin(), data.end(),
+      [](double a, double b) { return std::abs(a) < std::abs(b); })));
+
+  if (0.0 == max_value) {
+    *norm = 0.0;
+    return true;
+  }
+
+  if (p_ == std::numeric_limits<double>::infinity()) {
+    *norm = max_value;
+    return true;
+  }
+
+  const double* input(&(data[0]));
+  double sum(0.0);
+  for (int i(0); i <= num_order_; ++i) {
+    sum += std::pow(std::abs(input[i]) / max_value, p_);
+  }
+  *norm = max_value * std::pow(sum, 1.0 / p_);
+
+  return true;
+}
+
+}  // namespace sptk