Mutsel implementation

.github/workflows/ci.yaml

@@ -45,6 +45,13 @@ jobs:
           sudo apt-get update -q -y
           sudo apt-get install -q -y clang-${{ matrix.clang-version }} lld-${{ matrix.clang-version }} libomp-${{ matrix.clang-version }}-dev libeigen3-dev libboost-dev
           sudo ln -sf $(which ld.lld-${{ matrix.clang-version }}) /usr/local/bin/ld.lld
+      # Install Rust
+      - name: Install Rust
+        uses: actions-rs/toolchain@v1
+        with:
+          toolchain: stable
+          override: true
+          profile: minimal
 
       - name: Setup cmake
         uses: jwlawson/actions-setup-cmake@v2
@@ -55,7 +62,7 @@ jobs:
         run: |
           mkdir build
           cd build
-          cmake .. -DCMAKE_POLICY_VERSION_MINIMUM=3.5 -DIQTREE_FLAGS=static -DCMAKE_C_COMPILER=$(which clang-${{ matrix.clang-version }}) -DCMAKE_CXX_COMPILER=$(which clang++-${{ matrix.clang-version }})
+          cmake .. -DUSE_MUTSEL=ON -DCMAKE_POLICY_VERSION_MINIMUM=3.5 -DIQTREE_FLAGS=static -DCMAKE_C_COMPILER=$(which clang-${{ matrix.clang-version }}) -DCMAKE_CXX_COMPILER=$(which clang++-${{ matrix.clang-version }})
           make -j
           make package
           file iqtree3 | grep x86-64
@@ -109,11 +116,18 @@ jobs:
           sudo apt-get update -q -y
           sudo apt-get install -q -y gcc g++ libeigen3-dev libboost-dev
 
+      # Install Rust
+      - name: Install Rust
+        uses: actions-rs/toolchain@v1
+        with:
+          toolchain: stable
+          override: true
+          profile: minimal
       - name: Build
         run: |
           mkdir build
           cd build
-          cmake .. -DIQTREE_FLAGS=static
+          cmake .. -DIQTREE_FLAGS=static -DUSE_MUTSEL=ON 
           make -j
           make package
           file iqtree3 | grep aarch64
@@ -199,14 +213,21 @@ jobs:
       - name: Install dependencies
         run: brew install make eigen boost libomp
 
+      # Install Rust
+      - name: Install Rust
+        uses: actions-rs/toolchain@v1
+        with:
+          toolchain: stable
+          override: true
+          profile: minimal
       - name: Build
         run: |
           set -x
           mkdir build
           cd build
           export CPPFLAGS="-I/usr/local/opt/libomp/include"
           export CXXFLAGS="-I/usr/local/opt/libomp/include"
-          cmake .. -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ -DEIGEN3_INCLUDE_DIR=$(brew --prefix eigen)/include/eigen3
+          cmake .. -DUSE_MUTSEL=ON -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ -DEIGEN3_INCLUDE_DIR=$(brew --prefix eigen)/include/eigen3
           make -j
           make package
           file iqtree3 | grep x86_64
@@ -256,14 +277,21 @@ jobs:
       - name: Install dependencies
         run: brew install make eigen boost libomp
 
+      # Install Rust
+      - name: Install Rust
+        uses: actions-rs/toolchain@v1
+        with:
+          toolchain: stable
+          override: true
+          profile: minimal
       - name: Build
         run: |
           set -x
           mkdir build
           cd build
           export CPPFLAGS="-I/opt/homebrew/opt/libomp/include"
           export CXXFLAGS="-I/opt/homebrew/opt/libomp/include"
-          cmake .. -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ -DEIGEN3_INCLUDE_DIR=$(brew --prefix eigen)/include/eigen3
+          cmake .. -DUSE_MUTSEL=ON -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ -DEIGEN3_INCLUDE_DIR=$(brew --prefix eigen)/include/eigen3
           make -j
           make package
           file iqtree3 | grep arm64
@@ -378,21 +406,34 @@ jobs:
         shell: cmd
         run: choco install eigen
 
+      # Install Rust
+      - name: Install Rust
+        uses: actions-rs/toolchain@v1
+        with:
+          toolchain: stable
+          override: true
+          profile: minimal
+          target: x86_64-pc-windows-gnu
       - name: Compile
         shell: cmd
         run: |
           if exist build rd /s /q build
           mkdir build
           cd build
           cmake -G "MinGW Makefiles" ^
+            -DUSE_MUTSEL=ON ^
             -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ ^
