diff --git a/ChangeLog b/ChangeLog index 24e28841..d73006e9 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,38 @@ +2025-07-02 Dirk Eddelbuettel + + * DESCRIPTION (Version, Date): RcppArmadillo 14.6.0-1 + * inst/NEWS.Rd: Idem + * configure.ac: Idem + * configure: Idem + + * inst/include/armadillo/: Armadillo 14.6.0 + * inst/include/armadillo_bits/: Idem + +2025-07-01 Dirk Eddelbuettel + + * DESCRIPTION (Version, Date): RcppArmadillo 14.5.92-1 + * inst/NEWS.Rd: Idem + * configure.ac: Idem + * configure: Idem + + * inst/include/armadillo/: Armadillo 14.5.92 as RC2 for 14.6.0 + * inst/include/armadillo_bits/: Idem + +2025-06-30 Dirk Eddelbuettel + + * DESCRIPTION (Version, Date): RcppArmadillo 14.5.91-1 + * inst/NEWS.Rd: Idem + * configure.ac: Idem + * configure: Idem + + * inst/include/armadillo/: Armadillo 14.5.91 as RC1 for 14.6.0 + * inst/include/armadillo_bits/: Idem + +2025-06-21 Dirk Eddelbuettel + + * man/fastLm.Rd: Update discussion of rank-deficient systems to note + that newer Armadillo versions warn and fall back to approx solution + 2025-05-21 Dirk Eddelbuettel * DESCRIPTION (Version, Date): RcppArmadillo 14.4.3-1 diff --git a/DESCRIPTION b/DESCRIPTION index a930ea60..5cd56ea7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: RcppArmadillo Type: Package Title: 'Rcpp' Integration for the 'Armadillo' Templated Linear Algebra Library -Version: 14.4.3-1 -Date: 2025-05-21 +Version: 14.6.0-1 +Date: 2025-07-02 Authors@R: c(person("Dirk", "Eddelbuettel", role = c("aut", "cre"), email = "edd@debian.org", comment = c(ORCID = "0000-0001-6419-907X")), person("Romain", "Francois", role = "aut", diff --git a/configure b/configure index e2b9c49d..a0061227 100755 --- a/configure +++ b/configure @@ -1,6 +1,6 @@ #! /bin/sh # Guess values for system-dependent variables and create Makefiles. -# Generated by GNU Autoconf 2.72 for RcppArmadillo 14.4.3-1. +# Generated by GNU Autoconf 2.72 for RcppArmadillo 14.6.0-1. # # Report bugs to . # @@ -603,8 +603,8 @@ MAKEFLAGS= # Identity of this package. PACKAGE_NAME='RcppArmadillo' PACKAGE_TARNAME='rcpparmadillo' -PACKAGE_VERSION='14.4.3-1' -PACKAGE_STRING='RcppArmadillo 14.4.3-1' +PACKAGE_VERSION='14.6.0-1' +PACKAGE_STRING='RcppArmadillo 14.6.0-1' PACKAGE_BUGREPORT='edd@debian.org' PACKAGE_URL='' @@ -1222,7 +1222,7 @@ if test "$ac_init_help" = "long"; then # Omit some internal or obsolete options to make the list less imposing. # This message is too long to be a string in the A/UX 3.1 sh. cat <<_ACEOF -'configure' configures RcppArmadillo 14.4.3-1 to adapt to many kinds of systems. +'configure' configures RcppArmadillo 14.6.0-1 to adapt to many kinds of systems. Usage: $0 [OPTION]... [VAR=VALUE]... @@ -1284,7 +1284,7 @@ fi if test -n "$ac_init_help"; then case $ac_init_help in - short | recursive ) echo "Configuration of RcppArmadillo 14.4.3-1:";; + short | recursive ) echo "Configuration of RcppArmadillo 14.6.0-1:";; esac cat <<\_ACEOF @@ -1365,7 +1365,7 @@ fi test -n "$ac_init_help" && exit $ac_status if $ac_init_version; then cat <<\_ACEOF -RcppArmadillo configure 14.4.3-1 +RcppArmadillo configure 14.6.0-1 generated by GNU Autoconf 2.72 Copyright (C) 2023 Free Software Foundation, Inc. @@ -1481,7 +1481,7 @@ cat >config.log <<_ACEOF This file contains any messages produced by compilers while running configure, to aid debugging if configure makes a mistake. -It was created by RcppArmadillo $as_me 14.4.3-1, which was +It was created by RcppArmadillo $as_me 14.6.0-1, which was generated by GNU Autoconf 2.72. Invocation command line was $ $0$ac_configure_args_raw @@ -3980,7 +3980,7 @@ cat >>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1 # report actual input values of CONFIG_FILES etc. instead of their # values after options handling. ac_log=" -This file was extended by RcppArmadillo $as_me 14.4.3-1, which was +This file was extended by RcppArmadillo $as_me 14.6.0-1, which was generated by GNU Autoconf 2.72. Invocation command line was CONFIG_FILES = $CONFIG_FILES @@ -4035,7 +4035,7 @@ ac_cs_config_escaped=`printf "%s\n" "$ac_cs_config" | sed "s/^ //; s/'/'\\\\\\\\ cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1 ac_cs_config='$ac_cs_config_escaped' ac_cs_version="\\ -RcppArmadillo config.status 14.4.3-1 +RcppArmadillo config.status 14.6.0-1 configured by $0, generated by GNU Autoconf 2.72, with options \\"\$ac_cs_config\\" diff --git a/configure.ac b/configure.ac index 77da8d17..9776afcf 100644 --- a/configure.ac +++ b/configure.ac @@ -11,7 +11,7 @@ AC_PREREQ([2.69]) ## Process this file with autoconf to produce a configure script. -AC_INIT([RcppArmadillo],[14.4.3-1],[edd@debian.org]) +AC_INIT([RcppArmadillo],[14.6.0-1],[edd@debian.org]) ## Set R_HOME, respecting an environment variable if one is set : ${R_HOME=$(R RHOME)} diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd index caa763f4..03014978 100644 --- a/inst/NEWS.Rd +++ b/inst/NEWS.Rd @@ -3,6 +3,22 @@ \newcommand{\ghpr}{\href{https://github.com/RcppCore/RcppArmadillo/pull/#1}{##1}} \newcommand{\ghit}{\href{https://github.com/RcppCore/RcppArmadillo/issues/#1}{##1}} +\section{Changes in RcppArmadillo version 14.6.0-1 (2025-07-02)}{ + \itemize{ + \item Upgraded to Armadillo release 14.6.0 (Caffe Mocha) + \itemize{ + \item Added \code{balance()} to transform matrices so that column and + row norms are roughly the same + \item Added \code{omit_nan()} and \code{omit_nonfinite()} to extract + elements while omitting NaN and non-finite values + \item Added \code{find_nonnan()} for finding indices of non-NaN elements + \item Added standalone \code{replace()} function + } + \item The \code{fastLm()} help page now mentions that options to + \code{solve()} can control its behavior. + } +} + \section{Changes in RcppArmadillo version 14.4.3-1 (2025-05-21)}{ \itemize{ \item Upgraded to Armadillo release 14.4.3 (Filtered Espresso) @@ -20,9 +36,7 @@ \item Fix for \code{expmat()} \item Workaround for bugs in clang 20 compiler } - \itemize{ - \item Micro-cleanup in one test file - } + \item Micro-cleanup in one test file } } diff --git a/inst/include/armadillo b/inst/include/armadillo index da24cdf2..497ceec1 100644 --- a/inst/include/armadillo +++ b/inst/include/armadillo @@ -304,6 +304,7 @@ namespace arma #include "armadillo_bits/op_clamp_bones.hpp" #include "armadillo_bits/op_expmat_bones.hpp" #include "armadillo_bits/op_nonzeros_bones.hpp" + #include "armadillo_bits/op_omit_bones.hpp" #include "armadillo_bits/op_diff_bones.hpp" #include "armadillo_bits/op_norm_bones.hpp" #include "armadillo_bits/op_vecnorm_bones.hpp" @@ -377,6 +378,7 @@ namespace arma #include "armadillo_bits/spop_norm_bones.hpp" #include "armadillo_bits/spop_shift_bones.hpp" #include "armadillo_bits/spop_relational_bones.hpp" + #include "armadillo_bits/spop_omit_bones.hpp" #include "armadillo_bits/spglue_plus_bones.hpp" #include "armadillo_bits/spglue_minus_bones.hpp" @@ -564,6 +566,7 @@ namespace arma #include "armadillo_bits/fn_clamp.hpp" #include "armadillo_bits/fn_expmat.hpp" #include "armadillo_bits/fn_nonzeros.hpp" + #include "armadillo_bits/fn_omit.hpp" #include "armadillo_bits/fn_interp1.hpp" #include "armadillo_bits/fn_interp2.hpp" #include "armadillo_bits/fn_qz.hpp" @@ -590,6 +593,7 @@ namespace arma #include "armadillo_bits/fn_powmat.hpp" #include "armadillo_bits/fn_powext.hpp" #include "armadillo_bits/fn_diags_spdiags.hpp" + #include "armadillo_bits/fn_balance.hpp" #include "armadillo_bits/fn_speye.hpp" #include "armadillo_bits/fn_spones.hpp" @@ -758,6 +762,7 @@ namespace arma #include "armadillo_bits/op_clamp_meat.hpp" #include "armadillo_bits/op_expmat_meat.hpp" #include "armadillo_bits/op_nonzeros_meat.hpp" + #include "armadillo_bits/op_omit_meat.hpp" #include "armadillo_bits/op_diff_meat.hpp" #include "armadillo_bits/op_norm_meat.hpp" #include "armadillo_bits/op_vecnorm_meat.hpp" @@ -831,6 +836,7 @@ namespace arma #include "armadillo_bits/spop_norm_meat.hpp" #include "armadillo_bits/spop_shift_meat.hpp" #include "armadillo_bits/spop_relational_meat.hpp" + #include "armadillo_bits/spop_omit_meat.hpp" #include "armadillo_bits/spglue_plus_meat.hpp" #include "armadillo_bits/spglue_minus_meat.hpp" diff --git a/inst/include/armadillo_bits/BaseCube_meat.hpp b/inst/include/armadillo_bits/BaseCube_meat.hpp index 7fbf7f67..a106f8e3 100644 --- a/inst/include/armadillo_bits/BaseCube_meat.hpp +++ b/inst/include/armadillo_bits/BaseCube_meat.hpp @@ -328,7 +328,7 @@ BaseCube::is_finite() const for(uword c=0; c::has_nonfinite() const for(uword c=0; c::is_finite() const for(uword i=0; i::is_finite() const for(uword col=0; col::has_nonfinite() const for(uword i=0; i::has_nonfinite() const for(uword col=0; col::operator=(const Cube& x) arrayops::copy( memptr(), x.mem, n_elem ); } + else + { + arma_debug_print("Cube::operator=(): copy omitted"); + } return *this; } diff --git a/inst/include/armadillo_bits/MapMat_meat.hpp b/inst/include/armadillo_bits/MapMat_meat.hpp index 5a99578c..8a767138 100644 --- a/inst/include/armadillo_bits/MapMat_meat.hpp +++ b/inst/include/armadillo_bits/MapMat_meat.hpp @@ -106,7 +106,12 @@ MapMat::operator=(const MapMat& x) { arma_debug_sigprint(); - if(this == &x) { return; } + if(this == &x) + { + arma_debug_print("MapMat::operator=(): copy omitted"); + + return; + } access::rw(n_rows) = x.n_rows; access::rw(n_cols) = x.n_cols; diff --git a/inst/include/armadillo_bits/Mat_meat.hpp b/inst/include/armadillo_bits/Mat_meat.hpp index fac13bbc..5dc3dbce 100644 --- a/inst/include/armadillo_bits/Mat_meat.hpp +++ b/inst/include/armadillo_bits/Mat_meat.hpp @@ -915,6 +915,10 @@ Mat::operator=(const Mat& in_mat) arrayops::copy( memptr(), in_mat.mem, in_mat.n_elem ); } + else + { + arma_debug_print("Mat::operator=(): copy omitted"); + } return *this; } diff --git a/inst/include/armadillo_bits/SpBase_meat.hpp b/inst/include/armadillo_bits/SpBase_meat.hpp index d96cd051..25034952 100644 --- a/inst/include/armadillo_bits/SpBase_meat.hpp +++ b/inst/include/armadillo_bits/SpBase_meat.hpp @@ -748,7 +748,7 @@ SpBase::is_finite() const while(it != it_end) { - if(arma_isfinite(*it) == false) { return false; } + if(arma_isnonfinite(*it)) { return false; } ++it; } } @@ -850,7 +850,7 @@ SpBase::has_nonfinite() const while(it != it_end) { - if(arma_isfinite(*it) == false) { return true; } + if(arma_isnonfinite(*it)) { return true; } ++it; } } diff --git a/inst/include/armadillo_bits/SpMat_meat.hpp b/inst/include/armadillo_bits/SpMat_meat.hpp index 34d92024..9f381b27 100644 --- a/inst/include/armadillo_bits/SpMat_meat.hpp +++ b/inst/include/armadillo_bits/SpMat_meat.hpp @@ -5124,7 +5124,12 @@ SpMat::init(const SpMat& x) { arma_debug_sigprint(); - if(this == &x) { return; } + if(this == &x) + { + arma_debug_print("SpMat::init(): copy omitted"); + + return; + } bool init_done = false; diff --git a/inst/include/armadillo_bits/arma_cmath.hpp b/inst/include/armadillo_bits/arma_cmath.hpp index 22df4bf0..7ac283b4 100644 --- a/inst/include/armadillo_bits/arma_cmath.hpp +++ b/inst/include/armadillo_bits/arma_cmath.hpp @@ -65,6 +65,48 @@ arma_isfinite(const std::complex& x) } +// + + +template +inline +bool +arma_isnonfinite(eT) + { + return false; + } + + + +template<> +inline +bool +arma_isnonfinite(float x) + { + return (std::isfinite(x) == false); + } + + + +template<> +inline +bool +arma_isnonfinite(double x) + { + return (std::isfinite(x) == false); + } + + + +template +inline +bool +arma_isnonfinite(const std::complex& x) + { + return ( (std::isfinite(x.real()) == false) || (std::isfinite(x.imag()) == false) ); + } + + // // wrappers for isinf diff --git a/inst/include/armadillo_bits/arma_forward.hpp b/inst/include/armadillo_bits/arma_forward.hpp index 0eb857ee..20ee5032 100644 --- a/inst/include/armadillo_bits/arma_forward.hpp +++ b/inst/include/armadillo_bits/arma_forward.hpp @@ -93,6 +93,7 @@ class op_vectorise_row; class op_vectorise_col; class op_symmatu; class op_symmatl; +class op_omit; class op_row_as_mat; class op_col_as_mat; @@ -361,6 +362,7 @@ struct arma_zeros_indicator : public arma_initmode_indicator {}; struct arma_nozeros_indicator : public arma_initmode_indicator {}; + //! \addtogroup injector //! @{ diff --git a/inst/include/armadillo_bits/arma_ostream_meat.hpp b/inst/include/armadillo_bits/arma_ostream_meat.hpp index bf79a1b7..f6acf572 100644 --- a/inst/include/armadillo_bits/arma_ostream_meat.hpp +++ b/inst/include/armadillo_bits/arma_ostream_meat.hpp @@ -70,7 +70,7 @@ arma_ostream::modify_stream(std::ostream& o, const eT* data, const uword n_elem) { const eT val = data[i]; - if(arma_isfinite(val) == false) { continue; } + if(arma_isnonfinite(val)) { continue; } if( ( cond_rel< (sizeof(eT) > 4) && (is_same_type::yes || is_same_type::yes) >::geq(val, eT(+10000000000)) ) @@ -208,7 +208,7 @@ arma_ostream::modify_stream(std::ostream& o, typename SpMat::const_iterator { const eT val = (*it); - if(arma_isfinite(val) == false) { continue; } + if(arma_isnonfinite(val)) { continue; } if( val >= eT(+100) || diff --git a/inst/include/armadillo_bits/arma_version.hpp b/inst/include/armadillo_bits/arma_version.hpp index bff877e8..75c7e99c 100644 --- a/inst/include/armadillo_bits/arma_version.hpp +++ b/inst/include/armadillo_bits/arma_version.hpp @@ -22,9 +22,9 @@ #define ARMA_VERSION_MAJOR 14 -#define ARMA_VERSION_MINOR 4 -#define ARMA_VERSION_PATCH 3 -#define ARMA_VERSION_NAME "Filtered Espresso" +#define ARMA_VERSION_MINOR 6 +#define ARMA_VERSION_PATCH 0 +#define ARMA_VERSION_NAME "Caffe Mocha" diff --git a/inst/include/armadillo_bits/arrayops_meat.hpp b/inst/include/armadillo_bits/arrayops_meat.hpp index 825a292b..d20f958b 100644 --- a/inst/include/armadillo_bits/arrayops_meat.hpp +++ b/inst/include/armadillo_bits/arrayops_meat.hpp @@ -1050,13 +1050,13 @@ arrayops::is_finite(const eT* src, const uword n_elem) const eT val_i = (*src); src++; const eT val_j = (*src); src++; - if(arma_isfinite(val_i) == false) { return false; } - if(arma_isfinite(val_j) == false) { return false; } + if(arma_isnonfinite(val_i)) { return false; } + if(arma_isnonfinite(val_j)) { return false; } } if((j-1) < n_elem) { - if(arma_isfinite(*src) == false) { return false; } + if(arma_isnonfinite(*src)) { return false; } } return true; diff --git a/inst/include/armadillo_bits/auxlib_bones.hpp b/inst/include/armadillo_bits/auxlib_bones.hpp index f68ba1f5..aa073e3c 100644 --- a/inst/include/armadillo_bits/auxlib_bones.hpp +++ b/inst/include/armadillo_bits/auxlib_bones.hpp @@ -391,7 +391,7 @@ class auxlib // solve the Sylvester equation AX + XB = C template - inline static bool syl(Mat& X, const Mat& A, const Mat& B, const Mat& C); + inline static bool sylvester(Mat& X, const Mat& A, const Mat& B, const Mat& C); // @@ -404,6 +404,13 @@ class auxlib inline static bool qz(Mat< std::complex >& A, Mat< std::complex >& B, Mat< std::complex >& vsl, Mat< std::complex >& vsr, const Base< std::complex, T1 >& X_expr, const Base< std::complex, T2 >& Y_expr, const char mode); + // + // matrix balance + + template + inline static bool balance(Col::result>& S, Col& P, Mat& A, const bool calc_SP, const bool do_scal, const bool do_perm); + + // // rcond diff --git a/inst/include/armadillo_bits/auxlib_meat.hpp b/inst/include/armadillo_bits/auxlib_meat.hpp index eb33318f..38613ff4 100644 --- a/inst/include/armadillo_bits/auxlib_meat.hpp +++ b/inst/include/armadillo_bits/auxlib_meat.hpp @@ -6486,15 +6486,15 @@ auxlib::schur(Mat< std::complex >& U, Mat< std::complex >& S, const bool c template inline bool -auxlib::syl(Mat& X, const Mat& A, const Mat& B, const Mat& C) +auxlib::sylvester(Mat& X, const Mat& A, const Mat& B, const Mat& C) { arma_debug_sigprint(); #if defined(ARMA_USE_LAPACK) { - arma_conform_check( (A.is_square() == false) || (B.is_square() == false), "syl(): given matrices must be square sized" ); + arma_conform_check( (A.is_square() == false) || (B.is_square() == false), "sylvester(): given matrices must be square sized" ); - arma_conform_check( (C.n_rows != A.n_rows) || (C.n_cols != B.n_cols), "syl(): matrices are not conformant" ); + arma_conform_check( (C.n_rows != A.n_rows) || (C.n_cols != B.n_cols), "sylvester(): matrices are not conformant" ); if(A.is_empty() || B.is_empty() || C.is_empty()) { X.reset(); return true; } @@ -6534,7 +6534,7 @@ auxlib::syl(Mat& X, const Mat& A, const Mat& B, const Mat& C) arma_ignore(A); arma_ignore(B); arma_ignore(C); - arma_stop_logic_error("syl(): use of LAPACK must be enabled"); + arma_stop_logic_error("sylvester(): use of LAPACK must be enabled"); return false; } #endif @@ -6714,6 +6714,76 @@ auxlib::qz(Mat< std::complex >& A, Mat< std::complex >& B, Mat< std::compl +template +inline +bool +auxlib::balance(Col::result>& S, Col& P, Mat& A, const bool calc_SP, const bool do_scal, const bool do_perm) + { + arma_debug_sigprint(); + + #if defined(ARMA_USE_LAPACK) + { + typedef typename get_pod_type::result T; + + // assuming given matrix is square-sized + + if(A.n_elem == 0) { S.reset(); P.reset(); return true; } + + const char job = (do_scal && do_perm) ? 'B' : ((do_scal) ? 'S' : ((do_perm) ? 'P' : 'N')); + + blas_int n = blas_int(A.n_rows); + blas_int lda = blas_int(A.n_rows); + blas_int ilo = blas_int(0); + blas_int ihi = blas_int(0); + blas_int info = blas_int(0); + + podarray scale(A.n_rows); scale.zeros(); + + arma_debug_print("lapack::gebal()"); + lapack::gebal(&job, &n, A.memptr(), &lda, &ilo, &ihi, scale.memptr(), &info); + + if(info != blas_int(0)) { return false; } + + if(calc_SP == false) { return true; } + + const uword N = A.n_rows; + + // sanity check + if( (ilo < 1) || (uword(ihi) > N) ) { arma_debug_print("ilo and/or ihi out of bounds"); return false; } + + S.zeros(N); + P.zeros(N); + + T* S_mem = S.memptr(); + uword* P_mem = P.memptr(); + + const T* scale_mem = scale.memptr(); + + for(uword i = 0; i < uword(ilo)-1; ++i) { S_mem[i] = T(1); } + for(uword i = uword(ilo)-1; i < uword(ihi); ++i) { S_mem[i] = scale_mem[i]; } + for(uword i = uword(ihi); i < N; ++i) { S_mem[i] = T(1); } + + for(uword i=0; i < N; ++i) { P_mem[i] = i; } + + for(uword i=N-1; i >= uword(ihi) ; --i) { const uword j = uword(scale_mem[i]) - 1; std::swap(P_mem[i], P_mem[j]); } + for(uword i=0; i < uword(ilo)-1; ++i) { const uword j = uword(scale_mem[i]) - 1; std::swap(P_mem[i], P_mem[j]); } + + return true; + } + #else + { + arma_ignore(S); + arma_ignore(P); + arma_ignore(A); + arma_ignore(do_scal); + arma_ignore(do_perm); + return false; + } + #endif + } + + + template inline eT diff --git a/inst/include/armadillo_bits/compiler_setup.hpp b/inst/include/armadillo_bits/compiler_setup.hpp index a5034d0e..71bdb25a 100644 --- a/inst/include/armadillo_bits/compiler_setup.hpp +++ b/inst/include/armadillo_bits/compiler_setup.hpp @@ -169,7 +169,7 @@ // https://gcc.gnu.org/bugzilla/show_bug.cgi?id=57580 #if (ARMA_GCC_VERSION < 80100) - #pragma message("INFO: support for GCC versions older than 8.1 is deprecated" + #pragma message("INFO: support for GCC versions older than 8.1 is deprecated") #endif #if (ARMA_GCC_VERSION >= 170000) diff --git a/inst/include/armadillo_bits/def_lapack.hpp b/inst/include/armadillo_bits/def_lapack.hpp index 1f918b00..dcf8f133 100644 --- a/inst/include/armadillo_bits/def_lapack.hpp +++ b/inst/include/armadillo_bits/def_lapack.hpp @@ -293,6 +293,11 @@ #define arma_checon checon #define arma_zhecon zhecon + #define arma_sgebal sgebal + #define arma_dgebal dgebal + #define arma_cgebal cgebal + #define arma_zgebal zgebal + #else #define arma_sgetrf SGETRF @@ -553,6 +558,11 @@ #define arma_checon CHECON #define arma_zhecon ZHECON + #define arma_sgebal SGEBAL + #define arma_dgebal DGEBAL + #define arma_cgebal CGEBAL + #define arma_zgebal ZGEBAL + #endif @@ -926,6 +936,12 @@ extern "C" void arma_fortran(arma_checon)(const char* uplo, const blas_int* n, const blas_cxf* a, const blas_int* lda, const blas_int* ipiv, const float* anorm, float* rcond, blas_cxf* work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT; void arma_fortran(arma_zhecon)(const char* uplo, const blas_int* n, const blas_cxd* a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, blas_cxd* work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT; + // matrix balance + void arma_fortran(arma_sgebal)(const char* job, const blas_int* n, float* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, float* scale, blas_int* info, blas_len job_len) ARMA_NOEXCEPT; + void arma_fortran(arma_dgebal)(const char* job, const blas_int* n, double* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, double* scale, blas_int* info, blas_len job_len) ARMA_NOEXCEPT; + void arma_fortran(arma_cgebal)(const char* job, const blas_int* n, blas_cxf* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, float* scale, blas_int* info, blas_len job_len) ARMA_NOEXCEPT; + void arma_fortran(arma_zgebal)(const char* job, const blas_int* n, blas_cxd* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, double* scale, blas_int* info, blas_len job_len) ARMA_NOEXCEPT; + #else // prototypes without hidden arguments @@ -1282,6 +1298,12 @@ extern "C" void arma_fortran(arma_checon)(const char* uplo, const blas_int* n, const blas_cxf* a, const blas_int* lda, const blas_int* ipiv, const float* anorm, float* rcond, blas_cxf* work, blas_int* info) ARMA_NOEXCEPT; void arma_fortran(arma_zhecon)(const char* uplo, const blas_int* n, const blas_cxd* a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, blas_cxd* work, blas_int* info) ARMA_NOEXCEPT; + // matrix balance + void arma_fortran(arma_sgebal)(const char* job, const blas_int* n, float* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, float* scale, blas_int* info) ARMA_NOEXCEPT; + void arma_fortran(arma_dgebal)(const char* job, const blas_int* n, double* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, double* scale, blas_int* info) ARMA_NOEXCEPT; + void arma_fortran(arma_cgebal)(const char* job, const blas_int* n, blas_cxf* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, float* scale, blas_int* info) ARMA_NOEXCEPT; + void arma_fortran(arma_zgebal)(const char* job, const blas_int* n, blas_cxd* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, double* scale, blas_int* info) ARMA_NOEXCEPT; + #endif } diff --git a/inst/include/armadillo_bits/diagview_meat.hpp b/inst/include/armadillo_bits/diagview_meat.hpp index d9c08179..9c65bb10 100644 --- a/inst/include/armadillo_bits/diagview_meat.hpp +++ b/inst/include/armadillo_bits/diagview_meat.hpp @@ -995,9 +995,13 @@ diagview::randu() const uword local_n_elem = n_elem; + Col tmp(local_n_elem, arma_nozeros_indicator()); + + tmp.randu(); + for(uword ii=0; ii < local_n_elem; ++ii) { - x.at(ii+row_offset, ii+col_offset) = eT(arma_rng::randu()); + x.at(ii+row_offset, ii+col_offset) = tmp[ii]; } } @@ -1014,9 +1018,13 @@ diagview::randn() const uword local_n_elem = n_elem; + Col tmp(local_n_elem, arma_nozeros_indicator()); + + tmp.randn(); + for(uword ii=0; ii < local_n_elem; ++ii) { - x.at(ii+row_offset, ii+col_offset) = eT(arma_rng::randn()); + x.at(ii+row_offset, ii+col_offset) = tmp[ii]; } } diff --git a/inst/include/armadillo_bits/field_meat.hpp b/inst/include/armadillo_bits/field_meat.hpp index de613be9..6f7903d0 100644 --- a/inst/include/armadillo_bits/field_meat.hpp +++ b/inst/include/armadillo_bits/field_meat.hpp @@ -2120,7 +2120,12 @@ field::init(const field& x) { arma_debug_sigprint(); - if(this == &x) { return; } + if(this == &x) + { + arma_debug_print("field::init(): copy omitted"); + + return; + } field& t = (*this); diff --git a/inst/include/armadillo_bits/fn_accu.hpp b/inst/include/armadillo_bits/fn_accu.hpp index 14a9b428..250f2785 100644 --- a/inst/include/armadillo_bits/fn_accu.hpp +++ b/inst/include/armadillo_bits/fn_accu.hpp @@ -272,6 +272,99 @@ accu(const T1& X) +template +inline +typename T1::elem_type +accu_op_omit_helper(const Proxy& P, functor is_omitted) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + constexpr eT eT_zero = eT(0); + + eT acc = eT(0); + + if(Proxy::use_at) + { + const uword n_rows = P.get_n_rows(); + const uword n_cols = P.get_n_cols(); + + for(uword c=0; c < n_cols; ++c) + for(uword r=0; r < n_rows; ++r) + { + const eT val = P.at(r,c); + + acc += is_omitted(val) ? eT_zero : val; + } + } + else + { + typename Proxy::ea_type Pea = P.get_ea(); + + const uword n_elem = P.get_n_elem(); + + eT val1 = eT(0); + eT val2 = eT(0); + + uword i,j; + for(i=0, j=1; j < n_elem; i+=2, j+=2) + { + const eT tmp_i = Pea[i]; + const eT tmp_j = Pea[j]; + + val1 += is_omitted(tmp_i) ? eT_zero : tmp_i; + val2 += is_omitted(tmp_j) ? eT_zero : tmp_j; + } + + if(i < n_elem) + { + const eT tmp_i = Pea[i]; + + val1 += is_omitted(tmp_i) ? eT_zero : tmp_i; + } + + acc = val1 + val2; + } + + return acc; + } + + + +template +arma_warn_unused +inline +typename T1::elem_type +accu(const Op& in) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + const uword omit_mode = in.aux_uword_a; + + if(arma_config::fast_math_warn) + { + if(omit_mode == 1) { arma_warn(1, "omit_nan(): detection of NaN is not reliable in fast math mode"); } + if(omit_mode == 2) { arma_warn(1, "omit_nonfinite(): detection of non-finite values is not reliable in fast math mode"); } + } + + auto is_omitted_1 = [](const eT& x) -> bool { return arma_isnan(x); }; + auto is_omitted_2 = [](const eT& x) -> bool { return arma_isnonfinite(x); }; + + const Proxy P(in.m); + + eT acc = eT(0); + + if(omit_mode == 1) { acc = accu_op_omit_helper(P, is_omitted_1); } + if(omit_mode == 2) { acc = accu_op_omit_helper(P, is_omitted_2); } + + return acc; + } + + + template arma_warn_unused inline @@ -1005,6 +1098,101 @@ accu(const eGlueCube& expr) +template +inline +typename T1::elem_type +accu_cube_omit_helper(const ProxyCube& P, functor is_omitted) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + constexpr eT eT_zero = eT(0); + + eT acc = eT(0); + + if(ProxyCube::use_at) + { + const uword n_r = P.get_n_rows(); + const uword n_c = P.get_n_cols(); + const uword n_s = P.get_n_slices(); + + for(uword s=0; s < n_s; ++s) + for(uword c=0; c < n_c; ++c) + for(uword r=0; r < n_r; ++r) + { + const eT val = P.at(r,c,s); + + acc += is_omitted(val) ? eT_zero : val; + } + } + else + { + typename ProxyCube::ea_type Pea = P.get_ea(); + + const uword n_elem = P.get_n_elem(); + + eT val1 = eT(0); + eT val2 = eT(0); + + uword i,j; + for(i=0, j=1; j < n_elem; i+=2, j+=2) + { + const eT tmp_i = Pea[i]; + const eT tmp_j = Pea[j]; + + val1 += is_omitted(tmp_i) ? eT_zero : tmp_i; + val2 += is_omitted(tmp_j) ? eT_zero : tmp_j; + } + + if(i < n_elem) + { + const eT tmp_i = Pea[i]; + + val1 += is_omitted(tmp_i) ? eT_zero : tmp_i; + } + + acc = val1 + val2; + } + + return acc; + } + + + +template +arma_warn_unused +inline +typename T1::elem_type +accu(const CubeToMatOp& in) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + const ProxyCube P(in.m); + + const uword omit_mode = in.aux_uword; + + if(arma_config::fast_math_warn) + { + if(omit_mode == 1) { arma_warn(1, "omit_nan(): detection of NaN is not reliable in fast math mode"); } + if(omit_mode == 2) { arma_warn(1, "omit_nonfinite(): detection of non-finite values is not reliable in fast math mode"); } + } + + auto is_omitted_1 = [](const eT& x) -> bool { return arma_isnan(x); }; + auto is_omitted_2 = [](const eT& x) -> bool { return arma_isnonfinite(x); }; + + eT acc = eT(0); + + if(omit_mode == 1) { acc = accu_cube_omit_helper(P, is_omitted_1); } + if(omit_mode == 2) { acc = accu_cube_omit_helper(P, is_omitted_2); } + + return acc; + } + + + // @@ -1199,33 +1387,172 @@ accu(const SpOp& expr) if(is_vectorise) { return accu(expr.m); } - if(is_same_type::yes) + const SpMat tmp = expr; + + return accu(tmp); + } + + + +template +arma_warn_unused +inline +typename T1::elem_type +accu(const SpOp& expr) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + if(is_SpSubview_col::value) + { + const SpSubview_col& svcol = reinterpret_cast&>(expr.m); + + if(svcol.n_nonzero == 0) { return eT(0); } + + if(svcol.n_rows == svcol.m.n_rows) + { + arma_debug_print("accu(): SpSubview_col spop_square optimisation"); + + const SpMat& m = svcol.m; + const uword col = svcol.aux_col1; + + const eT* ptr = &(m.values[ m.col_ptrs[col] ]); + + return op_dot::direct_dot(svcol.n_nonzero, ptr, ptr); + } + } + + const SpProxy P(expr.m); + + const uword N = P.get_n_nonzero(); + + if(N == 0) { return eT(0); } + + if(SpProxy::use_iterator == false) { - const SpProxy P(expr.m); + return op_dot::direct_dot(N, P.get_values(), P.get_values()); + } + else + { + typename SpProxy::const_iterator_type it = P.begin(); - const uword N = P.get_n_nonzero(); + eT acc = eT(0); - if(N == 0) { return eT(0); } + for(uword i=0; i < N; ++i) { const eT tmp = (*it); acc += (tmp*tmp); ++it; } - if(SpProxy::use_iterator == false) + return acc; + } + } + + + +template +inline +typename T1::elem_type +accu_spop_omit_helper(const T1& expr, functor is_omitted) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + constexpr eT eT_zero = eT(0); + + if(is_SpSubview_col::value) + { + const SpSubview_col& svcol = reinterpret_cast&>(expr); + + if(svcol.n_nonzero == 0) { return eT(0); } + + if(svcol.n_rows == svcol.m.n_rows) { - return op_dot::direct_dot(N, P.get_values(), P.get_values()); + arma_debug_print("accu_spop_omit_helper(): SpSubview_col optimisation"); + + const SpMat& m = svcol.m; + const uword col = svcol.aux_col1; + + const eT* vals = &(m.values[ m.col_ptrs[col] ]); + + const uword N = svcol.n_nonzero; + + eT acc = eT(0); + + for(uword i=0; i < N; ++i) + { + const eT tmp = vals[i]; + + acc += is_omitted(tmp) ? eT_zero : tmp; + } + + return acc; } - else + } + + const SpProxy P(expr); + + const uword N = P.get_n_nonzero(); + + if(N == 0) { return eT(0); } + + eT acc = eT(0); + + if(SpProxy::use_iterator == false) + { + const eT* vals = P.get_values(); + + for(uword i=0; i < N; ++i) { - typename SpProxy::const_iterator_type it = P.begin(); + const eT tmp = vals[i]; - eT val = eT(0); + acc += is_omitted(tmp) ? eT_zero : tmp; + } + } + else + { + typename SpProxy::const_iterator_type it = P.begin(); + + for(uword i=0; i < N; ++i) + { + const eT tmp = (*it); - for(uword i=0; i < N; ++i) { const eT tmp = (*it); val += (tmp*tmp); ++it; } + acc += is_omitted(tmp) ? eT_zero : tmp; - return val; + ++it; } } - const SpMat tmp = expr; + return acc; + } + + + +template +arma_warn_unused +inline +typename T1::elem_type +accu(const SpOp& expr) + { + arma_debug_sigprint(); - return accu(tmp); + typedef typename T1::elem_type eT; + + const uword omit_mode = expr.aux_uword_a; + + if(arma_config::fast_math_warn) + { + if(omit_mode == 1) { arma_warn(1, "omit_nan(): detection of NaN is not reliable in fast math mode"); } + if(omit_mode == 2) { arma_warn(1, "omit_nonfinite(): detection of non-finite values is not reliable in fast math mode"); } + } + + auto is_omitted_1 = [](const eT& x) -> bool { return arma_isnan(x); }; + auto is_omitted_2 = [](const eT& x) -> bool { return arma_isnonfinite(x); }; + + eT acc = eT(0); + + if(omit_mode == 1) { acc = accu_spop_omit_helper(expr.m, is_omitted_1); } + if(omit_mode == 2) { acc = accu_spop_omit_helper(expr.m, is_omitted_2); } + + return acc; } diff --git a/inst/include/armadillo_bits/fn_balance.hpp b/inst/include/armadillo_bits/fn_balance.hpp new file mode 100644 index 00000000..142c6f75 --- /dev/null +++ b/inst/include/armadillo_bits/fn_balance.hpp @@ -0,0 +1,135 @@ +// SPDX-License-Identifier: Apache-2.0 +// +// Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au) +// Copyright 2008-2016 National ICT Australia (NICTA) +// +// 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. +// ------------------------------------------------------------------------ + + +//! \addtogroup fn_balance +//! @{ + + + +template +inline +typename enable_if2< is_supported_blas_type::value, bool >::result +balance(Col& S, Col& P, Mat& B, const Base& A, const char* method = "both") + { + arma_debug_sigprint(); + + arma_conform_check( (void_ptr(&S) == void_ptr(&B)), "eig_gen(): parameter 'S' is an alias of parameter 'B'" ); + + const char sig = (method != nullptr) ? method[0] : char(0); + + if( (sig != 'b') && (sig != 's') && (sig != 'p') ) { arma_stop_logic_error("balance(): unsupported method"); } + + const bool do_scale = (sig == 'b') || (sig == 's'); + const bool do_perm = (sig == 'b') || (sig == 'p'); + + const bool calc_SP = true; + + B = A.get_ref(); + + if(B.is_square() == false) + { + B.soft_reset(); + + arma_stop_logic_error("balance(): given matrix must be square sized"); + + return false; + } + + const bool status = auxlib::balance(S, P, B, calc_SP, do_scale, do_perm); + + if(status == false) + { + S.soft_reset(); + P.soft_reset(); + B.soft_reset(); + + arma_warn(3, "balance(): transformation failed"); + } + + return status; + } + + + +template +inline +typename enable_if2< is_supported_blas_type::value, bool >::result +balance(Mat& B, const Base& A, const char* method = "both") + { + arma_debug_sigprint(); + + typedef typename T1::pod_type T; + + const char sig = (method != nullptr) ? method[0] : char(0); + + if( (sig != 'b') && (sig != 's') && (sig != 'p') ) { arma_stop_logic_error("balance(): unsupported method"); } + + const bool do_scale = (sig == 'b') || (sig == 's'); + const bool do_perm = (sig == 'b') || (sig == 'p'); + + const bool calc_SP = false; + + B = A.get_ref(); + + if(B.is_square() == false) + { + B.soft_reset(); + + arma_stop_logic_error("balance(): given matrix must be square sized"); + + return false; + } + + Col S; + Col P; + + const bool status = auxlib::balance(S, P, B, calc_SP, do_scale, do_perm); + + if(status == false) + { + B.soft_reset(); + + arma_warn(3, "balance(): transformation failed"); + } + + return status; + } + + + +template +inline +typename enable_if2< is_supported_blas_type::value, Mat >::result +balance(const Base& A, const char* method = "both") + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + Mat B; + + const bool status = balance(B, A, method); + + if(status == false) { arma_stop_runtime_error("balance(): transformation failed"); } + + return B; + } + + + +//! @} diff --git a/inst/include/armadillo_bits/fn_elem.hpp b/inst/include/armadillo_bits/fn_elem.hpp index ad4b0d11..1031f50f 100644 --- a/inst/include/armadillo_bits/fn_elem.hpp +++ b/inst/include/armadillo_bits/fn_elem.hpp @@ -545,6 +545,45 @@ arg(const T1& X) +template +arma_warn_unused +inline +typename enable_if2< is_arma_type::value, const mtOp >::result +replace(const T1& X, typename T1::elem_type old_val, typename T1::elem_type new_val) + { + arma_debug_sigprint(); + + return mtOp(mtOp_dual_aux_indicator(), X, old_val, new_val); + } + + + +template +arma_warn_unused +inline +const mtOpCube +replace(const BaseCube& X, typename T1::elem_type old_val, typename T1::elem_type new_val) + { + arma_debug_sigprint(); + + return mtOpCube(mtOpCube_dual_aux_indicator(), X.get_ref(), old_val, new_val); + } + + + +template +arma_warn_unused +inline +typename enable_if2< is_arma_sparse_type::value, const mtSpOp >::result +replace(const T1& X, typename T1::elem_type old_val, typename T1::elem_type new_val) + { + arma_debug_sigprint(); + + return mtSpOp(mtSpOp_dual_aux_indicator(), X, old_val, new_val); + } + + + // // square diff --git a/inst/include/armadillo_bits/fn_find.hpp b/inst/include/armadillo_bits/fn_find.hpp index 5d5d46be..588389cd 100644 --- a/inst/include/armadillo_bits/fn_find.hpp +++ b/inst/include/armadillo_bits/fn_find.hpp @@ -269,6 +269,23 @@ find_nan(const T1& X) +template +arma_warn_unused +inline +typename enable_if2 + < + is_arma_type::value, + const mtOp + >::result +find_nonnan(const T1& X) + { + arma_debug_sigprint(); + + return mtOp(X); + } + + + // @@ -330,6 +347,25 @@ find_nan(const BaseCube& X) +template +arma_warn_unused +inline +uvec +find_nonnan(const BaseCube& X) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + const unwrap_cube tmp(X.get_ref()); + + const Mat R( const_cast< eT* >(tmp.M.memptr()), tmp.M.n_elem, 1, false ); + + return find_nonnan(R); + } + + + // @@ -401,7 +437,7 @@ find_nonfinite(const SpBase& X) for(uword i=0; i& X) +template +arma_warn_unused +inline +Col +find_nonnan(const SpBase& X) + { + arma_debug_sigprint(); + + const SpProxy P(X.get_ref()); + + const uword n_rows = P.get_n_rows(); + const uword n_nz = P.get_n_nonzero(); + + Mat tmp(n_nz, 1, arma_nozeros_indicator()); + + uword* tmp_mem = tmp.memptr(); + + typename SpProxy::const_iterator_type it = P.begin(); + + uword count = 0; + + for(uword i=0; i out; + + if(count > 0) { out.steal_mem_col(tmp, count); } + + return out; + } + + + //! @} diff --git a/inst/include/armadillo_bits/fn_misc.hpp b/inst/include/armadillo_bits/fn_misc.hpp index fcef65a4..c20c5f79 100644 --- a/inst/include/armadillo_bits/fn_misc.hpp +++ b/inst/include/armadillo_bits/fn_misc.hpp @@ -557,7 +557,7 @@ namespace priv const eT negdelta = log_b - log_a; - if( (negdelta < Datum::log_min) || (arma_isfinite(negdelta) == false) ) + if( (negdelta < Datum::log_min) || arma_isnonfinite(negdelta) ) { return log_a; } diff --git a/inst/include/armadillo_bits/fn_omit.hpp b/inst/include/armadillo_bits/fn_omit.hpp new file mode 100644 index 00000000..036ccaa2 --- /dev/null +++ b/inst/include/armadillo_bits/fn_omit.hpp @@ -0,0 +1,104 @@ +// SPDX-License-Identifier: Apache-2.0 +// +// Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au) +// Copyright 2008-2016 National ICT Australia (NICTA) +// +// 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. +// ------------------------------------------------------------------------ + + +//! \addtogroup fn_omit +//! @{ + + + +template +arma_warn_unused +inline +typename enable_if2< is_arma_type::value, const Op >::result +omit_nan(const T1& X) + { + arma_debug_sigprint(); + + return Op(X, 1, 0); + } + + + +template +arma_warn_unused +inline +typename enable_if2< is_arma_type::value, const Op >::result +omit_nonfinite(const T1& X) + { + arma_debug_sigprint(); + + return Op(X, 2, 0); + } + + + +template +arma_warn_unused +inline +CubeToMatOp +omit_nan(const BaseCube& X) + { + arma_debug_sigprint(); + + return CubeToMatOp(X.get_ref(), 1); + } + + + +template +arma_warn_unused +inline +CubeToMatOp +omit_nonfinite(const BaseCube& X) + { + arma_debug_sigprint(); + + return CubeToMatOp(X.get_ref(), 2); + } + + + +template +arma_warn_unused +inline +typename +enable_if2< is_arma_sparse_type::value, const SpOp >::result +omit_nan(const T1& X) + { + arma_debug_sigprint(); + + return SpOp(X, 1, 0); + } + + + +template +arma_warn_unused +inline +typename +enable_if2< is_arma_sparse_type::value, const SpOp >::result +omit_nonfinite(const T1& X) + { + arma_debug_sigprint(); + + return SpOp(X, 2, 0); + } + + + +//! @} diff --git a/inst/include/armadillo_bits/fn_powext.hpp b/inst/include/armadillo_bits/fn_powext.hpp index ae629c64..6538fe13 100644 --- a/inst/include/armadillo_bits/fn_powext.hpp +++ b/inst/include/armadillo_bits/fn_powext.hpp @@ -44,7 +44,7 @@ pow template -arma_deprecated +arma_frown("refactor your code to use pow() in conjunction with repmat()") inline Mat pow @@ -121,7 +121,7 @@ pow template -arma_deprecated +arma_frown("refactor your code to use pow() in conjunction with repmat()") inline typename enable_if2 diff --git a/inst/include/armadillo_bits/fn_randperm.hpp b/inst/include/armadillo_bits/fn_randperm.hpp index 5fb55bb7..6713bfd0 100644 --- a/inst/include/armadillo_bits/fn_randperm.hpp +++ b/inst/include/armadillo_bits/fn_randperm.hpp @@ -37,9 +37,18 @@ internal_randperm_helper(obj_type& x, const uword N, const uword N_keep) std::vector packet_vec(N); + podarray tmp(N); + + int* tmp_mem = tmp.memptr(); + + const int a = 0; + const int b = arma_rng::randi::max_val(); + + arma_rng::randi::fill(tmp_mem, N, a, b); + for(uword i=0; i < N; ++i) { - packet_vec[i].val = int(arma_rng::randi()); + packet_vec[i].val = tmp_mem[i]; packet_vec[i].index = i; } diff --git a/inst/include/armadillo_bits/fn_shift.hpp b/inst/include/armadillo_bits/fn_shift.hpp index 3ee13a21..a7f6b489 100644 --- a/inst/include/armadillo_bits/fn_shift.hpp +++ b/inst/include/armadillo_bits/fn_shift.hpp @@ -21,6 +21,9 @@ //! @{ +// TODO: deprecate shift() in favour of circshift() + + template arma_warn_unused arma_inline diff --git a/inst/include/armadillo_bits/fn_sort_index.hpp b/inst/include/armadillo_bits/fn_sort_index.hpp index 6796a9b6..8d53effe 100644 --- a/inst/include/armadillo_bits/fn_sort_index.hpp +++ b/inst/include/armadillo_bits/fn_sort_index.hpp @@ -67,10 +67,11 @@ sort_index +// DO NOT USE: kept only for compatibility with old user code template arma_warn_unused arma_inline -const mtOp +const mtOp stable_sort_index ( const Base& X @@ -78,11 +79,12 @@ stable_sort_index { arma_debug_sigprint(); - return mtOp(X.get_ref(), uword(0), uword(0)); + return mtOp(X.get_ref(), uword(0), uword(0)); } +// DO NOT USE: kept only for compatibility with old user code template arma_warn_unused inline @@ -90,7 +92,7 @@ typename enable_if2 < ( (is_arma_type::value) && (is_same_type::value) ), - const mtOp + const mtOp >::result stable_sort_index ( @@ -104,7 +106,7 @@ stable_sort_index arma_conform_check( ((sig != 'a') && (sig != 'd')), "stable_sort_index(): unknown sort direction" ); - return mtOp(X, ((sig == 'a') ? uword(0) : uword(1)), uword(0)); + return mtOp(X, ((sig == 'a') ? uword(0) : uword(1)), uword(0)); } diff --git a/inst/include/armadillo_bits/fn_sylvester.hpp b/inst/include/armadillo_bits/fn_sylvester.hpp index 323c9e86..f180549f 100644 --- a/inst/include/armadillo_bits/fn_sylvester.hpp +++ b/inst/include/armadillo_bits/fn_sylvester.hpp @@ -24,7 +24,7 @@ template inline bool -syl +sylvester ( Mat & out, const Base& in_A, @@ -46,12 +46,12 @@ syl const Mat& B = tmp_B.M; const Mat& C = tmp_C.M; - const bool status = auxlib::syl(out, A, B, C); + const bool status = auxlib::sylvester(out, A, B, C); if(status == false) { out.soft_reset(); - arma_warn(3, "syl(): solution not found"); + arma_warn(3, "sylvester(): solution not found"); } return status; @@ -59,10 +59,12 @@ syl +// kept for compatibility with old user code template +arma_frown("use sylvester() instead") inline bool -sylvester +syl ( Mat & out, const Base& in_A, @@ -72,7 +74,8 @@ sylvester ) { arma_ignore(junk); - return syl(out, in_A, in_B, in_C); + + return sylvester(out, in_A, in_B, in_C); } @@ -81,7 +84,7 @@ template arma_warn_unused inline Mat -syl +sylvester ( const Base& in_A, const Base& in_B, @@ -104,12 +107,12 @@ syl Mat out; - const bool status = auxlib::syl(out, A, B, C); + const bool status = auxlib::sylvester(out, A, B, C); if(status == false) { out.soft_reset(); - arma_stop_runtime_error("syl(): solution not found"); + arma_stop_runtime_error("sylvester(): solution not found"); } return out; @@ -117,11 +120,12 @@ syl +// kept for compatibility with old user code template -arma_warn_unused +arma_frown("use sylvester() instead") inline Mat -sylvester +syl ( const Base& in_A, const Base& in_B, @@ -130,7 +134,8 @@ sylvester ) { arma_ignore(junk); - return syl(in_A, in_B, in_C); + + return sylvester(in_A, in_B, in_C); } diff --git a/inst/include/armadillo_bits/glue_mvnrnd_meat.hpp b/inst/include/armadillo_bits/glue_mvnrnd_meat.hpp index 4b61820b..62f15f6f 100644 --- a/inst/include/armadillo_bits/glue_mvnrnd_meat.hpp +++ b/inst/include/armadillo_bits/glue_mvnrnd_meat.hpp @@ -139,13 +139,13 @@ glue_mvnrnd::apply_noalias(Mat& out, const Mat& M, const Mat& C, con const eT tol = eT(-100) * Datum::eps * norm(C, "fro"); - if(arma_isfinite(tol) == false) { return false; } + if(arma_isnonfinite(tol)) { return false; } for(uword i=0; i::em_iterate(const Mat& X, const uword max_iter, const eT var_fl get_cout_stream().flush(); } - if(arma_isfinite(new_avg_log_p) == false) { return false; } + if(arma_isnonfinite(new_avg_log_p)) { return false; } if(std::abs(old_avg_log_p - new_avg_log_p) <= Datum::eps) { break; } @@ -2473,7 +2473,7 @@ gmm_diag::em_update_params { const eT acc_norm_lhood = (std::max)( final_acc_norm_lhoods[g], std::numeric_limits::min() ); - if(arma_isfinite(acc_norm_lhood) == false) { continue; } + if(arma_isnonfinite(acc_norm_lhood)) { continue; } eT* acc_mean_mem = final_acc_means.colptr(g); eT* acc_dcov_mem = final_acc_dcovs.colptr(g); @@ -2488,7 +2488,7 @@ gmm_diag::em_update_params acc_mean_mem[d] = tmp1; acc_dcov_mem[d] = tmp2; - if(arma_isfinite(tmp2) == false) { ok = false; } + if(arma_isnonfinite(tmp2)) { ok = false; } } diff --git a/inst/include/armadillo_bits/gmm_full_meat.hpp b/inst/include/armadillo_bits/gmm_full_meat.hpp index 7ce1a009..736d0bd3 100644 --- a/inst/include/armadillo_bits/gmm_full_meat.hpp +++ b/inst/include/armadillo_bits/gmm_full_meat.hpp @@ -2392,7 +2392,7 @@ gmm_full::em_iterate(const Mat& X, const uword max_iter, const eT var_fl get_cout_stream().flush(); } - if(arma_isfinite(new_avg_log_p) == false) { return false; } + if(arma_isnonfinite(new_avg_log_p)) { return false; } if(std::abs(old_avg_log_p - new_avg_log_p) <= Datum::eps) { break; } @@ -2516,7 +2516,7 @@ gmm_full::em_update_params { const eT acc_norm_lhood = (std::max)( final_acc_norm_lhoods[g], std::numeric_limits::min() ); - if(arma_isfinite(acc_norm_lhood) == false) { continue; } + if(arma_isnonfinite(acc_norm_lhood)) { continue; } eT* acc_mean_mem = final_acc_means.colptr(g); diff --git a/inst/include/armadillo_bits/mtSpOp_bones.hpp b/inst/include/armadillo_bits/mtSpOp_bones.hpp index 335b25d8..88a7378e 100644 --- a/inst/include/armadillo_bits/mtSpOp_bones.hpp +++ b/inst/include/armadillo_bits/mtSpOp_bones.hpp @@ -21,6 +21,9 @@ +struct mtSpOp_dual_aux_indicator {}; + + template class mtSpOp : public SpBase< out_eT, mtSpOp > { @@ -39,6 +42,7 @@ class mtSpOp : public SpBase< out_eT, mtSpOp > inline mtSpOp(const T1& in_m, const in_eT in_aux); inline mtSpOp(const T1& in_m, const uword aux_uword_a, const uword aux_uword_b); inline mtSpOp(const char junk, const T1& in_m, const out_eT in_aux); + inline mtSpOp(const mtSpOp_dual_aux_indicator&, const T1& in_m, const in_eT in_aux_a, const out_eT in_aux_b); inline ~mtSpOp(); template diff --git a/inst/include/armadillo_bits/mtSpOp_meat.hpp b/inst/include/armadillo_bits/mtSpOp_meat.hpp index 59d2f3eb..49263f18 100644 --- a/inst/include/armadillo_bits/mtSpOp_meat.hpp +++ b/inst/include/armadillo_bits/mtSpOp_meat.hpp @@ -67,6 +67,18 @@ mtSpOp::mtSpOp(const char junk, const T1& in_m, const out_e +template +inline +mtSpOp::mtSpOp(const mtSpOp_dual_aux_indicator&, const T1& in_m, const typename T1::elem_type in_aux_a, const out_eT in_aux_b) + : m (in_m ) + , aux (in_aux_a) + , aux_out_eT(in_aux_b) + { + arma_debug_sigprint(); + } + + + template inline mtSpOp::~mtSpOp() diff --git a/inst/include/armadillo_bits/op_expmat_meat.hpp b/inst/include/armadillo_bits/op_expmat_meat.hpp index 9377d8e1..578fc1ef 100644 --- a/inst/include/armadillo_bits/op_expmat_meat.hpp +++ b/inst/include/armadillo_bits/op_expmat_meat.hpp @@ -120,7 +120,7 @@ op_expmat::apply_direct(Mat& out, const Base + inline static void apply(Mat& out, const mtOp& X); + }; + + + //! @} diff --git a/inst/include/armadillo_bits/op_find_meat.hpp b/inst/include/armadillo_bits/op_find_meat.hpp index 6654add2..e64df527 100644 --- a/inst/include/armadillo_bits/op_find_meat.hpp +++ b/inst/include/armadillo_bits/op_find_meat.hpp @@ -585,7 +585,7 @@ op_find_nonfinite::apply(Mat& out, const mtOp& out, const mtOp& out, const mtOp& X) +template +inline +void +op_find_nonnan::apply(Mat& out, const mtOp& X) + { + arma_debug_sigprint(); + + if(arma_config::fast_math_warn) { arma_warn(1, "find_nonnan(): detection of non-finite values is not reliable in fast math mode"); } + + const Proxy P(X.m); + + const uword n_elem = P.get_n_elem(); + + Mat indices(n_elem, 1, arma_nozeros_indicator()); + + uword* indices_mem = indices.memptr(); + uword count = 0; + + if(Proxy::use_at == false) + { + const typename Proxy::ea_type Pea = P.get_ea(); + + for(uword i=0; i < n_elem; ++i) + { + if( arma_isnan(Pea[i]) == false ) { indices_mem[count] = i; ++count; } + } + } + else + { + const uword n_rows = P.get_n_rows(); + const uword n_cols = P.get_n_cols(); + + uword i = 0; + + for(uword col=0; col < n_cols; ++col) + for(uword row=0; row < n_rows; ++row) + { + if( arma_isnan(P.at(row,col)) == false ) { indices_mem[count] = i; ++count; } + + ++i; + } + } + + out.steal_mem_col(indices, count); + } + + + //! @} diff --git a/inst/include/armadillo_bits/op_hist_meat.hpp b/inst/include/armadillo_bits/op_hist_meat.hpp index 27cffa83..675aaa92 100644 --- a/inst/include/armadillo_bits/op_hist_meat.hpp +++ b/inst/include/armadillo_bits/op_hist_meat.hpp @@ -66,8 +66,8 @@ op_hist::apply_noalias(Mat& out, const Mat& A, const uword n_bins, co max_val += (n_bins/2); } - if(arma_isfinite(min_val) == false) { min_val = priv::most_neg(); } - if(arma_isfinite(max_val) == false) { max_val = priv::most_pos(); } + if(arma_isnonfinite(min_val)) { min_val = priv::most_neg(); } + if(arma_isnonfinite(max_val)) { max_val = priv::most_pos(); } Col c(n_bins, arma_nozeros_indicator()); eT* c_mem = c.memptr(); diff --git a/inst/include/armadillo_bits/op_mean_bones.hpp b/inst/include/armadillo_bits/op_mean_bones.hpp index 20a86ae2..764b7976 100644 --- a/inst/include/armadillo_bits/op_mean_bones.hpp +++ b/inst/include/armadillo_bits/op_mean_bones.hpp @@ -31,14 +31,8 @@ class op_mean template inline static void apply(Mat& out, const Op& in); - template - inline static void apply_noalias(Mat& out, const Proxy& P, const uword dim); - - template - inline static void apply_noalias_unwrap(Mat& out, const Proxy& P, const uword dim); - - template - inline static void apply_noalias_proxy(Mat& out, const Proxy& P, const uword dim); + template + inline static void apply_noalias(Mat& out, const Mat& X, const uword dim); // cubes @@ -46,60 +40,28 @@ class op_mean template inline static void apply(Cube& out, const OpCube& in); - template - inline static void apply_noalias(Cube& out, const ProxyCube& P, const uword dim); - - template - inline static void apply_noalias_unwrap(Cube& out, const ProxyCube& P, const uword dim); - - template - inline static void apply_noalias_proxy(Cube& out, const ProxyCube& P, const uword dim); - - - // - - template - inline static eT direct_mean(const eT* const X, const uword N); - template - inline static eT direct_mean_robust(const eT* const X, const uword N); + inline static void apply_noalias(Cube& out, const Cube& X, const uword dim); - - // - - template - inline static eT direct_mean(const Mat& X, const uword row); - - template - inline static eT direct_mean_robust(const Mat& X, const uword row); - - // template - inline static eT mean_all(const subview& X); + inline static eT direct_mean(const eT* X_mem, const uword N); template - inline static eT mean_all_robust(const subview& X); - - - // - - template - inline static eT mean_all(const diagview& X); - - template - inline static eT mean_all_robust(const diagview& X); + inline static eT direct_mean_robust(const eT old_mean, const eT* X_mem, const uword N); // template - inline static typename T1::elem_type mean_all(const Op& X); + inline static typename T1::elem_type mean_all(const T1& X); template - inline static typename T1::elem_type mean_all(const Base& X); + inline static typename T1::elem_type mean_all(const Op& X); + template + inline static eT mean_all_omit(const eT* X_mem, const uword N, functor is_omitted); // @@ -111,5 +73,4 @@ class op_mean }; - //! @} diff --git a/inst/include/armadillo_bits/op_mean_meat.hpp b/inst/include/armadillo_bits/op_mean_meat.hpp index 1bb0c586..354ffc73 100644 --- a/inst/include/armadillo_bits/op_mean_meat.hpp +++ b/inst/include/armadillo_bits/op_mean_meat.hpp @@ -31,60 +31,35 @@ op_mean::apply(Mat& out, const Op& in) typedef typename T1::elem_type eT; const uword dim = in.aux_uword_a; + arma_conform_check( (dim > 1), "mean(): parameter 'dim' must be 0 or 1" ); - const Proxy P(in.m); + const quasi_unwrap U(in.m); - if(P.is_alias(out) == false) - { - op_mean::apply_noalias(out, P, dim); - } - else + if(U.is_alias(out)) { Mat tmp; - op_mean::apply_noalias(tmp, P, dim); + op_mean::apply_noalias(tmp, U.M, dim); out.steal_mem(tmp); } - } - - - -template -inline -void -op_mean::apply_noalias(Mat& out, const Proxy& P, const uword dim) - { - arma_debug_sigprint(); - - if((is_Mat::stored_type>::value) || (arma_config::openmp && Proxy::use_mp)) - { - op_mean::apply_noalias_unwrap(out, P, dim); - } else { - op_mean::apply_noalias_proxy(out, P, dim); + op_mean::apply_noalias(out, U.M, dim); } } -template +template inline void -op_mean::apply_noalias_unwrap(Mat& out, const Proxy& P, const uword dim) +op_mean::apply_noalias(Mat& out, const Mat& X, const uword dim) { arma_debug_sigprint(); - typedef typename T1::elem_type eT; - typedef typename get_pod_type::result T; - - typedef typename Proxy::stored_type P_stored_type; - - const unwrap tmp(P.Q); - - const typename unwrap::stored_type& X = tmp.M; + typedef typename get_pod_type::result T; const uword X_n_rows = X.n_rows; const uword X_n_cols = X.n_cols; @@ -115,96 +90,33 @@ op_mean::apply_noalias_unwrap(Mat& out, const Proxy& { const eT* col_mem = X.colptr(col); - for(uword row=0; row < X_n_rows; ++row) - { - out_mem[row] += col_mem[row]; - } + for(uword row=0; row < X_n_rows; ++row) { out_mem[row] += col_mem[row]; } } out /= T(X_n_cols); - for(uword row=0; row < X_n_rows; ++row) - { - if(arma_isfinite(out_mem[row]) == false) - { - out_mem[row] = op_mean::direct_mean_robust( X, row ); - } - } - } - } - - - -template -inline -void -op_mean::apply_noalias_proxy(Mat& out, const Proxy& P, const uword dim) - { - arma_debug_sigprint(); - - typedef typename T1::elem_type eT; - typedef typename get_pod_type::result T; - - const uword P_n_rows = P.get_n_rows(); - const uword P_n_cols = P.get_n_cols(); - - if(dim == 0) - { - out.set_size((P_n_rows > 0) ? 1 : 0, P_n_cols); - - if(P_n_rows == 0) { return; } - - eT* out_mem = out.memptr(); - - for(uword col=0; col < P_n_cols; ++col) + if(out.internal_has_nonfinite()) { - eT val1 = eT(0); - eT val2 = eT(0); + podarray tmp; - uword i,j; - for(i=0, j=1; j < P_n_rows; i+=2, j+=2) - { - val1 += P.at(i,col); - val2 += P.at(j,col); - } - - if(i < P_n_rows) + for(uword row=0; row < X_n_rows; ++row) { - val1 += P.at(i,col); + const eT old_mean = out_mem[row]; + + if(arma_isnonfinite(old_mean)) + { + tmp.copy_row(X, row); + + out_mem[row] = op_mean::direct_mean_robust(old_mean, tmp.memptr(), tmp.n_elem); + } } - - out_mem[col] = (val1 + val2) / T(P_n_rows); } } - else - if(dim == 1) - { - out.zeros(P_n_rows, (P_n_cols > 0) ? 1 : 0); - - if(P_n_cols == 0) { return; } - - eT* out_mem = out.memptr(); - - for(uword col=0; col < P_n_cols; ++col) - for(uword row=0; row < P_n_rows; ++row) - { - out_mem[row] += P.at(row,col); - } - - out /= T(P_n_cols); - } - - if(out.internal_has_nonfinite()) - { - // TODO: replace with dedicated handling to avoid unwrapping - op_mean::apply_noalias_unwrap(out, P, dim); - } } // -// cubes @@ -218,60 +130,35 @@ op_mean::apply(Cube& out, const OpCube& in) typedef typename T1::elem_type eT; const uword dim = in.aux_uword_a; + arma_conform_check( (dim > 2), "mean(): parameter 'dim' must be 0 or 1 or 2" ); - const ProxyCube P(in.m); + const unwrap_cube U(in.m); - if(P.is_alias(out) == false) - { - op_mean::apply_noalias(out, P, dim); - } - else + if(U.is_alias(out)) { Cube tmp; - op_mean::apply_noalias(tmp, P, dim); + op_mean::apply_noalias(tmp, U.M, dim); out.steal_mem(tmp); } - } - - - -template -inline -void -op_mean::apply_noalias(Cube& out, const ProxyCube& P, const uword dim) - { - arma_debug_sigprint(); - - if((is_Cube::stored_type>::value) || (arma_config::openmp && ProxyCube::use_mp)) - { - op_mean::apply_noalias_unwrap(out, P, dim); - } else { - op_mean::apply_noalias_proxy(out, P, dim); + op_mean::apply_noalias(out, U.M, dim); } } -template +template inline void -op_mean::apply_noalias_unwrap(Cube& out, const ProxyCube& P, const uword dim) +op_mean::apply_noalias(Cube& out, const Cube& X, const uword dim) { arma_debug_sigprint(); - typedef typename T1::elem_type eT; - typedef typename get_pod_type::result T; - - typedef typename ProxyCube::stored_type P_stored_type; - - const unwrap_cube U(P.Q); - - const Cube& X = U.M; + typedef typename get_pod_type::result T; const uword X_n_rows = X.n_rows; const uword X_n_cols = X.n_cols; @@ -308,21 +195,27 @@ op_mean::apply_noalias_unwrap(Cube& out, const ProxyCube { const eT* col_mem = X.slice_colptr(slice,col); - for(uword row=0; row < X_n_rows; ++row) - { - out_mem[row] += col_mem[row]; - } + for(uword row=0; row < X_n_rows; ++row) { out_mem[row] += col_mem[row]; } } - const Mat tmp('j', X.slice_memptr(slice), X_n_rows, X_n_cols); + for(uword row=0; row < X_n_rows; ++row) { out_mem[row] /= T(X_n_cols); } - for(uword row=0; row < X_n_rows; ++row) + if(arrayops::is_finite(out_mem, X_n_rows) == false) { - out_mem[row] /= T(X_n_cols); + const Mat tmp_mat('j', X.slice_memptr(slice), X_n_rows, X_n_cols); + + podarray tmp_vec; - if(arma_isfinite(out_mem[row]) == false) + for(uword row=0; row < X_n_rows; ++row) { - out_mem[row] = op_mean::direct_mean_robust( tmp, row ); + const eT old_mean = out_mem[row]; + + if(arma_isnonfinite(old_mean)) + { + tmp_vec.copy_row(tmp_mat, row); + + out_mem[row] = op_mean::direct_mean_robust(old_mean, tmp_vec.memptr(), tmp_vec.n_elem); + } } } } @@ -343,19 +236,21 @@ op_mean::apply_noalias_unwrap(Cube& out, const ProxyCube out /= T(X_n_slices); - podarray tmp(X_n_slices); - - for(uword col=0; col < X_n_cols; ++col) - for(uword row=0; row < X_n_rows; ++row) + if(out.internal_has_nonfinite()) { - if(arma_isfinite(out.at(row,col,0)) == false) + podarray tmp(X_n_slices); + + for(uword col=0; col < X_n_cols; ++col) + for(uword row=0; row < X_n_rows; ++row) { - for(uword slice=0; slice < X_n_slices; ++slice) + const eT old_mean = out.at(row,col,0); + + if(arma_isnonfinite(old_mean)) { - tmp[slice] = X.at(row,col,slice); + for(uword slice=0; slice < X_n_slices; ++slice) { tmp[slice] = X.at(row,col,slice); } + + out.at(row,col,0) = op_mean::direct_mean_robust(old_mean, tmp.memptr(), tmp.n_elem); } - - out.at(row,col,0) = op_mean::direct_mean_robust(tmp.memptr(), X_n_slices); } } } @@ -363,21 +258,6 @@ op_mean::apply_noalias_unwrap(Cube& out, const ProxyCube -template -inline -void -op_mean::apply_noalias_proxy(Cube& out, const ProxyCube& P, const uword dim) - { - arma_debug_sigprint(); - - op_mean::apply_noalias_unwrap(out, P, dim); - - // TODO: implement specialised handling - } - - - - // @@ -385,15 +265,15 @@ op_mean::apply_noalias_proxy(Cube& out, const ProxyCube< template inline eT -op_mean::direct_mean(const eT* const X, const uword n_elem) +op_mean::direct_mean(const eT* X_mem, const uword N) { arma_debug_sigprint(); typedef typename get_pod_type::result T; - const eT result = arrayops::accumulate(X, n_elem) / T(n_elem); + const eT mean = arrayops::accumulate(X_mem, N) / T(N); - return arma_isfinite(result) ? result : op_mean::direct_mean_robust(X, n_elem); + return arma_isfinite(mean) ? mean : op_mean::direct_mean_robust(mean, X_mem, N); } @@ -401,7 +281,7 @@ op_mean::direct_mean(const eT* const X, const uword n_elem) template inline eT -op_mean::direct_mean_robust(const eT* const X, const uword n_elem) +op_mean::direct_mean_robust(const eT old_mean, const eT* X_mem, const uword N) { arma_debug_sigprint(); @@ -409,25 +289,13 @@ op_mean::direct_mean_robust(const eT* const X, const uword n_elem) typedef typename get_pod_type::result T; - uword i,j; + if(arrayops::is_finite(X_mem, N) == false) { return old_mean; } eT r_mean = eT(0); - for(i=0, j=1; j -inline -eT -op_mean::direct_mean(const Mat& X, const uword row) - { - arma_debug_sigprint(); - - typedef typename get_pod_type::result T; - - const uword X_n_cols = X.n_cols; - - eT val = eT(0); - - uword i,j; - for(i=0, j=1; j < X_n_cols; i+=2, j+=2) - { - val += X.at(row,i); - val += X.at(row,j); - } - - if(i < X_n_cols) - { - val += X.at(row,i); - } - - const eT result = val / T(X_n_cols); - - return arma_isfinite(result) ? result : op_mean::direct_mean_robust(X, row); - } - - - -template -inline -eT -op_mean::direct_mean_robust(const Mat& X, const uword row) - { - arma_debug_sigprint(); - - typedef typename get_pod_type::result T; - - const uword X_n_cols = X.n_cols; - - eT r_mean = eT(0); - - for(uword col=0; col < X_n_cols; ++col) - { - r_mean = r_mean + (X.at(row,col) - r_mean)/T(col+1); - } - - return r_mean; - } +// -template +template inline -eT -op_mean::mean_all(const subview& X) +typename T1::elem_type +op_mean::mean_all(const T1& X) { arma_debug_sigprint(); - typedef typename get_pod_type::result T; + typedef typename T1::elem_type eT; - const uword X_n_rows = X.n_rows; - const uword X_n_cols = X.n_cols; - const uword X_n_elem = X.n_elem; + const quasi_unwrap U(X); - if(X_n_elem == 0) + if(U.M.n_elem == 0) { arma_conform_check(true, "mean(): object has no elements"); return Datum::nan; } - eT val = eT(0); - - if(X_n_rows == 1) - { - const Mat& A = X.m; - - const uword start_row = X.aux_row1; - const uword start_col = X.aux_col1; - - const uword end_col_p1 = start_col + X_n_cols; - - uword i,j; - for(i=start_col, j=start_col+1; j < end_col_p1; i+=2, j+=2) - { - val += A.at(start_row, i); - val += A.at(start_row, j); - } - - if(i < end_col_p1) - { - val += A.at(start_row, i); - } - } - else - { - for(uword col=0; col < X_n_cols; ++col) - { - val += arrayops::accumulate(X.colptr(col), X_n_rows); - } - } - - const eT result = val / T(X_n_elem); - - return arma_isfinite(result) ? result : op_mean::mean_all_robust(X); + return op_mean::direct_mean(U.M.memptr(), U.M.n_elem); } -template -inline -eT -op_mean::mean_all_robust(const subview& X) +template +inline +typename T1::elem_type +op_mean::mean_all(const Op& in) { arma_debug_sigprint(); - typedef typename get_pod_type::result T; - - const uword X_n_rows = X.n_rows; - const uword X_n_cols = X.n_cols; - - const uword start_row = X.aux_row1; - const uword start_col = X.aux_col1; - - const uword end_row_p1 = start_row + X_n_rows; - const uword end_col_p1 = start_col + X_n_cols; - - const Mat& A = X.m; + typedef typename T1::elem_type eT; + const uword omit_mode = in.aux_uword_a; - eT r_mean = eT(0); - - if(X_n_rows == 1) + if(arma_config::fast_math_warn) { - uword i=0; - - for(uword col = start_col; col < end_col_p1; ++col, ++i) - { - r_mean = r_mean + (A.at(start_row,col) - r_mean)/T(i+1); - } - } - else - { - uword i=0; - - for(uword col = start_col; col < end_col_p1; ++col) - for(uword row = start_row; row < end_row_p1; ++row, ++i) - { - r_mean = r_mean + (A.at(row,col) - r_mean)/T(i+1); - } + if(omit_mode == 1) { arma_warn(1, "omit_nan(): detection of NaN is not reliable in fast math mode"); } + if(omit_mode == 2) { arma_warn(1, "omit_nonfinite(): detection of non-finite values is not reliable in fast math mode"); } } - return r_mean; - } - - - -template -inline -eT -op_mean::mean_all(const diagview& X) - { - arma_debug_sigprint(); - - typedef typename get_pod_type::result T; - - const uword X_n_elem = X.n_elem; + const quasi_unwrap U(in.m); - if(X_n_elem == 0) + if(U.M.n_elem == 0) { arma_conform_check(true, "mean(): object has no elements"); return Datum::nan; } - eT val = eT(0); + auto is_omitted_1 = [](const eT& x) -> bool { return arma_isnan(x); }; + auto is_omitted_2 = [](const eT& x) -> bool { return arma_isnonfinite(x); }; - for(uword i=0; i -inline +template +inline eT -op_mean::mean_all_robust(const diagview& X) +op_mean::mean_all_omit(const eT* X_mem, const uword N, functor is_omitted) { arma_debug_sigprint(); typedef typename get_pod_type::result T; - const uword X_n_elem = X.n_elem; - - eT r_mean = eT(0); + uword count = 0; + eT acc = eT(0); - for(uword i=0; i -inline -typename T1::elem_type -op_mean::mean_all(const Op& X) - { - arma_debug_sigprint(); + acc /= T(count); - return op_mean::mean_all(X.m); - } - - - -template -inline -typename T1::elem_type -op_mean::mean_all(const Base& X) - { - arma_debug_sigprint(); + if(arma_isfinite(acc)) { return acc; } - typedef typename T1::elem_type eT; + // handle possible overflow - const quasi_unwrap tmp(X.get_ref()); - const Mat& A = tmp.M; + eT r_mean = eT(0); - const uword A_n_elem = A.n_elem; + count = 0; - if(A_n_elem == 0) + for(uword i=0; i < N; ++i) { - arma_conform_check(true, "mean(): object has no elements"); + const eT val = X_mem[i]; - return Datum::nan; + if(is_omitted(val) == false) + { + r_mean = r_mean + (val - r_mean) / T(count+1); // kept as count+1 to use same algorithm as op_mean::direct_mean_robust() + + ++count; + } } - return op_mean::direct_mean(A.memptr(), A_n_elem); + return r_mean; } +// + + + template arma_inline eT op_mean::robust_mean(const eT A, const eT B) { - return A + (B - A)/eT(2); + return (arma_isfinite(A) && arma_isfinite(B)) ? eT( A + (B - A)/eT(2) ) : eT( (A+B)/eT(2) ); } @@ -704,10 +434,11 @@ arma_inline std::complex op_mean::robust_mean(const std::complex& A, const std::complex& B) { - return A + (B - A)/T(2); + typedef typename std::complex eT; + + return (arma_isfinite(A) && arma_isfinite(B)) ? eT( A + (B - A)/T(2) ) : eT( (A+B)/T(2) ); } //! @} - diff --git a/inst/include/armadillo_bits/op_median_meat.hpp b/inst/include/armadillo_bits/op_median_meat.hpp index a2efd84d..681889b4 100644 --- a/inst/include/armadillo_bits/op_median_meat.hpp +++ b/inst/include/armadillo_bits/op_median_meat.hpp @@ -41,7 +41,7 @@ op_median::apply(Mat& out, const Op& expr) { Mat tmp; - op_median::apply_noalias(out, U.M, dim); + op_median::apply_noalias(tmp, U.M, dim); out.steal_mem(tmp); } diff --git a/inst/include/armadillo_bits/op_misc_bones.hpp b/inst/include/armadillo_bits/op_misc_bones.hpp index 5fd15714..39ac853b 100644 --- a/inst/include/armadillo_bits/op_misc_bones.hpp +++ b/inst/include/armadillo_bits/op_misc_bones.hpp @@ -77,4 +77,16 @@ class op_arg +class op_replace + : public traits_op_passthru + { + public: + + template inline static void apply(Mat& out, const mtOp& in); + + template inline static void apply(Cube& out, const mtOpCube& in); + }; + + + //! @} diff --git a/inst/include/armadillo_bits/op_misc_meat.hpp b/inst/include/armadillo_bits/op_misc_meat.hpp index 69c837ca..f03346d3 100644 --- a/inst/include/armadillo_bits/op_misc_meat.hpp +++ b/inst/include/armadillo_bits/op_misc_meat.hpp @@ -411,4 +411,38 @@ op_arg::apply( Cube& out, const mtOpCube +inline +void +op_replace::apply(Mat& out, const mtOp& in) + { + arma_debug_sigprint(); + + const eT old_val = in.aux; + const eT new_val = in.aux_out_eT; + + out = in.m; + + out.replace(old_val, new_val); + } + + + +template +inline +void +op_replace::apply(Cube& out, const mtOpCube& in) + { + arma_debug_sigprint(); + + const eT old_val = in.aux; + const eT new_val = in.aux_out_eT; + + out = in.m; + + out.replace(old_val, new_val); + } + + + //! @} diff --git a/inst/include/armadillo_bits/op_nonzeros_meat.hpp b/inst/include/armadillo_bits/op_nonzeros_meat.hpp index cdc0b932..0b97a451 100644 --- a/inst/include/armadillo_bits/op_nonzeros_meat.hpp +++ b/inst/include/armadillo_bits/op_nonzeros_meat.hpp @@ -31,6 +31,8 @@ op_nonzeros::apply_noalias(Mat& out, const Proxy& P) typedef typename T1::elem_type eT; + constexpr eT eT_zero = eT(0); + const uword N_max = P.get_n_elem(); Mat tmp(N_max, 1, arma_nozeros_indicator()); @@ -47,7 +49,7 @@ op_nonzeros::apply_noalias(Mat& out, const Proxy& P) { const eT val = Pea[i]; - if(val != eT(0)) { tmp_mem[N_nz] = val; ++N_nz; } + if(val != eT_zero) { tmp_mem[N_nz] = val; ++N_nz; } } } else @@ -60,7 +62,7 @@ op_nonzeros::apply_noalias(Mat& out, const Proxy& P) { const eT val = P.at(row,col); - if(val != eT(0)) { tmp_mem[N_nz] = val; ++N_nz; } + if(val != eT_zero) { tmp_mem[N_nz] = val; ++N_nz; } } } @@ -84,11 +86,11 @@ op_nonzeros::apply(Mat& out, const Op& if(P.is_alias(out)) { - Mat out2; + Mat tmp; - op_nonzeros::apply_noalias(out2, P); + op_nonzeros::apply_noalias(tmp, P); - out.steal_mem(out2); + out.steal_mem(tmp); } else { diff --git a/inst/include/armadillo_bits/op_norm2est_meat.hpp b/inst/include/armadillo_bits/op_norm2est_meat.hpp index bf7d4497..953cb168 100644 --- a/inst/include/armadillo_bits/op_norm2est_meat.hpp +++ b/inst/include/armadillo_bits/op_norm2est_meat.hpp @@ -138,7 +138,7 @@ op_norm2est::norm2est T x_norm = op_norm::vec_norm_2( Proxy< Col >(x) ); - if(x_norm == T(0) || (arma_isfinite(x_norm) == false) || (x.internal_has_nonfinite())) + if( (x_norm == T(0)) || arma_isnonfinite(x_norm) || x.internal_has_nonfinite() ) { randu_filler.fill(x.memptr(), x.n_elem); @@ -155,7 +155,7 @@ op_norm2est::norm2est arma_debug_print(arma_str::format("norm2est(): est_old: %e") % est_old); arma_debug_print(arma_str::format("norm2est(): est_cur: %e") % est_cur); - if(arma_isfinite(est_cur) == false) { return est_old; } + if(arma_isnonfinite(est_cur)) { return est_old; } if( ((std::abs)(est_cur - est_old)) <= (tol * (std::max)(est_cur,est_old)) ) { break; } } @@ -218,7 +218,7 @@ op_norm2est::norm2est T x_norm = op_norm::vec_norm_2( Proxy< Mat >(x) ); - if(x_norm == T(0) || (arma_isfinite(x_norm) == false) || (x.internal_has_nonfinite())) + if( (x_norm == T(0)) || arma_isnonfinite(x_norm) || x.internal_has_nonfinite() ) { randu_filler.fill(x.memptr(), x.n_elem); @@ -235,7 +235,7 @@ op_norm2est::norm2est arma_debug_print(arma_str::format("norm2est(): est_old: %e") % est_old); arma_debug_print(arma_str::format("norm2est(): est_cur: %e") % est_cur); - if(arma_isfinite(est_cur) == false) { return est_old; } + if(arma_isnonfinite(est_cur)) { return est_old; } if( ((std::abs)(est_cur - est_old)) <= (tol * (std::max)(est_cur,est_old)) ) { break; } } diff --git a/inst/include/armadillo_bits/op_omit_bones.hpp b/inst/include/armadillo_bits/op_omit_bones.hpp new file mode 100644 index 00000000..f8dc7911 --- /dev/null +++ b/inst/include/armadillo_bits/op_omit_bones.hpp @@ -0,0 +1,50 @@ +// SPDX-License-Identifier: Apache-2.0 +// +// Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au) +// Copyright 2008-2016 National ICT Australia (NICTA) +// +// 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. +// ------------------------------------------------------------------------ + + +//! \addtogroup op_omit +//! @{ + + +class op_omit + : public traits_op_col + { + public: + + template + inline static void apply(Mat& out, const Op& in); + + template + inline static void apply(Mat& out, const T1& X, functor is_omitted); + }; + + + +class op_omit_cube + : public traits_op_col + { + public: + + template inline static void apply(Mat& out, const CubeToMatOp& in); + + template + inline static void apply(Mat& out, const T1& X, functor is_omitted); + }; + + + +//! @} diff --git a/inst/include/armadillo_bits/op_omit_meat.hpp b/inst/include/armadillo_bits/op_omit_meat.hpp new file mode 100644 index 00000000..16c6cbaa --- /dev/null +++ b/inst/include/armadillo_bits/op_omit_meat.hpp @@ -0,0 +1,230 @@ +// SPDX-License-Identifier: Apache-2.0 +// +// Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au) +// Copyright 2008-2016 National ICT Australia (NICTA) +// +// 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. +// ------------------------------------------------------------------------ + + +//! \addtogroup op_omit +//! @{ + + + +template +inline +void +op_omit::apply(Mat& out, const Op& in) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + const uword omit_mode = in.aux_uword_a; + + if(arma_config::fast_math_warn) + { + if(omit_mode == 1) { arma_warn(1, "omit_nan(): detection of NaN is not reliable in fast math mode"); } + if(omit_mode == 2) { arma_warn(1, "omit_nonfinite(): detection of non-finite values is not reliable in fast math mode"); } + } + + auto is_omitted_1 = [](const eT& x) -> bool { return arma_isnan(x); }; + auto is_omitted_2 = [](const eT& x) -> bool { return arma_isnonfinite(x); }; + + if(omit_mode == 1) { op_omit::apply(out, in.m, is_omitted_1); } + if(omit_mode == 2) { op_omit::apply(out, in.m, is_omitted_2); } + } + + + +template +inline +void +op_omit::apply(Mat& out, const T1& X, functor is_omitted) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + if(is_Mat::value || is_subview_col::value || is_Mat::stored_type>::value || Proxy::use_mp) + { + const quasi_unwrap U(X); + + const eT* X_mem = U.M.memptr(); + const uword N = U.M.n_elem; + + Mat Y(N, 1, arma_nozeros_indicator()); + + eT* Y_mem = Y.memptr(); + + uword count = 0; + + for(uword i=0; i < N; ++i) + { + const eT val = X_mem[i]; + + if(is_omitted(val) == false) { Y_mem[count] = val; ++count; } + } + + out.steal_mem_col(Y, count); + } + else + { + const Proxy P(X); + + const uword N = P.get_n_elem(); + + Mat Y(N, 1, arma_nozeros_indicator()); + + eT* Y_mem = Y.memptr(); + + uword count = 0; + + if(Proxy::use_at == false) + { + const typename Proxy::ea_type Pea = P.get_ea(); + + for(uword i=0; i < N; ++i) + { + const eT val = Pea[i]; + + if(is_omitted(val) == false) { Y_mem[count] = val; ++count; } + } + } + else + { + const uword n_rows = P.get_n_rows(); + const uword n_cols = P.get_n_cols(); + + for(uword c=0; c < n_cols; ++c) + for(uword r=0; r < n_rows; ++r) + { + const eT val = P.at(r,c); + + if(is_omitted(val) == false) { Y_mem[count] = val; ++count; } + } + } + + out.steal_mem_col(Y, count); + } + } + + + +// + + + +template +inline +void +op_omit_cube::apply(Mat& out, const CubeToMatOp& in) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + const uword omit_mode = in.aux_uword; + + if(arma_config::fast_math_warn) + { + if(omit_mode == 1) { arma_warn(1, "omit_nan(): detection of NaN is not reliable in fast math mode"); } + if(omit_mode == 2) { arma_warn(1, "omit_nonfinite(): detection of non-finite values is not reliable in fast math mode"); } + } + + auto is_omitted_1 = [](const eT& x) -> bool { return arma_isnan(x); }; + auto is_omitted_2 = [](const eT& x) -> bool { return arma_isnonfinite(x); }; + + if(omit_mode == 1) { op_omit_cube::apply(out, in.m, is_omitted_1); } + if(omit_mode == 2) { op_omit_cube::apply(out, in.m, is_omitted_2); } + } + + + +template +inline +void +op_omit_cube::apply(Mat& out, const T1& X, functor is_omitted) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + if(is_Cube::value || is_Cube::stored_type>::value || ProxyCube::use_mp) + { + const unwrap_cube U(X); + + const eT* X_mem = U.M.memptr(); + const uword N = U.M.n_elem; + + Mat Y(N, 1, arma_nozeros_indicator()); + + eT* Y_mem = Y.memptr(); + + uword count = 0; + + for(uword i=0; i < N; ++i) + { + const eT val = X_mem[i]; + + if(is_omitted(val) == false) { Y_mem[count] = val; ++count; } + } + + out.steal_mem_col(Y, count); + } + else + { + const ProxyCube P(X); + + const uword N = P.get_n_elem(); + + Mat Y(N, 1, arma_nozeros_indicator()); + + eT* Y_mem = Y.memptr(); + + uword count = 0; + + if(ProxyCube::use_at == false) + { + const typename ProxyCube::ea_type Pea = P.get_ea(); + + for(uword i=0; i < N; ++i) + { + const eT val = Pea[i]; + + if(is_omitted(val) == false) { Y_mem[count] = val; ++count; } + } + } + else + { + const uword n_r = P.get_n_rows(); + const uword n_c = P.get_n_cols(); + const uword n_s = P.get_n_slices(); + + for(uword s=0; s < n_s; ++s) + for(uword c=0; c < n_c; ++c) + for(uword r=0; r < n_r; ++r) + { + const eT val = P.at(r,c,s); + + if(is_omitted(val) == false) { Y_mem[count] = val; ++count; } + } + } + + out.steal_mem_col(Y, count); + } + } + + + +//! @} diff --git a/inst/include/armadillo_bits/op_shuffle_meat.hpp b/inst/include/armadillo_bits/op_shuffle_meat.hpp index d1228ad5..f3092dd8 100644 --- a/inst/include/armadillo_bits/op_shuffle_meat.hpp +++ b/inst/include/armadillo_bits/op_shuffle_meat.hpp @@ -40,9 +40,18 @@ op_shuffle::apply_direct(Mat& out, const Mat& X, const uword dim) std::vector packet_vec(N); + podarray tmp(N); + + int* tmp_mem = tmp.memptr(); + + const int a = 0; + const int b = arma_rng::randi::max_val(); + + arma_rng::randi::fill(tmp_mem, N, a, b); + for(uword i=0; i()); + packet_vec[i].val = tmp_mem[i]; packet_vec[i].index = i; } diff --git a/inst/include/armadillo_bits/op_sort_index_bones.hpp b/inst/include/armadillo_bits/op_sort_index_bones.hpp index 2eb9cd46..513d565b 100644 --- a/inst/include/armadillo_bits/op_sort_index_bones.hpp +++ b/inst/include/armadillo_bits/op_sort_index_bones.hpp @@ -27,10 +27,7 @@ class op_sort_index public: template - static inline bool apply_noalias_proxy(Mat& out, const Proxy& P, const uword sort_mode); - - template - static inline void apply_noalias_mat(Mat& out, const Mat& X, const uword sort_mode); + static inline bool apply_helper(Mat& out, const Proxy& P, const uword sort_mode); template static inline void apply(Mat& out, const mtOp& in); @@ -38,20 +35,6 @@ class op_sort_index -class op_stable_sort_index - : public traits_op_col - { - public: - - template - static inline bool apply_noalias(Mat& out, const Proxy& P, const uword sort_mode); - - template - static inline void apply(Mat& out, const mtOp& in); - }; - - - template struct arma_sort_index_packet { @@ -98,16 +81,6 @@ struct arma_sort_index_helper_ascend< std::complex > { return (std::abs(A.val) < std::abs(B.val)); } - - // inline - // bool - // operator() (const arma_sort_index_packet& A, const arma_sort_index_packet& B) const - // { - // const T abs_A_val = std::abs(A.val); - // const T abs_B_val = std::abs(B.val); - // - // return ( (abs_A_val != abs_B_val) ? (abs_A_val < abs_B_val) : (std::arg(A.val) < std::arg(B.val)) ); - // } }; @@ -123,20 +96,14 @@ struct arma_sort_index_helper_descend< std::complex > { return (std::abs(A.val) > std::abs(B.val)); } - - // inline - // bool - // operator() (const arma_sort_index_packet& A, const arma_sort_index_packet& B) const - // { - // const T abs_A_val = std::abs(A.val); - // const T abs_B_val = std::abs(B.val); - // - // return ( (abs_A_val != abs_B_val) ? (abs_A_val > abs_B_val) : (std::arg(A.val) > std::arg(B.val)) ); - // } }; +// + + + template struct arma_sort_index_helper_prepare { diff --git a/inst/include/armadillo_bits/op_sort_index_meat.hpp b/inst/include/armadillo_bits/op_sort_index_meat.hpp index 356dbd3b..2e229516 100644 --- a/inst/include/armadillo_bits/op_sort_index_meat.hpp +++ b/inst/include/armadillo_bits/op_sort_index_meat.hpp @@ -21,10 +21,10 @@ -template +template inline bool -arma_sort_index_helper(Mat& out, const Proxy& P, const uword sort_mode) +op_sort_index::apply_helper(Mat& out, const Proxy& P, const uword sort_mode) { arma_debug_sigprint(); @@ -41,11 +41,13 @@ arma_sort_index_helper(Mat& out, const Proxy& P, const uword sort_mod if(Proxy::use_at == false) { - for(uword i=0; i::ea_type Pea = P.get_ea(); + + for(uword i=0; i < n_elem; ++i) { - const eT val = P[i]; + const eT val = Pea[i]; - if(arma_isnan(val)) { out.soft_reset(); return false; } + if(arma_isnan(val)) { return false; } packet_vec[i].val = prepare(val); packet_vec[i].index = i; @@ -63,7 +65,7 @@ arma_sort_index_helper(Mat& out, const Proxy& P, const uword sort_mod { const eT val = P.at(row,col); - if(arma_isnan(val)) { out.soft_reset(); return false; } + if(arma_isnan(val)) { return false; } packet_vec[i].val = prepare(val); packet_vec[i].index = i; @@ -79,14 +81,7 @@ arma_sort_index_helper(Mat& out, const Proxy& P, const uword sort_mod arma_sort_index_helper_ascend comparator; - if(sort_stable == false) - { - std::sort( packet_vec.begin(), packet_vec.end(), comparator ); - } - else - { - std::stable_sort( packet_vec.begin(), packet_vec.end(), comparator ); - } + std::stable_sort( packet_vec.begin(), packet_vec.end(), comparator ); } else { @@ -94,14 +89,7 @@ arma_sort_index_helper(Mat& out, const Proxy& P, const uword sort_mod arma_sort_index_helper_descend comparator; - if(sort_stable == false) - { - std::sort( packet_vec.begin(), packet_vec.end(), comparator ); - } - else - { - std::stable_sort( packet_vec.begin(), packet_vec.end(), comparator ); - } + std::stable_sort( packet_vec.begin(), packet_vec.end(), comparator ); } uword* out_mem = out.memptr(); @@ -116,40 +104,6 @@ arma_sort_index_helper(Mat& out, const Proxy& P, const uword sort_mod -// - - - -template -inline -bool -op_sort_index::apply_noalias_proxy(Mat& out, const Proxy& P, const uword sort_mode) - { - arma_debug_sigprint(); - - return arma_sort_index_helper(out, P, sort_mode); - } - - - -template -inline -void -op_sort_index::apply_noalias_mat(Mat& out, const Mat& X, const uword sort_mode) - { - arma_debug_sigprint(); - - if(X.n_elem == 0) { out.set_size(0,1); return; } - - const Proxy< Mat > P(X); - - const bool all_non_nan = op_sort_index::apply_noalias_proxy(out, P, sort_mode); - - arma_conform_check( (all_non_nan == false), "sort_index(): detected NaN" ); - } - - - template inline void @@ -169,65 +123,18 @@ op_sort_index::apply(Mat& out, const mtOp& in) { Mat tmp; - all_non_nan = op_sort_index::apply_noalias_proxy(tmp, P, sort_mode); + all_non_nan = op_sort_index::apply_helper(tmp, P, sort_mode); out.steal_mem(tmp); } else { - all_non_nan = op_sort_index::apply_noalias_proxy(out, P, sort_mode); + all_non_nan = op_sort_index::apply_helper(out, P, sort_mode); } - arma_conform_check( (all_non_nan == false), "sort_index(): detected NaN" ); - } - - - -// - - - -template -inline -bool -op_stable_sort_index::apply_noalias(Mat& out, const Proxy& P, const uword sort_mode) - { - arma_debug_sigprint(); - - return arma_sort_index_helper(out, P, sort_mode); - } - - - -template -inline -void -op_stable_sort_index::apply(Mat& out, const mtOp& in) - { - arma_debug_sigprint(); + if(all_non_nan == false) { out.soft_reset(); } - const Proxy P(in.m); - - if(P.get_n_elem() == 0) { out.set_size(0,1); return; } - - const uword sort_mode = in.aux_uword_a; - - bool all_non_nan = false; - - if(P.is_alias(out)) - { - Mat tmp; - - all_non_nan = op_stable_sort_index::apply_noalias(tmp, P, sort_mode); - - out.steal_mem(tmp); - } - else - { - all_non_nan = op_stable_sort_index::apply_noalias(out, P, sort_mode); - } - - arma_conform_check( (all_non_nan == false), "stable_sort_index(): detected NaN" ); + arma_conform_check( (all_non_nan == false), "sort_index(): detected NaN" ); } diff --git a/inst/include/armadillo_bits/op_sort_meat.hpp b/inst/include/armadillo_bits/op_sort_meat.hpp index df3c40a5..0305c67f 100644 --- a/inst/include/armadillo_bits/op_sort_meat.hpp +++ b/inst/include/armadillo_bits/op_sort_meat.hpp @@ -86,7 +86,18 @@ op_sort::apply_noalias(Mat& out, const Mat& X, const uword sort_mode, co { const Col X_col( const_cast(X.colptr(col)), n_rows, false ); - op_sort_index::apply_noalias_mat(indices, X_col, sort_mode); + const Proxy< Col > P(X_col); + + const bool all_non_nan = op_sort_index::apply_helper(indices, P, sort_mode); + + if(all_non_nan == false) + { + out.soft_reset(); + + arma_conform_check(true, "sort(): detected NaN"); + + return; + } const uword* indices_mem = indices.memptr(); const eT* X_col_mem = X_col.memptr(); @@ -115,7 +126,18 @@ op_sort::apply_noalias(Mat& out, const Mat& X, const uword sort_mode, co { const Col Y_col( const_cast(Y.colptr(col)), n_rows, false ); - op_sort_index::apply_noalias_mat(indices, Y_col, sort_mode); + const Proxy< Col > P(Y_col); + + const bool all_non_nan = op_sort_index::apply_helper(indices, P, sort_mode); + + if(all_non_nan == false) + { + out.soft_reset(); + + arma_conform_check(true, "sort(): detected NaN"); + + return; + } const uword* indices_mem = indices.memptr(); const eT* Y_col_mem = Y_col.memptr(); @@ -195,9 +217,8 @@ op_sort::apply(Mat& out, const Op& in) const uword sort_mode = in.aux_uword_a; const uword dim = in.aux_uword_b; - arma_conform_check( (sort_mode > 1), "sort(): parameter 'sort_mode' must be 0 or 1" ); - arma_conform_check( (dim > 1), "sort(): parameter 'dim' must be 0 or 1" ); - arma_conform_check( (X.internal_has_nan()), "sort(): detected NaN" ); + arma_conform_check( (sort_mode > 1), "sort(): parameter 'sort_mode' must be 0 or 1" ); + arma_conform_check( (dim > 1), "sort(): parameter 'dim' must be 0 or 1" ); if(U.is_alias(out)) { @@ -229,8 +250,7 @@ op_sort_vec::apply(Mat& out, const Op& i const uword sort_mode = in.aux_uword_a; - arma_conform_check( (sort_mode > 1), "sort(): parameter 'sort_mode' must be 0 or 1" ); - arma_conform_check( (X.internal_has_nan()), "sort(): detected NaN" ); + arma_conform_check( (sort_mode > 1), "sort(): parameter 'sort_mode' must be 0 or 1" ); if(X.n_elem <= 1) { out = X; return; } @@ -238,7 +258,18 @@ op_sort_vec::apply(Mat& out, const Op& i { uvec indices; - op_sort_index::apply_noalias_mat(indices, X, sort_mode); + const Proxy< Mat > P(X); + + const bool all_non_nan = op_sort_index::apply_helper(indices, P, sort_mode); + + if(all_non_nan == false) + { + out.soft_reset(); + + arma_conform_check(true, "sort(): detected NaN"); + + return; + } const uword N = indices.n_elem; diff --git a/inst/include/armadillo_bits/op_sp_nonzeros_meat.hpp b/inst/include/armadillo_bits/op_sp_nonzeros_meat.hpp index 33a8e1c3..01a9a553 100644 --- a/inst/include/armadillo_bits/op_sp_nonzeros_meat.hpp +++ b/inst/include/armadillo_bits/op_sp_nonzeros_meat.hpp @@ -31,31 +31,30 @@ op_sp_nonzeros::apply(Mat& out, const SpToDOp P(X.m); - - const uword N = P.get_n_nonzero(); - - out.set_size(N,1); - - if(N == 0) { return; } - - if(is_SpMat::stored_type>::value) + if(is_SpMat::value || is_SpMat::stored_type>::value) { - const unwrap_spmat::stored_type> U(P.Q); + const unwrap_spmat U(X.m); - arrayops::copy(out.memptr(), U.M.values, N); + out.set_size(U.M.n_nonzero,1); + + arrayops::copy(out.memptr(), U.M.values, U.M.n_nonzero); return; } - if(is_SpSubview::stored_type>::value) + if(is_SpSubview::value) { - const SpSubview& sv = reinterpret_cast< const SpSubview& >(P.Q); + const SpSubview& sv = reinterpret_cast< const SpSubview& >(X.m); if(sv.n_rows == sv.m.n_rows) { + arma_debug_print("op_sp_nonzeros::apply(): SpSubview optimisation"); + const SpMat& m = sv.m; const uword col = sv.aux_col1; + const uword N = sv.n_nonzero; + + out.set_size(N, 1); arrayops::copy(out.memptr(), &(m.values[ m.col_ptrs[col] ]), N); @@ -63,6 +62,14 @@ op_sp_nonzeros::apply(Mat& out, const SpToDOp P(X.m); + + const uword N = P.get_n_nonzero(); + + out.set_size(N,1); + + if(N == 0) { return; } + eT* out_mem = out.memptr(); typename SpProxy::const_iterator_type it = P.begin(); diff --git a/inst/include/armadillo_bits/op_stddev_bones.hpp b/inst/include/armadillo_bits/op_stddev_bones.hpp index fbef5d35..08370f38 100644 --- a/inst/include/armadillo_bits/op_stddev_bones.hpp +++ b/inst/include/armadillo_bits/op_stddev_bones.hpp @@ -32,16 +32,8 @@ class op_stddev template inline static void apply_noalias(Mat::result>& out, const Mat& X, const uword norm_type, const uword dim); - // - - template - inline static typename get_pod_type::result stddev_vec(const subview_col& X, const uword norm_type = 0); - - template - inline static typename get_pod_type::result stddev_vec(const subview_row& X, const uword norm_type = 0); - template - inline static typename T1::pod_type stddev_vec(const Base& X, const uword norm_type = 0); + inline static typename T1::pod_type stddev_vec(const T1& X, const uword norm_type = 0); }; diff --git a/inst/include/armadillo_bits/op_stddev_meat.hpp b/inst/include/armadillo_bits/op_stddev_meat.hpp index 9bf9e307..b54a3a4f 100644 --- a/inst/include/armadillo_bits/op_stddev_meat.hpp +++ b/inst/include/armadillo_bits/op_stddev_meat.hpp @@ -76,7 +76,7 @@ op_stddev::apply_noalias(Mat::result>& out, const M { out_eT* out_mem = out.memptr(); - for(uword col=0; col::result>& out, const M if(X_n_cols > 0) { - podarray dat(X_n_cols); - - in_eT* dat_mem = dat.memptr(); out_eT* out_mem = out.memptr(); - for(uword row=0; row tmp; + + for(uword row=0; row < X_n_rows; ++row) { - dat.copy_row(X, row); + tmp.copy_row(X, row); - out_mem[row] = std::sqrt( op_var::direct_var( dat_mem, X_n_cols, norm_type) ); + out_mem[row] = std::sqrt( op_var::direct_var( tmp.memptr(), tmp.n_elem, norm_type) ); } } } @@ -111,7 +110,7 @@ op_stddev::apply_noalias(Mat::result>& out, const M template inline typename T1::pod_type -op_stddev::stddev_vec(const Base& X, const uword norm_type) +op_stddev::stddev_vec(const T1& X, const uword norm_type) { arma_debug_sigprint(); @@ -119,7 +118,7 @@ op_stddev::stddev_vec(const Base& X, const uword nor arma_conform_check( (norm_type > 1), "stddev(): parameter 'norm_type' must be 0 or 1" ); - const quasi_unwrap U(X.get_ref()); + const quasi_unwrap U(X); if(U.M.n_elem == 0) { @@ -133,67 +132,4 @@ op_stddev::stddev_vec(const Base& X, const uword nor -template -inline -typename get_pod_type::result -op_stddev::stddev_vec(const subview_col& X, const uword norm_type) - { - arma_debug_sigprint(); - - typedef typename get_pod_type::result T; - - arma_conform_check( (norm_type > 1), "stddev(): parameter 'norm_type' must be 0 or 1" ); - - if(X.n_elem == 0) - { - arma_conform_check(true, "stddev(): object has no elements"); - - return Datum::nan; - } - - return std::sqrt( op_var::direct_var(X.colptr(0), X.n_rows, norm_type) ); - } - - - - -template -inline -typename get_pod_type::result -op_stddev::stddev_vec(const subview_row& X, const uword norm_type) - { - arma_debug_sigprint(); - - typedef typename get_pod_type::result T; - - arma_conform_check( (norm_type > 1), "stddev(): parameter 'norm_type' must be 0 or 1" ); - - if(X.n_elem == 0) - { - arma_conform_check(true, "stddev(): object has no elements"); - - return Datum::nan; - } - - const Mat& A = X.m; - - const uword start_row = X.aux_row1; - const uword start_col = X.aux_col1; - - const uword end_col_p1 = start_col + X.n_cols; - - podarray tmp(X.n_elem); - eT* tmp_mem = tmp.memptr(); - - for(uword i=0, col=start_col; col < end_col_p1; ++col, ++i) - { - tmp_mem[i] = A.at(start_row, col); - } - - return std::sqrt( op_var::direct_var(tmp.memptr(), tmp.n_elem, norm_type) ); - } - - - //! @} - diff --git a/inst/include/armadillo_bits/op_var_bones.hpp b/inst/include/armadillo_bits/op_var_bones.hpp index ee13bd48..690a71cb 100644 --- a/inst/include/armadillo_bits/op_var_bones.hpp +++ b/inst/include/armadillo_bits/op_var_bones.hpp @@ -34,23 +34,17 @@ class op_var // - template - inline static typename get_pod_type::result var_vec(const subview_col& X, const uword norm_type = 0); - - template - inline static typename get_pod_type::result var_vec(const subview_row& X, const uword norm_type = 0); - template - inline static typename T1::pod_type var_vec(const Base& X, const uword norm_type = 0); + inline static typename T1::pod_type var_vec(const T1& X, const uword norm_type = 0); // template - inline static eT direct_var(const eT* const X, const uword N, const uword norm_type = 0); + inline static eT direct_var(const eT* X, const uword N, const uword norm_type = 0); template - inline static eT direct_var_robust(const eT* const X, const uword N, const uword norm_type = 0); + inline static eT direct_var_robust(const eT* X, const uword N, const uword norm_type = 0); // diff --git a/inst/include/armadillo_bits/op_var_meat.hpp b/inst/include/armadillo_bits/op_var_meat.hpp index 64d82d04..bdf5cdfb 100644 --- a/inst/include/armadillo_bits/op_var_meat.hpp +++ b/inst/include/armadillo_bits/op_var_meat.hpp @@ -76,7 +76,7 @@ op_var::apply_noalias(Mat::result>& out, const Mat< { out_eT* out_mem = out.memptr(); - for(uword col=0; col::result>& out, const Mat< if(X_n_cols > 0) { - podarray dat(X_n_cols); - - in_eT* dat_mem = dat.memptr(); out_eT* out_mem = out.memptr(); - for(uword row=0; row tmp; + + for(uword row=0; row < X_n_rows; ++row) { - dat.copy_row(X, row); + tmp.copy_row(X, row); - out_mem[row] = op_var::direct_var( dat_mem, X_n_cols, norm_type ); + out_mem[row] = op_var::direct_var( tmp.memptr(), tmp.n_elem, norm_type ); } } } @@ -111,7 +110,7 @@ op_var::apply_noalias(Mat::result>& out, const Mat< template inline typename T1::pod_type -op_var::var_vec(const Base& X, const uword norm_type) +op_var::var_vec(const T1& X, const uword norm_type) { arma_debug_sigprint(); @@ -119,7 +118,7 @@ op_var::var_vec(const Base& X, const uword norm_type arma_conform_check( (norm_type > 1), "var(): parameter 'norm_type' must be 0 or 1" ); - const quasi_unwrap U(X.get_ref()); + const quasi_unwrap U(X); if(U.M.n_elem == 0) { @@ -133,73 +132,11 @@ op_var::var_vec(const Base& X, const uword norm_type -template -inline -typename get_pod_type::result -op_var::var_vec(const subview_col& X, const uword norm_type) - { - arma_debug_sigprint(); - - typedef typename get_pod_type::result T; - - arma_conform_check( (norm_type > 1), "var(): parameter 'norm_type' must be 0 or 1" ); - - if(X.n_elem == 0) - { - arma_conform_check(true, "var(): object has no elements"); - - return Datum::nan; - } - - return op_var::direct_var(X.colptr(0), X.n_rows, norm_type); - } - - - - -template -inline -typename get_pod_type::result -op_var::var_vec(const subview_row& X, const uword norm_type) - { - arma_debug_sigprint(); - - typedef typename get_pod_type::result T; - - arma_conform_check( (norm_type > 1), "var(): parameter 'norm_type' must be 0 or 1" ); - - if(X.n_elem == 0) - { - arma_conform_check(true, "var(): object has no elements"); - - return Datum::nan; - } - - const Mat& A = X.m; - - const uword start_row = X.aux_row1; - const uword start_col = X.aux_col1; - - const uword end_col_p1 = start_col + X.n_cols; - - podarray tmp(X.n_elem); - eT* tmp_mem = tmp.memptr(); - - for(uword i=0, col=start_col; col < end_col_p1; ++col, ++i) - { - tmp_mem[i] = A.at(start_row, col); - } - - return op_var::direct_var(tmp.memptr(), tmp.n_elem, norm_type); - } - - - //! find the variance of an array template inline eT -op_var::direct_var(const eT* const X, const uword n_elem, const uword norm_type) +op_var::direct_var(const eT* X, const uword n_elem, const uword norm_type) { arma_debug_sigprint(); @@ -207,6 +144,8 @@ op_var::direct_var(const eT* const X, const uword n_elem, const uword norm_type) { const eT acc1 = op_mean::direct_mean(X, n_elem); + if(arma_isnonfinite(acc1)) { return Datum::nan; } + eT acc2 = eT(0); eT acc3 = eT(0); @@ -251,7 +190,7 @@ op_var::direct_var(const eT* const X, const uword n_elem, const uword norm_type) template inline eT -op_var::direct_var_robust(const eT* const X, const uword n_elem, const uword norm_type) +op_var::direct_var_robust(const eT* X, const uword n_elem, const uword norm_type) { arma_debug_sigprint(); @@ -294,6 +233,8 @@ op_var::direct_var(const std::complex* const X, const uword n_elem, const uwo { const eT acc1 = op_mean::direct_mean(X, n_elem); + if(arma_isnonfinite(acc1)) { return Datum::nan; } + T acc2 = T(0); eT acc3 = eT(0); @@ -354,4 +295,3 @@ op_var::direct_var_robust(const std::complex* const X, const uword n_elem, co //! @} - diff --git a/inst/include/armadillo_bits/podarray_meat.hpp b/inst/include/armadillo_bits/podarray_meat.hpp index 17be3216..09a2bef6 100644 --- a/inst/include/armadillo_bits/podarray_meat.hpp +++ b/inst/include/armadillo_bits/podarray_meat.hpp @@ -37,8 +37,8 @@ podarray::~podarray() template inline podarray::podarray() - : n_elem(0) - , mem (0) + : n_elem(0 ) + , mem (nullptr) { arma_debug_sigprint_this(this); } @@ -128,7 +128,7 @@ arma_inline eT& podarray::operator[] (const uword i) { - return access::rw(mem[i]); + return mem[i]; } @@ -152,7 +152,7 @@ podarray::operator() (const uword i) { arma_conform_check_bounds( (i >= n_elem), "podarray::operator(): index out of bounds" ); - return access::rw(mem[i]); + return mem[i]; } @@ -258,11 +258,11 @@ podarray::copy_row(const Mat& A, const uword row) { arma_debug_sigprint(); - // note: this function assumes that the podarray has been set to the correct size beforehand - const uword n_rows = A.n_rows; const uword n_cols = A.n_cols; + init_warm(n_cols); + const eT* A_mem = &(A.at(row,0)); eT* out_mem = memptr(); diff --git a/inst/include/armadillo_bits/running_stat_meat.hpp b/inst/include/armadillo_bits/running_stat_meat.hpp index 00f73deb..9c3997f9 100644 --- a/inst/include/armadillo_bits/running_stat_meat.hpp +++ b/inst/include/armadillo_bits/running_stat_meat.hpp @@ -162,9 +162,9 @@ running_stat::operator() (const typename running_stat::T sample) { arma_debug_sigprint(); - if( arma_isfinite(sample) == false ) + if(arma_isnonfinite(sample)) { - arma_warn(3, "running_stat: sample ignored as it is non-finite" ); + arma_warn(3, "running_stat: non-finite sample ignored" ); return; } @@ -181,9 +181,9 @@ running_stat::operator() (const std::complex< typename running_stat::T > { arma_debug_sigprint(); - if( arma_isfinite(sample) == false ) + if(arma_isnonfinite(sample)) { - arma_warn(3, "running_stat: sample ignored as it is non-finite" ); + arma_warn(3, "running_stat: non-finite sample ignored" ); return; } diff --git a/inst/include/armadillo_bits/running_stat_vec_meat.hpp b/inst/include/armadillo_bits/running_stat_vec_meat.hpp index d949ac81..b40272e7 100644 --- a/inst/include/armadillo_bits/running_stat_vec_meat.hpp +++ b/inst/include/armadillo_bits/running_stat_vec_meat.hpp @@ -100,7 +100,7 @@ running_stat_vec::operator() (const Base::operator() (const Base< std::complex + inline static void apply(SpMat& out, const mtSpOp& in); + }; + + + //! @} diff --git a/inst/include/armadillo_bits/spop_misc_meat.hpp b/inst/include/armadillo_bits/spop_misc_meat.hpp index 5da73142..e1b9c4d4 100644 --- a/inst/include/armadillo_bits/spop_misc_meat.hpp +++ b/inst/include/armadillo_bits/spop_misc_meat.hpp @@ -547,4 +547,21 @@ spop_fliplr::apply(SpMat& out, const SpOp +inline +void +spop_replace::apply(SpMat& out, const mtSpOp& in) + { + arma_debug_sigprint(); + + const eT old_val = in.aux; + const eT new_val = in.aux_out_eT; + + out = in.m; + + out.replace(old_val, new_val); + } + + + //! @} diff --git a/inst/include/armadillo_bits/spop_omit_bones.hpp b/inst/include/armadillo_bits/spop_omit_bones.hpp new file mode 100644 index 00000000..f798273c --- /dev/null +++ b/inst/include/armadillo_bits/spop_omit_bones.hpp @@ -0,0 +1,35 @@ +// SPDX-License-Identifier: Apache-2.0 +// +// Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au) +// Copyright 2008-2016 National ICT Australia (NICTA) +// +// 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. +// ------------------------------------------------------------------------ + + + +//! \addtogroup spop_omit +//! @{ + + +class spop_omit + : public traits_op_col + { + public: + + template inline static void apply(SpMat& out, const SpOp& in); + + template inline static void apply_noalias(SpMat& out, const SpProxy& P, functor is_omitted); + }; + + +//! @} diff --git a/inst/include/armadillo_bits/spop_omit_meat.hpp b/inst/include/armadillo_bits/spop_omit_meat.hpp new file mode 100644 index 00000000..e19c6abd --- /dev/null +++ b/inst/include/armadillo_bits/spop_omit_meat.hpp @@ -0,0 +1,119 @@ +// SPDX-License-Identifier: Apache-2.0 +// +// Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au) +// Copyright 2008-2016 National ICT Australia (NICTA) +// +// 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. +// ------------------------------------------------------------------------ + + + +//! \addtogroup spop_omit +//! @{ + + + +template +inline +void +spop_omit::apply(SpMat& out, const SpOp& in) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + const uword omit_mode = in.aux_uword_a; + + if(arma_config::fast_math_warn) + { + if(omit_mode == 1) { arma_warn(1, "omit_nan(): detection of NaN is not reliable in fast math mode"); } + if(omit_mode == 2) { arma_warn(1, "omit_nonfinite(): detection of non-finite values is not reliable in fast math mode"); } + } + + auto is_omitted_1 = [](const eT& x) -> bool { return arma_isnan(x); }; + auto is_omitted_2 = [](const eT& x) -> bool { return arma_isnonfinite(x); }; + + const SpProxy P(in.m); + + if(P.is_alias(out)) + { + SpMat tmp; + + if(omit_mode == 1) { spop_omit::apply_noalias(tmp, P, is_omitted_1); } + if(omit_mode == 2) { spop_omit::apply_noalias(tmp, P, is_omitted_2); } + + out.steal_mem(tmp); + } + else + { + if(omit_mode == 1) { spop_omit::apply_noalias(out, P, is_omitted_1); } + if(omit_mode == 2) { spop_omit::apply_noalias(out, P, is_omitted_2); } + } + } + + + +template +inline +void +spop_omit::apply_noalias(SpMat& out, const SpProxy& P, functor is_omitted) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + const uword n_rows = P.get_n_rows(); + const uword max_n_nonzero = P.get_n_nonzero(); + + if(max_n_nonzero == 0) { out.reset(); return; } + + out.reserve(P.get_n_elem(), 1, max_n_nonzero); + + typename SpProxy::const_iterator_type it = P.begin(); + typename SpProxy::const_iterator_type it_end = P.end(); + + uword count = 0; + + for(; it != it_end; ++it) + { + const eT val = (*it); + + if(is_omitted(val) == false) + { + const uword index = it.row() + it.col()*n_rows; + + access::rw(out.values[count]) = val; + access::rw(out.row_indices[count]) = index; + access::rw(out.col_ptrs[1])++; + ++count; + } + } + + if(count < max_n_nonzero) + { + if(count <= (max_n_nonzero/2)) + { + out.mem_resize(count); + } + else + { + // quick resize without reallocating memory and copying data + access::rw( out.n_nonzero) = count; + access::rw( out.values[count]) = eT(0); + access::rw(out.row_indices[count]) = uword(0); + } + } + } + + + +//! @} diff --git a/inst/include/armadillo_bits/subview_elem1_meat.hpp b/inst/include/armadillo_bits/subview_elem1_meat.hpp index 3f62f104..ee3ebf52 100644 --- a/inst/include/armadillo_bits/subview_elem1_meat.hpp +++ b/inst/include/armadillo_bits/subview_elem1_meat.hpp @@ -474,8 +474,8 @@ subview_elem1::randu() eT* m_mem = m_local.memptr(); const uword m_n_elem = m_local.n_elem; - const unwrap_check_mixed tmp(a.get_ref(), m_local); - const umat& aa = tmp.M; + const unwrap_check_mixed U(a.get_ref(), m_local); + const umat& aa = U.M; if(resolves_to_vector::no) { @@ -485,28 +485,19 @@ subview_elem1::randu() const uword* aa_mem = aa.memptr(); const uword aa_n_elem = aa.n_elem; - uword iq,jq; - for(iq=0, jq=1; jq < aa_n_elem; iq+=2, jq+=2) - { - const uword ii = aa_mem[iq]; - const uword jj = aa_mem[jq]; - - arma_conform_check_bounds( ( (ii >= m_n_elem) || (jj >= m_n_elem) ), "Mat::elem(): index out of bounds" ); - - const eT val1 = eT(arma_rng::randu()); - const eT val2 = eT(arma_rng::randu()); - - m_mem[ii] = val1; - m_mem[jj] = val2; - } + podarray tmp(aa_n_elem); - if(iq < aa_n_elem) + eT* tmp_mem = tmp.memptr(); + + arma_rng::randu::fill(tmp_mem, aa_n_elem); + + for(uword iq=0; iq < aa_n_elem; ++iq) { const uword ii = aa_mem[iq]; - arma_conform_check_bounds( (ii >= m_n_elem) , "Mat::elem(): index out of bounds" ); + arma_conform_check_bounds( (ii >= m_n_elem), "Mat::elem(): index out of bounds" ); - m_mem[ii] = eT(arma_rng::randu()); + m_mem[ii] = tmp_mem[iq]; } } @@ -524,8 +515,8 @@ subview_elem1::randn() eT* m_mem = m_local.memptr(); const uword m_n_elem = m_local.n_elem; - const unwrap_check_mixed tmp(a.get_ref(), m_local); - const umat& aa = tmp.M; + const unwrap_check_mixed U(a.get_ref(), m_local); + const umat& aa = U.M; if(resolves_to_vector::no) { @@ -535,24 +526,19 @@ subview_elem1::randn() const uword* aa_mem = aa.memptr(); const uword aa_n_elem = aa.n_elem; - uword iq,jq; - for(iq=0, jq=1; jq < aa_n_elem; iq+=2, jq+=2) - { - const uword ii = aa_mem[iq]; - const uword jj = aa_mem[jq]; - - arma_conform_check_bounds( ( (ii >= m_n_elem) || (jj >= m_n_elem) ), "Mat::elem(): index out of bounds" ); - - arma_rng::randn::dual_val( m_mem[ii], m_mem[jj] ); - } + podarray tmp(aa_n_elem); - if(iq < aa_n_elem) + eT* tmp_mem = tmp.memptr(); + + arma_rng::randn::fill(tmp_mem, aa_n_elem); + + for(uword iq=0; iq < aa_n_elem; ++iq) { const uword ii = aa_mem[iq]; - arma_conform_check_bounds( (ii >= m_n_elem) , "Mat::elem(): index out of bounds" ); + arma_conform_check_bounds( (ii >= m_n_elem), "Mat::elem(): index out of bounds" ); - m_mem[ii] = eT(arma_rng::randn()); + m_mem[ii] = tmp_mem[iq]; } } diff --git a/inst/include/armadillo_bits/subview_elem2_bones.hpp b/inst/include/armadillo_bits/subview_elem2_bones.hpp index bec561d7..4c5fe55f 100644 --- a/inst/include/armadillo_bits/subview_elem2_bones.hpp +++ b/inst/include/armadillo_bits/subview_elem2_bones.hpp @@ -67,6 +67,8 @@ class subview_elem2 : public Base< eT, subview_elem2 > inline void fill(const eT val); inline void zeros(); inline void ones(); + inline void randu(); + inline void randn(); inline void operator+= (const eT val); inline void operator-= (const eT val); diff --git a/inst/include/armadillo_bits/subview_elem2_meat.hpp b/inst/include/armadillo_bits/subview_elem2_meat.hpp index 4a160642..6d25c7b8 100644 --- a/inst/include/armadillo_bits/subview_elem2_meat.hpp +++ b/inst/include/armadillo_bits/subview_elem2_meat.hpp @@ -299,6 +299,240 @@ subview_elem2::inplace_op(const Base& x) +template +inline +void +subview_elem2::randu() + { + arma_debug_sigprint(); + + Mat& m_local = const_cast< Mat& >(m); + + const uword m_n_rows = m_local.n_rows; + const uword m_n_cols = m_local.n_cols; + + if( (all_rows == false) && (all_cols == false) ) + { + const unwrap_check_mixed U1(base_ri.get_ref(), m_local); + const unwrap_check_mixed U2(base_ci.get_ref(), m_local); + + const umat& ri = U1.M; + const umat& ci = U2.M; + + arma_conform_check + ( + ( ((ri.is_vec() == false) && (ri.is_empty() == false)) || ((ci.is_vec() == false) && (ci.is_empty() == false)) ), + "Mat::elem(): given object must be a vector" + ); + + const uword* ri_mem = ri.memptr(); + const uword ri_n_elem = ri.n_elem; + + const uword* ci_mem = ci.memptr(); + const uword ci_n_elem = ci.n_elem; + + podarray tmp(ri_n_elem); + + eT* tmp_mem = tmp.memptr(); + + for(uword ci_count=0; ci_count < ci_n_elem; ++ci_count) + { + const uword col = ci_mem[ci_count]; + + arma_conform_check_bounds( (col >= m_n_cols), "Mat::elem(): index out of bounds" ); + + arma_rng::randu::fill(tmp_mem, ri_n_elem); + + for(uword ri_count=0; ri_count < ri_n_elem; ++ri_count) + { + const uword row = ri_mem[ri_count]; + + arma_conform_check_bounds( (row >= m_n_rows), "Mat::elem(): index out of bounds" ); + + m_local.at(row,col) = tmp_mem[ri_count]; + } + } + } + else + if( (all_rows == true) && (all_cols == false) ) + { + const unwrap_check_mixed U2(base_ci.get_ref(), m_local); + + const umat& ci = U2.M; + + arma_conform_check + ( + ( (ci.is_vec() == false) && (ci.is_empty() == false) ), + "Mat::elem(): given object must be a vector" + ); + + const uword* ci_mem = ci.memptr(); + const uword ci_n_elem = ci.n_elem; + + for(uword ci_count=0; ci_count < ci_n_elem; ++ci_count) + { + const uword col = ci_mem[ci_count]; + + arma_conform_check_bounds( (col >= m_n_cols), "Mat::elem(): index out of bounds" ); + + arma_rng::randu::fill(m_local.colptr(col), m_n_rows); + } + } + else + if( (all_rows == false) && (all_cols == true) ) + { + const unwrap_check_mixed U1(base_ri.get_ref(), m_local); + + const umat& ri = U1.M; + + arma_conform_check + ( + ( (ri.is_vec() == false) && (ri.is_empty() == false) ), + "Mat::elem(): given object must be a vector" + ); + + const uword* ri_mem = ri.memptr(); + const uword ri_n_elem = ri.n_elem; + + podarray tmp(ri_n_elem); + + eT* tmp_mem = tmp.memptr(); + + for(uword col=0; col < m_n_cols; ++col) + { + arma_rng::randu::fill(tmp_mem, ri_n_elem); + + for(uword ri_count=0; ri_count < ri_n_elem; ++ri_count) + { + const uword row = ri_mem[ri_count]; + + arma_conform_check_bounds( (row >= m_n_rows), "Mat::elem(): index out of bounds" ); + + m_local.at(row,col) = tmp_mem[ri_count]; + } + } + } + } + + + +template +inline +void +subview_elem2::randn() + { + arma_debug_sigprint(); + + Mat& m_local = const_cast< Mat& >(m); + + const uword m_n_rows = m_local.n_rows; + const uword m_n_cols = m_local.n_cols; + + if( (all_rows == false) && (all_cols == false) ) + { + const unwrap_check_mixed U1(base_ri.get_ref(), m_local); + const unwrap_check_mixed U2(base_ci.get_ref(), m_local); + + const umat& ri = U1.M; + const umat& ci = U2.M; + + arma_conform_check + ( + ( ((ri.is_vec() == false) && (ri.is_empty() == false)) || ((ci.is_vec() == false) && (ci.is_empty() == false)) ), + "Mat::elem(): given object must be a vector" + ); + + const uword* ri_mem = ri.memptr(); + const uword ri_n_elem = ri.n_elem; + + const uword* ci_mem = ci.memptr(); + const uword ci_n_elem = ci.n_elem; + + podarray tmp(ri_n_elem); + + eT* tmp_mem = tmp.memptr(); + + for(uword ci_count=0; ci_count < ci_n_elem; ++ci_count) + { + const uword col = ci_mem[ci_count]; + + arma_conform_check_bounds( (col >= m_n_cols), "Mat::elem(): index out of bounds" ); + + arma_rng::randn::fill(tmp_mem, ri_n_elem); + + for(uword ri_count=0; ri_count < ri_n_elem; ++ri_count) + { + const uword row = ri_mem[ri_count]; + + arma_conform_check_bounds( (row >= m_n_rows), "Mat::elem(): index out of bounds" ); + + m_local.at(row,col) = tmp_mem[ri_count]; + } + } + } + else + if( (all_rows == true) && (all_cols == false) ) + { + const unwrap_check_mixed U2(base_ci.get_ref(), m_local); + + const umat& ci = U2.M; + + arma_conform_check + ( + ( (ci.is_vec() == false) && (ci.is_empty() == false) ), + "Mat::elem(): given object must be a vector" + ); + + const uword* ci_mem = ci.memptr(); + const uword ci_n_elem = ci.n_elem; + + for(uword ci_count=0; ci_count < ci_n_elem; ++ci_count) + { + const uword col = ci_mem[ci_count]; + + arma_conform_check_bounds( (col >= m_n_cols), "Mat::elem(): index out of bounds" ); + + arma_rng::randn::fill(m_local.colptr(col), m_n_rows); + } + } + else + if( (all_rows == false) && (all_cols == true) ) + { + const unwrap_check_mixed U1(base_ri.get_ref(), m_local); + + const umat& ri = U1.M; + + arma_conform_check + ( + ( (ri.is_vec() == false) && (ri.is_empty() == false) ), + "Mat::elem(): given object must be a vector" + ); + + const uword* ri_mem = ri.memptr(); + const uword ri_n_elem = ri.n_elem; + + podarray tmp(ri_n_elem); + + eT* tmp_mem = tmp.memptr(); + + for(uword col=0; col < m_n_cols; ++col) + { + arma_rng::randn::fill(tmp_mem, ri_n_elem); + + for(uword ri_count=0; ri_count < ri_n_elem; ++ri_count) + { + const uword row = ri_mem[ri_count]; + + arma_conform_check_bounds( (row >= m_n_rows), "Mat::elem(): index out of bounds" ); + + m_local.at(row,col) = tmp_mem[ri_count]; + } + } + } + } + + + // // diff --git a/inst/include/armadillo_bits/sym_helper.hpp b/inst/include/armadillo_bits/sym_helper.hpp index 2134885d..6f27530a 100644 --- a/inst/include/armadillo_bits/sym_helper.hpp +++ b/inst/include/armadillo_bits/sym_helper.hpp @@ -58,8 +58,8 @@ guess_sympd_worker(const Mat& A) { const eT A_jj = A_col[j]; - if( A_jj <= eT(0)) { return false; } - if(arma_isfinite(A_jj) == false) { return false; } + if( A_jj <= eT(0)) { return false; } + if(arma_isnonfinite(A_jj)) { return false; } if(A_jj >= tol) { diag_below_tol = false; } @@ -147,8 +147,8 @@ guess_sympd_worker(const Mat& A) const T A_jj_rabs = std::abs(A_jj_r); const T A_jj_iabs = std::abs(A_jj_i); - if( A_jj_r <= T(0) ) { return false; } // real should be positive - if(arma_isfinite(A_jj_r) == false) { return false; } + if( A_jj_r <= T(0) ) { return false; } // real should be positive + if(arma_isnonfinite(A_jj_r)) { return false; } if(A_jj_iabs > tol ) { return false; } // imag should be approx zero if(A_jj_iabs > A_jj_rabs) { return false; } // corner case: real and imag are close to zero, and imag is dominant @@ -164,7 +164,7 @@ guess_sympd_worker(const Mat& A) const T square_max_diag = max_diag * max_diag; - if(arma_isfinite(square_max_diag) == false) { return false; } + if(arma_isnonfinite(square_max_diag)) { return false; } A_col = A_mem; @@ -188,7 +188,7 @@ guess_sympd_worker(const Mat& A) // avoid using std::abs(), as that is time consuming due to division and std::sqrt() const T square_A_ij_abs = (A_ij_real * A_ij_real) + (A_ij_imag * A_ij_imag); - if(arma_isfinite(square_A_ij_abs) == false) { return false; } + if(arma_isnonfinite(square_A_ij_abs)) { return false; } if(square_A_ij_abs >= square_max_diag) { return false; } @@ -285,7 +285,7 @@ is_approx_sym_worker(const Mat& A) { const eT A_jj = A_col[j]; - if(arma_isfinite(A_jj) == false) { return false; } + if(arma_isnonfinite(A_jj)) { return false; } if(std::abs(A_jj) >= tol) { diag_below_tol = false; } @@ -359,7 +359,7 @@ is_approx_sym_worker(const Mat& A) if(A_jj_iabs > tol ) { return false; } // imag should be approx zero if(A_jj_iabs > A_jj_rabs) { return false; } // corner case: real and imag are close to zero, and imag is dominant - if(arma_isfinite(A_jj_r) == false) { return false; } + if(arma_isnonfinite(A_jj_r)) { return false; } if(A_jj_rabs >= tol) { diag_below_tol = false; } diff --git a/inst/include/armadillo_bits/translate_lapack.hpp b/inst/include/armadillo_bits/translate_lapack.hpp index ddf5e0a1..afec7fed 100644 --- a/inst/include/armadillo_bits/translate_lapack.hpp +++ b/inst/include/armadillo_bits/translate_lapack.hpp @@ -1485,6 +1485,28 @@ namespace lapack #endif } + + + template + inline + void + gebal(const char* job, const blas_int* n, eT* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, typename get_pod_type::result* scale, blas_int* info) + { + arma_type_check(( is_supported_blas_type::value == false )); + + #if defined(ARMA_USE_FORTRAN_HIDDEN_ARGS) + if( is_float::value) { typedef float pod_T; typedef float T; arma_fortran(arma_sgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info, 1); } + else if( is_double::value) { typedef double pod_T; typedef double T; arma_fortran(arma_dgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info, 1); } + else if( is_cx_float::value) { typedef float pod_T; typedef blas_cxf T; arma_fortran(arma_cgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info, 1); } + else if(is_cx_double::value) { typedef double pod_T; typedef blas_cxd T; arma_fortran(arma_zgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info, 1); } + #else + if( is_float::value) { typedef float pod_T; typedef float T; arma_fortran(arma_sgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info); } + else if( is_double::value) { typedef double pod_T; typedef double T; arma_fortran(arma_dgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info); } + else if( is_cx_float::value) { typedef float pod_T; typedef blas_cxf T; arma_fortran(arma_cgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info); } + else if(is_cx_double::value) { typedef double pod_T; typedef blas_cxd T; arma_fortran(arma_zgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info); } + #endif + } + } diff --git a/man/fastLm.Rd b/man/fastLm.Rd index 0f807557..5ad1d256 100644 --- a/man/fastLm.Rd +++ b/man/fastLm.Rd @@ -38,13 +38,14 @@ fastLm(X, \dots) \code{\link{lm}} and friends versus this approach based on \code{Armadillo}. The reason that \code{Armadillo} can do something like \code{\link{lm.fit}} faster than the functions in the stats - package is because \code{Armadillo} uses the Lapack version of the QR - decomposition while the stats package uses a \emph{modified} Linpack - version. Hence \code{Armadillo} uses level-3 BLAS code whereas the - stats package uses level-1 BLAS. However, \code{Armadillo} will + package is because \code{Armadillo} can use different solvers, including + fast / approximative ones. Older versions of Armadillo could therefore either fail or, worse, produce completely incorrect answers on rank-deficient model matrices whereas the functions from the stats package will handle them properly due to the modified Linpack code. + Newer Armadillo version pivot (with warning) to an approximate solutions. + This behavior can be controlled with options to the \code{solve} function, + see the Armadillo documentation. An example of the type of situation requiring extra care in checking for rank deficiency is a two-way layout with missing cells (see the @@ -60,7 +61,7 @@ fastLm(X, \dots) \code{fastLm} returns a richer object which also includes the residuals, fitted values and call argument similar to the \code{\link{lm}} or - \code{\link[MASS]{rlm}} functions.. + \code{\link[MASS]{rlm}} functions. } \seealso{\code{\link{lm}}, \code{\link{lm.fit}}} \references{Armadillo project: \url{https://arma.sourceforge.net/}} @@ -93,7 +94,7 @@ fastLm(X, \dots) set.seed(1) dd$y <- mm \%*\% seq_len(ncol(mm)) + rnorm(nrow(mm), sd = 0.1) summary(lm(y ~ f1 * f2, dd)) # detects rank deficiency - summary(fastLm(y ~ f1 * f2, dd)) # some huge coefficients + summary(fastLm(y ~ f1 * f2, dd)) # fits all coefficients via approx solution \dontshow{armadillo_reset_cores()} }