diff --git a/CHANGELOG.md b/CHANGELOG.md index 588ffde9c..095ee35bb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,32 @@ Below, you'll find a list of notable changes for each version of the ViennaRNA P ## Version 2.4.x -### [Unreleased](https://github.com/ViennaRNA/ViennaRNA/compare/v2.4.10...HEAD) +### [Unreleased](https://github.com/ViennaRNA/ViennaRNA/compare/v2.4.11...HEAD) + +### [v2.4.11](https://github.com/ViennaRNA/ViennaRNA/compare/v2.4.10...v2.4.11) (2018-12-17) + +#### Programs + * Add `--commands` option to `RNAsubopt` + * Add non-redundant Boltzmann sampling mode for `RNAsubopt` + +#### Library + * API: Fix wrong access to base pair soft constraints in equilibrium probability computations + * API: Fix behavior of `vrna_nucleotide_encode()` with lowercase characters in sequence + * API: Fix behavior of `encode_char()` with lowercase characters in sequence + * API: Fix forbidden `GU` pairs behavior in pscore computation for comparative folding + * API: Fix potential errors due to uninitialized `next` pointers in `vrna_move_t` of `vrna_eval_move_shift_pt` + * API: Add `AVX 512` optimized version of MFE multibranch loop decomposition + * API: Add functions for CPU SIMD feature detection + * API: Add dispatcher to automatically delegate exterior-/multibranch loop MFE decomposition to supported SIMD optimized implementation + * API: Add function `vrna_dist_mountain()` to compute mountain distance between two structures + * API: Add function `vrna_ensemble_defect()` to compute ensemble defect given a target structure + * API: Add non-redundant Boltzmann sampling + * API: Change behavior of `vrna_cstr_free()` and `vrna_cstr_close()` to always flush output before unregistering the stream + * SWIG: Add interface for `vrna_loopidx_from_ptable()` + +#### Package + * Activate compilation for compile-time supported SIMD optimized implementations by default + * Replace `--enable-sse` configure script option with `--disable-simd` ### [v2.4.10](https://github.com/ViennaRNA/ViennaRNA/compare/v2.4.9...v2.4.10) (2018-09-26) diff --git a/README.md b/README.md index bf28b944d..e37af1bd5 100644 --- a/README.md +++ b/README.md @@ -70,8 +70,8 @@ and installation process.* Usually you'll simply unpack the distribution tarball, configure and make: ``` -tar -zxvf ViennaRNA-2.4.10.tar.gz -cd ViennaRNA-2.4.10 +tar -zxvf ViennaRNA-2.4.11.tar.gz +cd ViennaRNA-2.4.11 ./configure make sudo make install @@ -102,7 +102,7 @@ script: program `RNALfold`: ``` cd src -tar -xzf libsvm-3.22.tar.gz +tar -xzf libsvm-3.23.tar.gz cd .. ``` @@ -184,7 +184,7 @@ will turn-off the `Perl 5` and `Python 2` interfaces. Disabling the entire scripting language support alltogether can be accomplished with the `--without-swig` switch. -#### Streaming SIMD Extension (SSE) support +#### Streaming SIMD Extension support Our latest version contains code that implements a faster multibranch loop decomposition in global MFE predictions, as used e.g. in `RNAfold`. This implementation makes use of modern processors capability to execute particular @@ -192,11 +192,13 @@ instructions on multiple data simultaneously (SIMD - single instruction multiple data, thanks to W. B. Langdon for providing the modified code). Consequently, the time required to assess the minimum of all multibranch loop decompositions is reduced up to about one half compared to the runtime of the original -implementation. To make use of this piece of code you need a CPU capable to -handle SSE4.1 instructions and enable the feature at compile-time using the -following configure flag: +implementation. This feature is enabled by default since version 2.4.11 and a +dispatcher ensures that the correct implementation will be selected at runtime. +If for any reason you want to disable this feature at compile-time use the +following configure flag + ``` -./configure --enable-sse +./configure --disable-simd ``` #### Link Time Optimization (LTO) diff --git a/RNAlib2.pc.in b/RNAlib2.pc.in index 2119ac2c1..1fd8197bf 100644 --- a/RNAlib2.pc.in +++ b/RNAlib2.pc.in @@ -7,6 +7,6 @@ Name: @PACKAGE_NAME@ Description: ViennaRNA Package 2 - Core library. Version: @PACKAGE_VERSION@ URL: @PACKAGE_URL@ -Libs: -L${libdir} -lRNA @LIBGOMPFLAG@ @GSL_LIBS@ @PTHREAD_LIBS@ +Libs: -L${libdir} -lRNA @LIBGOMPFLAG@ @GSL_LIBS@ @PTHREAD_LIBS@ @MPFR_LIBS@ Libs.private: -lm Cflags: -I${includedir} -I${includedir}/ViennaRNA @FLOAT_PF_FLAG@ @DEPRECATION_WARNING@ @DISABLE_C11_FEATURES@ @PTHREAD_CFLAGS@ diff --git a/THANKS b/THANKS index 9c7190401..f3c27b6e8 100644 --- a/THANKS +++ b/THANKS @@ -2,12 +2,13 @@ Over the past decades since the ViennaRNA Package first sprang to life of Ivo Hofackers PhD project, several different authors contributed more and more algorithm implementations. In 2008, Ronny Lorenz took over the extensive task to harmonize and simplify the already existing implementations for the sake of easier -feature addition. This eventually lead to version 2.0 of the ViennaRNA PAckage. +feature addition. This eventually lead to version 2.0 of the ViennaRNA Package. Since then, he (re-)implemented a large portion of the currently existing library features, such as the new, generalized constraints framework, RNA folding grammar domain extensions, and the major part of the scripting language interface. Below is a list of most people -who contributed more or less large portions of the implementations: +who contributed larger parts of the implementations: +- Juraj Michalik (non-redundant Boltzmann sampling) - Gregor Entzian (neighbor, walk) - Mario Koestl (worked on SWIG interface and related unit testing) - Dominik Luntzer (pertubation fold) diff --git a/configure.ac b/configure.ac index 51595a4f3..8fc0c14f6 100644 --- a/configure.ac +++ b/configure.ac @@ -1,7 +1,7 @@ dnl Process this file with autoconf to produce a configure script. m4_define([vrna_version_major], [2]) m4_define([vrna_version_minor], [4]) -m4_define([vrna_version_patch], [10]) +m4_define([vrna_version_patch], [11]) dnl Every other copy of the package version number gets its value from here AC_INIT([ViennaRNA],vrna_version_major.vrna_version_minor.vrna_version_patch,[rna@tbi.univie.ac.at],[ViennaRNA], [http://www.tbi.univie.ac.at/RNA]) diff --git a/doc/refman.include/install.dox b/doc/refman.include/install.dox index 8a59af5a6..1f7605ebe 100644 --- a/doc/refman.include/install.dox +++ b/doc/refman.include/install.dox @@ -17,8 +17,8 @@ as well as for Windows users are available from Download section of the @@ -121,6 +126,25 @@ extern char *pbacktrack(char *sequence); { return vrna_pbacktrack5($self, length); } + + std::vector + pbacktrack_nr(int num_samples) + { + std::vector str_vec; + + char **ptr, **output = vrna_pbacktrack_nr($self, num_samples); + + if (output) { + for (ptr = output; *ptr != NULL; ptr++) { + str_vec.push_back(std::string(*ptr)); + free(*ptr); + } + + free(output); + } + + return str_vec; + } } %include diff --git a/interfaces/structure_utils.i b/interfaces/structure_utils.i index 18427c17a..387206715 100644 --- a/interfaces/structure_utils.i +++ b/interfaces/structure_utils.i @@ -198,6 +198,37 @@ std::vector my_ptable_pk(std::string str); char *my_db_from_ptable(std::vector pt); +%rename (loopidx_from_ptable) my_loopidx_from_ptable; + +%{ + + std::vector + my_loopidx_from_ptable(std::vector pt) + { + int i, *idx; + std::vector vc; + std::vector v_idx; + + transform(pt.begin(), pt.end(), back_inserter(vc), convert_vecint2vecshort); + + idx = vrna_loopidx_from_ptable((short *)&vc[0]); + + v_idx.assign(idx, idx + pt.size()); + + free(idx); + + return v_idx; + } +%} + +#ifdef SWIGPYTHON +%feature("autodoc") my_loopidx_from_ptable; +%feature("kwargs") my_loopidx_from_ptable; +#endif + +std::vector my_loopidx_from_ptable(std::vector pt); + + %rename (bp_distance) my_bp_distance; %{ int my_bp_distance(const char *str1, const char *str2){ @@ -212,6 +243,24 @@ char *my_db_from_ptable(std::vector pt); int my_bp_distance(const char *str1, const char *str2); +%rename (dist_mountain) my_dist_mountain; + +%{ + double + my_dist_mountain( std::string str1, + std::string str2, + unsigned int p = 1) + { + return vrna_dist_mountain(str1.c_str(), str2.c_str(), p); + } + +%} + +double +my_dist_mountain(std::string str1, + std::string str2, + unsigned int p = 1); + /* * Rename 'struct vrna_ep_t' to target language class 'ep' diff --git a/m4/ac_rna.m4 b/m4/ac_rna.m4 index e4b7dc927..0305f0de3 100644 --- a/m4/ac_rna.m4 +++ b/m4/ac_rna.m4 @@ -52,12 +52,14 @@ RNA_ENABLE_GSL RNA_ENABLE_OPENMP RNA_ENABLE_PTHREADS RNA_ENABLE_BOUSTROPHEDON +RNA_ENABLE_NR_SAMPLE_HASH RNA_ENABLE_FLOATPF RNA_ENABLE_DEPRECATION_WARNINGS RNA_ENABLE_COLORED_TTY RNA_ENABLE_STATIC_BIN -RNA_ENABLE_SSE +RNA_ENABLE_SIMD RNA_ENABLE_VECTORIZE +RNA_ENABLE_MPFR ## Set post conditions for feature ## settings @@ -191,6 +193,16 @@ AS_IF([test "x$with_tutorial" != "xno"],[ _pdfdir="" ]) + +AS_IF([test "x$ac_rna_warning" != "x"],[ + ac_rna_warning=" +================================================== +Warning: +$ac_rna_warning +================================================== +" +]) + # Notify the user AC_MSG_NOTICE([ @@ -211,11 +223,13 @@ Extra Libraries --------------- * Support Vector Machine : ${with_svm:-no} * GNU Scientific Library : ${with_gsl:-no} + * GNU MPFR : ${enable_mpfr:-no} * JSON : ${with_json:-no} Features -------- * Boustrophedon : ${enable_boustrophedon:-no} + * Use hash for NR Sampling : ${enable_NRhash:-no} * C11 features : ${enable_c11:-no} * TTY colors : ${enable_tty_colors:-no} * Float Precision(PF} : ${enable_floatpf:-no} @@ -224,7 +238,7 @@ Features Optimizations ------------- * Auto Vectorization : ${enable_vectorize:-no} - * Streaming SIMD Extension : ${enable_sse:-no} + * Explicit SIMD Extension : ${enable_simd:-no} ${simd_failed} * Link Time Optimization : ${enable_lto:-no} * POSIX Threads : ${enable_pthreads:-no} * OpenMP : ${enable_openmp:-no} @@ -275,9 +289,8 @@ Install Directories * Python3 Interface : $_python3_install (binaries) : $_python3_arch_dir (scripts) : $_python3_lib_dir - -You can run 'make', 'make check' and 'make install' now! -]) +$ac_rna_warning +You can run 'make', 'make check' and 'make install' now!]) ]) diff --git a/m4/ac_rna_features.m4 b/m4/ac_rna_features.m4 index b94664e80..d74cedff4 100644 --- a/m4/ac_rna_features.m4 +++ b/m4/ac_rna_features.m4 @@ -99,6 +99,75 @@ AC_DEFUN([RNA_ENABLE_BOUSTROPHEDON],[ AC_SUBST(CONFIG_BOUSTROPHEDON) ]) +# +# Use hash for non-redundant sampling +# + +AC_DEFUN([RNA_ENABLE_NR_SAMPLE_HASH],[ + + RNA_ADD_FEATURE([NRhash], + [Hash for non-redundant sampling datas structure], + [no]) + + ## Add preprocessor define statement for Boustrophedon scheme in stochastic backtracking in part_func.c + RNA_FEATURE_IF_ENABLED([NRhash],[ + AC_DEFINE([VRNA_NR_SAMPLING_HASH], [1], [Use Hash for non-redundant sampling data structure]) + CONFIG_NR_SAMPLING="#define VRNA_NR_SAMPLING_HASH" + ]) + + AC_SUBST(CONFIG_NR_SAMPLING) +]) + + +AC_DEFUN([RNA_ENABLE_MPFR], [ + + RNA_ADD_FEATURE([mpfr], + [Use MPFR library for aribtrary precision computations in non-redundant sampling], + [yes]) + + RNA_FEATURE_IF_ENABLED([mpfr],[ + ## Check for mpfr.h header first + AC_CHECK_HEADER([mpfr.h], [ + ## now, check if we can compile a program + AC_MSG_CHECKING([whether we can compile programs with mpfr support]) + ac_save_LIBS="$LIBS" + LIBS="$ac_save_LIBS -lmpfr -lgmp" + + AC_LANG_PUSH([C]) + + AC_COMPILE_IFELSE([ + AC_LANG_PROGRAM( + [[#include + #include + ]], + [[ printf ("MPFR library: %-12s\nMPFR header: %s (based on %d.%d.%d)\n", + mpfr_get_version (), MPFR_VERSION_STRING, MPFR_VERSION_MAJOR, + MPFR_VERSION_MINOR, MPFR_VERSION_PATCHLEVEL); + return 0; + ]]) + ],[ + MPFR_LIBS="-lmpfr -lgmp" + AC_DEFINE([VRNA_NR_SAMPLING_MPFR], [1], [Use MPFR for non-redundant sampling data structure operations]) + ],[ + enable_mpfr=no + ]) + AC_LANG_POP([C]) + LIBS="$ac_save_LIBS" + AC_MSG_RESULT([$enable_mpfr]) + ], [ + AC_MSG_WARN([ +========================== +Failed to find mpfr.h! + +You probably need to install the mpfr-devel package or similar +========================== + ]) + enable_mpfr=no]) + ]) + + AC_SUBST(MPFR_LIBS) + AM_CONDITIONAL(VRNA_AM_SWITCH_MPFR, test "x$enable_mpfr" = "xyes") +]) # # OpenMP support @@ -355,41 +424,75 @@ into the executables are present! # -# SSE implementations +# SIMD optimizations # -AC_DEFUN([RNA_ENABLE_SSE],[ +AC_DEFUN([RNA_ENABLE_SIMD],[ + + RNA_ADD_FEATURE([simd], + [Speed-up MFE computations using explicit SIMD instructions.], + [yes]) RNA_ADD_FEATURE([sse], - [Speed-up MFE computations using SSE 4.1 implementations], + [Deprecated switch for SIMD optimizations. Use --enable-simd/--disable-simd instead], [no]) - ## Add preprocessor define statement for Boustrophedon scheme in stochastic backtracking in part_func.c - RNA_FEATURE_IF_ENABLED([sse],[ - if test "x$SIND_CFLAGS" = x; then - case $ax_cv_c_compiler_vendor in - gnu) - AC_LANG_PUSH([C]) - AX_CHECK_COMPILE_FLAG([-msse4.1], [ac_sse41_supported=yes],[ac_sse41_supported=no],[],[]) - AC_LANG_POP([C]) - ;; - intel) - ;; - *) - ;; - esac - fi + AS_IF([test "x$enable_sse" != "xno"],[ + AC_MSG_WARN([[ - if test $ac_sse41_supported = no; then - enable_sse=no; - fi +############################################ +Option --enable-sse is deprecated! + +Please consider using the successor option --enable-simd instead. +############################################ + ]]) + enable_simd="yes" + AC_RNA_ADD_WARNING([Deprecated option --enable-simd detected => Please use --enable-simd instead!]) ]) - RNA_FEATURE_IF_ENABLED([sse],[ - AC_MSG_CHECKING([compiler support for SSE4.1 min function]) + AS_IF([test "x$enable_simd" != "xno"],[ + ## Check for all supported SIMD features first + AC_MSG_CHECKING([compiler support for AVX 512 instructions]) + ac_save_CFLAGS="$CFLAGS" - CFLAGS="$ac_save_CFLAGS -msse4.1" - AC_LANG(C) + CFLAGS="$ac_save_CFLAGS -Werror -mavx512f" + AC_LANG_PUSH([C]) + + AC_COMPILE_IFELSE( + [ + AC_LANG_PROGRAM([[ + #include + #include + ]], + [[__m512i a = _mm512_set1_epi32(INT_MAX); + __m512i b = _mm512_set1_epi32(INT_MIN); + __m512i c = _mm512_set1_epi32(INT_MIN); + __mmask16 mask = _kand_mask16(_mm512_cmplt_epi32_mask(a, c), + _mm512_cmplt_epi32_mask(b, c)); + + b = _mm512_min_epi32(a, b); + int e = _mm512_mask_reduce_min_epi32(mask, b); + ]]) + ], + [ + AC_MSG_RESULT([yes]) + AC_DEFINE([VRNA_WITH_SIMD_AVX512], [1], [use AVX 512 implementations]) + ac_simd_capability_avx512f=yes + SIMD_AVX512_FLAGS="-mavx512f" + ], + [ + AC_MSG_RESULT([no]) + ]) + + AC_LANG_POP([C]) + CFLAGS="$ac_save_CFLAGS" + + AC_MSG_CHECKING([compiler support for SSE 4.1 instructions]) + + ac_save_CFLAGS="$CFLAGS" + CFLAGS="$ac_save_CFLAGS -Werror -msse4.1" + AC_LANG_PUSH([C]) + AC_COMPILE_IFELSE( [ AC_LANG_PROGRAM([[ @@ -403,19 +506,22 @@ AC_DEFUN([RNA_ENABLE_SSE],[ ], [ AC_MSG_RESULT([yes]) - SIMD_CFLAGS="${SIMD_CFLAGS} -msse4.1" - AC_DEFINE([VRNA_WITH_SSE_IMPLEMENTATION], [1], [use SSE implementations]) - CONFIG_SSE_IMPLEMENTATION="#define VRNA_WITH_SSE_IMPLEMENTATION" + AC_DEFINE([VRNA_WITH_SIMD_SSE41], [1], [use SSE 4.1 implementations]) + ac_simd_capability_sse41=yes + SIMD_SSE41_FLAGS="-msse4.1" ], [ - enable_sse=no - AC_MSG_RESULT([no; using default implementation]) + AC_MSG_RESULT([no]) ]) + + AC_LANG_POP([C]) CFLAGS="$ac_save_CFLAGS" ]) - AC_SUBST(SIMD_CFLAGS) - AC_SUBST(CONFIG_SSE_IMPLEMENTATION) + AC_SUBST(SIMD_AVX512_FLAGS) + AC_SUBST(SIMD_SSE41_FLAGS) + AM_CONDITIONAL(VRNA_AM_SWITCH_SIMD_AVX512, test "x$ac_simd_capability_avx512f" = "xyes") + AM_CONDITIONAL(VRNA_AM_SWITCH_SIMD_SSE41, test "x$ac_simd_capability_sse41" = "xyes") ]) @@ -426,7 +532,7 @@ AC_DEFUN([RNA_ENABLE_SSE],[ AC_DEFUN([RNA_ENABLE_VECTORIZE],[ RNA_ADD_FEATURE([vectorize], - [Apply SIMD vectorization to optimize execution speed], + [Apply automatic SIMD vectorization to optimize execution speed], [yes]) ## Add preprocessor define statement for Boustrophedon scheme in stochastic backtracking in part_func.c diff --git a/m4/ac_rna_macros.m4 b/m4/ac_rna_macros.m4 index 8c080f41e..3eaca7f92 100644 --- a/m4/ac_rna_macros.m4 +++ b/m4/ac_rna_macros.m4 @@ -5,6 +5,15 @@ ## Several macros for AC_RNA_INIT ## ##--------------------------------## +AC_DEFUN([AC_RNA_ADD_WARNING],[ + if test "x$ac_rna_warning" != "x"; then + ac_rna_warning="${ac_rna_warning} +$1" + else + ac_rna_warning="$1" + fi +]) + AC_DEFUN([AC_RNA_APPEND_VAR_COMMA],[ if test "x$$1" != "x" ; then $1="${$1}, " diff --git a/packaging/PKGBUILD.in b/packaging/PKGBUILD.in index 69a533e50..40bf98396 100644 --- a/packaging/PKGBUILD.in +++ b/packaging/PKGBUILD.in @@ -13,7 +13,9 @@ makedepends=( 'perl' 'python2' 'python' 'libtool' - 'check') + 'check' + 'mpfr' + 'gsl') source=(http://www.tbi.univie.ac.at/RNA/packages/source/ViennaRNA-${pkgver}.tar.gz) options=('staticlibs' '!strip') @@ -32,9 +34,10 @@ check() { } package_viennarna() { - depends=('perl') - optdepends=('gsl: use a variety of optimization methods in RNApvmin' - 'gri: create 2D plots of secondary structure landscape with RNA2Dfold') + depends=('perl' + 'mpfr' + 'gsl') + optdepends=('gri: create 2D plots of secondary structure landscape with RNA2Dfold') provides=('Kinfold=1.3' 'RNAforester=2.0' "viennarna2=${pkgver}" "libRNA=${pkgver}") cd "${srcdir}/ViennaRNA-${pkgver}" diff --git a/packaging/debian/changelog b/packaging/debian/changelog index b1c9433bf..b91bbe058 100644 --- a/packaging/debian/changelog +++ b/packaging/debian/changelog @@ -1,3 +1,9 @@ +viennarna (2.4.11-1) UNRELEASED; urgency=low + + * Bump to upstream release 2.4.11 + + -- Ronny Lorenz Mon, 17 Dec 2018 10:00:00 +0100 + viennarna (2.4.10-1) UNRELEASED; urgency=low * Bump to upstream release 2.4.10 diff --git a/packaging/debian/control b/packaging/debian/control index a2a07e51b..cdcd1fd66 100644 --- a/packaging/debian/control +++ b/packaging/debian/control @@ -4,11 +4,11 @@ Section: science Priority: optional Standards-Version: 3.9.5 X-Python-Version: >= 2.7 -Build-Depends: debhelper (>= 9), pkg-config, libgsl0-dev | libgsl-dev, gawk, perl, libperl-dev, python-all, python-all-dev, python3-all, python3-all-dev, lsb-release +Build-Depends: debhelper (>= 9), pkg-config, libgsl0-dev | libgsl-dev, libmpfr-dev, gawk, perl, libperl-dev, python-all, python-all-dev, python3-all, python3-all-dev, lsb-release Package: viennarna Architecture: any -Depends: libgomp1 (>= 4.4.5-1), libgsl0ldbl | libgsl2 | libgsl23, ${shlibs:Depends}, ${misc:Depends} +Depends: libgomp1 (>= 4.4.5-1), libgsl0ldbl | libgsl2 | libgsl23, libmpfr4, ${shlibs:Depends}, ${misc:Depends} Replaces: vienna-rna (<< 2.2.1-2) Breaks: vienna-rna (<< 2.2.1-2) Recommends: python-rna, perl-rna diff --git a/packaging/make_windows_installer.sh b/packaging/make_windows_installer.sh index 2bb745e38..37fa8dd9e 100755 --- a/packaging/make_windows_installer.sh +++ b/packaging/make_windows_installer.sh @@ -12,16 +12,18 @@ CONFIGURE_OPTIONS=" --without-swig \ --with-cluster \ --with-kinwalker \ --disable-pthreads \ - --disable-tty-colors" + --disable-tty-colors \ + --disable-lto" case "$1" in arch64) echo "Making Windows Installer using Arch Linux mingw-w64 installation" echo -ne "...making 32bit version..." cd .. + >packaging/${WIN_INSTALLER_LOG} ./configure --host=i686-w64-mingw32 ${CONFIGURE_OPTIONS} >> packaging/${WIN_INSTALLER_LOG} 2>&1 make clean >> packaging/${WIN_INSTALLER_LOG} 2>&1 - make >>${WIN_INSTALLER_LOG} 2>&1 + make >> packaging/${WIN_INSTALLER_LOG} 2>&1 cd packaging makensis win_installer_archlinux_i686.nsi >>${WIN_INSTALLER_LOG} 2>&1 echo -ne " done\n" @@ -29,7 +31,7 @@ arch64) echo "Making Windows Installer using Arch Linux mingw-w64 installation" cd .. ./configure --host=x86_64-w64-mingw32 ${CONFIGURE_OPTIONS} >> packaging/${WIN_INSTALLER_LOG} 2>&1 make clean >> packaging/${WIN_INSTALLER_LOG} 2>&1 - make >>${WIN_INSTALLER_LOG} 2>&1 + make >> packaging/${WIN_INSTALLER_LOG} 2>&1 cd packaging makensis win_installer_archlinux_x86_64.nsi >>${WIN_INSTALLER_LOG} 2>&1 echo -ne " done\n" diff --git a/packaging/viennarna.spec.in b/packaging/viennarna.spec.in index 232ca2188..24e684d08 100644 --- a/packaging/viennarna.spec.in +++ b/packaging/viennarna.spec.in @@ -48,11 +48,12 @@ BuildRequires: rpm-devel BuildRequires: libstdc++-devel BuildRequires: gcc gcc-c++ glibc-devel info -BuildRequires: gsl-devel +BuildRequires: gsl-devel +BuildRequires: mpfr-devel BuildRequires: perl(ExtUtils::Embed) BuildRequires: perl-Test-Simple -Requires: perl glibc info gsl +Requires: perl glibc info gsl mpfr Provides: Kinfold %if 0%{?suse_version} diff --git a/packaging/win_installer_archlinux_i686.nsi.in b/packaging/win_installer_archlinux_i686.nsi.in index 5a273beef..f1559895a 100644 --- a/packaging/win_installer_archlinux_i686.nsi.in +++ b/packaging/win_installer_archlinux_i686.nsi.in @@ -120,7 +120,10 @@ section "Core" sec_core File "/usr/i686-w64-mingw32/bin/libgcc_s_sjlj-1.dll" File "/usr/i686-w64-mingw32/bin/libgsl-23.dll" File "/usr/i686-w64-mingw32/bin/libgslcblas-0.dll" - + File "/usr/i686-w64-mingw32/bin/libmpfr-6.dll" + File "/usr/i686-w64-mingw32/bin/libgmp-10.dll" + File "/usr/i686-w64-mingw32/bin/libgmpxx-4.dll" + File ".local" # we also want an uninstall to be installed @@ -236,6 +239,9 @@ section "Uninstall" delete $INSTDIR\libgcc_s_sjlj-1.dll delete $INSTDIR\libgsl-23.dll delete $INSTDIR\libgslcblas-0.dll + delete $INSTDIR\libmpfr-6.dll + delete $INSTDIR\libgmp-10.dll + delete $INSTDIR\libgmpxx-4.dll delete $INSTDIR\.local delete "$INSTDIR\Uninstall-${PACKAGE}.exe" diff --git a/packaging/win_installer_archlinux_x86_64.nsi.in b/packaging/win_installer_archlinux_x86_64.nsi.in index 6c6da0f80..2ccb1412a 100644 --- a/packaging/win_installer_archlinux_x86_64.nsi.in +++ b/packaging/win_installer_archlinux_x86_64.nsi.in @@ -120,7 +120,10 @@ section "Core" sec_core File "/usr/x86_64-w64-mingw32/bin/libgcc_s_seh-1.dll" File "/usr/x86_64-w64-mingw32/bin/libgsl-23.dll" File "/usr/x86_64-w64-mingw32/bin/libgslcblas-0.dll" - + File "/usr/x86_64-w64-mingw32/bin/libmpfr-6.dll" + File "/usr/x86_64-w64-mingw32/bin/libgmp-10.dll" + File "/usr/x86_64-w64-mingw32/bin/libgmpxx-4.dll" + File ".local" # we also want an uninstall to be installed @@ -236,6 +239,9 @@ section "Uninstall" delete $INSTDIR\libgcc_s_seh-1.dll delete $INSTDIR\libgsl-23.dll delete $INSTDIR\libgslcblas-0.dll + delete $INSTDIR\libmpfr-6.dll + delete $INSTDIR\libgmp-10.dll + delete $INSTDIR\libgmpxx-4.dll delete $INSTDIR\.local delete "$INSTDIR\Uninstall-${PACKAGE}.exe" diff --git a/packaging/win_installer_fedora_i686.nsi.in b/packaging/win_installer_fedora_i686.nsi.in index baf910a10..ac1ce6b39 100644 --- a/packaging/win_installer_fedora_i686.nsi.in +++ b/packaging/win_installer_fedora_i686.nsi.in @@ -120,7 +120,10 @@ section "Core" sec_core File "/usr/i686-w64-mingw32/sys-root/mingw/bin/libgcc_s_sjlj-1.dll" File "/usr/i686-w64-mingw32/sys-root/mingw/bin/libgsl-0.dll" File "/usr/i686-w64-mingw32/sys-root/mingw/bin/libgslcblas-0.dll" - + File "/usr/i686-w64-mingw32/sys-root/mingw/bin/libmpfr-4.dll" + File "/usr/i686-w64-mingw32/sys-root/mingw/bin/libgmp-10.dll" + File "/usr/i686-w64-mingw32/sys-root/mingw/bin/libgmpxx-4.dll" + File ".local" # we also want an uninstall to be installed @@ -236,6 +239,9 @@ section "Uninstall" delete $INSTDIR\libgcc_s_sjlj-1.dll delete $INSTDIR\libgsl-0.dll delete $INSTDIR\libgslcblas-0.dll + delete $INSTDIR\libmpfr-4.dll + delete $INSTDIR\libgmp-10.dll + delete $INSTDIR\libgmpxx-4.dll delete $INSTDIR\.local delete "$INSTDIR\Uninstall-${PACKAGE}.exe" diff --git a/packaging/win_installer_fedora_x86_64.nsi.in b/packaging/win_installer_fedora_x86_64.nsi.in index 76f903014..c7cb2f9f7 100644 --- a/packaging/win_installer_fedora_x86_64.nsi.in +++ b/packaging/win_installer_fedora_x86_64.nsi.in @@ -68,7 +68,7 @@ outfile "Install-ViennaRNA-${VERSION}_64bit.exe" !insertmacro MUI_UNPAGE_FINISH !insertmacro MUI_LANGUAGE "English" - + # default section section "Core" sec_core SectionIn RO @@ -120,6 +120,9 @@ section "Core" sec_core File "/usr/x86_64-w64-mingw32/sys-root/mingw/bin/libgcc_s_seh-1.dll" File "/usr/x86_64-w64-mingw32/sys-root/mingw/bin/libgsl-0.dll" File "/usr/x86_64-w64-mingw32/sys-root/mingw/bin/libgslcblas-0.dll" + File "/usr/x86_64-w64-mingw32/sys-root/mingw/bin/libmpfr-4.dll" + File "/usr/x86_64-w64-mingw32/sys-root/mingw/bin/libgmp-10.dll" + File "/usr/x86_64-w64-mingw32/sys-root/mingw/bin/libgmpxx-4.dll" File ".local" @@ -236,6 +239,9 @@ section "Uninstall" delete $INSTDIR\libgcc_s_seh-1.dll delete $INSTDIR\libgsl-0.dll delete $INSTDIR\libgslcblas-0.dll + delete $INSTDIR\libmpfr-4.dll + delete $INSTDIR\libgmp-10.dll + delete $INSTDIR\libgmpxx-4.dll delete $INSTDIR\.local delete "$INSTDIR\Uninstall-${PACKAGE}.exe" diff --git a/src/ViennaRNA/Makefile.am b/src/ViennaRNA/Makefile.am index 3109bc284..894184ef5 100644 --- a/src/ViennaRNA/Makefile.am +++ b/src/ViennaRNA/Makefile.am @@ -1,7 +1,7 @@ AUTOMAKE_OPTIONS = subdir-objects -AM_CFLAGS = $(RNA_CFLAGS) $(SIMD_CFLAGS) $(PTHREAD_CFLAGS) -AM_CXXFLAGS = $(RNA_CXXFLAGS) $(SIMD_CFLAGS) $(PTHREAD_CFLAGS) +AM_CFLAGS = $(RNA_CFLAGS) $(PTHREAD_CFLAGS) +AM_CXXFLAGS = $(RNA_CXXFLAGS) $(PTHREAD_CFLAGS) AM_CPPFLAGS = $(RNA_CPPFLAGS) ${SVM_INC} -I$(top_srcdir)/src ${JSON_INC} AM_LDFLAGS = $(RNA_LDFLAGS) $(PTHREAD_LIBS) @@ -47,6 +47,19 @@ libRNA_conv_la_LIBADD = \ libRNA_eval.la \ libRNA_loops.la + +if VRNA_AM_SWITCH_SIMD_SSE41 +noinst_LTLIBRARIES += libRNA_utils_sse41.la +libRNA_conv_la_LIBADD += libRNA_utils_sse41.la +libRNA_utils_sse41_la_CFLAGS = $(SIMD_SSE41_FLAGS) +endif + +if VRNA_AM_SWITCH_SIMD_AVX512 +noinst_LTLIBRARIES += libRNA_utils_avx512.la +libRNA_conv_la_LIBADD += libRNA_utils_avx512.la +libRNA_utils_avx512_la_CFLAGS = $(SIMD_AVX512_FLAGS) +endif + # Dummy C++ source to cause C++ linking. if VRNA_AM_SWITCH_SVM nodist_EXTRA_libRNA_la_SOURCES = dummy.cxx @@ -181,6 +194,8 @@ vrna_utils_HEADERS = \ utils/strings.h \ utils/structures.h \ utils/alignments.h \ + utils/higher_order_functions.h \ + utils/cpu.h \ ${SVM_UTILS_H} @@ -292,6 +307,8 @@ libRNA_utils_la_SOURCES = \ utils/structure_utils.c \ utils/structure_tree.c \ utils/msa_utils.c \ + utils/higher_order_functions.c \ + utils/cpu.c \ io/io_utils.c \ io/file_formats.c \ io/file_formats_msa.c \ @@ -301,6 +318,16 @@ libRNA_utils_la_SOURCES = \ combinatorics.c \ ${SVM_UTILS} +if VRNA_AM_SWITCH_SIMD_SSE41 +libRNA_utils_sse41_la_SOURCES = \ + utils/higher_order_functions_sse41.c +endif + +if VRNA_AM_SWITCH_SIMD_AVX512 +libRNA_utils_avx512_la_SOURCES = \ + utils/higher_order_functions_avx512.c +endif + libRNA_plotting_la_SOURCES = \ plotting/alignments.c \ plotting/layouts.c \ @@ -382,6 +409,7 @@ EXTRA_DIST = $(pkginclude_HEADERS) \ loops/multibranch_sc_pf.inc \ params/svm_model_avg.inc \ params/svm_model_sd.inc \ + data_structures_nonred.inc \ ${SVM_H} \ ${JSON_H} \ color_output.inc \ diff --git a/src/ViennaRNA/alphabet.c b/src/ViennaRNA/alphabet.c index e2ed7782b..a006c764d 100644 --- a/src/ViennaRNA/alphabet.c +++ b/src/ViennaRNA/alphabet.c @@ -297,7 +297,7 @@ vrna_seq_encode_simple(const char *sequence, S = (short *)vrna_alloc(sizeof(short) * (l + 2)); for (i = 1; i <= l; i++) /* make numerical encoding of sequence */ - S[i] = (short)vrna_nucleotide_encode(toupper(sequence[i - 1]), md); + S[i] = (short)vrna_nucleotide_encode(sequence[i - 1], md); S[l + 1] = S[1]; S[0] = (short)l; @@ -314,6 +314,8 @@ vrna_nucleotide_encode(char c, /* return numerical representation of nucleotide used e.g. in vrna_md_t.pair[][] */ int code = -1; + c = toupper(c); + if (md) { if (md->energy_set > 0) { code = (int)(c - 'A') + 1; diff --git a/src/ViennaRNA/boltzmann_sampling.c b/src/ViennaRNA/boltzmann_sampling.c index c3aeedd4e..57dfc902d 100644 --- a/src/ViennaRNA/boltzmann_sampling.c +++ b/src/ViennaRNA/boltzmann_sampling.c @@ -12,6 +12,7 @@ #include #include #include +#include #include #include "ViennaRNA/utils/basic.h" @@ -24,6 +25,33 @@ #include "ViennaRNA/alphabet.h" #include "ViennaRNA/boltzmann_sampling.h" + +#include "ViennaRNA/data_structures_nonred.inc" + +/* + ################################# + # PREPROCESSOR DEFININTIONS # + ################################# + */ + +#ifdef VRNA_NR_SAMPLING_HASH +# define NR_NODE tr_node +# define NR_TOTAL_WEIGHT(a) total_weight_par(a) +# define NR_TOTAL_WEIGHT_TYPE(a, b) total_weight_par_type(a, b) +# define NR_GET_WEIGHT(a, b, c, d, e) tr_node_weight(a, c, d, e) +#else +# define NR_NODE tllr_node +# define NR_TOTAL_WEIGHT(a) get_weight_all(a) +# define NR_TOTAL_WEIGHT_TYPE(a, b) get_weight_type_spec(a, b) +# define NR_GET_WEIGHT(a, b, c, d, e) get_weight(b, c, d, e) +#endif + + +struct nr_structure_list { + unsigned int num; + char **list; +}; + /* ################################# # GLOBAL VARIABLES # @@ -46,45 +74,96 @@ PRIVATE char *info_set_uniq_ml = # PRIVATE FUNCTION DECLARATIONS # ################################# */ -PRIVATE void backtrack(int i, - int j, - char *pstruc, - vrna_fold_compound_t *vc); -PRIVATE void backtrack_qm(int i, - int j, - char *pstruc, - vrna_fold_compound_t *vc); +PRIVATE void +save_nr_samples(const char *structure, + void *data); -PRIVATE void backtrack_qm1(int i, - int j, - char *pstruc, - vrna_fold_compound_t *vc); +/* In the following: + * - q_remain is a pointer to value of sum of Boltzmann factors of still accessible solutions at that point + * - current_node is a double pointer to current node in datastructure memorizing the solutions and paths taken */ +PRIVATE char * +pbacktrack5_gen(vrna_fold_compound_t *vc, + int length, + double *q_remain, + NR_NODE **current_node, + struct nr_memory **memory_dat); + + +PRIVATE int +backtrack(int i, + int j, + char *pstruc, + vrna_fold_compound_t *vc, + double *q_remain, + NR_NODE **current_node, + struct nr_memory **memory_dat); -PRIVATE void backtrack_qm2(int u, - int n, - char *pstruc, - vrna_fold_compound_t *vc); +PRIVATE int +backtrack_ext_loop(int init_val, + char *pstruc, + vrna_fold_compound_t *vc, + int length, + double *q_remain, + NR_NODE **current_node, + struct nr_memory **memory_dat); -PRIVATE char *wrap_pbacktrack_circ(vrna_fold_compound_t *vc); +PRIVATE int +backtrack_qm(int i, + int j, + char *pstruc, + vrna_fold_compound_t *vc); -PRIVATE void backtrack_comparative(vrna_fold_compound_t *vc, - char *pstruc, - int i, - int j, - double *prob); +PRIVATE int +backtrack_qm_nr(int i, + int j, + char *pstruc, + vrna_fold_compound_t *vc, + double *q_remain, + NR_NODE **current_node, + struct nr_memory **memory_dat); -PRIVATE void backtrack_qm1_comparative(vrna_fold_compound_t *vc, - char *pstruc, - int i, - int j, - double *prob); +PRIVATE int +backtrack_qm1(int i, + int j, + char *pstruc, + vrna_fold_compound_t *vc, + double *q_remain, + NR_NODE **current_node, + struct nr_memory **memory_dat); + + +PRIVATE void +backtrack_qm2(int u, + int n, + char *pstruc, + vrna_fold_compound_t *vc); + + +PRIVATE char * +wrap_pbacktrack_circ(vrna_fold_compound_t *vc); + + +PRIVATE void +backtrack_comparative(vrna_fold_compound_t *vc, + char *pstruc, + int i, + int j, + double *prob); + + +PRIVATE void +backtrack_qm1_comparative(vrna_fold_compound_t *vc, + char *pstruc, + int i, + int j, + double *prob); /* @@ -98,8 +177,9 @@ PRIVATE void backtrack_qm1_comparative(vrna_fold_compound_t *vc, * @param prob to be described (berni) * @return A sampled consensus secondary structure in dot-bracket notation */ -PRIVATE char *pbacktrack_comparative(vrna_fold_compound_t *vc, - double *prob); +PRIVATE char * +pbacktrack_comparative(vrna_fold_compound_t *vc, + double *prob); /* @@ -153,13 +233,167 @@ vrna_pbacktrack(vrna_fold_compound_t *vc) } +/* adapter function for more general expression using non-redundant sampling */ PUBLIC char * vrna_pbacktrack5(vrna_fold_compound_t *vc, int length) { - FLT_OR_DBL r, qt, q_temp, qkl; - int i, j, ij, n, k, u, type; - char *pstruc; + return pbacktrack5_gen(vc, length, NULL, NULL, NULL); +} + + +PUBLIC char ** +vrna_pbacktrack_nr(vrna_fold_compound_t *vc, + int num_samples) +{ + struct nr_structure_list data; + + data.num = 0; + data.list = (char **)vrna_alloc(sizeof(char *) * num_samples); + data.list[0] = NULL; + vrna_pbacktrack_nr_cb(vc, num_samples, &save_nr_samples, (void *)&data); + + /* re-allocate memory */ + data.list = (char **)vrna_realloc(data.list, sizeof(char *) * (data.num + 1)); + data.list[data.num] = NULL; + + return data.list; +} + + +PUBLIC void +vrna_pbacktrack_nr_cb(vrna_fold_compound_t *vc, + int num_samples, + vrna_boltzmann_sampling_callback *bs_cb, + void *data) +{ + if (vc) { + if (!vc->exp_params) { + vrna_message_warning("vrna_pbacktrack_nr_cb: DP matrices are missing! Call vrna_pf() first!"); + } else if (!vc->exp_params->model_details.uniq_ML) { + vrna_message_warning("vrna_pbacktrack_nr_cb: Unique multiloop decomposition is unset!"); + vrna_message_info(stderr, info_set_uniq_ml); + } else if (vc->type != VRNA_FC_TYPE_SINGLE) { + vrna_message_warning( + "vrna_pbacktrack_nr_cb: No implementation for comparative structure prediction available yet!"); + } else if (vc->exp_params->model_details.circ) { + vrna_message_warning( + "vrna_pbacktrack_nr_cb: No implementation for circular RNAs available yet!"); + } else { + int i, is_dup, pf_overflow; + double q_remain, part_fci; + size_t block_size; + NR_NODE *current_node; + NR_NODE *root_node; + struct nr_memory *memory_dat; + + memory_dat = NULL; + block_size = 5000 * sizeof(NR_NODE); + q_remain = 0; + pf_overflow = 0; + part_fci = vc->exp_matrices->q[vc->iindx[1] - vc->length]; + +#ifdef VRNA_NR_SAMPLING_HASH + root_node = create_root(vc->length, part_fci); +#else + memory_dat = create_nr_memory(sizeof(NR_NODE), block_size, NULL); // memory pre-allocation + root_node = create_ll_root(&memory_dat, part_fci); + +#endif + current_node = root_node; + + for (i = 0; i < num_samples; i++) { + is_dup = 1; + q_remain = vc->exp_matrices->q[vc->iindx[1] - vc->length]; + char *ss = pbacktrack5_gen(vc, vc->length, &q_remain, ¤t_node, &memory_dat); + +#ifdef VRNA_NR_SAMPLING_HASH + current_node = traceback_to_root(current_node, q_remain, &is_dup, &pf_overflow); +#else + current_node = traceback_to_ll_root(current_node, q_remain, &is_dup, &pf_overflow); +#endif + + if (pf_overflow) { + vrna_message_warning("vrna_pbacktrack_nr*(): " + "Partition function overflow detected for forbidden structures, presumably due to numerical instabilities."); + free(ss); + break; + } + + if (is_dup) { + vrna_message_warning("vrna_pbacktrack_nr*(): " + "Duplicate structures detected, presumably due to numerical instabilities"); + free(ss); + break; + } + + bs_cb(ss, data); + free(ss); + /* finish if no more structures available */ + if (!ss) + break; + } + /* print warning if we've aborted backtracking too early */ + if ((i > 0) && (i < num_samples)) { + vrna_message_warning("vrna_pbacktrack_nr*(): " + "Stopped backtracking after %d samples due to numeric instabilities!\n" + "Coverage of partition function so far: %.6f%%", + i, + 100. * return_node_weight(root_node) / part_fci); + } + +#ifdef VRNA_NR_SAMPLING_HASH + free_all_nr(current_node); +#else + free_all_nrll(&memory_dat); +#endif + } + } +} + + +/* general expr of vrna5_pbacktrack with possibility of non-redundant sampling */ +PRIVATE char * +pbacktrack5_gen(vrna_fold_compound_t *vc, + int length, + double *q_remain, + NR_NODE **current_node, + struct nr_memory **memory_dat) +{ + int ret, i; + char *pstruc; + + pstruc = vrna_alloc((length + 1) * sizeof(char)); + + for (i = 0; i < length; i++) + pstruc[i] = '.'; + +#ifdef VRNA_WITH_BOUSTROPHEDON + ret = backtrack_ext_loop(length, pstruc, vc, length, q_remain, current_node, memory_dat); +#else + ret = backtrack_ext_loop(1, pstruc, vc, length, q_remain, current_node, memory_dat); +#endif + + if (ret > 0) + return pstruc; + + free(pstruc); + return NULL; +} + + +/* backtrack one external */ +PRIVATE int +backtrack_ext_loop(int init_val, + char *pstruc, + vrna_fold_compound_t *vc, + int length, + double *q_remain, + NR_NODE **current_node, + struct nr_memory **memory_dat) +{ + FLT_OR_DBL r, fbd, fbds, qt, q_temp, qkl; + int ret, i, j, ij, n, k, u, type; int *my_iindx, hc_decompose, *hc_up_ext; FLT_OR_DBL *q, *qb, *q1k, *qln, *scale; unsigned char *hard_constraints; @@ -170,6 +404,15 @@ vrna_pbacktrack5(vrna_fold_compound_t *vc, vrna_sc_t *sc; vrna_exp_param_t *pf_params; +#ifndef VRNA_NR_SAMPLING_HASH + /* non-redundant data-structure memorization nodes */ + NR_NODE *memorized_node_prev = NULL; /* remembers previous-to-current node in linked list */ + NR_NODE *memorized_node_cur = NULL; /* remembers actual node in linked list */ +#endif + + fbd = 0.; /* stores weight of forbidden terms for given q[ij]*/ + fbds = 0.; /* stores weight of forbidden term for given motif */ + n = vc->length; pf_params = vc->exp_params; @@ -187,30 +430,28 @@ vrna_pbacktrack5(vrna_fold_compound_t *vc, if (length > n) { vrna_message_warning("vrna_pbacktrack5: 3'-end exceeds sequence length"); - return NULL; + return -1; } else if (length < 1) { vrna_message_warning("vrna_pbacktrack5: 3'-end too small"); - return NULL; + return -1; } else if ((!matrices) || (!matrices->q) || (!matrices->qb) || (!matrices->qm) || (!pf_params)) { vrna_message_warning("vrna_pbacktrack5: DP matrices are missing! Call vrna_pf() first!"); - return NULL; + return -1; } else if ((!vc->exp_params->model_details.uniq_ML) || (!matrices->qm1)) { vrna_message_warning("vrna_pbacktrack5: Unique multiloop decomposition is unset!"); vrna_message_info(stderr, info_set_uniq_ml); - return NULL; + return -1; } + /* assume successful backtracing by default */ + ret = 1; + q = matrices->q; qb = matrices->qb; q1k = matrices->q1k; qln = matrices->qln; scale = matrices->scale; - pstruc = vrna_alloc((length + 1) * sizeof(char)); - - for (i = 0; i < length; i++) - pstruc[i] = '.'; - if (!(q1k && qln)) { matrices->q1k = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (n + 1)); matrices->qln = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (n + 2)); @@ -224,14 +465,38 @@ vrna_pbacktrack5(vrna_fold_compound_t *vc, qln[n + 1] = 1.0; } +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) { + memorized_node_prev = NULL; + memorized_node_cur = (*current_node)->head; + } + +#endif + + q_temp = 0.; + #ifdef VRNA_WITH_BOUSTROPHEDON - j = length; - while (j > 1) { + j = init_val; + if (j > 1) { /* find j position of first pair */ for (; j > 1; j--) { if (hc_up_ext[j]) { - r = vrna_urn() * q1k[j]; - q_temp = q[my_iindx[1] - j + 1] * scale[1]; + if (current_node) { + fbd = NR_TOTAL_WEIGHT(*current_node) * + q1k[j] / + (*q_remain); + +#ifdef USE_FLOAT_PF + if (fabsf(NR_TOTAL_WEIGHT(*current_node) - (*q_remain)) / (*q_remain) <= FLT_EPSILON) +#else + if (fabs(NR_TOTAL_WEIGHT(*current_node) - (*q_remain)) / (*q_remain) <= DBL_EPSILON) +#endif + /* exhausted ensemble */ + return 0; + } + + r = vrna_urn() * (q1k[j] - fbd); + q_temp = q1k[j - 1] * scale[1]; if (sc) { if (sc->exp_energy_up) @@ -241,16 +506,51 @@ vrna_pbacktrack5(vrna_fold_compound_t *vc, q_temp *= sc->exp_f(1, j, 1, j - 1, VRNA_DECOMP_EXT_EXT, sc->data); } - if (r > q_temp) + if (current_node) { + fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_UNPAIRED_SG, j - 1, j) * + q1k[j] / + (*q_remain); + } + + if (r > (q_temp - fbds)) { break; /* j is paired */ + } else if (current_node) { + /* j is unpaired */ + *q_remain *= q_temp / q1k[j]; +#ifdef VRNA_NR_SAMPLING_HASH + *current_node = add_if_nexists(NRT_UNPAIRED_SG, j - 1, j, *current_node, *q_remain); +#else + *current_node = add_if_nexists_ll(memory_dat, + NRT_UNPAIRED_SG, + j - 1, + j, + memorized_node_prev, + memorized_node_cur, + *current_node, + *q_remain); + reset_cursor(&memorized_node_prev, &memorized_node_cur, *current_node); /* resets cursor */ +#endif + } } } if (j <= md->min_loop_size + 1) - break; /* no more pairs */ + return ret; /* no more pairs, but still successful */ + +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) + advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_UNPAIRED_SG, j - 1, j); +#endif /* now find the pairing partner i */ - r = vrna_urn() * (q1k[j] - q_temp); + if (current_node) { + fbd = NR_TOTAL_WEIGHT_TYPE(NRT_EXT_LOOP, *current_node) * + q1k[j] / + (*q_remain); + } + + r = vrna_urn() * (q1k[j] - q_temp - fbd); u = j - 1; + i = 2; for (qt = 0, k = 1; k < j; k++) { /* apply alternating boustrophedon scheme to variable i */ @@ -276,24 +576,71 @@ vrna_pbacktrack5(vrna_fold_compound_t *vc, qkl *= sc->exp_f(i, j, i, j, VRNA_DECOMP_EXT_STEM, sc->data); } - qt += qkl; - if (qt > r) + if (current_node) { + fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_EXT_LOOP, i, j) * + q1k[j] / + (*q_remain); + qt += qkl - fbds; + } else { + qt += qkl; + } + + if (qt > r) { + if (current_node) { + *q_remain *= qkl / q1k[j]; +#ifdef VRNA_NR_SAMPLING_HASH + *current_node = add_if_nexists(NRT_EXT_LOOP, i, j, *current_node, *q_remain); +#else + *current_node = add_if_nexists_ll(memory_dat, + NRT_EXT_LOOP, + i, + j, + memorized_node_prev, + memorized_node_cur, + *current_node, + *q_remain); +#endif + } + break; /* j is paired */ + } + +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) + advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_EXT_LOOP, i, j); + +#endif + } + } + if (k == j) { + if (current_node) { + /* exhausted ensemble */ + return 0; + } else { + vrna_message_warning("backtracking failed in ext loop"); + /* error */ + return -1; } } - if (k == j) - vrna_message_error("backtracking failed in ext loop"); - backtrack(i, j, pstruc, vc); - j = i - 1; + backtrack(i, j, pstruc, vc, q_remain, current_node, memory_dat); + j = i - 1; + ret = backtrack_ext_loop(j, pstruc, vc, length, q_remain, current_node, memory_dat); } + #else - start = 1; - while (start < length) { + int start = init_val; + if (start < length) { /* find i position of first pair */ for (i = start; i < length; i++) { if (hc_up_ext[i]) { - r = vrna_urn() * qln[i]; + if (current_node) { + fbd = NR_TOTAL_WEIGHT(*current_node) * + qln[i] / + (*q_remain); + } + + r = vrna_urn() * (qln[i] - fbd); q_temp = qln[i + 1] * scale[1]; if (sc) { @@ -304,18 +651,52 @@ vrna_pbacktrack5(vrna_fold_compound_t *vc, q_temp *= sc->exp_f(i, length, i + 1, length, VRNA_DECOMP_EXT_EXT, sc->data); } - if (r > q_temp) + if (current_node) { + fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_UNPAIRED_SG, i, i + 1) * + qln[i] / + (*q_remain); + } + + if (r > (q_temp - fbds)) { break; /* i is paired */ + } else if (current_node) { + *q_remain *= q_temp / qln[i]; +#ifdef VRNA_NR_SAMPLING_HASH + *current_node = add_if_nexists(NRT_UNPAIRED_SG, i, i + 1, *current_node, *q_remain); +#else + *current_node = add_if_nexists_ll(memory_dat, + NRT_UNPAIRED_SG, + i, + i + 1, + memorized_node_prev, + memorized_node_cur, + *current_node, + *q_remain); + reset_cursor(&memorized_node_prev, &memorized_node_cur, *current_node); /* resets cursor */ +#endif + } } } if (i >= length) - break; /* no more pairs */ + return ret; /* no more pairs, but still successful */ + +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) + advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_UNPAIRED_SG, i, i + 1); + +#endif /* now find the pairing partner j */ - r = vrna_urn() * (qln[i] - q_temp); + if (current_node) { + fbd = NR_TOTAL_WEIGHT_TYPE(NRT_EXT_LOOP, *current_node) * + qln[i] / + (*q_remain); + } + + r = vrna_urn() * (qln[i] - q_temp - fbd); for (qt = 0, j = i + 1; j <= length; j++) { ij = my_iindx[i] - j; - type = vrna_get_ptype(jindx[j] + i, ptype); + type = vrna_get_ptype_md(S2[i], S2[j], md); hc_decompose = hard_constraints[n * i + j]; if (hc_decompose & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) { qkl = qb[ij] * exp_E_ExtLoop(type, @@ -334,23 +715,79 @@ vrna_pbacktrack5(vrna_fold_compound_t *vc, qkl *= sc->exp_f(i, j, i, j, VRNA_DECOMP_EXT_STEM, sc->data); } - qt += qkl; - if (qt > r) - break; /* j is paired */ + if (current_node) { + fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_EXT_LOOP, i, j) * + qln[i] / + (*q_remain); + qt += qkl - fbds; + } else { + qt += qkl; + } + + if (qt > r) { + if (current_node) { + *q_remain *= qkl / qln[i]; +#ifdef VRNA_NR_SAMPLING_HASH + *current_node = add_if_nexists(NRT_EXT_LOOP, i, j, *current_node, *q_remain); +#else + *current_node = add_if_nexists_ll(memory_dat, + NRT_EXT_LOOP, + i, + j, + memorized_node_prev, + memorized_node_cur, + *current_node, + *q_remain); +#endif + } + + break; /* j is paired */ + } + +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) + advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_EXT_LOOP, i, j); + +#endif + } + } + if (j == length + 1) { + if (current_node) { + /* exhausted ensemble */ + return 0; + } else { + vrna_message_warning("backtracking failed in ext loop"); + /* error */ + return -1; } } - if (j == length + 1) - vrna_message_error("backtracking failed in ext loop"); start = j + 1; - backtrack(i, j, pstruc, vc); + backtrack(i, j, pstruc, vc, q_remain, current_node, memory_dat); + + ret = backtrack_ext_loop(start, pstruc, vc, length, q_remain, current_node, memory_dat); } + #endif - return pstruc; + + return ret; } PRIVATE void +save_nr_samples(const char *structure, + void *data) +{ + struct nr_structure_list *d = (struct nr_structure_list *)data; + + if (structure) + d->list[d->num++] = strdup(structure); + else + d->list[d->num++] = NULL; +} + + +PRIVATE int backtrack_qm(int i, int j, char *pstruc, @@ -360,12 +797,13 @@ backtrack_qm(int i, FLT_OR_DBL qmt, r, q_temp; int k, u, cnt, span, turn; FLT_OR_DBL *qm, *qm1, *expMLbase; - int *my_iindx, *jindx, *hc_up_ml; + int *my_iindx, *jindx, *hc_up_ml, ret; vrna_sc_t *sc; vrna_hc_t *hc; vrna_mx_pf_t *matrices = vc->exp_matrices; + ret = 1; my_iindx = vc->iindx; jindx = vc->jindx; @@ -423,10 +861,15 @@ backtrack_qm(int i, } } - if (cnt > j) + if (cnt > j) { vrna_message_error("backtrack failed in qm"); + return 0; + } - backtrack_qm1(k, j, pstruc, vc); + ret = backtrack_qm1(k, j, pstruc, vc, NULL, NULL, NULL); + + if (ret == 0) + return ret; if (k < i + turn) break; /* no more pairs */ @@ -451,19 +894,235 @@ backtrack_qm(int i, j = k - 1; } + + return ret; } -PRIVATE void +/* non redundant version of function bactrack_qm */ +PRIVATE int +backtrack_qm_nr(int i, + int j, + char *pstruc, + vrna_fold_compound_t *vc, + double *q_remain, + NR_NODE **current_node, + struct nr_memory **memory_dat) +{ + /* divide multiloop into qm and qm1 */ + FLT_OR_DBL qmt, fbd, fbds, r, q_temp; + int k, u, cnt, span, turn; + int is_unpaired; /* 1 if [i ... k-1] is unpaired */ + FLT_OR_DBL *qm, *qm1, *expMLbase; + int *my_iindx, *jindx, *hc_up_ml, ret; + vrna_sc_t *sc; + vrna_hc_t *hc; + +#ifndef VRNA_NR_SAMPLING_HASH + /* non-redundant data-structure memorization nodes */ + NR_NODE *memorized_node_prev = NULL; /* remembers previous-to-current node in linked list */ + NR_NODE *memorized_node_cur = NULL; /* remembers actual node in linked list */ +#endif + + ret = 1; + fbd = 0.; /* stores weight of forbidden terms for given q[ij]*/ + fbds = 0.; /* stores weight of forbidden term for given motif */ + + is_unpaired = 0; + + vrna_mx_pf_t *matrices = vc->exp_matrices; + + my_iindx = vc->iindx; + jindx = vc->jindx; + + hc = vc->hc; + sc = vc->sc; + hc_up_ml = hc->up_ml; + + qm = matrices->qm; + qm1 = matrices->qm1; + expMLbase = matrices->expMLbase; + + turn = vc->exp_params->model_details.min_loop_size; + +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) { + memorized_node_prev = NULL; + memorized_node_cur = (*current_node)->head; + } + +#endif + + if (j > i) { + /* now backtrack [i ... j] in qm[] */ + + if (current_node) { + fbd = NR_TOTAL_WEIGHT(*current_node) * + qm[my_iindx[i] - j] / + (*q_remain); + } + + r = vrna_urn() * (qm[my_iindx[i] - j] - fbd); + if (current_node) { + fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_QM_UNPAIR, i, 0) * + qm[my_iindx[i] - j] / + (*q_remain); + qmt = qm1[jindx[j] + i] - fbds; + } else { + qmt = qm1[jindx[j] + i]; + } + + k = cnt = i; + q_temp = qm1[jindx[j] + i]; + + if (qmt < r) { +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) + advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_QM_UNPAIR, i, 0); + +#endif + + for (span = j - i, cnt = i + 1; cnt <= j; cnt++) { +#ifdef VRNA_WITH_BOUSTROPHEDON + k = (int)(i + 1 + span * ((cnt - i - 1) % 2)) + + (int)((1 - (2 * ((cnt - i - 1) % 2))) * ((cnt - i) / 2)); +#else + k = cnt; +#endif + q_temp = 0.; + u = k - i; + /* [i...k] is unpaired */ + if (hc_up_ml[i] >= u) { + q_temp += expMLbase[u] * qm1[jindx[j] + k]; + + if (sc) { + if (sc->exp_energy_up) + q_temp *= sc->exp_energy_up[i][u]; + + if (sc->exp_f) + q_temp *= sc->exp_f(i, j, k, j, VRNA_DECOMP_ML_ML, sc->data); + } + + if (current_node) { + fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_QM_UNPAIR, k, 0) * + qm[my_iindx[i] - j] / + (*q_remain); + qmt += q_temp - fbds; + } else { + qmt += q_temp; + } + } + + if (qmt >= r) { + /* we have chosen unpaired version */ + is_unpaired = 1; + break; + } + +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) + advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_QM_UNPAIR, k, 0); + +#endif + + /* split between k-1, k */ + q_temp = qm[my_iindx[i] - (k - 1)] * + qm1[jindx[j] + k]; + + if (sc) + if (sc->exp_f) + q_temp *= sc->exp_f(i, j, k - 1, k, VRNA_DECOMP_ML_ML_ML, sc->data); + + if (current_node) { + fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_QM_PAIR, k, 0) * + qm[my_iindx[i] - j] / + (*q_remain); + qmt += q_temp - fbds; + } else { + qmt += q_temp; + } + + if (qmt >= r) + break; + +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) + advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_QM_PAIR, k, 0); + +#endif + } + } else { + is_unpaired = 1; + } + + if (current_node) { + *q_remain *= q_temp / qm[my_iindx[i] - j]; +#ifdef VRNA_NR_SAMPLING_HASH + if (is_unpaired) + *current_node = add_if_nexists(NRT_QM_UNPAIR, k, 0, *current_node, *q_remain); + else + *current_node = add_if_nexists(NRT_QM_PAIR, k, 0, *current_node, *q_remain); + +#else + if (is_unpaired) + *current_node = add_if_nexists_ll(memory_dat, + NRT_QM_UNPAIR, + k, + 0, + memorized_node_prev, + memorized_node_cur, + *current_node, + *q_remain); + else + *current_node = add_if_nexists_ll(memory_dat, + NRT_QM_PAIR, + k, + 0, + memorized_node_prev, + memorized_node_cur, + *current_node, + *q_remain); + +#endif + } + + if (cnt > j) + return 0; + + ret = backtrack_qm1(k, j, pstruc, vc, q_remain, current_node, memory_dat); + + if (ret == 0) + return ret; + + if (k < i + turn) + return ret; /* no more pairs */ + + if (!is_unpaired) { + /* if we've chosen creating a branch in [i..k-1] */ + ret = backtrack_qm_nr(i, k - 1, pstruc, vc, q_remain, current_node, memory_dat); + + if (ret == 0) + return ret; + } + } + + return ret; +} + + +PRIVATE int backtrack_qm1(int i, int j, char *pstruc, - vrna_fold_compound_t *vc) + vrna_fold_compound_t *vc, + double *q_remain, + NR_NODE **current_node, + struct nr_memory **memory_dat) { /* i is paired to l, ilength; + fbd = 0.; + fbds = 0.; pf_params = vc->exp_params; my_iindx = vc->iindx; jindx = vc->jindx; - - ptype = vc->ptype; + ptype = vc->ptype; sc = vc->sc; hc = vc->hc; @@ -494,7 +1160,21 @@ backtrack_qm1(int i, turn = pf_params->model_details.min_loop_size; - r = vrna_urn() * qm1[jindx[j] + i]; +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) { + memorized_node_prev = NULL; + memorized_node_cur = (*current_node)->head; + } + +#endif + + if (current_node) { + fbd = NR_TOTAL_WEIGHT(*current_node) * + qm1[jindx[j] + i] / + (*q_remain); + } + + r = vrna_urn() * (qm1[jindx[j] + i] - fbd); ii = my_iindx[i]; for (qt = 0., l = j; l > i + turn; l--) { il = jindx[l] + i; @@ -514,19 +1194,56 @@ backtrack_qm1(int i, q_temp *= sc->exp_f(i, j, i, l, VRNA_DECOMP_ML_STEM, sc->data); } - qt += q_temp; - if (qt >= r) + if (current_node) { + fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_QM1_BRANCH, i, l) * + qm1[jindx[j] + i] / + (*q_remain); + qt += q_temp - fbds; + } else { + qt += q_temp; + } + + if (qt >= r) { + if (current_node) { + *q_remain *= q_temp / qm1[jindx[j] + i]; +#ifdef VRNA_NR_SAMPLING_HASH + *current_node = add_if_nexists(NRT_QM1_BRANCH, i, l, *current_node, *q_remain); +#else + *current_node = add_if_nexists_ll(memory_dat, + NRT_QM1_BRANCH, + i, + l, + memorized_node_prev, + memorized_node_cur, + *current_node, + *q_remain); +#endif + } + break; + } + +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) + advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_QM1_BRANCH, i, l); + +#endif } else { l = i + turn; break; } } } - if (l < i + turn + 1) - vrna_message_error("backtrack failed in qm1"); + if (l < i + turn + 1) { + if (current_node) { + return 0; + } else { + vrna_message_error("backtrack failed in qm1"); + return 0; + } + } - backtrack(i, l, pstruc, vc); + return backtrack(i, l, pstruc, vc, q_remain, current_node, memory_dat); } @@ -552,9 +1269,9 @@ backtrack_qm2(int k, /* we have to search for our barrier u between qm1 and qm1 */ if ((sc) && (sc->exp_f)) { for (qom2t = 0., u = k + turn + 1; u < n - turn - 1; u++) { - qom2t += qm1[jindx[u] + k] * - qm1[jindx[n] + (u + 1)] * - sc->exp_f(k, n, u, u + 1, VRNA_DECOMP_ML_ML_ML, sc->data); + qom2t += qm1[jindx[u] + k] * + qm1[jindx[n] + (u + 1)] * + sc->exp_f(k, n, u, u + 1, VRNA_DECOMP_ML_ML_ML, sc->data); if (qom2t > r) break; @@ -570,29 +1287,44 @@ backtrack_qm2(int k, if (u == n - turn) vrna_message_error("backtrack failed in qm2"); - backtrack_qm1(k, u, pstruc, vc); - backtrack_qm1(u + 1, n, pstruc, vc); + backtrack_qm1(k, u, pstruc, vc, NULL, NULL, NULL); + backtrack_qm1(u + 1, n, pstruc, vc, NULL, NULL, NULL); } -PRIVATE void +PRIVATE int backtrack(int i, int j, char *pstruc, - vrna_fold_compound_t *vc) + vrna_fold_compound_t *vc, + double *q_remain, + NR_NODE **current_node, + struct nr_memory **memory_dat) { char *ptype; unsigned char *hard_constraints, hc_decompose; vrna_exp_param_t *pf_params; FLT_OR_DBL *qb, *qm, *qm1, *scale; - FLT_OR_DBL r, qbt1, qt, q_temp; + FLT_OR_DBL r, fbd, fbds, qbt1, qbr, qt, q_temp; /* qbr stores qb used for generating r */ vrna_mx_pf_t *matrices; unsigned int n; - int *my_iindx, *jindx, *hc_up_int; + int *my_iindx, *jindx, *hc_up_int, ret; vrna_sc_t *sc; vrna_hc_t *hc; short *S1; +#ifndef VRNA_NR_SAMPLING_HASH + /* non-redundant data-structure memorization nodes */ + NR_NODE *memorized_node_prev = NULL; /* remembers previous-to-current node in linked list */ + NR_NODE *memorized_node_cur = NULL; /* remembers actual node in linked list */ +#endif + + ret = 1; /* default is success */ + fbd = 0.; /* stores weight of forbidden terms for given q[ij] */ + fbds = 0.; /* stores weight of forbidden term for given motif */ + qbt1 = 0.; + q_temp = 0.; + n = vc->length; pf_params = vc->exp_params; ptype = vc->ptype; @@ -611,10 +1343,18 @@ backtrack(int i, qm1 = matrices->qm1; scale = matrices->scale; - int turn = pf_params->model_details.min_loop_size; - int *rtype = &(pf_params->model_details.rtype[0]); + int turn = pf_params->model_details.min_loop_size; + int *rtype = &(pf_params->model_details.rtype[0]); + +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) { + memorized_node_prev = NULL; + memorized_node_cur = (*current_node)->head; + } + +#endif - qbt1 = 0.; + hc_decompose = hard_constraints[n * j + i]; do { int k, l, kl, u1, u2, max_k, min_l; @@ -622,18 +1362,61 @@ backtrack(int i, k = i; l = j; + qbr = qb[my_iindx[i] - j]; + + if (current_node) + fbd = NR_TOTAL_WEIGHT(*current_node) * qbr / (*q_remain); + pstruc[i - 1] = '('; pstruc[j - 1] = ')'; - r = vrna_urn() * qb[my_iindx[i] - j]; + r = vrna_urn() * (qb[my_iindx[i] - j] - fbd); + type = vrna_get_ptype(jindx[j] + i, ptype); + qbt1 = 0.; + + r = vrna_urn() * (qb[my_iindx[i] - j] - fbd); + qbr = qb[my_iindx[i] - j]; type = vrna_get_ptype(jindx[j] + i, ptype); hc_decompose = hard_constraints[n * i + j]; /* hairpin contribution */ - qbt1 = vrna_exp_E_hp_loop(vc, i, j); + q_temp = vrna_exp_E_hp_loop(vc, i, j); + + if (current_node) { + fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_HAIRPIN, 0, 0) * + qbr / + (*q_remain); + qbt1 += (q_temp - fbds); + } else { + qbt1 += q_temp; + } + + if (qbt1 >= r) { + /* found the hairpin we're done */ + if (current_node) { + *q_remain *= q_temp / qbr; +#ifdef VRNA_NR_SAMPLING_HASH + *current_node = add_if_nexists(NRT_HAIRPIN, 0, 0, *current_node, *q_remain); +#else + *current_node = add_if_nexists_ll(memory_dat, + NRT_HAIRPIN, + 0, + 0, + memorized_node_prev, + memorized_node_cur, + *current_node, + *q_remain); +#endif + } + + return ret; + } - if (qbt1 >= r) - return; /* found the hairpin we're done */ +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) + advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_HAIRPIN, 0, 0); + +#endif if (hc_decompose & VRNA_CONSTRAINT_CONTEXT_INT_LOOP) { /* interior loop contributions */ @@ -685,17 +1468,46 @@ backtrack(int i, q_temp *= sc->exp_f(i, j, k, l, VRNA_DECOMP_PAIR_IL, sc->data); } - qbt1 += q_temp; + if (current_node) { + fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_IT_LOOP, k, l) * + qbr / + (*q_remain); + qbt1 += q_temp - fbds; + } else { + qbt1 += q_temp; + } + if (qbt1 >= r) break; + +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) + advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_IT_LOOP, k, l); + +#endif } } if (qbt1 >= r) break; } if (k <= max_k) { - i = k; - j = l; + if (current_node) { + *q_remain *= q_temp / qbr; +#ifdef VRNA_NR_SAMPLING_HASH + *current_node = add_if_nexists(NRT_IT_LOOP, k, l, *current_node, *q_remain); +#else + *current_node = add_if_nexists_ll(memory_dat, + NRT_IT_LOOP, + k, + l, + memorized_node_prev, + memorized_node_cur, + *current_node, + *q_remain); +#endif + } + + return backtrack(k, l, pstruc, vc, q_remain, current_node, memory_dat); /* found the interior loop, repeat for inside */ } else { /* interior loop contributions did not exceed threshold, so we break */ break; @@ -707,7 +1519,7 @@ backtrack(int i, } while (1); /* backtrack in multi-loop */ - { + if (hard_constraints[n * j + i] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) { int k, ii, jj, tt; FLT_OR_DBL closingPair; tt = rtype[vrna_get_ptype(jindx[j] + i, ptype)]; @@ -716,7 +1528,7 @@ backtrack(int i, * scale[2]; if (sc) { if (sc->exp_energy_bp) - q_temp *= sc->exp_energy_bp[jindx[j] + i]; + closingPair *= sc->exp_energy_bp[jindx[j] + i]; if (sc->exp_f) closingPair *= sc->exp_f(i, j, i, j, VRNA_DECOMP_PAIR_ML, sc->data); @@ -730,38 +1542,101 @@ backtrack(int i, if ((sc) && (sc->exp_f)) { for (qt = qbt1, k = i + 1; k < j; k++) { - q_temp = qm[ii - (k - 1)] * - qm1[jj + k] * - closingPair * - sc->exp_f(i, j, k - 1, k, VRNA_DECOMP_ML_ML_ML, sc->data); - - qt += q_temp; - qbt1 += q_temp; + q_temp = qm[ii - (k - 1)] * + qm1[jj + k] * + closingPair * + sc->exp_f(i, + j, + k - 1, + k, + VRNA_DECOMP_ML_ML_ML, + sc->data); + + + if (current_node) { + fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_MT_LOOP, k, 0) * + qbr / + (*q_remain); + qt += q_temp - fbds; + qbt1 += q_temp - fbds; + } else { + qt += q_temp; + qbt1 += q_temp; + } if (qt >= r) break; + +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) + advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_MT_LOOP, k, 0); + +#endif } } else { for (qt = qbt1, k = i + 1; k < j; k++) { - q_temp = qm[ii - (k - 1)] * - qm1[jj + k] * - closingPair; - - qt += q_temp; - qbt1 += q_temp; + q_temp = qm[ii - (k - 1)] * + qm1[jj + k] * + closingPair; + + if (current_node) { + fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_MT_LOOP, k, 0) * + qbr / + (*q_remain); + qt += q_temp - fbds; + qbt1 += q_temp - fbds; + } else { + qt += q_temp; + qbt1 += q_temp; + } if (qt >= r) break; + +#ifndef VRNA_NR_SAMPLING_HASH + if (current_node) + advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_MT_LOOP, k, 0); + +#endif } } - if (k >= j) - vrna_message_error("backtrack failed, can't find split index "); - backtrack_qm1(k, j, pstruc, vc); + if (k >= j) { + if (current_node) + return 0; /* backtrack failed for non-redundant mode most likely due to numerical instabilities */ + else + vrna_message_error("backtrack failed, can't find split index "); + } + + if (current_node) { + *q_remain *= q_temp / qbr; +#ifdef VRNA_NR_SAMPLING_HASH + *current_node = add_if_nexists(NRT_MT_LOOP, k, 0, *current_node, *q_remain); +#else + *current_node = add_if_nexists_ll(memory_dat, + NRT_MT_LOOP, + k, + 0, + memorized_node_prev, + memorized_node_cur, + *current_node, + *q_remain); +#endif + } + + ret = backtrack_qm1(k, j, pstruc, vc, q_remain, current_node, memory_dat); + + if (ret == 0) + return ret; j = k - 1; - backtrack_qm(i, j, pstruc, vc); + + ret = (current_node) ? + backtrack_qm_nr(i, j, pstruc, vc, q_remain, current_node, memory_dat) : + backtrack_qm(i, j, pstruc, vc); } + + return ret; } @@ -822,7 +1697,7 @@ wrap_pbacktrack_circ(vrna_fold_compound_t *vc) for (i = 0; i < n; i++) pstruc[i] = '.'; - qt = 1.0 * scale[n]; + qt = 1.0 * scale[n]; /* add soft constraints for open chain configuration */ if (sc) { if (sc->exp_energy_up) @@ -832,7 +1707,7 @@ wrap_pbacktrack_circ(vrna_fold_compound_t *vc) qt *= sc->exp_f(1, n, 1, n, VRNA_DECOMP_EXT_UP, sc->data); } - r = vrna_urn() * qo; + r = vrna_urn() * qo; /* open chain? */ if (qt > r) @@ -863,7 +1738,7 @@ wrap_pbacktrack_circ(vrna_fold_compound_t *vc) /* found a hairpin? so backtrack in the enclosed part and we're done */ if (qt > r) { - backtrack(i, j, pstruc, vc); + backtrack(i, j, pstruc, vc, NULL, NULL, NULL); return pstruc; } @@ -921,12 +1796,12 @@ wrap_pbacktrack_circ(vrna_fold_compound_t *vc) q_temp *= sc->exp_f(i, j, k, l, VRNA_DECOMP_PAIR_IL, sc->data); } - qt += q_temp; + qt += q_temp; /* found an exterior interior loop? also this time, we can go straight */ /* forward and backtracking the both enclosed parts and we're done */ if (qt > r) { - backtrack(i, j, pstruc, vc); - backtrack(k, l, pstruc, vc); + backtrack(i, j, pstruc, vc, NULL, NULL, NULL); + backtrack(k, l, pstruc, vc, NULL, NULL, NULL); return pstruc; } } @@ -942,7 +1817,13 @@ wrap_pbacktrack_circ(vrna_fold_compound_t *vc) qt += qm[my_iindx[1] - k] * qm2[k + 1] * expMLclosing * - sc->exp_f(1, n, k, k + 1, VRNA_DECOMP_ML_ML_ML, sc->data); + sc->exp_f(1, + n, + k, + k + 1, + VRNA_DECOMP_ML_ML_ML, + sc->data); + /* backtrack in qm and qm2 if we've found a valid barrier k */ if (qt > r) { diff --git a/src/ViennaRNA/boltzmann_sampling.h b/src/ViennaRNA/boltzmann_sampling.h index 42025974c..f101de169 100644 --- a/src/ViennaRNA/boltzmann_sampling.h +++ b/src/ViennaRNA/boltzmann_sampling.h @@ -18,6 +18,9 @@ * @brief Functions to draw random structure samples from the ensemble according to their * equilibrium probability */ +typedef void (vrna_boltzmann_sampling_callback)(const char *stucture, + void *data); + /** * @brief Sample a secondary structure of a subsequence from the Boltzmann ensemble according its probability @@ -35,6 +38,33 @@ char *vrna_pbacktrack5(vrna_fold_compound_t *vc, int length); +/** + * @brief Samples multiple secondary structures non-redundantly from the Boltzmann ensemble according its probability + * + * @ingroup subopt_stochbt + * @pre The fold compound has to be obtained using the #VRNA_OPTION_HYBRID option in vrna_fold_compound() + * @pre vrna_pf() has to be called first to fill the partition function matrices + * + * @note In some cases, this function does not return the number of requested samples but a smaller number. + * This may happen if a) the number of requested structures is larger than the total number of structures + * in the ensemble, or b) numeric instabilities prevent the backtracking function to enumerate structures + * with very high free energies. + * + * @param vc The fold compound data structure + * @param num_samples The number of desired non-redundant samples + * @return A list of sampled secondary structures in dot-bracket notation, terminated by @em NULL + */ +char **vrna_pbacktrack_nr(vrna_fold_compound_t *vc, + int num_samples); + + +void +vrna_pbacktrack_nr_cb(vrna_fold_compound_t *vc, + int num_samples, + vrna_boltzmann_sampling_callback *cb, + void *data); + + /** * @brief Sample a secondary structure (consensus structure) from the Boltzmann ensemble according its probability * diff --git a/src/ViennaRNA/data_structures_nonred.inc b/src/ViennaRNA/data_structures_nonred.inc new file mode 100644 index 000000000..6f546c4f8 --- /dev/null +++ b/src/ViennaRNA/data_structures_nonred.inc @@ -0,0 +1,1120 @@ +#ifndef DATA_STRUCTURES_NONRED_H_ +#define DATA_STRUCTURES_NONRED_H_ + +#ifdef VRNA_NR_SAMPLING_MPFR +#include +#endif + +#define DEBUG 0 + +/* General mpfr funtions */ + +#ifdef VRNA_NR_SAMPLING_MPFR +PRIVATE int precision(); // returns precision +PRIVATE mpfr_rnd_t default_rnd(); // returns default rounding mode (to nearest) +#endif + +/* defines types for different sequence stetches/loops */ +typedef enum _type Type; +enum stetch_type { + NRT_NONE_TYPE = 0, /* nonetype: both 0 (root) */ + NRT_HAIRPIN = 1, /* hairpin: both 0 (unused) */ + NRT_IT_LOOP = 2, /* internal loop: (i', j') - closed pair */ + NRT_MT_LOOP = 3, /* multiloop: K (second is 0 - unused) */ + NRT_EXT_LOOP = 4, /* external loop: (i, j) */ + NRT_UNPAIRED_SG = 5, /* unpaired positions: (first unpaired position - end known of stretch) */ + NRT_QM1_BRANCH = 6, /* branch in multiloop in qm1: (i,j) */ + NRT_QM_PAIR = 7, /* QM split on K with pairing within [i..k-1] : (k, 0) */ + NRT_QM_UNPAIR = 8, /* QM split on K with no pairing within [i..k-1] : (k, 0) */ +}; + +#ifdef VRNA_NR_SAMPLING_HASH +/******************************************/ +/* version with hash table */ +/******************************************/ + +/* explanation: each node (tr_node) holds a pointer to a start of a linked + * list to the next node, link to previous node, Boltzmann's factor + * for given flat structure and itentifier, which corresponds to an id of + * chosen flat structure + */ + +/* define node of a tree */ +typedef struct tr_node tr_node; + + +/* define elements of hashtable */ +typedef struct hash_node hash_node; /* list of arrays sorting triplets for every hashvalue*/ + +struct tr_node { + int type; + int loop_spec_1; + int loop_spec_2; + int seqlen; + tr_node *parent; + tr_node *child; /* used when this node has at most one child */ + hash_node *hash_tab; /* used when there is at least two children */ + tr_node *next_in_hash; /* used for node that is next in hash table */ + int hash_value; /* hash value computed for this node */ +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_t weight; + mpfr_t max_weight; /* maximum allowed weight (maximum of partition function) */ +#else + double weight; + double max_weight; /* maximum allowed weight (maximum of partition function) */ +#endif + int created_recently; /* 1 if was created during last iteration, otherwise 0 */ +}; + +struct hash_node { + tr_node **array; + int hash_elemts; + int hash_size; +}; + +/* tree + hash functions */ +/** @brief creates a root of datastructure tree (hash table version) **/ +PRIVATE tr_node *create_root(int seqlen, + double max_weight); + + +/** @brief returns a weigh of node (type, loop_spec_1, loop_spec_2) if child of last_node, otherwise returns 0.0 **/ +PRIVATE double tr_node_weight(tr_node *last_node, + int type, + int loop_spec_1, + int loop_spec_2); + + +/** @brief sums weight of all children of par_node and returns it **/ +PRIVATE double total_weight_par(tr_node *par_node); + + +/** @brief sums weight of all children of par_node with certain type and returns it **/ +PRIVATE double total_weight_par_type(int type, + tr_node *par_node); + + +/** @brief creates node (type, loop_spec_1, loop_spec_2) if not existing and returns pointer to it, + * or returns pointer to exisiting case **/ +PRIVATE tr_node *add_if_nexists(int type, + int loop_spec_1, + int loop_spec_2, + tr_node *par_node, + double max_weight); + + + +/** @brief traces back from leaf to root while updating weights of leaf to all nodes in path, + * returns pointer to root **/ +PRIVATE tr_node *traceback_to_root(tr_node *leaf, + double struct_weight, + int *is_dup, + int *pf_overflow); + + +/** @brief destructor **/ +PRIVATE void free_all_nr(tr_node *root); + + +#else +/******************************************/ +/* version with linked lists */ +/******************************************/ + +typedef struct tllr_node tllr_node; + +struct tllr_node { + int type; + int loop_spec_1; + int loop_spec_2; + tllr_node *parent; /* vertical chaining - ancestor */ + tllr_node *head; /* vertical chaining - successor */ + tllr_node *next_node; /*horizontal chaining - linked list */ +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_t weight; + mpfr_t max_weight; /* maximum allowed weight (maximum of partition function) */ +#else + double weight; + double max_weight; /* maximum allowed weight (maximum of partition function) */ +#endif + int created_recently; /* 1 if was created during last iteration, otherwise 0 */ +}; + + +/* memory object for non-redundant sampling approach using linked lists*/ +typedef struct nr_memory nr_memory; + +struct nr_memory { + tllr_node *nr_memory_allocated; + int memory_index; + size_t tllr_node_size; + size_t block_size; /* block size */ + nr_memory *prev_block; /* stores previous block */ +}; + +/* creates an object nr_memory that pre-allocates a block of memory for memory allocation */ +PRIVATE nr_memory * create_nr_memory(size_t node_size, size_t block_size, nr_memory * prev_memory); + + +/* tree + linked list functions */ +/** @brief creates a root of datastructure tree (linked list version) **/ +PRIVATE tllr_node *create_ll_root(struct nr_memory **memory_dat, + double max_weight); + + +/** resets cursor to current_node and start of linked list **/ +PRIVATE void reset_cursor(tllr_node **memorized_node_prev, + tllr_node **memorized_node_cur, + tllr_node *current_node); + + +/** @brief moves cursor to next node if current_node is identical to one in loop, otherwise does nothing **/ +PRIVATE void advance_cursor(tllr_node **memorized_node_prev, + tllr_node **memorized_node_cur, + int type, + int loop_spec_1, + int loop_spec_2); + + +/** @brief returns a weigh of node (type, loop_spec_1, loop_spec_2) if child of last_node, otherwise returns 0.0 **/ +PRIVATE double get_weight(tllr_node *memorized_node_cur, + int type, + int loop_spec_1, + int loop_spec_2); + + +/** @brief sums weight of all children of par_node and returns it **/ +PRIVATE double get_weight_all(tllr_node *par_node); + + +/** @brief sums weight of all children of par_node with certain type and returns it **/ +PRIVATE double get_weight_type_spec(int type, + tllr_node *par_node); + + +/** @brief creates node (type, loop_spec_1, loop_spec_2) if not existing and returns pointer to it, + * or returns pointer to exisiting case **/ +PRIVATE tllr_node *add_if_nexists_ll(struct nr_memory **memory_dat, + int type, + int loop_spec_1, + int loop_spec_2, + tllr_node *memorized_node_prev, + tllr_node *memorized_node_cur, + tllr_node *parent_node, + double max_weight); + + +/** @brief traces back from leaf to root while updating weights of leaf to all nodes in path, + * returns pointer to root **/ +PRIVATE tllr_node *traceback_to_ll_root(tllr_node *leaf, + double weight, + int *is_dup, + int *pf_overflow); + + +/** @brief destructor **/ +PRIVATE void free_all_nrll(nr_memory **memory_dat); + +#endif + + +/**************************************************/ +/* mpfr related functions */ +/**************************************************/ + +#ifdef VRNA_NR_SAMPLING_MPFR +PRIVATE int precision(){ // returns precision + return 128; +} + + +PRIVATE mpfr_rnd_t default_rnd(){ // returns default rounding mode (to nearest) + return mpfr_get_default_rounding_mode(); +} +#endif + +/**************************************************/ +/* node weight getter */ +/**************************************************/ + +#ifdef VRNA_NR_SAMPLING_HASH +PRIVATE double +return_node_weight(tr_node * node){ +#ifdef VRNA_NR_SAMPLING_MPFR + return mpfr_get_d(node->weight, default_rnd()); +#else + return node->weight; +#endif +} +#else +PRIVATE double +return_node_weight(tllr_node * node){ +#ifdef VRNA_NR_SAMPLING_MPFR + return mpfr_get_d(node->weight, default_rnd()); +#else + return node->weight; +#endif +} +#endif + +/**************************************************/ +/* data structure (tree + hash table) */ +/**************************************************/ +/* Contains + * + * + * + * Creates a structure memorizing the weights of already used solutions for the + * non-redundant sampling. It uses tree structure where each node contains a hash + * table. + */ + + +#ifdef VRNA_NR_SAMPLING_HASH + + +PRIVATE int +pow_of_2_modulus(int numerator, + int denominator) +{ + return numerator & (denominator - 1); +} + + +/* creates new tree node with undefined weight with parent node */ +PRIVATE tr_node * +create_tr_node(int type, + int loop_spec_1, + int loop_spec_2, + int seqlen, + tr_node *parent, + double max_weight) +{ + tr_node *new_tr_node = (tr_node *)vrna_alloc(sizeof(tr_node)); + new_tr_node->type = type; + /* Types and properties specific to loops: + * type 0 : nonetype: both 0 (root) + * type 1 : hairpin: both 0 (unused) + * type 2 : internal loop: (i', j') - closed pair + * type 3 : multiloop: K (second is 0 - unused) + * type 4 : external loop: (i, j) + * type 5 : unpaired positions: (first unpaired position - end known of stretch) + * type 6 : special type qm1: (i, j) (used in multiloops) + * type 7 : branch in multiloop in qm1: (i,j) + */ + new_tr_node->loop_spec_1 = loop_spec_1; + new_tr_node->loop_spec_2 = loop_spec_2; + new_tr_node->seqlen = seqlen; + new_tr_node->parent = parent; + new_tr_node->child = NULL; + new_tr_node->hash_tab = NULL; + new_tr_node->next_in_hash = NULL; + new_tr_node->hash_value = 0; +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_init2(new_tr_node->weight, precision()); + mpfr_set_d(new_tr_node->weight, 0., default_rnd()); + mpfr_init2(new_tr_node->max_weight, precision()); + mpfr_set_d(new_tr_node->max_weight, max_weight, default_rnd()); +#else + new_tr_node->weight = 0; + new_tr_node->max_weight = max_weight; +#endif + new_tr_node->created_recently = 1; + return new_tr_node; +} + + +/* compares hash children values with actual value in parent */ +PRIVATE void +compare_parent_children_weight(tr_node *parent){ + tr_node *node_t; + int i; + +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_t total; + mpfr_init2(total, precision()); + mpfr_set_d(total, 0., default_rnd()); +#else + double total = 0.; +#endif + + if(!parent->hash_tab){ + if(parent->child){ + +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_add(total, total, parent->child->weight, default_rnd()); +#else + total += parent->child->weight; +#endif + } + } + else{ + for(i = 0; ihash_tab->hash_size; i++){ + node_t = parent->hash_tab->array[i]; + while(node_t){ +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_add(total, total, node_t->weight, default_rnd()); +#else + total += node_t->weight; +#endif + node_t = node_t->next_in_hash; + } + } + } +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_clear(total); +#endif +} + + +/* creates hash entry point (its start) */ +PRIVATE hash_node * +init_hash() +{ + int i; + hash_node *new_hash = (hash_node *)vrna_alloc(sizeof(hash_node)); + + new_hash->hash_elemts = 0; + new_hash->hash_size = 4; /* will be expanded when necessary */ + new_hash->array = (tr_node **)vrna_alloc(sizeof(tr_node *) * (new_hash->hash_size)); + for (i = 0; i < new_hash->hash_size; i++) + new_hash->array[i] = NULL; + return new_hash; +} + + +/* creates root (start of a tree) */ +PRIVATE tr_node * +create_root(int seqlen, + double max_weight) +{ + tr_node *root = create_tr_node(NRT_NONE_TYPE, 0, 0, seqlen, NULL, max_weight); + + return root; +} + + +/* returns hash key for given triple of values */ +PRIVATE int +hash_fci(int type, + int loop_spec_1, + int loop_spec_2, + int seqlen) +{ + int hash_f = 3 * type + 11 * loop_spec_1 + 50021 * loop_spec_2; + + //int hash_f = type + 9*loop_spec_1 + 9*(seqlen+1)*loop_spec_2; + return hash_f; +} + + +/* Expands hash when its load becomes higher than 50% number and relinks elements. + * New size is lowest prime number higher than twice the size. + * */ +PRIVATE void +extend_hash(hash_node *hash_to_res) +{ + /* extends table */ + tr_node *next; + tr_node *tmp; + int hash_val, hash_val_mod; + int hash_size_tmp = hash_to_res->hash_size; + + hash_to_res->hash_size *= 2; + hash_to_res->array = (tr_node **)realloc(hash_to_res->array, sizeof(tr_node *) * (hash_to_res->hash_size)); + for (int i = hash_size_tmp; i < hash_to_res->hash_size; i++) + hash_to_res->array[i] = NULL; + /* relinks all elements within table + * this table points to beginning of each chain within hash table */ + tr_node **hash_tmp = (tr_node **)vrna_alloc(sizeof(tr_node *) * (hash_to_res->hash_size)); + for (int i = 0; i < hash_to_res->hash_size; i++) + hash_tmp[i] = NULL; + + /* this table points at the current end of eaxh hash table */ + tr_node **hash_cur_end = (tr_node **)vrna_alloc(sizeof(tr_node *) * (hash_to_res->hash_size)); + for (int i = 0; i < hash_to_res->hash_size; i++) + hash_cur_end[i] = NULL; + for (int i = 0; i < hash_size_tmp; i++) { + next = hash_to_res->array[i]; + while (next) { + hash_val = hash_fci(next->type, next->loop_spec_1, next->loop_spec_2, next->seqlen); + hash_val_mod = pow_of_2_modulus(hash_val, hash_to_res->hash_size); + if (!hash_tmp[hash_val_mod]) { + hash_tmp[hash_val_mod] = next; + hash_cur_end[hash_val_mod] = next; + } else { + hash_cur_end[hash_val_mod]->next_in_hash = next; + hash_cur_end[hash_val_mod] = next; + } + + tmp = next; + next = next->next_in_hash; + tmp->next_in_hash = NULL; + } + } + for (int i = 0; i < hash_to_res->hash_size; i++) + hash_to_res->array[i] = hash_tmp[i]; + free(hash_tmp); + free(hash_cur_end); +} + + +/* inserting new value into corresponding hash */ +PRIVATE void +insert_val(hash_node *hash_t, + int type, + int loop_spec_1, + int loop_spec_2, + tr_node *inserted_node) +{ + int hash_val_mod; + + tr_node *ptr; + + /* where to insert */ + inserted_node->hash_value = hash_fci(type, loop_spec_1, loop_spec_2, inserted_node->seqlen); + hash_val_mod = pow_of_2_modulus(inserted_node->hash_value, hash_t->hash_size); + + + if (!hash_t->array[hash_val_mod]) { + /* head case */ + hash_t->array[hash_val_mod] = inserted_node; + return; + } + ptr = hash_t->array[hash_val_mod]; + while (ptr->next_in_hash) /* inside list */ + ptr = ptr->next_in_hash; + ptr->next_in_hash = inserted_node; +} + + +/* creates and inserts a new tr_node (along with creating hash) into tree */ +PRIVATE tr_node * +insert_tr_node(int type, + int loop_spec_1, + int loop_spec_2, + tr_node *last_node, + double max_weight) +{ + /* forge new tr_node */ + tr_node *new_tr_node = + create_tr_node(type, loop_spec_1, loop_spec_2, last_node->seqlen, last_node, max_weight); + + if (!last_node->child) { + /* empty so no child yet - insert it directly */ + last_node->child = new_tr_node; + } else { + if (!last_node->hash_tab) { + /* hash not defined yet (we have exactly two values) */ + last_node->hash_tab = init_hash(last_node->seqlen); + insert_val(last_node->hash_tab, + last_node->child->type, + last_node->child->loop_spec_1, + last_node->child->loop_spec_2, + last_node->child); + last_node->hash_tab->hash_elemts++; + } + + /* insert it into hashtable of previous node */ + insert_val(last_node->hash_tab, type, loop_spec_1, loop_spec_2, new_tr_node); + last_node->hash_tab->hash_elemts++; + if (((double)last_node->hash_tab->hash_elemts / last_node->hash_tab->hash_size) > 0.5) + extend_hash(last_node->hash_tab); + } + + return new_tr_node; +} + + +/* tracks whether there is node in hash with given values */ +PRIVATE tr_node * +access_val(tr_node *parent, + int type, + int loop_spec_1, + int loop_spec_2) +{ + tr_node *ptr = NULL; + + if (!parent->hash_tab) { + /* hash does not exist (at most one child) */ + if (parent->child) { + /* child exists (exactly one child) */ + if ((parent->child->type == type) && (parent->child->loop_spec_1 == loop_spec_1) + && (parent->child->loop_spec_2 == loop_spec_2)) + return parent->child; + } + } else { + /* hash exists -> 2+ children */ + int hash_f = hash_fci(type, loop_spec_1, loop_spec_2, parent->seqlen); + int hash_f_mod = pow_of_2_modulus(hash_f, parent->hash_tab->hash_size); + ptr = parent->hash_tab->array[hash_f_mod]; + while (ptr) { + if (ptr->hash_value == hash_f) + break; + + ptr = ptr->next_in_hash; + } + } + + if (ptr) + return ptr; + + return NULL; +} + + +/* creates node with properties that don't already exist for given node */ +PRIVATE tr_node * +add_if_nexists(int type, + int loop_spec_1, + int loop_spec_2, + tr_node *last_node, + double max_weight) +{ + tr_node *new_tr_node = access_val(last_node, type, loop_spec_1, loop_spec_2); + + if (!new_tr_node) + new_tr_node = insert_tr_node(type, loop_spec_1, loop_spec_2, last_node, max_weight); + + return new_tr_node; +} + + +/* returns weight of node if it exists; otherwise returns 0 */ +PRIVATE double +tr_node_weight(tr_node *last_node, + int type, + int loop_spec_1, + int loop_spec_2) +{ + tr_node *next_node = access_val(last_node, type, loop_spec_1, loop_spec_2); + + if (next_node){ +#ifdef VRNA_NR_SAMPLING_MPFR + return mpfr_get_d(next_node->weight, default_rnd()); +#else + return next_node->weight; +#endif + } + return 0.; +} + + +/* returns sum of all weights for given parent */ +PRIVATE double +total_weight_par(tr_node *last_node) +{ + if ((!last_node->hash_tab) && (!last_node->child)) + return 0.; +#ifdef VRNA_NR_SAMPLING_MPFR + return mpfr_get_d(last_node->weight, default_rnd()); +#else + return last_node->weight; +#endif +} + + +/* returns sum of all weights for given parent and for given type*/ +PRIVATE double +total_weight_par_type(int type, + tr_node *last_node) +{ + tr_node *ptr = NULL; + int i; + +#ifdef VRNA_NR_SAMPLING_MPFR + double weight_sum_d; + mpfr_t weight_sum; + mpfr_init2(weight_sum, precision()); + mpfr_set_d(weight_sum, 0., default_rnd()); +#else + double weight_sum = 0.; +#endif + + + if (!last_node->hash_tab) { + if (last_node->child) + if (last_node->child->type == type){ +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_clear(weight_sum); + return mpfr_get_d(last_node->child->weight, default_rnd()); +#else + return last_node->child->weight; +#endif + } + } else { + for (i = 0; i < last_node->hash_tab->hash_size; i++) { + ptr = last_node->hash_tab->array[i]; + while (ptr) { + if (ptr->type == type){ +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_add(weight_sum, weight_sum, ptr->weight, default_rnd()); +#else + weight_sum += ptr->weight; +#endif + } + ptr = ptr->next_in_hash; + } + } + } + +#ifdef VRNA_NR_SAMPLING_MPFR + weight_sum_d = mpfr_get_d(weight_sum, default_rnd()); + mpfr_clear(weight_sum); + return weight_sum_d; +#else + return weight_sum; +#endif +} + + +/* updates weight of a given node */ +PRIVATE int +update_weight(tr_node *node, + double weight) +{ +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_t intermediate; + mpfr_init2(intermediate, precision()); + mpfr_add_d(intermediate, node->weight, weight, default_rnd()); + mpfr_sub(intermediate, node->max_weight, intermediate, default_rnd()); + + if(mpfr_cmp_d(intermediate, -1E-14)<0){ + mpfr_clear(intermediate); + return 1; + }else{ + mpfr_clear(intermediate); + mpfr_add_d(node->weight, node->weight, weight, default_rnd()); + } +#else + if(node->max_weight - (node->weight + weight) < -(1E-14)){ + return 1; + }else{ + node->weight += weight; + } +#endif + return 0; +} + +/* tracebacks to root while updating values for each node passed through + * - also verifies unicity (at least one node differs) */ +PRIVATE tr_node * +traceback_to_root(tr_node *leaf, + double struct_weight, + int *is_dup, + int *pf_overflow) +{ + *pf_overflow = update_weight(leaf, struct_weight); + if(leaf->created_recently){ + leaf->created_recently = 0; + *is_dup = 0; + } + while (leaf->parent) { + *pf_overflow = update_weight(leaf->parent, struct_weight); + if(leaf->parent->created_recently){ + leaf->parent->created_recently = 0; + *is_dup = 0; + } + leaf = leaf->parent; + } + return leaf; +} + + +/* destructor */ +PRIVATE void +free_all_nr(tr_node *root) +{ + if (!root->hash_tab) { + if (root->child){ + free_all_nr(root->child); + } + } else { + for (int i = 0; i < root->hash_tab->hash_size; i++) { + tr_node *ptr = root->hash_tab->array[i]; + tr_node *tmp; + while (ptr) { + tmp = ptr->next_in_hash; + free_all_nr(ptr); + ptr = tmp; + } + } + free(root->hash_tab->array); + free(root->hash_tab); + } +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_clear(root->weight); + mpfr_clear(root->max_weight); +#endif + free(root); +} +#endif + + +/*********************************************************/ +/* data structure (linked_list + hash table) */ +/*********************************************************/ + +#ifndef VRNA_NR_SAMPLING_HASH +/* allocates a nr_memory object - pre-allocator for tllr_nodes */ +PRIVATE nr_memory * +create_nr_memory(size_t tllr_node_size, size_t block_size, nr_memory * prev_block) +{ + struct nr_memory * memory_dat = vrna_alloc(sizeof(nr_memory)); + memory_dat->nr_memory_allocated = vrna_alloc(block_size); + memory_dat->memory_index = 0; + memory_dat->tllr_node_size = tllr_node_size; + memory_dat->block_size = block_size; + memory_dat->prev_block = prev_block; + + return memory_dat; +} + + +/* This creates structure that uses linked list instead of hash. The thought behind this is + * the order of investigated nodes is always the same so we can add them to specific place. + * It is thus a bit faster. + */ +PRIVATE tllr_node * +create_tllr_node(struct nr_memory **memory_dat, + int type, + int loop_spec_1, + int loop_spec_2, + tllr_node *parent, + double max_weight) +{ + tllr_node *new_tllr_node; + if((*memory_dat)->tllr_node_size * ((*memory_dat)->memory_index + 1) > (*memory_dat)->block_size){ + struct nr_memory * memory_dat_tmp = create_nr_memory((*memory_dat)->tllr_node_size, (*memory_dat)->block_size, *memory_dat); + *memory_dat = memory_dat_tmp; + new_tllr_node = &((*memory_dat)->nr_memory_allocated[(*memory_dat)->memory_index]); + }else{ + new_tllr_node = &((*memory_dat)->nr_memory_allocated[(*memory_dat)->memory_index]); + } + new_tllr_node->type = type; + /* Types and properties specific to loops: + * type 0 : nonetype: both 0 (root) + * type 1 : hairpin: both 0 (unused) + * type 2 : internal loop: (i', j') - closed pair + * type 3 : multiloop: K (second is 0 - unused) + * type 4 : external loop: (i, j) + * type 5 : unpaired positions: (first unpaired position - end known of stretch) + * type 6 : special type qm1: (i, j) (used in multiloops) + * type 7 : branch in multiloop in qm1: (i,j) + */ + new_tllr_node->loop_spec_1 = loop_spec_1; + new_tllr_node->loop_spec_2 = loop_spec_2; + new_tllr_node->parent = parent; + new_tllr_node->next_node = NULL; + new_tllr_node->head = NULL; +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_init2(new_tllr_node->weight, precision()); + mpfr_set_d(new_tllr_node->weight, 0., default_rnd()); + mpfr_init2(new_tllr_node->max_weight, precision()); + mpfr_set_d(new_tllr_node->max_weight, max_weight, default_rnd()); +#else + new_tllr_node->weight = 0; + new_tllr_node->max_weight = max_weight; +#endif + new_tllr_node->created_recently = 1; + + (*memory_dat)->memory_index++; + return new_tllr_node; +} + + +/* compares hash children values with actual value in parent */ +#if DEBUG +PRIVATE void +compare_parent_children_weight_tr(tllr_node *parent){ + tllr_node *node_t; + +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_t total; + mpfr_init2(total, precision()); + mpfr_set_d(total, 0., default_rnd()); +#else + double total = 0.; +#endif + + node_t = parent->head; + while(node_t){ +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_add(total, total, node_t->weight, default_rnd()); +#else + total += node_t->weight; +#endif + node_t = node_t->next_node; + } +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_clear(total); +#endif +} +#endif + +/* creates root (start of a tree) */ +PRIVATE tllr_node * +create_ll_root(struct nr_memory **memory_dat, + double max_weight) +{ + tllr_node *root = create_tllr_node(memory_dat, NRT_NONE_TYPE, 0, 0, NULL, max_weight); + + return root; +} + + +/* inserts a tllr_node before 'next_node' and after previous node + * Node cursor hols previous and current ll_node */ +PRIVATE tllr_node * +insert_tllr_node(struct nr_memory **memory_dat, + tllr_node *memorized_node_prev, + tllr_node *memorized_node_cur, + int type, + int loop_spec_1, + int loop_spec_2, + tllr_node *parent_node, + double max_weight) +{ + tllr_node *new_node = create_tllr_node(memory_dat, type, loop_spec_1, loop_spec_2, parent_node, max_weight); + + if (!memorized_node_prev) /* first node to be inserted */ + parent_node->head = new_node; + else + memorized_node_prev->next_node = new_node; + + new_node->next_node = memorized_node_cur; + return new_node; +} + + +/* resets cursor to beginning of loop*/ +PRIVATE void +reset_cursor(tllr_node **memorized_node_prev, + tllr_node **memorized_node_cur, + tllr_node *current_node) +{ + (*memorized_node_prev) = NULL; + (*memorized_node_cur) = current_node->head; +} + + +/* advances pointer in loop if the identifier coincide with current pointer and returns weight */ +PRIVATE inline void +advance_cursor(tllr_node **memorized_node_prev, + tllr_node **memorized_node_cur, + int type, + int loop_spec_1, + int loop_spec_2) +{ + if (*memorized_node_cur) { + if ((*memorized_node_cur)->type == type + && (*memorized_node_cur)->loop_spec_1 == loop_spec_1 + && (*memorized_node_cur)->loop_spec_2 == loop_spec_2) { + (*memorized_node_prev) = (*memorized_node_cur); + (*memorized_node_cur) = (*memorized_node_cur)->next_node; + } + } +} + + +/* gets weight of actual node */ +PRIVATE inline double +get_weight(tllr_node *memorized_node_cur, + int type, + int loop_spec_1, + int loop_spec_2) +{ + double weight = 0; + + if (memorized_node_cur) { + if (memorized_node_cur->type == type + && memorized_node_cur->loop_spec_1 == loop_spec_1 + && memorized_node_cur->loop_spec_2 == loop_spec_2) +#ifdef VRNA_NR_SAMPLING_MPFR + weight = mpfr_get_d(memorized_node_cur->weight, default_rnd()); +#else + weight = memorized_node_cur->weight; +#endif + } + + return weight; +} + + +/* get weight of all child nodes */ +PRIVATE double +get_weight_all(tllr_node *last_node) +{ + if (!last_node->head) + return 0; + +#ifdef VRNA_NR_SAMPLING_MPFR + return mpfr_get_d(last_node->weight, default_rnd()); +#else + return last_node->weight; +#endif +} + + +/* get weight of all child nodes of certain type */ +PRIVATE double +get_weight_type_spec(int type, + tllr_node *last_node) +{ + //double weight_total = 0; + +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_t weight_total; + double weight_total_d; + mpfr_init2(weight_total, precision()); + mpfr_set_d(weight_total, 0., default_rnd()); +#else + double weight_total = 0.; +#endif + + tllr_node *ptr = last_node->head; + + while (ptr) { + if (ptr->type == type){ +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_add(weight_total, weight_total, ptr->weight, default_rnd()); +#else + weight_total += ptr->weight; +#endif + } + + ptr = ptr->next_node; + } + +#ifdef VRNA_NR_SAMPLING_MPFR + weight_total_d = mpfr_get_d(weight_total, default_rnd()); + mpfr_clear(weight_total); + return weight_total_d; +#else + return weight_total; +#endif +} + + +/* adds node if the current one isn't the one we want */ +PRIVATE inline tllr_node * +add_if_nexists_ll(struct nr_memory **memory_dat, + int type, + int loop_spec_1, + int loop_spec_2, + tllr_node *memorized_node_prev, + tllr_node *memorized_node_cur, + tllr_node *parent_node, + double max_weight) +{ + tllr_node *returned_node; + + if (memorized_node_cur) { + if (memorized_node_cur->type == type + && memorized_node_cur->loop_spec_1 == loop_spec_1 + && memorized_node_cur->loop_spec_2 == loop_spec_2) + returned_node = memorized_node_cur; + else + returned_node = insert_tllr_node(memory_dat, + memorized_node_prev, + memorized_node_cur, + type, + loop_spec_1, + loop_spec_2, + parent_node, + max_weight); + } else { + returned_node = insert_tllr_node(memory_dat, + memorized_node_prev, + memorized_node_cur, + type, + loop_spec_1, + loop_spec_2, + parent_node, + max_weight); + } + + return returned_node; +} + + +/* updates weight of a given node */ +PRIVATE int +update_weight_ll(tllr_node *node, + double weight) +{ + +#ifdef VRNA_NR_SAMPLING_MPFR + mpfr_t intermediate; + mpfr_init2(intermediate, precision()); + mpfr_add_d(intermediate, node->weight, weight, default_rnd()); + mpfr_sub(intermediate, node->max_weight, intermediate, default_rnd()); + + //if(node->max_weight - (node->weight + weight) < -(1E-14)){ + if(mpfr_cmp_d(intermediate, -1E-14)<0){ + mpfr_clear(intermediate); + return 1; + }else{ + mpfr_clear(intermediate); + mpfr_add_d(node->weight, node->weight, weight, default_rnd()); + } +#else + if(node->max_weight - (node->weight + weight) < -(1E-14)){ + return 1; + }else{ + node->weight += weight; + } +#endif + return 0; +} + + +/* tracebacks to root while updating values for each node passed through + * - also verifies unicity (at least one node differs) */ +PRIVATE tllr_node * +traceback_to_ll_root(tllr_node *leaf, + double weight, + int *is_dup, + int *pf_overflow) +{ + *pf_overflow = update_weight_ll(leaf, weight); + if(leaf->created_recently){ /* check whether the last sequence is not a duplicate */ + leaf->created_recently = 0; + *is_dup = 0; + } + while (leaf->parent){ + *pf_overflow = update_weight_ll(leaf->parent, weight); + if(leaf->parent->created_recently){ + leaf->parent->created_recently = 0; + *is_dup = 0; + } + leaf = leaf->parent; + } + return leaf; +} + + +/* destructor */ +PRIVATE void +free_all_nrll(struct nr_memory **memory_dat) +{ + int i; + nr_memory *memory_block = *memory_dat; + nr_memory *memory_block_next; + while(memory_block){ + memory_block_next = memory_block->prev_block; +#ifdef VRNA_NR_SAMPLING_MPFR + for(i = 0; i < memory_block->memory_index; i++){ + mpfr_clear(memory_block->nr_memory_allocated[i].weight); + mpfr_clear(memory_block->nr_memory_allocated[i].max_weight); + } +#endif + free(memory_block->nr_memory_allocated); + free(memory_block); + memory_block = memory_block_next; + } +} +#endif + + +#endif diff --git a/src/ViennaRNA/datastructures/char_stream.c b/src/ViennaRNA/datastructures/char_stream.c index 9338a3af7..533a7358c 100644 --- a/src/ViennaRNA/datastructures/char_stream.c +++ b/src/ViennaRNA/datastructures/char_stream.c @@ -69,6 +69,7 @@ void vrna_cstr_free(struct vrna_cstr_s *buf) { if (buf) { + vrna_cstr_fflush(buf); free(buf->string); free(buf); } @@ -79,9 +80,11 @@ void vrna_cstr_close(struct vrna_cstr_s *buf) { if (buf) { + vrna_cstr_fflush(buf); + free(buf->string); - if (buf->output != stdout) + if ((buf->output != stdout) && (buf->output != stderr)) fclose(buf->output); free(buf); diff --git a/src/ViennaRNA/datastructures/char_stream.h b/src/ViennaRNA/datastructures/char_stream.h index 1af87259f..e106412db 100644 --- a/src/ViennaRNA/datastructures/char_stream.h +++ b/src/ViennaRNA/datastructures/char_stream.h @@ -18,19 +18,61 @@ /* below is our own implementation of a dynamic char * stream */ typedef struct vrna_cstr_s *vrna_cstr_t; +/** + * @brief Create a dynamic char * stream data structure + * + * @see vrna_cstr_free(), vrna_cstr_close(), vrna_cstr_fflush(), vrna_cstr_printf() + * + * @param size The initial size of the buffer in characters + * @param output An optional output file stream handle that is used to write the collected data to (defaults to @em stdout if @em NULL) + */ vrna_cstr_t vrna_cstr(size_t size, FILE *output); +/** + * @brief Free the memory occupied by a dynamic char * stream data structure + * + * This function first flushes any remaining character data within the stream + * and then free's the memory occupied by the data structure. + * + * @see vrna_cstr_close(), vrna_cstr_fflush(), vrna_cstr() + * + * @param buf The dynamic char * stream data structure to free + */ void vrna_cstr_free(vrna_cstr_t buf); +/** + * @brief Free the memory occupied by a dynamic char * stream and close the output stream + * + * This function first flushes any remaining character data within the stream + * then closes the attached output file stream (if any), and finally free's the + * memory occupied by the data structure. + * + * @see vrna_cstr_free(), vrna_cstr_fflush(), vrna_cstr() + * + * @param buf The dynamic char * stream data structure to free + */ void vrna_cstr_close(vrna_cstr_t buf); +/** + * @brief Flush the dynamic char * output stream + * + * This function flushes the collected char * stream, either by writing + * to the attached file handle, or simply by writing to @em stdout if + * no file handle has been attached upon construction using vrna_cstr(). + * + * @post The stream buffer is empty after execution of this function + * + * @see vrna_cstr(), vrna_cstr_close(), vrna_cstr_free() + * + * @param buf The dynamic char * stream data structure to flush + */ void vrna_cstr_fflush(struct vrna_cstr_s *buf); diff --git a/src/ViennaRNA/equilibrium_probs.c b/src/ViennaRNA/equilibrium_probs.c index e9101d14a..5a3366b46 100644 --- a/src/ViennaRNA/equilibrium_probs.c +++ b/src/ViennaRNA/equilibrium_probs.c @@ -334,6 +334,54 @@ vrna_mean_bp_distance(vrna_fold_compound_t *vc) } +PUBLIC double +vrna_ensemble_defect(vrna_fold_compound_t *fc, + const char *structure) +{ + unsigned int i, j, n; + int ii; + double ed = -1.; + + if ((fc) && + (structure) && + (strlen(structure) == fc->length) && + (fc->exp_matrices) && + (fc->exp_matrices->probs)) { + n = fc->length; + + short *pt = vrna_ptable(structure); + FLT_OR_DBL *probs = fc->exp_matrices->probs; + int *idx = fc->iindx; + ed = 0.; + + for (i = 1; i < n; i++) { + ii = idx[i]; + double pi; + + /* compute probability to be paired */ + for (pi = 0., j = 1; j < i; j++) + pi += probs[idx[j] - i]; + + for (j = i + 1; j <= n; j++) + pi += probs[ii - j]; + + if (pt[i] == 0) + ed += pi; + else if (pt[i] > i) + ed += 1 - probs[ii - pt[i]]; + else + ed += 1 - probs[idx[pt[i]] - i]; + } + + ed /= (double)n; + + free(pt); + } + + return ed; +} + + PUBLIC vrna_ep_t * vrna_stack_prob(vrna_fold_compound_t *vc, double cutoff) @@ -997,7 +1045,8 @@ compute_bpp_internal(vrna_fold_compound_t *fc, break; if (hc->mx[i * n + j] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP) { - type = vrna_get_ptype(jindx[j] + i, ptype); + int jij = jindx[j] + i; + type = vrna_get_ptype(jij, ptype); if ((sn[k] == sn[i]) && (sn[j] == sn[l])) { @@ -1019,7 +1068,7 @@ compute_bpp_internal(vrna_fold_compound_t *fc, * sc->exp_energy_up[l + 1][u2]; if (sc->exp_energy_bp) - tmp2 *= sc->exp_energy_bp[ij]; + tmp2 *= sc->exp_energy_bp[jij]; if (sc->exp_energy_stack) { if ((i + 1 == k) && (j - 1 == l)) { @@ -1328,7 +1377,7 @@ compute_bpp_multibranch(vrna_fold_compound_t *fc, if (sc) { if (sc->exp_energy_bp) - ppp *= sc->exp_energy_bp[ij]; + ppp *= sc->exp_energy_bp[jindx[j] + i]; /* * if(sc->exp_f) @@ -1355,7 +1404,7 @@ compute_bpp_multibranch(vrna_fold_compound_t *fc, if (sc) { /* which decompositions are covered here? => (i, l+1) -> enclosing pair */ if (sc->exp_energy_bp) - prmt1 *= sc->exp_energy_bp[ii - (l + 1)]; + prmt1 *= sc->exp_energy_bp[jindx[l + 1] + i]; /* * if(sc->exp_f) @@ -2846,7 +2895,7 @@ ud_outside_mb_loops(vrna_fold_compound_t *vc) if (sc) if (sc->exp_energy_bp) - qqq *= sc->exp_energy_bp[kl]; + qqq *= sc->exp_energy_bp[jkl]; temp += qqq; } @@ -2873,7 +2922,8 @@ ud_outside_mb_loops(vrna_fold_compound_t *vc) kl = my_iindx[k] - l; if ((hc[jindx[l] + k] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) && (probs[kl] > 0.) && (hc_up[k + 1] >= up)) { - tt = rtype[vrna_get_ptype(jindx[l] + k, ptype)]; + int jkl = jindx[l] + k; + tt = rtype[vrna_get_ptype(jkl, ptype)]; temp = probs[kl] * expMLbase[up] * exp_E_MLstem(tt, S[l - 1], S[k + 1], pf_params) @@ -2881,7 +2931,7 @@ ud_outside_mb_loops(vrna_fold_compound_t *vc) * scale[2]; if (sc) { if (sc->exp_energy_bp) - temp *= sc->exp_energy_bp[kl]; + temp *= sc->exp_energy_bp[jkl]; if (sc->exp_energy_up) temp *= sc->exp_energy_up[k + 1][up]; @@ -2998,8 +3048,9 @@ ud_outside_mb_loops(vrna_fold_compound_t *vc) for (l = j + 1; l <= n; l++) { kl = my_iindx[k] - l; if (hc[jindx[l] + k] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) { - int up; - tt = rtype[vrna_get_ptype(jindx[l] + k, ptype)]; + int up, jkl; + jkl = jindx[l] + k; + tt = rtype[vrna_get_ptype(jkl, ptype)]; up = l - j - 1; if (hc_up[j + 1] >= up) { temp = probs[kl] @@ -3010,7 +3061,7 @@ ud_outside_mb_loops(vrna_fold_compound_t *vc) if (sc) { if (sc->exp_energy_bp) - temp *= sc->exp_energy_bp[kl]; + temp *= sc->exp_energy_bp[jkl]; if (sc->exp_energy_up) temp *= sc->exp_energy_up[j + 1][up]; @@ -3156,7 +3207,7 @@ ud_outside_mb_loops2(vrna_fold_compound_t *vc) if (sc) if (sc->exp_energy_bp) - qqq *= sc->exp_energy_bp[kl]; + qqq *= sc->exp_energy_bp[jkl]; temp += qqq; } @@ -3183,7 +3234,8 @@ ud_outside_mb_loops2(vrna_fold_compound_t *vc) kl = my_iindx[k] - l; if ((hc[l * n + k] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) && (probs[kl] > 0.) && (hc_up[k + 1] >= up)) { - tt = rtype[vrna_get_ptype(jindx[l] + k, ptype)]; + int jkl = jindx[l] + k; + tt = rtype[vrna_get_ptype(jkl, ptype)]; temp = probs[kl] * expMLbase[up] * exp_E_MLstem(tt, S[l - 1], S[k + 1], pf_params) @@ -3191,7 +3243,7 @@ ud_outside_mb_loops2(vrna_fold_compound_t *vc) * scale[2]; if (sc) { if (sc->exp_energy_bp) - temp *= sc->exp_energy_bp[kl]; + temp *= sc->exp_energy_bp[jkl]; if (sc->exp_energy_up) temp *= sc->exp_energy_up[k + 1][up]; @@ -3304,8 +3356,9 @@ ud_outside_mb_loops2(vrna_fold_compound_t *vc) for (l = j + 1; l <= n; l++) { kl = my_iindx[k] - l; if (hc[k * n + l] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) { - int up; - tt = rtype[vrna_get_ptype(jindx[l] + k, ptype)]; + int up, jkl; + jkl = jindx[l] + k; + tt = rtype[vrna_get_ptype(jkl, ptype)]; up = l - j - 1; if (hc_up[j + 1] >= up) { temp = probs[kl] @@ -3316,7 +3369,7 @@ ud_outside_mb_loops2(vrna_fold_compound_t *vc) if (sc) { if (sc->exp_energy_bp) - temp *= sc->exp_energy_bp[kl]; + temp *= sc->exp_energy_bp[jkl]; if (sc->exp_energy_up) temp *= sc->exp_energy_up[j + 1][up]; diff --git a/src/ViennaRNA/equilibrium_probs.h b/src/ViennaRNA/equilibrium_probs.h index 05b3a4b06..a1e8b4c1c 100644 --- a/src/ViennaRNA/equilibrium_probs.h +++ b/src/ViennaRNA/equilibrium_probs.h @@ -58,6 +58,35 @@ double vrna_mean_bp_distance_pr(int length, FLT_OR_DBL *pr); */ double vrna_mean_bp_distance(vrna_fold_compound_t *vc); +/** + * @brief Compute the Ensemble Defect for a given target structure + * + * Given a target structure @f$s@f$, compute the average dissimilarity of a randomly + * drawn structure from the ensemble, i.e.: + * @f[ + * ED(s) = 1 - \frac{1}{n} \sum_{ij, (i,j) \in s} p_{ij} - \frac{1}{n} \sum_{i}(1 - s_i)q_i + * @f] + * with sequence length @f$n@f$, the probability @f$p_{ij}@f$ of a base pair @f$(i,j)@f$, + * the probability @f$q_i = 1 - \sum_j p_{ij}@f$ of nucleotide @f$i@f$ being unpaired, and + * the indicator variable @f$s_i = 1@f$ if @f$\exists (i,j) \in s@f$, and @f$s_i = 0@f$ otherwise. + * + * @pre The #vrna_fold_compound_t input parameter @p fc must contain a valid base pair + * probability matrix. This means that partition function and base pair probabilities + * must have been computed using @p fc before execution of this function! + * + * @ingroup part_func_global + * + * @see vrna_pf(), vrna_pairing_probs() + * + * @param fc A fold_compound with pre-computed base pair probabilities + * @param structure A target structure in dot-bracket notation + * @return The ensemble defect with respect to the target structure, or -1. upon failure, e.g. pre-conditions are not met + */ +double +vrna_ensemble_defect(vrna_fold_compound_t *fc, + const char *structure); + + /** * @brief Compute stacking probabilities * diff --git a/src/ViennaRNA/eval_wrappers.c b/src/ViennaRNA/eval_wrappers.c index 68b593dcb..3deb93fef 100644 --- a/src/ViennaRNA/eval_wrappers.c +++ b/src/ViennaRNA/eval_wrappers.c @@ -407,24 +407,18 @@ vrna_eval_move_shift_pt(vrna_fold_compound_t *fc, int d1 = -structure[unchangedPosition]; int d2 = -unchangedPosition; vrna_move_t deletion; - if (d1 < d2) { - deletion.pos_5 = d2; - deletion.pos_3 = d1; - } else { - deletion.pos_5 = d1; - deletion.pos_3 = d2; - } + if (d1 < d2) + deletion = vrna_move_init(d2, d1); + else + deletion = vrna_move_init(d1, d2); int i1 = unchangedPosition; int i2 = insertedPosition; vrna_move_t insertion; - if (i1 > i2) { - insertion.pos_5 = i2; - insertion.pos_3 = i1; - } else { - insertion.pos_5 = i1; - insertion.pos_3 = i2; - } + if (i1 > i2) + insertion = vrna_move_init(i2, i1); + else + insertion = vrna_move_init(i1, i2); int energy = vrna_eval_move_pt(fc, structure, deletion.pos_5, deletion.pos_3); short *tmpS = vrna_ptable_copy(structure); diff --git a/src/ViennaRNA/fold_compound.c b/src/ViennaRNA/fold_compound.c index 1541a1f15..f895b36fe 100644 --- a/src/ViennaRNA/fold_compound.c +++ b/src/ViennaRNA/fold_compound.c @@ -770,10 +770,13 @@ make_pscores(vrna_fold_compound_t *fc) if (S[s][i] == 0 && S[s][j] == 0) { type = 7; /* gap-gap */ } else { - if ((AS[s][i] == '~') || (AS[s][j] == '~')) + if ((AS[s][i] == '~') || (AS[s][j] == '~')) { type = 7; - else + } else { type = md->pair[S[s][i]][S[s][j]]; + if ((md->noGU) && ((type == 3) || (type == 4))) + type = 0; + } } pfreq[type]++; diff --git a/src/ViennaRNA/grammar.c b/src/ViennaRNA/grammar.c index 7b7736544..c6c300f4f 100644 --- a/src/ViennaRNA/grammar.c +++ b/src/ViennaRNA/grammar.c @@ -167,7 +167,7 @@ vrna_gr_set_aux_exp_m1(vrna_fold_compound_t *fc, PUBLIC int vrna_gr_set_aux(vrna_fold_compound_t *fc, - vrna_callback_gr_rule *cb) + vrna_callback_gr_rule_aux *cb) { int ret = 0; @@ -186,7 +186,7 @@ vrna_gr_set_aux(vrna_fold_compound_t *fc, PUBLIC int vrna_gr_set_aux_exp(vrna_fold_compound_t *fc, - vrna_callback_gr_rule_exp *cb) + vrna_callback_gr_rule_aux_exp *cb) { int ret = 0; diff --git a/src/ViennaRNA/grammar.h b/src/ViennaRNA/grammar.h index 7517c3cb6..91a4a9298 100644 --- a/src/ViennaRNA/grammar.h +++ b/src/ViennaRNA/grammar.h @@ -21,12 +21,24 @@ typedef int (vrna_callback_gr_rule)(vrna_fold_compound_t *vc, void *data); +typedef void (vrna_callback_gr_rule_aux)(vrna_fold_compound_t *vc, + int i, + int j, + void *data); + + typedef FLT_OR_DBL (vrna_callback_gr_rule_exp)(vrna_fold_compound_t *vc, int i, int j, void *data); +typedef void (vrna_callback_gr_rule_aux_exp)(vrna_fold_compound_t *vc, + int i, + int j, + void *data); + + typedef void (vrna_callback_gr_cond)(vrna_fold_compound_t *fc, unsigned char stage, void *data); @@ -45,13 +57,13 @@ struct vrna_gr_aux_s { vrna_callback_gr_rule *cb_aux_c; vrna_callback_gr_rule *cb_aux_m; vrna_callback_gr_rule *cb_aux_m1; - vrna_callback_gr_rule *cb_aux; + vrna_callback_gr_rule_aux *cb_aux; vrna_callback_gr_rule_exp *cb_aux_exp_f; vrna_callback_gr_rule_exp *cb_aux_exp_c; - vrna_callback_gr_rule_exp *cb_aux_exp_m; - vrna_callback_gr_rule_exp *cb_aux_exp_m1; - vrna_callback_gr_rule_exp *cb_aux_exp; + vrna_callback_gr_rule_exp *cb_aux_exp_m; + vrna_callback_gr_rule_exp *cb_aux_exp_m1; + vrna_callback_gr_rule_aux_exp *cb_aux_exp; void *data; vrna_callback_gr_free_data *free_data; @@ -100,12 +112,12 @@ vrna_gr_set_aux_exp_m1(vrna_fold_compound_t *fc, int vrna_gr_set_aux(vrna_fold_compound_t *fc, - vrna_callback_gr_rule *cb); + vrna_callback_gr_rule_aux *cb); int vrna_gr_set_aux_exp(vrna_fold_compound_t *fc, - vrna_callback_gr_rule_exp *cb); + vrna_callback_gr_rule_aux_exp *cb); int diff --git a/src/ViennaRNA/loops/external.c b/src/ViennaRNA/loops/external.c index abea7362c..167cc4231 100644 --- a/src/ViennaRNA/loops/external.c +++ b/src/ViennaRNA/loops/external.c @@ -19,6 +19,7 @@ #include "ViennaRNA/structured_domains.h" #include "ViennaRNA/unstructured_domains.h" #include "ViennaRNA/loops/external.h" +#include "ViennaRNA/utils/higher_order_functions.h" #ifdef __GNUC__ # define INLINE inline @@ -2110,72 +2111,20 @@ add_f3_gquad(vrna_fold_compound_t *fc, } -#ifdef VRNA_WITH_SSE_IMPLEMENTATION -/* SSE modular decomposition -------------------------------*/ -#include -#include - -PRIVATE INLINE int -horizontal_min_Vec4i(__m128i x) -{ - __m128i min1 = _mm_shuffle_epi32(x, _MM_SHUFFLE(0, 0, 3, 2)); - __m128i min2 = _mm_min_epi32(x, min1); - __m128i min3 = _mm_shuffle_epi32(min2, _MM_SHUFFLE(0, 0, 0, 1)); - __m128i min4 = _mm_min_epi32(min2, min3); - - return _mm_cvtsi128_si32(min4); -} - - -#endif - - PRIVATE INLINE int decompose_f5_ext_stem(vrna_fold_compound_t *fc, int j, int *stems) { - int e, i, *f5, turn; + int e, *f5, turn; f5 = fc->matrices->f5; turn = fc->params->model_details.min_loop_size; e = INF; - /* modular decomposition */ -#if VRNA_WITH_SSE_IMPLEMENTATION - __m128i inf = _mm_set1_epi32(INF); - - const int end = j - turn; - - for (i = 2; i < end - 3; i += 4) { - __m128i a = _mm_loadu_si128((__m128i *)&f5[i - 1]); - __m128i b = _mm_loadu_si128((__m128i *)&stems[i]); - __m128i c = _mm_add_epi32(a, b); - __m128i mask1 = _mm_cmplt_epi32(a, inf); - __m128i mask2 = _mm_cmplt_epi32(b, inf); - __m128i res = _mm_or_si128(_mm_and_si128(mask1, c), - _mm_andnot_si128(mask1, a)); - - res = _mm_or_si128(_mm_and_si128(mask2, res), - _mm_andnot_si128(mask2, b)); - const int en = horizontal_min_Vec4i(res); - e = MIN2(e, en); - } - - for (; i < end; i++) { - if ((f5[i - 1] != INF) && (stems[i] != INF)) { - const int en = f5[i - 1] + stems[i]; - e = MIN2(e, en); - } - } -#else - for (i = 2; i < j - turn; i++) - if ((f5[i - 1] != INF) && (stems[i] != INF)) { - const int en = f5[i - 1] + stems[i]; - e = MIN2(e, en); - } + const int count = j - turn; -#endif + e = vrna_fun_zip_add_min(f5 + 1, stems + 2, count - 2); return e; } @@ -2187,45 +2136,15 @@ decompose_f3_ext_stem(vrna_fold_compound_t *fc, int max_j, int *stems) { - int e, j, *f3, turn; + int e, *f3, turn, count; f3 = fc->matrices->f3_local; turn = fc->params->model_details.min_loop_size; + count = max_j - i - turn; e = INF; /* modular decomposition */ -#if VRNA_WITH_SSE_IMPLEMENTATION - __m128i inf = _mm_set1_epi32(INF); - - for (j = i + turn + 1; j < max_j - 3; j += 4) { - __m128i a = _mm_loadu_si128((__m128i *)&f3[j + 1]); - __m128i b = _mm_loadu_si128((__m128i *)&stems[j]); - __m128i c = _mm_add_epi32(a, b); - __m128i mask1 = _mm_cmplt_epi32(a, inf); - __m128i mask2 = _mm_cmplt_epi32(b, inf); - __m128i res = _mm_or_si128(_mm_and_si128(mask1, c), - _mm_andnot_si128(mask1, a)); - - res = _mm_or_si128(_mm_and_si128(mask2, res), - _mm_andnot_si128(mask2, b)); - const int en = horizontal_min_Vec4i(res); - e = MIN2(e, en); - } - - for (; j <= max_j; j++) { - if ((f3[j + 1] != INF) && (stems[j] != INF)) { - const int en = f3[j + 1] + stems[j]; - e = MIN2(e, en); - } - } -#else - for (j = i + turn + 1; j <= max_j; j++) - if ((f3[j + 1] != INF) && (stems[j] != INF)) { - const int en = stems[j] + f3[j + 1]; - e = MIN2(e, en); - } - -#endif + e = vrna_fun_zip_add_min(stems + i + turn + 1, f3 + i + turn + 2, count); return e; } diff --git a/src/ViennaRNA/loops/external.h b/src/ViennaRNA/loops/external.h index f9193bfb6..6de43a9aa 100644 --- a/src/ViennaRNA/loops/external.h +++ b/src/ViennaRNA/loops/external.h @@ -252,6 +252,10 @@ vrna_BT_ext_loop_f3_pp(vrna_fold_compound_t *fc, * @see E_ExtLoop() * @note This function is threadsafe * + * @deprecated Please use one of the functions vrna_E_ext_stem() and + * E_MLstem() instead! Use the former for cases where @p extLoop != 0 + * and the latter otherwise. + * * @param type The pair type of the first base pair un the stem * @param si1 The 5'-mismatching nucleotide * @param sj1 The 3'-mismatching nucleotide @@ -265,7 +269,7 @@ DEPRECATED(int E_Stem(int type, int sj1, int extLoop, vrna_param_t *P), - "This function is obsolete"); + "This function is obsolete. Use vrna_E_ext_stem() or E_MLstem() instead"); DEPRECATED(int E_ExtLoop(int type, diff --git a/src/ViennaRNA/loops/external_pf.c b/src/ViennaRNA/loops/external_pf.c index e17d56d0b..5838ce0e0 100644 --- a/src/ViennaRNA/loops/external_pf.c +++ b/src/ViennaRNA/loops/external_pf.c @@ -638,6 +638,10 @@ exp_E_ext_fast(vrna_fold_compound_t *fc, qbt1 += split_ext_fast(fc, i, j, aux_mx, evaluate, &hc_dat_local, &sc_wrapper); + /* apply auxiliary grammar rule for exterior loop case */ + if ((fc->aux_grammar) && (fc->aux_grammar->cb_aux_exp_f)) + qbt1 += fc->aux_grammar->cb_aux_exp_f(fc, i, j, fc->aux_grammar->data); + free_sc_wrapper_pf(&sc_wrapper); return qbt1; diff --git a/src/ViennaRNA/loops/hairpin_pf.c b/src/ViennaRNA/loops/hairpin_pf.c index 92a05c020..8dd064cbe 100644 --- a/src/ViennaRNA/loops/hairpin_pf.c +++ b/src/ViennaRNA/loops/hairpin_pf.c @@ -98,7 +98,7 @@ exp_eval_hp_loop_fake(vrna_fold_compound_t *fc, { short *S, *S2, s5, s3; unsigned int *sn, *ss, *se; - int u, type, *iidx; + int u, type, *iidx, *jidx; FLT_OR_DBL qq, temp, *q, *scale; vrna_exp_param_t *pf_params; vrna_sc_t *sc; @@ -106,6 +106,7 @@ exp_eval_hp_loop_fake(vrna_fold_compound_t *fc, vrna_ud_t *domains_up; iidx = fc->iindx; + jidx = fc->jindx; pf_params = fc->exp_params; md = &(pf_params->model_details); q = fc->exp_matrices->q; @@ -151,7 +152,7 @@ exp_eval_hp_loop_fake(vrna_fold_compound_t *fc, qq *= sc->exp_energy_up[i + 1][u]; if (sc->exp_energy_bp) - qq *= sc->exp_energy_bp[iidx[i] - j]; + qq *= sc->exp_energy_bp[jidx[j] + i]; if (sc->exp_f) qq *= sc->exp_f(i, j, i, j, VRNA_DECOMP_PAIR_HP, sc->data); diff --git a/src/ViennaRNA/loops/multibranch.c b/src/ViennaRNA/loops/multibranch.c index 6469c763d..eecb60fa6 100644 --- a/src/ViennaRNA/loops/multibranch.c +++ b/src/ViennaRNA/loops/multibranch.c @@ -1,3 +1,6 @@ +/*WBL 24 Aug 2018 Add AVX512 based on sources_034_578/modular_decomposition_id3.c */ +/*WBL 22 Aug 2018 by hand d3c17fd3e04e2419c147a1e097d3c4d2c5a6f11d lines 1355-1357*/ + #ifdef HAVE_CONFIG_H #include "config.h" #endif @@ -19,6 +22,7 @@ #include "ViennaRNA/structured_domains.h" #include "ViennaRNA/unstructured_domains.h" #include "ViennaRNA/loops/multibranch.h" +#include "ViennaRNA/utils/higher_order_functions.h" #ifdef __GNUC__ # define INLINE inline @@ -168,6 +172,11 @@ E_ml_rightmost_stem(int i, e = extend_fm_3p(i, j, fc->matrices->fM1, fc, evaluate, &hc_dat_local, &sc_wrapper); + if ((fc->aux_grammar) && (fc->aux_grammar->cb_aux_m1)) { + int ee = fc->aux_grammar->cb_aux_m1(fc, i, j, fc->aux_grammar->data); + e = MIN2(e, ee); + } + free_sc_wrapper(&sc_wrapper); } @@ -1021,29 +1030,6 @@ extend_fm_3p(int i, } -#ifdef VRNA_WITH_SSE_IMPLEMENTATION -#include -#include - -/* - * SSE minimum - * see also: http://stackoverflow.com/questions/9877700/getting-max-value-in-a-m128i-vector-with-sse - */ -int -horizontal_min_Vec4i(__m128i x) -{ - __m128i min1 = _mm_shuffle_epi32(x, _MM_SHUFFLE(0, 0, 3, 2)); - __m128i min2 = _mm_min_epi32(x, min1); - __m128i min3 = _mm_shuffle_epi32(min2, _MM_SHUFFLE(0, 0, 0, 1)); - __m128i min4 = _mm_min_epi32(min2, min3); - - return _mm_cvtsi128_si32(min4); -} - - -#endif - - PRIVATE int E_ml_stems_fast(vrna_fold_compound_t *fc, int i, @@ -1325,47 +1311,14 @@ E_ml_stems_fast(vrna_fold_compound_t *fc, if (last_nt < i) last_nt = i; /* do not start before i */ - const int stop = last_nt; -#ifdef VRNA_WITH_SSE_IMPLEMENTATION - const int end = 1 + stop - k; - int cnt; - __m128i inf = _mm_set1_epi32(INF); - - for (cnt = 0; cnt < end - 3; cnt += 4) { - __m128i a = _mm_loadu_si128((__m128i *)&fmi_tmp[k + cnt]); - __m128i b = _mm_loadu_si128((__m128i *)&fm[k1j + cnt]); - __m128i c = _mm_add_epi32(a, b); - __m128i mask1 = _mm_cmplt_epi32(a, inf); - __m128i mask2 = _mm_cmplt_epi32(b, inf); - __m128i res = _mm_or_si128(_mm_and_si128(mask1, c), - _mm_andnot_si128(mask1, a)); - - res = _mm_or_si128(_mm_and_si128(mask2, res), - _mm_andnot_si128(mask2, b)); - const int en = horizontal_min_Vec4i(res); - decomp = MIN2(decomp, en); - } + const int count = last_nt - k + 1; - for (; cnt < end; cnt++) { - if ((fmi[k + cnt] != INF) && (fm[k1j + cnt] != INF)) { - const int en = fmi[k + cnt] + fm[k1j + cnt]; - decomp = MIN2(decomp, en); - } - } + en = vrna_fun_zip_add_min(fmi_tmp + k, fm + k1j, count); + decomp = MIN2(decomp, en); - k += cnt; - k1j += cnt; -#else - for (; k <= stop; k++, k1j++) { - if ((fmi_tmp[k] != INF) && (fm[k1j] != INF)) { - en = fmi_tmp[k] + fm[k1j]; - decomp = MIN2(decomp, en); - } - } -#endif - - k++; - k1j++; + /* advance counters by processed subsegment and add 1 for the split point between strands */ + k += count + 1; + k1j += count + 1; if (k > j - turn - 2) break; @@ -1483,6 +1436,11 @@ E_ml_stems_fast(vrna_fold_compound_t *fc, e = MIN2(e, decomp); } + if ((fc->aux_grammar) && (fc->aux_grammar->cb_aux_m)) { + en = fc->aux_grammar->cb_aux_m(fc, i, j, fc->aux_grammar->data); + e = MIN2(e, en); + } + fmi[j] = e; free_sc_wrapper(&sc_wrapper); diff --git a/src/ViennaRNA/loops/multibranch_pf.c b/src/ViennaRNA/loops/multibranch_pf.c index d5cc1ba60..fb4d7ac57 100644 --- a/src/ViennaRNA/loops/multibranch_pf.c +++ b/src/ViennaRNA/loops/multibranch_pf.c @@ -647,6 +647,10 @@ exp_E_ml_fast(vrna_fold_compound_t *fc, free(qqm_tmp); } + /* apply auxiliary grammar rule for multibranch loop case */ + if ((fc->aux_grammar) && (fc->aux_grammar->cb_aux_exp_m)) + temp += fc->aux_grammar->cb_aux_exp_m(fc, i, j, fc->aux_grammar->data); + free_sc_wrapper(&sc_wrapper); return temp + qqm[i]; diff --git a/src/ViennaRNA/mfe.c b/src/ViennaRNA/mfe.c index 484af8895..3ab7da119 100644 --- a/src/ViennaRNA/mfe.c +++ b/src/ViennaRNA/mfe.c @@ -228,7 +228,7 @@ vrna_backtrack_from_intervals(vrna_fold_compound_t *fc, PRIVATE int fill_arrays(vrna_fold_compound_t *fc) { - int i, j, ij, length, turn, uniq_ML, e, *indx, *f5, *c, *fML, *fM1; + int i, j, ij, length, turn, uniq_ML, *indx, *f5, *c, *fML, *fM1; vrna_param_t *P; vrna_mx_mfe_t *matrices; vrna_ud_t *domains_up; @@ -281,26 +281,14 @@ fill_arrays(vrna_fold_compound_t *fc) c[ij] = decompose_pair(fc, i, j, helper_arrays); /* decompose subsegment [i, j] that is multibranch loop part with at least one branch */ - e = vrna_E_ml_stems_fast(fc, i, j, helper_arrays->Fmi, helper_arrays->DMLi); + fML[ij] = vrna_E_ml_stems_fast(fc, i, j, helper_arrays->Fmi, helper_arrays->DMLi); - if ((fc->aux_grammar) && (fc->aux_grammar->cb_aux_m)) { - int ee = fc->aux_grammar->cb_aux_m(fc, i, j, fc->aux_grammar->data); - e = MIN2(e, ee); - } - - fML[ij] = e; - - if (uniq_ML) { - /* decompose subsegment [i, j] that is multibranch loop part with exactly one branch */ - e = E_ml_rightmost_stem(i, j, fc); - - if ((fc->aux_grammar) && (fc->aux_grammar->cb_aux_m1)) { - int ee = fc->aux_grammar->cb_aux_m1(fc, i, j, fc->aux_grammar->data); - e = MIN2(e, ee); - } + /* decompose subsegment [i, j] that is multibranch loop part with exactly one branch */ + if (uniq_ML) + fM1[ij] = E_ml_rightmost_stem(i, j, fc); - fM1[ij] = e; - } + if ((fc->aux_grammar) && (fc->aux_grammar->cb_aux)) + fc->aux_grammar->cb_aux(fc, i, j, fc->aux_grammar->data); } /* end of j-loop */ rotate_aux_arrays(helper_arrays, length); diff --git a/src/ViennaRNA/pair_mat.h b/src/ViennaRNA/pair_mat.h index 33234bc23..8569366b1 100644 --- a/src/ViennaRNA/pair_mat.h +++ b/src/ViennaRNA/pair_mat.h @@ -50,6 +50,8 @@ encode_char(char c) /* return numerical representation of base used e.g. in pair[][] */ int code; + c = toupper(c); + if (energy_set > 0) { code = (int)(c - 'A') + 1; } else { @@ -164,14 +166,14 @@ encode_sequence(const char *sequence, /* standard encoding as always used for S */ case 0: for (i = 1; i <= l; i++) /* make numerical encoding of sequence */ - S[i] = (short)encode_char(toupper(sequence[i - 1])); + S[i] = (short)encode_char(sequence[i - 1]); S[l + 1] = S[1]; S[0] = (short)l; break; /* encoding for mismatches of nostandard bases (normally used for S1) */ case 1: for (i = 1; i <= l; i++) - S[i] = alias[(short)encode_char(toupper(sequence[i - 1]))]; + S[i] = alias[(short)encode_char(sequence[i - 1])]; S[l + 1] = S[1]; S[0] = S[l]; break; diff --git a/src/ViennaRNA/part_func.c b/src/ViennaRNA/part_func.c index 90687bb82..c4fa7b58d 100644 --- a/src/ViennaRNA/part_func.c +++ b/src/ViennaRNA/part_func.c @@ -413,13 +413,7 @@ fill_arrays(vrna_fold_compound_t *fc) qb[ij] = decompose_pair(fc, i, j, aux_mx_ml); /* Multibranch loop */ - temp = vrna_exp_E_ml_fast(fc, i, j, aux_mx_ml); - - /* apply auxiliary grammar rule for multibranch loop case */ - if ((fc->aux_grammar) && (fc->aux_grammar->cb_aux_exp_m)) - temp += fc->aux_grammar->cb_aux_exp_m(fc, i, j, fc->aux_grammar->data); - - qm[ij] = temp; + qm[ij] = vrna_exp_E_ml_fast(fc, i, j, aux_mx_ml); if (qm1) { temp = vrna_exp_E_ml_fast_qqm(aux_mx_ml)[i]; /* for stochastic backtracking and circfold */ @@ -432,21 +426,19 @@ fill_arrays(vrna_fold_compound_t *fc) } /* Exterior loop */ - temp = vrna_exp_E_ext_fast(fc, i, j, aux_mx_el); - - /* apply auxiliary grammar rule for exterior loop case */ - if ((fc->aux_grammar) && (fc->aux_grammar->cb_aux_exp_f)) - temp += fc->aux_grammar->cb_aux_exp_f(fc, i, j, fc->aux_grammar->data); + q[ij] = vrna_exp_E_ext_fast(fc, i, j, aux_mx_el); - q[ij] = temp; + /* apply auxiliary grammar rule (storage takes place in user-defined data structure */ + if ((fc->aux_grammar) && (fc->aux_grammar->cb_aux_exp)) + fc->aux_grammar->cb_aux_exp(fc, i, j, fc->aux_grammar->data); - if (temp > Qmax) { - Qmax = temp; + if (q[ij] > Qmax) { + Qmax = q[ij]; if (Qmax > max_real / 10.) - vrna_message_warning("Q close to overflow: %d %d %g", i, j, temp); + vrna_message_warning("Q close to overflow: %d %d %g", i, j, q[ij]); } - if (temp >= max_real) { + if (q[ij] >= max_real) { vrna_message_warning("overflow while computing partition function for segment q[%d,%d]\n" "use larger pf_scale", i, j); diff --git a/src/ViennaRNA/plotting/alignments.c b/src/ViennaRNA/plotting/alignments.c index e8bf47331..18af7d24f 100644 --- a/src/ViennaRNA/plotting/alignments.c +++ b/src/ViennaRNA/plotting/alignments.c @@ -294,13 +294,13 @@ vrna_file_PS_aln_sub(const char *filename, for (i = currPos; (i < currPos + columnWidth && i < length); i++) { match = 0; for (j = 0; j < N; j++) { - if (cons[i] == seqs[j][i]) + if (cons[i] == toupper(seqs[j][i])) match++; - if (cons[i] == 'U' && seqs[j][i] == 'T') + if (cons[i] == 'U' && toupper(seqs[j][i]) == 'T') match++; - if (cons[i] == 'T' && seqs[j][i] == 'U') + if (cons[i] == 'T' && toupper(seqs[j][i]) == 'U') match++; } score = (float)(match - 1) / (N - 1); diff --git a/src/ViennaRNA/plotting/plot_utils.c b/src/ViennaRNA/plotting/plot_utils.c index e356e4328..fc8501613 100644 --- a/src/ViennaRNA/plotting/plot_utils.c +++ b/src/ViennaRNA/plotting/plot_utils.c @@ -92,8 +92,8 @@ vrna_annotate_covar_struct(const char **alignment, continue; for (s = 0; alignment[s] != NULL; s++) { - a = vrna_nucleotide_encode(toupper(alignment[s][i - 1]), &md); - b = vrna_nucleotide_encode(toupper(alignment[s][j - 1]), &md); + a = vrna_nucleotide_encode(alignment[s][i - 1], &md); + b = vrna_nucleotide_encode(alignment[s][j - 1], &md); type = md.pair[a][b]; pfreq[type]++; if (type) { @@ -189,8 +189,8 @@ vrna_annotate_covar_pairs(const char **alignment, for (z = 0; z < 7; z++) pfreq[z] = 0; for (s = 0; s < n_seq; s++) { - a = vrna_nucleotide_encode(toupper(alignment[s][cp[c].i - 1]), &md); - b = vrna_nucleotide_encode(toupper(alignment[s][cp[c].j - 1]), &md); + a = vrna_nucleotide_encode(alignment[s][cp[c].i - 1], &md); + b = vrna_nucleotide_encode(alignment[s][cp[c].j - 1], &md); if ((alignment[s][cp[c].j - 1] == '~') || (alignment[s][cp[c].i - 1] == '~')) continue; diff --git a/src/ViennaRNA/utils/basic.h b/src/ViennaRNA/utils/basic.h index 597003dce..d6b6a83a8 100644 --- a/src/ViennaRNA/utils/basic.h +++ b/src/ViennaRNA/utils/basic.h @@ -30,7 +30,7 @@ #include /* two helper macros to indicate whether a function should be exported in -the library or stays hidden */ + * the library or stays hidden */ #define PUBLIC #define PRIVATE static @@ -51,7 +51,7 @@ the library or stays hidden */ * @brief Input/Output flag of get_input_line():\n * if used as input option this tells get_input_line() that the data to be read should comply * with the FASTA format - * + * * the function will return this flag if a fasta header was read */ #define VRNA_INPUT_FASTA_HEADER 8U @@ -59,13 +59,13 @@ the library or stays hidden */ /* * @brief Input flag for get_input_line():\n * Tell get_input_line() that we assume to read a nucleotide sequence - * + * */ #define VRNA_INPUT_SEQUENCE 16U /** @brief Input flag for get_input_line():\n * Tell get_input_line() that we assume to read a structure constraint - * + * */ #define VRNA_INPUT_CONSTRAINT 32U @@ -118,18 +118,18 @@ the library or stays hidden */ /** * @brief Get the minimum of three comparable values */ -#define MIN3(A, B, C) (MIN2( (MIN2((A),(B))) ,(C))) +#define MIN3(A, B, C) (MIN2((MIN2((A), (B))), (C))) /** * @brief Get the maximum of three comparable values */ -#define MAX3(A, B, C) (MAX2( (MAX2((A),(B))) ,(C))) +#define MAX3(A, B, C) (MAX2((MAX2((A), (B))), (C))) #ifdef WITH_DMALLOC /* use dmalloc library to check for memory management bugs */ #include "dmalloc.h" -#define vrna_alloc(S) calloc(1,(S)) +#define vrna_alloc(S) calloc(1, (S)) #define vrna_realloc(p, S) xrealloc(p, S) #else @@ -139,7 +139,8 @@ the library or stays hidden */ * @param size The size of the memory to be allocated in bytes * @return A pointer to the allocated memory */ -void *vrna_alloc(unsigned size); +void *vrna_alloc(unsigned size); + /** * @brief Reallocate space safely @@ -148,7 +149,9 @@ void *vrna_alloc(unsigned size); * @param size The size of the memory to be allocated in bytes * @return A pointer to the newly allocated memory */ -void *vrna_realloc(void *p, unsigned size); +void *vrna_realloc(void *p, + unsigned size); + #endif @@ -157,6 +160,7 @@ void *vrna_realloc(void *p, unsigned size); */ void vrna_init_rand(void); + /** * @brief Current 48 bit random number * @@ -176,6 +180,7 @@ extern unsigned short xsubi[3]; */ double vrna_urn(void); + /** * @brief Generates a pseudo random integer in a specified range * @@ -184,7 +189,9 @@ double vrna_urn(void); * @param to The last number in range * @return A pseudo random number in range [from, to] */ -int vrna_int_urn(int from, int to); +int vrna_int_urn(int from, + int to); + /** * @brief Get a timestamp @@ -194,7 +201,8 @@ int vrna_int_urn(int from, int to); * * @return A string containing the timestamp */ -char *vrna_time_stamp(void); +char *vrna_time_stamp(void); + /** * Retrieve a line from 'stdin' savely while skipping comment characters and @@ -203,21 +211,21 @@ char *vrna_time_stamp(void); * An option argument allows one to switch between different reading modes.\n * Currently available options are:\n * #VRNA_INPUT_NOPRINT_COMMENTS, #VRNA_INPUT_NOSKIP_COMMENTS, #VRNA_INPUT_NOELIM_WS_SUFFIX - * + * * pass a collection of options as one value like this: * @verbatim get_input_line(string, option_1 | option_2 | option_n) @endverbatim - * + * * If the function recognizes the type of input, it will report it in the return * value. It also reports if a user defined 'quit' command (@-sign on 'stdin') * was given. Possible return values are:\n * #VRNA_INPUT_FASTA_HEADER, #VRNA_INPUT_ERROR, #VRNA_INPUT_MISC, #VRNA_INPUT_QUIT - * + * * @param string A pointer to the character array that contains the line read * @param options A collection of options for switching the functions behavior * @return A flag with information about what has been read */ -unsigned int get_input_line(char **string, - unsigned int options); +unsigned int get_input_line(char **string, + unsigned int options); /** @@ -226,31 +234,33 @@ unsigned int get_input_line(char **string, * Access of a position "(i,j)" is then accomplished by using @verbatim (i,j) ~ iindx[i]-j @endverbatim * This function is necessary as most of the two-dimensional energy matrices are actually one-dimensional arrays throughout * the ViennaRNA Package - * + * * Consult the implemented code to find out about the mapping formula ;) - * + * * @see vrna_idx_col_wise() * @param length The length of the RNA sequence * @return The mapper array */ int *vrna_idx_row_wise(unsigned int length); + /** * @brief Get an index mapper array (indx) for accessing the energy matrices, e.g. in MFE related functions. * * Access of a position "(i,j)" is then accomplished by using @verbatim (i,j) ~ indx[j]+i @endverbatim * This function is necessary as most of the two-dimensional energy matrices are actually one-dimensional arrays throughout * the ViennaRNAPackage - * + * * Consult the implemented code to find out about the mapping formula ;) - * + * * @see vrna_idx_row_wise() * @param length The length of the RNA sequence * @return The mapper array - * + * */ int *vrna_idx_col_wise(unsigned int length); + /** * @} */ @@ -273,7 +283,8 @@ int *vrna_idx_col_wise(unsigned int length); * @param format The error message to be printed * @param ... Optional arguments for the formatted message string */ -void vrna_message_error(const char *format, ...); +void vrna_message_error(const char *format, + ...); /** @@ -288,7 +299,8 @@ void vrna_message_error(const char *format, ...); * @param format The error message to be printed * @param args The argument list for the formatted message string */ -void vrna_message_verror(const char *format, va_list args); +void vrna_message_verror(const char *format, + va_list args); /** @@ -302,7 +314,8 @@ void vrna_message_verror(const char *format, va_list args); * @param format The warning message to be printed * @param ... Optional arguments for the formatted message string */ -void vrna_message_warning(const char *format, ...); +void vrna_message_warning(const char *format, + ...); /** @@ -316,7 +329,8 @@ void vrna_message_warning(const char *format, ...); * @param format The warning message to be printed * @param args The argument list for the formatted message string */ -void vrna_message_vwarning(const char *format, va_list args); +void vrna_message_vwarning(const char *format, + va_list args); /** @@ -330,7 +344,9 @@ void vrna_message_vwarning(const char *format, va_list args); * @param format The warning message to be printed * @param ... Optional arguments for the formatted message string */ -void vrna_message_info(FILE *fp, const char *format, ...); +void vrna_message_info(FILE *fp, + const char *format, + ...); /** @@ -344,7 +360,9 @@ void vrna_message_info(FILE *fp, const char *format, ...); * @param format The info message to be printed * @param args The argument list for the formatted message string */ -void vrna_message_vinfo(FILE *fp, const char *format, va_list args); +void vrna_message_vinfo(FILE *fp, + const char *format, + va_list args); /** @@ -360,22 +378,24 @@ void vrna_message_input_seq_simple(void); * * (usually this is used to ask for user input) * There will also be a ruler (scale line) printed that helps orientation of the sequence positions - * + * * @param s A user defined string that will be printed to stdout */ void vrna_message_input_seq(const char *s); + void vrna_message_input_msa(const char *s); + /** * @} */ #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY -DEPRECATED(int *get_indx(unsigned int length), "Use vrna_idx_col_wise() instead"); +DEPRECATED(int *get_indx(unsigned int length), "Use vrna_idx_col_wise() instead"); -DEPRECATED(int *get_iindx(unsigned int length), "Use vrna_idx_row_wise() instead"); +DEPRECATED(int *get_iindx(unsigned int length), "Use vrna_idx_row_wise() instead"); /** * @brief Read a line of arbitrary length from a stream @@ -389,7 +409,7 @@ DEPRECATED(int *get_iindx(unsigned int length), "Use vrna_idx_row_wise() inste * @param fp A file pointer to the stream where the function should read from * @return A pointer to the resulting string */ -DEPRECATED(char *get_line(FILE *fp), "Use vrna_read_line() instead"); +DEPRECATED(char *get_line(FILE *fp), "Use vrna_read_line() instead"); /** * @brief Print a line to @e stdout that asks for an input sequence @@ -404,7 +424,7 @@ DEPRECATED(void print_tty_input_seq(void), "Use vrna_message_input_seq_simple() * * (usually this is used to ask for user input) * There will also be a ruler (scale line) printed that helps orientation of the sequence positions - * + * * @deprecated Use vrna_message_input_seq() instead! */ DEPRECATED(void print_tty_input_seq_str(const char *s), "Use vrna_message_input_seq() instead"); @@ -437,7 +457,8 @@ DEPRECATED(void *space(unsigned size), "Use vrna_alloc() instead"); * * @deprecated Use vrna_realloc() instead! */ -DEPRECATED(void *xrealloc(void *p, unsigned size), "Use vrna_realloc() instead"); +DEPRECATED(void *xrealloc(void *p, + unsigned size), "Use vrna_realloc() instead"); /** * @brief Make random number seeds @@ -457,14 +478,16 @@ DEPRECATED(double urn(void), "Use vrna_urn() instead"); * * @deprecated Use vrna_int_urn() instead! */ -DEPRECATED(int int_urn(int from, int to), "Use vrna_int_urn() instead()"); +DEPRECATED(int int_urn(int from, + int to), "Use vrna_int_urn() instead()"); /** * @brief Inefficient `cp` * * @deprecated Use vrna_file_copy() instead! */ -DEPRECATED(void filecopy(FILE *from, FILE *to), "Use vrna_file_copy() instead"); +DEPRECATED(void filecopy(FILE *from, + FILE *to), "Use vrna_file_copy() instead"); /** * @brief Get a timestamp diff --git a/src/ViennaRNA/utils/cpu.c b/src/ViennaRNA/utils/cpu.c new file mode 100644 index 000000000..3f7543d4e --- /dev/null +++ b/src/ViennaRNA/utils/cpu.c @@ -0,0 +1,171 @@ +/* + * A collection of useful functions to detect various CPU features + * + * (c) 2018 - Ronny Lorenz - ViennaRNA Package + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include +#include + +#include "ViennaRNA/utils/basic.h" +#include "ViennaRNA/utils/cpu.h" + +#ifdef __GNUC__ +# define INLINE inline +#else +# define INLINE +#endif + +/* + * cpuid register flags for (64bit) x86 CPUs + * See also https://en.wikipedia.org/wiki/CPUID + */ + +#define bit_SSE2 (1 << 26) /* stored in EDX after cpuid with EAX=1 */ +#define bit_SSE3 (1 << 0) /* stored in ECX after cpuid with EAX=1 */ +#define bit_SSE41 (1 << 19) /* stored in ECX after cpuid with EAX=1 */ +#define bit_SSE42 (1 << 20) /* stored in ECX after cpuid with EAX=1 */ +#define bit_AVX (1 << 28) /* stored in ECX after cpuid with EAX=1 */ +#define bit_AVX2 (1 << 5) /* stored in EBX after cpuid with EAX=7, ECX=0 */ +#define bit_AVX512F (1 << 16) /* stored in EBX after cpuid with EAX=7, ECX=0 */ + +/* + ################################# + # PRIVATE FUNCTION DECLARATIONS # + ################################# + */ +PRIVATE INLINE int +execute_cpuid(uint32_t *regs); + + +PRIVATE unsigned int +cpu_feature_bits(void); + + +PRIVATE unsigned int +cpu_extended_feature_bits(void); + + +/* + ################################# + # BEGIN OF FUNCTION DEFINITIONS # + ################################# + */ +PUBLIC char * +vrna_cpu_vendor_string(void) +{ + static char name[13] = { + 0 + }; + uint32_t regs[4] = { + 0 + }; + + if (execute_cpuid(®s[0])) { + memcpy(name + 0, ®s[1], 4); + memcpy(name + 4, ®s[3], 4); + memcpy(name + 8, ®s[2], 4); + name[12] = '\0'; + } + + return name; +} + + +PUBLIC unsigned int +vrna_cpu_simd_capabilities(void) +{ + unsigned int capabilities = VRNA_CPU_SIMD_NONE; + + capabilities |= cpu_feature_bits(); + capabilities |= cpu_extended_feature_bits(); + + return capabilities; +} + + +/* + ################################# + # STATIC helper functions below # + ################################# + */ + +/* + * execute 'cpuid' command with values stored in registers regs[0]..regs[3] + * that directly correspond to EAX, EBX, ECX, and EDX. The result of the + * 'cpuid' command is then returned in the same register array + */ +PRIVATE INLINE int +execute_cpuid(uint32_t *regs) +{ +#if defined(__x86_64__) || defined(_M_AMD64) || defined (_M_X64) + __asm__ __volatile__ ("cpuid" + : "=a" (regs[0]), + "=b" (regs[1]), + "=c" (regs[2]), + "=d" (regs[3]) + : "0" (regs[0]), + "2" (regs[2])); + + return 1; +#else + return 0; +#endif +} + + +PRIVATE unsigned int +cpu_feature_bits(void) +{ + unsigned int features = VRNA_CPU_SIMD_NONE; + + uint32_t regs[4] = { + 1, 0, 0, 0 + }; + + if (execute_cpuid(®s[0])) { + if (regs[3] & bit_SSE2) + features |= VRNA_CPU_SIMD_SSE2; + + if (regs[2] & bit_SSE3) + features |= VRNA_CPU_SIMD_SSE3; + + if (regs[2] & bit_SSE41) + features |= VRNA_CPU_SIMD_SSE41; + + if (regs[2] & bit_SSE42) + features |= VRNA_CPU_SIMD_SSE42; + + if (regs[2] & bit_AVX) + features |= VRNA_CPU_SIMD_AVX; + } + + return features; +} + + +PRIVATE unsigned int +cpu_extended_feature_bits(void) +{ + unsigned int features = VRNA_CPU_SIMD_NONE; + + uint32_t regs[4] = { + 7, 0, 0, 0 + }; + + if (execute_cpuid(®s[0])) { + if (regs[1] & bit_AVX2) + features |= VRNA_CPU_SIMD_AVX2; + + if (regs[1] & bit_AVX512F) + features |= VRNA_CPU_SIMD_AVX512F; + } + + return features; +} diff --git a/src/ViennaRNA/utils/cpu.h b/src/ViennaRNA/utils/cpu.h new file mode 100644 index 000000000..b21941171 --- /dev/null +++ b/src/ViennaRNA/utils/cpu.h @@ -0,0 +1,22 @@ +#ifndef VIENNA_RNA_PACKAGE_UTILS_CPU_H +#define VIENNA_RNA_PACKAGE_UTILS_CPU_H + +#define VRNA_CPU_SIMD_NONE 0U +#define VRNA_CPU_SIMD_SSE2 1U +#define VRNA_CPU_SIMD_SSE3 2U +#define VRNA_CPU_SIMD_SSE41 4U +#define VRNA_CPU_SIMD_SSE42 8U +#define VRNA_CPU_SIMD_AVX 16U +#define VRNA_CPU_SIMD_AVX2 32U +#define VRNA_CPU_SIMD_AVX512F 64U + + +char * +vrna_cpu_vendor_string(void); + + +unsigned int +vrna_cpu_simd_capabilities(void); + + +#endif diff --git a/src/ViennaRNA/utils/higher_order_functions.c b/src/ViennaRNA/utils/higher_order_functions.c new file mode 100644 index 000000000..6bbdf30d8 --- /dev/null +++ b/src/ViennaRNA/utils/higher_order_functions.c @@ -0,0 +1,138 @@ +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include + +#include "ViennaRNA/utils/basic.h" +#include "ViennaRNA/utils/cpu.h" + + +typedef int (proto_fun_zip_reduce)(const int *a, + const int *b, + int size); + + +/* + ################################# + # PRIVATE FUNCTION DECLARATIONS # + ################################# + */ +static int zip_add_min_dispatcher(const int *a, + const int *b, + int size); + + +static int +fun_zip_add_min_default(const int *e1, + const int *e2, + int count); + + +#if VRNA_WITH_SIMD_AVX512 +int +vrna_fun_zip_add_min_avx512(const int *e1, + const int *e2, + int count); + + +#endif + +#if VRNA_WITH_SIMD_SSE41 +int +vrna_fun_zip_add_min_sse41(const int *e1, + const int *e2, + int count); + + +#endif + + +static proto_fun_zip_reduce *fun_zip_add_min = &zip_add_min_dispatcher; + + +/* + ################################# + # BEGIN OF FUNCTION DEFINITIONS # + ################################# + */ +PUBLIC void +vrna_fun_dispatch_disable(void) +{ + fun_zip_add_min = &fun_zip_add_min_default; +} + + +PUBLIC void +vrna_fun_dispatch_enable(void) +{ + fun_zip_add_min = &zip_add_min_dispatcher; +} + + +PUBLIC int +vrna_fun_zip_add_min(const int *e1, + const int *e2, + int count) +{ + return (*fun_zip_add_min)(e1, e2, count); +} + + +/* + ################################# + # STATIC helper functions below # + ################################# + */ + +/* zip_add_min() dispatcher */ +static int +zip_add_min_dispatcher(const int *a, + const int *b, + int size) +{ + unsigned int features = vrna_cpu_simd_capabilities(); + +#if VRNA_WITH_SIMD_AVX512 + if (features & VRNA_CPU_SIMD_AVX512F) { + fun_zip_add_min = &vrna_fun_zip_add_min_avx512; + goto exec_fun_zip_add_min; + } + +#endif + +#if VRNA_WITH_SIMD_SSE41 + if (features & VRNA_CPU_SIMD_SSE41) { + fun_zip_add_min = &vrna_fun_zip_add_min_sse41; + goto exec_fun_zip_add_min; + } + +#endif + + fun_zip_add_min = &fun_zip_add_min_default; + +exec_fun_zip_add_min: + + return (*fun_zip_add_min)(a, b, size); +} + + +static int +fun_zip_add_min_default(const int *e1, + const int *e2, + int count) +{ + int i; + int decomp = INF; + + for (i = 0; i < count; i++) { + if ((e1[i] != INF) && (e2[i] != INF)) { + const int en = e1[i] + e2[i]; + decomp = MIN2(decomp, en); + } + } + + return decomp; +} diff --git a/src/ViennaRNA/utils/higher_order_functions.h b/src/ViennaRNA/utils/higher_order_functions.h new file mode 100644 index 000000000..ff032f258 --- /dev/null +++ b/src/ViennaRNA/utils/higher_order_functions.h @@ -0,0 +1,18 @@ +#ifndef VIENNA_RNA_PACKAGE_UTILS_FUN_H +#define VIENNA_RNA_PACKAGE_UTILS_FUN_H + +void +vrna_fun_dispatch_disable(void); + + +void +vrna_fun_dispatch_enable(void); + + +int +vrna_fun_zip_add_min(const int *e1, + const int *e2, + int count); + + +#endif diff --git a/src/ViennaRNA/utils/higher_order_functions_avx512.c b/src/ViennaRNA/utils/higher_order_functions_avx512.c new file mode 100644 index 000000000..1b19374a8 --- /dev/null +++ b/src/ViennaRNA/utils/higher_order_functions_avx512.c @@ -0,0 +1,50 @@ +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include + +#include "ViennaRNA/utils/basic.h" + +#include + + +PUBLIC int +vrna_fun_zip_add_min_avx512(const int *e1, + const int *e2, + int count) +{ + int i = 0; + int decomp = INF; + + __m512i inf = _mm512_set1_epi32(INF); + + /* WBL 21 Aug 2018 Add SSE512 code from sources_034_578/modular_decomposition_id3.c by hand */ + for (i = 0; i < count - 15; i += 16) { + __m512i a = _mm512_loadu_si512((__m512i *)&e1[i]); + __m512i b = _mm512_loadu_si512((__m512i *)&e2[i]); + + /* compute mask for entries where both, a and b, are less than INF */ + __mmask16 mask = _kand_mask16(_mm512_cmplt_epi32_mask(a, inf), + _mm512_cmplt_epi32_mask(b, inf)); + + /* add values */ + __m512i c = _mm512_add_epi32(a, b); + + /* reduce to minimum (only those where one of the source values was not INF before) */ + const int en = _mm512_mask_reduce_min_epi32(mask, c); + + decomp = MIN2(decomp, en); + } + + for (; i < count; i++) { + if ((e1[i] != INF) && (e2[i] != INF)) { + const int en = e1[i] + e2[i]; + decomp = MIN2(decomp, en); + } + } + + return decomp; +} diff --git a/src/ViennaRNA/utils/higher_order_functions_sse41.c b/src/ViennaRNA/utils/higher_order_functions_sse41.c new file mode 100644 index 000000000..ac17a085c --- /dev/null +++ b/src/ViennaRNA/utils/higher_order_functions_sse41.c @@ -0,0 +1,71 @@ +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include + +#include "ViennaRNA/utils/basic.h" + +#include +#include + +static int +horizontal_min_Vec4i(__m128i x); + + +PUBLIC int +vrna_fun_zip_add_min_sse41(const int *e1, + const int *e2, + int count) +{ + int i = 0; + int decomp = INF; + + __m128i inf = _mm_set1_epi32(INF); + + for (i = 0; i < count - 3; i += 4) { + __m128i a = _mm_loadu_si128((__m128i *)&e1[i]); + __m128i b = _mm_loadu_si128((__m128i *)&e2[i]); + __m128i c = _mm_add_epi32(a, b); + + /* create mask for non-INF values */ + __m128i mask = _mm_and_si128(_mm_cmplt_epi32(a, inf), + _mm_cmplt_epi32(b, inf)); + + /* delete results where a or b has been INF before */ + c = _mm_and_si128(mask, c); + + /* fill all values with INF if they've been INF in a or b before */ + __m128i res = _mm_or_si128(c, _mm_andnot_si128(mask, inf)); + const int en = horizontal_min_Vec4i(res); + + decomp = MIN2(decomp, en); + } + + for (; i < count; i++) { + if ((e1[i] != INF) && (e2[i] != INF)) { + const int en = e1[i] + e2[i]; + decomp = MIN2(decomp, en); + } + } + + return decomp; +} + + +/* + * SSE minimum + * see also: http://stackoverflow.com/questions/9877700/getting-max-value-in-a-m128i-vector-with-sse + */ +static int +horizontal_min_Vec4i(__m128i x) +{ + __m128i min1 = _mm_shuffle_epi32(x, _MM_SHUFFLE(0, 0, 3, 2)); + __m128i min2 = _mm_min_epi32(x, min1); + __m128i min3 = _mm_shuffle_epi32(min2, _MM_SHUFFLE(0, 0, 0, 1)); + __m128i min4 = _mm_min_epi32(min2, min3); + + return _mm_cvtsi128_si32(min4); +} diff --git a/src/ViennaRNA/utils/msa_utils.c b/src/ViennaRNA/utils/msa_utils.c index f219a55d7..b265930e4 100644 --- a/src/ViennaRNA/utils/msa_utils.c +++ b/src/ViennaRNA/utils/msa_utils.c @@ -450,8 +450,8 @@ vrna_aln_conservation_struct(const char **alignment, if (i < pt[i]) { j = pt[i]; for (s = 0; s < n_seq; s++) { - a = vrna_nucleotide_encode(toupper(alignment[s][i - 1]), &md); - b = vrna_nucleotide_encode(toupper(alignment[s][j - 1]), &md); + a = vrna_nucleotide_encode(alignment[s][i - 1], &md); + b = vrna_nucleotide_encode(alignment[s][j - 1], &md); if (md.pair[a][b]) { conservation[i] += 1.; conservation[j] += 1.; @@ -517,7 +517,7 @@ vrna_aln_conservation_col(const char **alignment, /* count frequencies of individual nucleotides */ for (s = 0; s < n_seq; s++) - nfreq[vrna_nucleotide_encode(toupper(alignment[s][i - 1]), &md)]++; + nfreq[vrna_nucleotide_encode(alignment[s][i - 1], &md)]++; if (options & VRNA_MEASURE_SHANNON_ENTROPY) { double sum = 0; @@ -572,7 +572,7 @@ vrna_aln_consensus_sequence(const char **alignment, else vrna_md_set_default(&md); - consensus = (char *)vrna_alloc((n + 1) * sizeof(char)); + consensus = (char *)vrna_alloc((n + 1) * sizeof(char)); for (i = 0; i < n; i++) { int c, fm; @@ -585,8 +585,8 @@ vrna_aln_consensus_sequence(const char **alignment, for (s = c = fm = 0; s < 8; s++) /* find the most frequent char */ if (freq[s] > fm) { - c = s; - fm = freq[c]; + c = s; + fm = freq[c]; } if (s > 4) diff --git a/src/ViennaRNA/utils/structure_utils.c b/src/ViennaRNA/utils/structure_utils.c index 0774bcf88..1482f12a5 100644 --- a/src/ViennaRNA/utils/structure_utils.c +++ b/src/ViennaRNA/utils/structure_utils.c @@ -442,6 +442,74 @@ vrna_bp_distance(const char *str1, } +PUBLIC double +vrna_dist_mountain(const char *str1, + const char *str2, + unsigned int p) +{ + short *pt1, *pt2; + unsigned int i, n; + double distance, w, *f1, *f2; + + distance = -1.; + f1 = NULL; + f2 = NULL; + + if ((str1) && (str2)) { + n = strlen(str1); + + if (n != strlen(str2)) { + vrna_message_warning("vrna_dist_mountain: input structures have unequal lengths!"); + return distance; + } + + pt1 = vrna_ptable(str1); + pt2 = vrna_ptable(str2); + f1 = (double *)vrna_alloc(sizeof(double) * (n + 1)); + f2 = (double *)vrna_alloc(sizeof(double) * (n + 1)); + + /* count (mountain)heights for positions 1 <= i <= n */ + for (w = 0., i = 1; i <= n; i++) { + if (pt1[i] == 0) + continue; + + if (pt1[i] > i) + w += 1. / (double)(pt1[i] - i); + else + w -= 1. / (double)(i - pt1[i]); + + f1[i] = w; + } + + for (w = 0., i = 1; i <= n; i++) { + if (pt2[i] == 0) + continue; + + if (pt2[i] > i) + w += 1. / (double)(pt2[i] - i); + else + w -= 1. / (double)(i - pt2[i]); + + f2[i] = w; + } + + + /* finally, compute L_p-norm */ + for (distance = 0., i = 1; i <= n; i++) + distance += pow(fabs(f1[i] - f2[i]), (double)p); + + distance = pow(distance, 1. / (double)p); + + free(pt1); + free(pt2); + free(f1); + free(f2); + } + + return distance; +} + + /* get a matrix containing the number of basepairs of a reference structure for each interval [i,j] with ierr); + /* flush/free errors first */ vrna_cstr_free(s->err); - /* flush data[k] */ - vrna_cstr_fflush(s->data); - /* free data[k] */ + /* flush/free data[k] */ vrna_cstr_free(s->data); free(s); diff --git a/src/bin/RNAcofold.c b/src/bin/RNAcofold.c index 9450a26d0..4911b97d5 100644 --- a/src/bin/RNAcofold.c +++ b/src/bin/RNAcofold.c @@ -237,13 +237,10 @@ flush_cstr_callback(void *auxdata, { struct output_stream *s = (struct output_stream *)data; - /* flush errors first */ - vrna_cstr_fflush(s->err); + /* flush/free errors first */ vrna_cstr_free(s->err); - /* flush data[k] */ - vrna_cstr_fflush(s->data); - /* free data[k] */ + /* flush/free data[k] */ vrna_cstr_free(s->data); free(s); diff --git a/src/bin/RNAeval.c b/src/bin/RNAeval.c index cee6948f5..e6ad58593 100644 --- a/src/bin/RNAeval.c +++ b/src/bin/RNAeval.c @@ -183,13 +183,10 @@ flush_cstr_callback(void *auxdata, { struct output_stream *s = (struct output_stream *)data; - /* flush errors first */ - vrna_cstr_fflush(s->err); + /* flush/free errors first */ vrna_cstr_free(s->err); - /* flush data[k] */ - vrna_cstr_fflush(s->data); - /* free data[k] */ + /* flush/free data[k] */ vrna_cstr_free(s->data); free(s); @@ -631,6 +628,8 @@ process_record(struct record_data *record) &(opt->md), VRNA_OPTION_MFE | VRNA_OPTION_EVAL_ONLY); + n = (int)vc->length; + if (opt->shape) { vrna_constraints_add_SHAPE(vc, opt->shape_file, @@ -725,8 +724,8 @@ static void process_alignment_record(struct record_data_msa *record) { char **alignment, *consensus_sequence, *structure; - unsigned int n, n_seq; - double energy, real_en, cov_en; + unsigned int n; + double real_en, cov_en; struct options *opt; vrna_fold_compound_t *vc; struct output_stream *o_stream; @@ -737,7 +736,6 @@ process_alignment_record(struct record_data_msa *record) vrna_message_error("structure missing for record %d\n", record->number); opt = record->options; - n_seq = record->n_seq; structure = vrna_db_from_WUSS(record->structure); /* diff --git a/src/bin/RNAfold.c b/src/bin/RNAfold.c index b3fcb96e7..b7f0be341 100644 --- a/src/bin/RNAfold.c +++ b/src/bin/RNAfold.c @@ -201,9 +201,7 @@ flush_cstr_callback(void *auxdata, struct output_stream *s = (struct output_stream *)data; if (s) { - /* flush data[k] */ - vrna_cstr_fflush(s->data); - /* free/close data[k] */ + /* flush/free/close data[k] */ if (s->individual) vrna_cstr_close(s->data); else diff --git a/src/bin/RNAplot.c b/src/bin/RNAplot.c index bb9702ad6..42592e092 100644 --- a/src/bin/RNAplot.c +++ b/src/bin/RNAplot.c @@ -641,14 +641,12 @@ process_alignment_record(struct record_data_msa *record) { char *consensus_sequence, *structure, *ffname, *tmp_string, **A, *pre, *post; - unsigned int n_seq; struct options *opt; if (!record->structure) vrna_message_error("structure missing for record %d\n", record->number); opt = record->options; - n_seq = record->n_seq; structure = vrna_db_from_WUSS(record->structure); pre = opt->pre; post = opt->post; diff --git a/src/bin/RNAsubopt.c b/src/bin/RNAsubopt.c index ab644c8a6..76b2c6acd 100644 --- a/src/bin/RNAsubopt.c +++ b/src/bin/RNAsubopt.c @@ -30,6 +30,7 @@ #include "ViennaRNA/constraints/SHAPE.h" #include "ViennaRNA/io/file_formats.h" #include "ViennaRNA/io/utils.h" +#include "ViennaRNA/commands.h" #include "RNAsubopt_cmdl.h" #include "gengetopt_helper.h" #include "input_id_helpers.h" @@ -40,6 +41,24 @@ PRIVATE void putoutzuker(FILE *output, vrna_subopt_solution_t *zukersolution); +struct nr_en_data { + FILE *output; + vrna_fold_compound_t *fc; + double kT; + double ens_en; +}; + + +PRIVATE void +print_nr_samples(const char *structure, + void *data); + + +PRIVATE void +print_nr_samples_en(const char *structure, + void *data); + + int main(int argc, char *argv[]) @@ -53,10 +72,11 @@ main(int argc, unsigned int rec_type, read_opt; int i, length, cl, istty, delta, n_back, noconv, dos, zuker, with_shapes, verbose, enforceConstraints, st_back_en, batch, - tofile, filename_full, canonicalBPonly; + tofile, filename_full, canonicalBPonly, nonRedundant; double deltap; vrna_md_t md; dataset_id id_control; + vrna_cmd_t commands; do_backtrack = 1; delta = 100; @@ -73,6 +93,8 @@ main(int argc, tofile = 0; filename_full = 0; canonicalBPonly = 0; + commands = NULL; + nonRedundant = 0; set_model_details(&md); @@ -210,6 +232,14 @@ main(int argc, if (args_info.filename_full_given) filename_full = 1; + /* non-redundant backtracing */ + if (args_info.nonRedundant_given) + nonRedundant = 1; + + if (args_info.commands_given) + commands = vrna_file_commands_read(args_info.commands_arg, + VRNA_CMD_PARSE_HC | VRNA_CMD_PARSE_SC); + /* free allocated memory of command line data structure */ RNAsubopt_cmdline_parser_free(&args_info); @@ -322,8 +352,7 @@ main(int argc, /* parse the rest of the current dataset to obtain a structure constraint */ if (fold_constrained) { if (constraints_file) { - vrna_constraints_add(vc, constraints_file, - VRNA_OPTION_MFE | ((n_back > 0) ? VRNA_OPTION_PF : 0)); + vrna_constraints_add(vc, constraints_file, VRNA_OPTION_DEFAULT); } else { cstruc = NULL; int cp = -1; @@ -369,6 +398,11 @@ main(int argc, VRNA_OPTION_MFE | ((n_back > 0) ? VRNA_OPTION_PF : 0)); } + if (commands) + vrna_commands_apply(vc, + commands, + VRNA_CMD_PARSE_HC | VRNA_CMD_PARSE_SC); + if (istty) { if (cut_point == -1) { vrna_message_info(stdout, "length = %d", length); @@ -411,19 +445,33 @@ main(int argc, ens_en = vrna_pf(vc, structure); kT = vc->exp_params->kT / 1000.; - for (i = 0; i < n_back; i++) { - char *s, *e_string = NULL; - s = vrna_pbacktrack(vc); + if (nonRedundant) { if (st_back_en) { - double e, prob; - e = vrna_eval_structure(vc, s); - prob = exp((ens_en - e) / kT); - e_string = vrna_strdup_printf(" %6.2f %6g", e, prob); + struct nr_en_data dat; + dat.output = output; + dat.fc = vc; + dat.kT = kT; + dat.ens_en = ens_en; + + vrna_pbacktrack_nr_cb(vc, n_back, &print_nr_samples_en, (void *)&dat); + } else { + vrna_pbacktrack_nr_cb(vc, n_back, &print_nr_samples, (void *)output); + } + } else { + for (i = 0; i < n_back; i++) { + char *s, *e_string = NULL; + s = vrna_pbacktrack(vc); + if (st_back_en) { + double e, prob; + e = vrna_eval_structure(vc, s); + prob = exp((ens_en - e) / kT); + e_string = vrna_strdup_printf(" %6.2f %6g", e, prob); + } + + print_structure(output, s, e_string); + free(s); + free(e_string); } - - print_structure(output, s, e_string); - free(s); - free(e_string); } } /* normal subopt */ @@ -523,6 +571,7 @@ main(int argc, free(shape_method); free(shape_conversion); free(filename_delim); + vrna_commands_free(commands); free_id_data(id_control); @@ -530,6 +579,37 @@ main(int argc, } +PRIVATE void +print_nr_samples(const char *structure, + void *data) +{ + if (structure) + print_structure((FILE *)data, structure, NULL); +} + + +PRIVATE void +print_nr_samples_en(const char *structure, + void *data) +{ + if (structure) { + struct nr_en_data *d = (struct nr_en_data *)data; + FILE *output = d->output; + vrna_fold_compound_t *fc = d->fc; + double kT = d->kT; + double ens_en = d->ens_en; + + double e = vrna_eval_structure(fc, structure); + double prob = exp((ens_en - e) / kT); + char *e_string = vrna_strdup_printf(" %6.2f %6g", e, prob); + + print_structure(output, structure, e_string); + + free(e_string); + } +} + + PRIVATE void putoutzuker(FILE *output, vrna_subopt_solution_t *zukersolution) diff --git a/src/bin/RNAsubopt.ggo b/src/bin/RNAsubopt.ggo index 2e88970ca..ceb9d2428 100644 --- a/src/bin/RNAsubopt.ggo +++ b/src/bin/RNAsubopt.ggo @@ -230,6 +230,17 @@ typestr="M/C/S/L/O + [optional parameters]" default="O" optional +option "commands" - +"Read additional commands from file\n" +details="Commands include hard and soft constraints, but also structure motifs in hairpin and interior\ + loops that need to be treeted differently. Furthermore, commands can be set for unstructured and structured\ + domains.\n\n" +string +typestr="" +optional +hidden + + section "Algorithms" sectiondesc="Select the algorithms which should be applied to the given RNA sequence.\n\n" @@ -270,6 +281,13 @@ int typestr="number" optional + +option "nonRedundant" N +"Enable non-redundant sampling strategy.\n\n" +flag +off + + option "pfScale" S "In the calculation of the pf use scale*mfe as an estimate for the ensemble free energy (used to avoid\ overflows). Needed by stochastic backtracking\n" diff --git a/tests/Makefile.am b/tests/Makefile.am index f1417e2d5..8510ba47d 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -11,6 +11,10 @@ AM_LDFLAGS = $(RNA_LDFLAGS) $(PTHREAD_LIBS) LDADD = $(top_builddir)/src/ViennaRNA/libRNA_conv.la @CHECK_LIBS@ +if VRNA_AM_SWITCH_MPFR +LDADD += $(MPFR_LIBS) +endif + SUFFIXES = .c .ts TEST_EXTENSIONS = .pl .t .py .py3 diff --git a/tests/RNAHelpers.pm.in b/tests/RNAHelpers.pm.in index 6bb48289f..b993835d1 100644 --- a/tests/RNAHelpers.pm.in +++ b/tests/RNAHelpers.pm.in @@ -6,11 +6,12 @@ use warnings; use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS); -$VERSION = 1.00; +$VERSION = "@VERSION@"; @ISA = qw(Exporter); @EXPORT = (); @EXPORT_OK = qw(MsgChecking getDataDirPath); -%EXPORT_TAGS = ( Messages => [qw(&MsgChecking)], +%EXPORT_TAGS = ( VERSION => $VERSION, + Messages => [qw(&MsgChecking)], Paths => [qw(&getDataDirPath)] ); sub MsgChecking diff --git a/tests/RNApath.py.in b/tests/RNApath.py.in index 245c6fecf..bbefdcfc5 100644 --- a/tests/RNApath.py.in +++ b/tests/RNApath.py.in @@ -1,5 +1,7 @@ import sys +VERSION_NUMBER = "@VERSION@" + def addSwigInterfacePath(version=2): if version == 3: sys.path.insert(0, '@top_builddir@/interfaces/Python3') diff --git a/tests/perl5/test-RNA.pl b/tests/perl5/test-RNA.pl index 5de75d8be..1b458f3ab 100755 --- a/tests/perl5/test-RNA.pl +++ b/tests/perl5/test-RNA.pl @@ -3,12 +3,13 @@ ######################### We start with some black magic to print on failure. # (It may become useful if the test is moved to ./t subdirectory.) use strict; -use Test::More tests => 55; +use Test::More tests => 56; use Data::Dumper; use FileHandle; use RNA; use warnings; +use RNAHelpers qw(:VERSION); ######################### End of black magic. @@ -22,6 +23,9 @@ my $struc2="((((((...))).)))"; +# Version number +is($RNA::VERSION, $RNAHelpers::VERSION); + # calculate a hamming distance (from util.c) is(RNA::hamming($seq1, $seq2), 16); diff --git a/tests/python/test-RNA.py b/tests/python/test-RNA.py index a506383c1..2e38f1b35 100644 --- a/tests/python/test-RNA.py +++ b/tests/python/test-RNA.py @@ -19,6 +19,10 @@ class GeneralTests(unittest.TestCase): + def test_version(self): + print "test version number" + self.assertEqual(RNA.__version__, RNApath.VERSION_NUMBER) + def test_hammingDistance(self): print "test_hammingDistance \t calculate a hamming distance" self.assertEqual(RNA.hamming(seq1,seq2),16) diff --git a/tests/python3/test-RNA.py3 b/tests/python3/test-RNA.py3 index 1335c6826..6df589c9d 100644 --- a/tests/python3/test-RNA.py3 +++ b/tests/python3/test-RNA.py3 @@ -23,6 +23,10 @@ align = [seq1,seq2,seq3] class GeneralTests(unittest.TestCase): + def test_version(self): + print("test version number") + self.assertEqual(RNA.__version__, RNApath.VERSION_NUMBER) + def test_hammingDistance(self): print("test_hammingDistance \t calculate a hamming distance") self.assertEqual(RNA.hamming(seq1,seq2),16)