-            -DCMAKE_C_FLAGS=--target=x86_64-pc-windows-gnu ^
-            -DCMAKE_CXX_FLAGS=--target=x86_64-pc-windows-gnu ^
+            -DCMAKE_C_FLAGS="--target=x86_64-pc-windows-gnu --sysroot=C:/mingw64 -B C:/mingw64/bin" ^
+            -DCMAKE_CXX_FLAGS="--target=x86_64-pc-windows-gnu --sysroot=C:/mingw64 -B C:/mingw64/bin" ^
             -DCMAKE_MAKE_PROGRAM=mingw32-make ^
             -DBoost_INCLUDE_DIR="C:/local/boost/include" ^
             -DBoost_LIBRARY_DIRS="C:/local/boost/lib" ^
             -DCMAKE_POLICY_VERSION_MINIMUM=3.5 ^
-            -DIQTREE_FLAGS="static cpp14" ..
+            -DIQTREE_FLAGS="static" ..
+
+          type CMakeFiles\iqtree3.dir\link.txt
+          type CMakeFiles\iqtree3.dir\linkLibs.rsp
+
           make -j
           make package
         env:

.gitignore

@@ -35,3 +35,4 @@ zlib-1.2.7/
 /.direnv/
 /.envrc
 softwipe_build
+mutsel_rust/target/

CMakeLists.txt

@@ -284,6 +284,43 @@ else()
     message("Integrate CMAPLE_AA: OFF")
 endif()
 
+##################################################################
+# Include Rust dependecies: phylo_grad, MUTSEL model
+##################################################################
+
+if (USE_MUTSEL STREQUAL "ON")
+    message("Compile Rust Dependencies: ON")
+    if (WIN32)
+        set(MUTSEL_RUST_LIB_FILE "x86_64-pc-windows-gnu/release/libmutsel_rust.a")
+    else()
+        set(MUTSEL_RUST_LIB_FILE "release/libmutsel_rust.a")
+    endif()
+    set(MUTSEL_RUST_LIB_PATH "${CMAKE_SOURCE_DIR}/mutsel_rust/target/${MUTSEL_RUST_LIB_FILE}")
+
+    if (WIN32)
+        set(CARGO_TARGET --target x86_64-pc-windows-gnu)
+    else()
+        set(CARGO_TARGET "")
+    endif()
+
+
+    include(ExternalProject)
+    ExternalProject_Add(
+        mutsel
+        SOURCE_DIR "${CMAKE_SOURCE_DIR}/mutsel_rust"  # Where your Rust source is
+        CONFIGURE_COMMAND ""                           # No configure step for Rust
+        BUILD_COMMAND cargo build --release --lib ${CARGO_TARGET}
+        BUILD_BYPRODUCTS "${MUTSEL_RUST_LIB_PATH}"
+        INSTALL_COMMAND ""  # No install step for Rust
+        BUILD_IN_SOURCE 1  # Build in source directory
+        LOG_BUILD ON
+        LOG_OUTPUT_ON_FAILURE ON
+    )
+    add_definitions(-DUSE_MUTSEL)
+else()
+    message("Compile Rust Dependencies: OFF")
+endif()
+
 ##################################################################
 # build IQTREE library
 ##################################################################
@@ -313,16 +350,10 @@ elseif (APPLE)
         # set the first version of mac 11 that supports ARM
         set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} --target=arm64-apple-macos11.0.1")
 	else()
-		# to be compatible back to Mac OS X 10.7
-    	if (IQTREE_FLAGS MATCHES "oldmac")
-        	add_definitions("-mmacosx-version-min=10.5")
-        	set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -mmacosx-version-min=10.5")
-    	else()
-            # supports the first version changed from OS X to macOS
-        	add_definitions("--target=x86_64-apple-macos10.12")
-            set(CMAKE_OSX_DEPLOYMENT_TARGET 10.12)
-            set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} --target=x86_64-apple-macos10.12")
-    	endif()
+        # supports the first version changed from OS X to macOS
+        add_definitions("--target=x86_64-apple-macos11.0.1")
+        set(CMAKE_OSX_DEPLOYMENT_TARGET 11.0.1)
+        set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} --target=x86_64-apple-macos11.0.1")
 	endif()
     SET(CMAKE_FIND_LIBRARY_SUFFIXES .a ${CMAKE_FIND_LIBRARY_SUFFIXES})
 elseif (UNIX)
@@ -1017,7 +1048,8 @@ endif()
 
 # basic linking librararies
 target_link_libraries(iqtree3 pll ncl nclextra utils pda lbfgsb whtest sprng vectorclass model
-    gsl alignment tree simulator terrace yaml-cpp phyloYAML main ${TARGET_CMAPLE} ${PLATFORM_LIB} ${STD_LIB} ${THREAD_LIB} ${ATOMIC_LIB})
+    gsl alignment tree simulator terrace yaml-cpp phyloYAML main ${TARGET_CMAPLE} ${PLATFORM_LIB} ${STD_LIB} ${THREAD_LIB} ${ATOMIC_LIB} ${CMAKE_DL_LIBS})
+
 
 # Link the bundled LLVM libomp.a even for GCC builds
 # dl is needed by libomp on Linux for dlsym

alignment/alignment.cpp

@@ -4102,6 +4102,15 @@ void Alignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern_freq
                         memcpy(state_freq, aln->ptn_state_freq[ptn], num_states*sizeof(double));
                     }
                     ptn_state_freq.push_back(state_freq);
+                    if (!aln->site_rate_matrices.empty()) {
+                        /* Minh/Thomas: Better change 190 to num_states*(numstates-1)/2 so that if you want to
+                         extend the model in the future, no change is needed here
+                         Also: This is only for reversible models. For non-rev models
+                         you need to store the full matrix, i.e., num_states*num_states entries
+                         */
+                        const double *rate_matrix = aln->site_rate_matrices.data() + ptn * 190;
+                        site_rate_matrices.insert(site_rate_matrices.end(), rate_matrix, rate_matrix + 190);
+                    }
                 }
                 if (pattern_freq) {
                     ((*pattern_freq)[ptn])++;
@@ -4186,6 +4195,9 @@ void Alignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern_freq
     		out_site += site_vec[part+1];
     	}
     }
+    if (!aln->site_rate_matrices.empty()) {
+        ASSERT(aln->site_rate_matrices.size() == aln->getNPattern() * 190);
+    }
     if (aln->isSSF()) {
         ASSERT(ptn_state_freq.size() == getNPattern());
     }

alignment/alignment.h

@@ -949,6 +949,10 @@ class Alignment : public vector<Pattern>, public CharSet, public StateSpace {
     /** pattern index to state frequency vector map */
     vector<double*> ptn_state_freq;
 
+    /** site to rate matrix. Stored in row-major order [num_sites, 190] used for MUTSEL */
+    vector<double> site_rate_matrices;
+    // Minh/Thomas: TODO rename this to ptn_rate_matrices to avoid confusions
+
     /**
      * @return true if data type is SEQ_CODON and state is a stop codon
      */

main/CMakeLists.txt

@@ -14,6 +14,7 @@ add_library(main
     terraceanalysis.h
 )
 
+
 if (USE_BOOSTER)
     target_link_libraries(main booster)
 endif()
@@ -58,4 +59,4 @@ if (USE_CMAPLE)
 
         target_link_libraries(main-aa pda whtest vectorclass terrace maple-aa)
     endif()
-endif()
+endif()

main/main.cpp

@@ -64,6 +64,7 @@
 #include "utils/operatingsystem.h" //for getOSName()
 #include <stdlib.h>
 #include "vectorclass/instrset.h"
+#include "main/outstreambuf.h"
 #include "alignment/alignment.h"
 
 #include "utils/MPIHelper.h"
@@ -1686,30 +1687,6 @@ void processNCBITree(Params &params) {
     }
 }
 
-/* write simultaneously to cout/cerr and a file */
-class outstreambuf : public streambuf {
-public:
-    outstreambuf* open( const char* name, ios::openmode mode = ios::out);
-    bool is_open();
-    outstreambuf* close();
-    ~outstreambuf() { close(); }
-    streambuf *get_fout_buf() {
-        return fout_buf;
-    }
-    streambuf *get_cout_buf() {
-        return cout_buf;
-    }
-    ofstream *get_fout() {
-        return &fout;
-    }
-
-protected:
-    ofstream fout;
-    streambuf *cout_buf;
-    streambuf *fout_buf;
-    virtual int     overflow( int c = EOF);
-    virtual int     sync();
-};
 
 outstreambuf* outstreambuf::open( const char* name, ios::openmode mode) {
     if (!(Params::getInstance().suppress_output_flags & OUT_LOG)) {

main/outstreambuf.h

@@ -0,0 +1,33 @@
+#ifndef OUTSTREAMBUF_H
+#define OUTSTREAMBUF_H
+
+#include <fstream>
+#include <streambuf>
+#include <iostream>
+
+/* write simultaneously to cout/cerr and a file */
+class outstreambuf : public std::streambuf {
+public:
+    outstreambuf* open( const char* name, std::ios::openmode mode = std::ios::out);
+    bool is_open();
+    outstreambuf* close();
+    ~outstreambuf() { close(); }
+    std::streambuf *get_fout_buf() {
+        return fout_buf;
+    }
+    std::streambuf *get_cout_buf() {
+        return cout_buf;
+    }
+    std::ofstream *get_fout() {
+        return &fout;
+    }
+
+protected:
+    std::ofstream fout;
+    std::streambuf *cout_buf;
+    std::streambuf *fout_buf;
+    virtual int     overflow( int c = EOF);
+    virtual int     sync();
+};
+
+#endif // OUTSTREAMBUF_H