Repository: CitrineInformatics/lolo (language: python)

Issue:
Removing scalactic
Removing scalactic library from lolo as it is not being used. This transitive dependency is being pulled into the platform-backend causing a version clash.


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 8a3093d6..5d98e9d1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -11,3 +11,6 @@ tests/
 doc/
 .DS_Store
 foobar
+python/lolopy.egg-info/
+python/lolopy/__pycache__/
+python/lolopy/jar/
diff --git a/.travis.yml b/.travis.yml
index 3445d270..eae4fb47 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,49 +1,67 @@
+language: scala
+scala: 2.12.4
+jdk: openjdk8
+
 cache:
   directories:
   - "~/.m2/repository"
   - "~/.cache/pip"
+  - $HOME/.ivy2/cache
+
+before_cache:
+  # Cleanup the cached directories to avoid unnecessary cache updates
+  - find $HOME/.ivy2/cache -name "ivydata-*.properties" -print -delete
+
+before_install:
+  - git fetch --tags
+
+stages:
+  - name: test
+  - name: release
+    if: ((branch = master AND type = push) OR (tag IS present)) AND NOT fork
+
 jobs:
   include:
   - stage: test
-    language: scala
-    scala: 2.12.4
-    jdk: openjdk8
-    install:
-    - mvn install -P !build-extras -DskipTests=true -Dmaven.javadoc.skip=true -B
     script:
-    - mvn test -P !build-extras -B
+    - sbt +test
   - stage: test
-    language: python
-    python:
-    - 3.6
+    before_install:
+      - - wget http://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh
+        - chmod +x miniconda.sh
+        - ./miniconda.sh -b
+        - export PATH=/home/travis/miniconda3/bin:$PATH
+        - conda update --yes conda
+        - conda create --yes -n condaenv python=3.6
+        - conda install --yes -n condaenv pip
+        - source activate condaenv
     install:
     - cd python
     - make
     - pip install --only-binary=numpy,scipy -r requirements.txt
+    - pip install -r test_requirements.txt
     script:
     - nosetests
-  - stage: deploy
+  - stage: release
     name: maven central
-    language: scala
-    scala: 2.12.4
-    jdk: openjdk8
-    script: true
-    deploy:
-      provider: script
-      script: bash cd/deploy.sh
-      on:
-        branch: master
-  - stage: deploy
+    script: sbt ci-release
+  - stage: release
     name: pypi
-    language: python
-    python: 3.6
     addons:
       apt:
         packages:
         - pandoc
     before_install:
-    - cd python
+      - wget http://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh
+      - chmod +x miniconda.sh
+      - ./miniconda.sh -b
+      - export PATH=/home/travis/miniconda3/bin:$PATH
+      - conda update --yes conda
+      - conda create --yes -n condaenv python=3.6
+      - conda install --yes -n condaenv pip
+      - source activate condaenv
     install:
+    - cd python
     - make
     - pip install pypandoc
     - pip install --only-binary=numpy,scipy -r requirements.txt
@@ -56,10 +74,4 @@ jobs:
         secure: fb/NiAH7E1zA2iRUT+NMywH7dh0HF3qLM/JhjWLQWQ8l/pxOy7HPgk1CDEkZip/4emH2ur5+bGAD0pm0ewfDWfvjTPA+woV5IEebK2KL3Gm4Oam4fZP27NZ5zfYs6Q1dx3YPjRUgI3gtJ+0Y1tVKvZrUaZHityDMJDuEsYjdmPnFsA6s6U/5GYeQfWaYjZcflqYd01H8K9kbCeBsRBlOUm1lCXmt13R4uBOoMNAkuR7knYUOE1VM6VJWTN2T0iKphS9agxPHh/9/B9gjCQExhynOWSs5E2WhhDnWUoAJgyaZdxZQOxT2jJKA4dJGMCtyQScWItjsglcpthd+DsNtx/0vF1fDV3tvdjRVfegnhcF9fMUeI30O5jGTqNhimkiQsV1L+Sn6LFVAJKZ2yGzjqyS+8wn/uGOLQ04R46aw7KcfDJthzPU+dx0W78scTnkOnOx0R6eLH/HBn619h415JSsYpQ/D4A2VzODhHfScgzEe9xoL61ArhYz/wAqEAvit8WbFs4D7ZVnBc/98wsBD2NKUlmTZfkjPDAZQr49S1Na6UKlK6p4gLkXE6hfwgkFwxLiaQxOMfg0u/7VbbJzfrB2568HerS/oilDq+Bp+jNKUwyo08P0qD+8a4Gif1uk1UcipixDWbCvNlfr1H7+2FrN0+b92Xmfjd+CTMKmfv7o=
       distributions: sdist bdist_wheel
       on:
-        branch: master
-env:
-  global:
-  - secure: KXYC2UsqWnM+6ghbVviBC4J7tin4RDH+Ylgihk107h59s4osg/dFmoeZs0qmfhcFhauqT/+wSSaSbEqv1XCg79yv4JuYvlnLosQI7kJsNh0sjqnwNpw9qvvwWGph9ZCKskUMU27f0J57bdBi37spXGApL6I7vYia/0e3grqERqWp+p68CbD7H72Zb4V0j3WkU+zURFZbjl8D5EmIgkb/My9olYq9QPblYj79o+5LfJRO8upVuGPO8uUsbxtx4x9T/blwKd4tXULFNiRF067qfTxuemKdbCuZ8fcg7VDJBeMj6wXVQRcY6Nvfcyi2mJh0RVYWHvdHabzsHWlwAJr0tQ0rv1V7jv1CH2dwyL6fZ+jWU9dB47fhr4tg3cpRqFC4a/rRNabbWvwk2+SuWC33fDIlq68fpqpOQ5113ePqTIrKrnYybyzpjEHEdBQgxcQC6A0Ni0K8BiVWVphTxznt1XA0IYCptuNrvCCJq2c5WFXodsl52elkwLnyTOJ91yzymPRQYr4c3OhXmVi7E456qCGZuZ8lHDlYHcftC5Fsft385CxurZ+RfL3xH+DP0Xj4m27RCQiiXllK6TGCXnc9L7NlcuIwx3zu7ysEBtWWRtSdCnvR550FLf+iUuBIgjm4uNTuaTjlD3ddK7rK/fU08PPpqq80pYfqDge6yTIFDJs=
-  - secure: ELoVdfDuQp+/M00FDy1MH6tDkpWCzEkx4fY3j7mP+DV3pBrgLlp84gm/v91At6hNrtLFzJTpBqmPecZuqdDYp9jLy2cSAEMM1dWAgSZ8/Xlz/aeGTFvM1jnXr423ErrnlCxyczxd7BNC1iulUDvqg+4z+C6nxM1a3CCNmly7kjFGo1Naq3YCWbkrjNMqZNW0Qg6nh0i1EpVvFGE2KJdoh78KEXvOCc8+AyLHKM1L8T41eucWYdxYXbSxMG9cJONDlO7jol4IyaCLqc7d1YMcF9LT7x0hHJL8Mr8I5CRRlI71FlytijrmMQgePQdYwPGoce3vzi2q8qNfOzM1AL46c9aHidz6tjYVxRNEko+KpZOe7bOyEYS1bITotEUO9TPiJX/puiac4frDfEWjafg63eWs6gsSDvmGlHkxxDmJY9i0sYYYKigakFEeRz2/vFlJ6bUTiKpiVYcaVMFrtaYyW+4LRnMybTZmYnc7aM30D6QnueXygyjLBPY2W3TKXiUV3tnfKwm+XPGxsJK8Rqiw4GUJbTP9ZREF4hjm0JtNYKqtUhnJslGIqlLpBXH/iu0Kfsfg7WfG8iEBUQ/v0u9zrMJGgqXl7I8uUBAtoYkph+o/ydVE8p8yO2KeNfRm6JVat9UN3FN/ubfD5qnNIRn0J1oIbWdhtr9fttN1ZznFA3I=
-  - secure: kDU4TblhVO2PmtZuI8jB/EHW6/mv5w5FenArY0AfhyoP5oF6YWA62oLTVvEvmtTaCZrt33RKeDyRaeJshCEMgt1kY1oWXpSzZLnBlnXDRj7yNzaVQcp9aSjZd/AijMwTuvG31tiqB+DXp4nsv7JP0pppQ0NGZAClrBkHIeJVWjeFxcacwffJLUxHjKvKJ39y12eraXN5ROYNxbIuv40n8hrtRgrxSNz0+1BBggT+AMgtFk6FhuEerulZ/SVTUgflB0MdvvjBLkq5Gzv9YqzOnEGE0iB5r+xIGzjIEQGtxXduQcj8lAlzTBWPcdDUKt1qrPJ6E6F5c+tNzyi1I6E4Lu+tW1z+ApNP7U0FuTYG0ddK/JLe0PCr1RMlaY4lUDBcQDRgTfXNNJ1uJTdbR6nCVUY19Diyu7Hg+idcMnRsV1B/GCoRjQT9HwuxzrSHplfoUNVVEU0toEYL3FOx+hF8wpBga3gixbRGVYtSUob89nu+0oAvvHAU1Mjz0VySPsnEUENXkL1qZqPvLIBJ1HNMb866BiboVkolx8dkIr4yqDPHc3Bm+MpC5rd9DTOF4aF5eNSnEBwCjI27EvV+2aQlEDleHqq0N2x6LuDpBYa2mWoLw91RiHW9+fYhmxxBGEppUKC+u/w8dCUqm9lHu6azKWw4Pw1AMKTJZ3/3gXEu+uI=
-  - secure: YUFtsHv7u1q/2Rt1ZYI7Pfog9O5itVdTuHcx6Mxo4NAok+Up0WLijmn4Bwr/VrXStkxxHsjALm6eG8cVO8HGtFQ1ANwsOq6Nbwm18xAqgRYvlJNKtEllSLJPvTuBM1uoCnZZWBXG2iAAq48KoPUiR7i08IUMvSdcTEtHYjJXLO8lRGdoaKlSNOYAqxP7WkOrCEoDNi8VfDKrzDeE1G+5YRPluRsEUk9FvsgQnbUvFsyjy5LHOZmXMsLFJ05+aRe9VN5DbzYGPCs2h47uJQaDZP/MvM/gPs5ypKt/glYSUOehbHeaxCr0qqCRPISDPxlSongDbo0odgd4dv+gruqubHQTBeHfnKwPNOVWvPMO+w787umkD1qP0C9AnuuaCaKMHN7q/j354i/7avjo2fBVcsQz+Mgg/7VBDKaeycBJM/i8+2lCcOvfk3bMXr6BXM1O9IA2TnuGuEP8crf56dGa8ntZYyFRgTEEUCdm+ffroFifhtRBchblX/1FSQIEl1brHOz5CIAhq3oPa+IUrtcWbKz9wkbO0qrsM7NxulhuBgEtqBSau8DWg/pc9yYCsQTr02KOawzC2brryat5nAykJqsgT95EYtca2gsSXg6VZsLZM+NDJDqYD9A0v6CQ5FmTCm64capxix+PpIfD3H2TcAzZlfNa8rfICp7VPBczoJ4=
+        tags: true
\ No newline at end of file
diff --git a/build.sbt b/build.sbt
new file mode 100644
index 00000000..00a99a7a
--- /dev/null
+++ b/build.sbt
@@ -0,0 +1,30 @@
+import Dependencies._
+
+name               := "lolo"
+scalaVersion       := "2.13.4"
+crossScalaVersions := List("2.13.4", "2.12.4")
+organization       := "io.citrine"
+organizationName   := "Citrine Informatics"
+homepage           := Some(url("https://github.com/CitrineInformatics/lolo"))
+developers         := List(Developer(id="maxhutch", name="Max Hutchinson", email="maxhutch@citrine.io", url=url("https://github.com/maxhutch")))
+description        := "A random forest-centered machine learning library in Scala."
+licenses           += "Apache-2.0" ->  url("http://www.apache.org/licenses/LICENSE-2.0.txt")
+scmInfo            := Some(ScmInfo(url("https://github.com/CitrineInformatics/lolo"), "scm:git@github.com:CitrineInformatics/lolo.git"))
+
+pomIncludeRepository := { _ => false }
+test in assembly := {}
+assemblyJarName in assembly := s"lolo-jar-with-dependencies.jar"
+
+libraryDependencies ++= loloDeps
+libraryDependencies ++= {
+  CrossVersion.partialVersion(scalaVersion.value) match {
+    case Some((2, major)) if major <= 12 =>
+      Seq()
+    case _ =>
+      Seq("org.scala-lang.modules" %% "scala-parallel-collections" % scalaParallelCollectionsVersion)
+  }
+}
+
+enablePlugins(SiteScaladocPlugin)
+enablePlugins(GhpagesPlugin)
+git.remoteRepo := "git@github.com:CitrineInformatics/lolo.git"
diff --git a/cd/codesigning.asc.enc b/cd/codesigning.asc.enc
deleted file mode 100644
index 7c72320a..00000000
Binary files a/cd/codesigning.asc.enc and /dev/null differ
diff --git a/cd/deploy.sh b/cd/deploy.sh
deleted file mode 100755
index 2a65987a..00000000
--- a/cd/deploy.sh
+++ /dev/null
@@ -1,6 +0,0 @@
-#!/usr/bin/env bash
-if [ "$TRAVIS_BRANCH" = 'master' ] && [ "$TRAVIS_PULL_REQUEST" == 'false' ]; then
-    openssl aes-256-cbc -K $encrypted_9c3484b2e90f_key -iv $encrypted_9c3484b2e90f_iv -in cd/codesigning.asc.enc -out cd/codesigning.asc -d
-    gpg --fast-import cd/codesigning.asc
-    mvn deploy -P sign,build-extras,lolopy --settings cd/mvnsettings.xml
-fi
diff --git a/cd/mvnsettings.xml b/cd/mvnsettings.xml
deleted file mode 100644
index 785a3d94..00000000
--- a/cd/mvnsettings.xml
+++ /dev/null
@@ -1,24 +0,0 @@
-<settings>
-  <servers>
-    <server>
-      <id>ossrh</id>
-      <username>${env.OSSRH_JIRA_USERNAME}</username>
-      <password>${env.OSSRH_JIRA_PASSWORD}</password>
-    </server>
-  </servers>
-  
-  <profiles>
-    <profile>
-      <id>ossrh</id>
-      <activation>
-        <activeByDefault>true</activeByDefault>
-      </activation>
-      <properties>
-        <gpg.executable>gpg</gpg.executable>
-        <gpg.keyname>${env.GPG_KEY_NAME}</gpg.keyname>
-        <gpg.passphrase>${env.GPG_PASSPHRASE}</gpg.passphrase>
-      </properties>
-
-    </profile>
-  </profiles>
-</settings>
diff --git a/pom.xml b/pom.xml
deleted file mode 100644
index 26f99d3c..00000000
--- a/pom.xml
+++ /dev/null
@@ -1,309 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<project xmlns="http://maven.apache.org/POM/4.0.0"
-         xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
-         xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/maven-v4_0_0.xsd">
-    <modelVersion>4.0.0</modelVersion>
-
-    <groupId>io.citrine</groupId>
-    <artifactId>lolo</artifactId>
-    <version>3.0.1</version>
-    <packaging>jar</packaging>
-
-    <name>${project.groupId}:${project.artifactId}</name>
-    <description>A random forest-centered machine learning library in Scala.</description>
-    <url>https://github.com/CitrineInformatics/lolo</url>
-
-    <licenses>
-        <license>
-            <name>The Apache License, Version 2.0</name>
-            <url>http://www.apache.org/licenses/LICENSE-2.0.txt</url>
-        </license>
-    </licenses>
-
-    <developers>
-        <developer>
-            <name>Max Hutchinson</name>
-            <email>maxhutch@citrine.io</email>
-            <organization>Citrine Informatics</organization>
-            <organizationUrl>https://www.citrine.io</organizationUrl>
-        </developer>
-    </developers>
-
-    <scm>
-        <connection>scm:git:git://github.com/CitrineInformatics/lolo.git</connection>
-        <developerConnection>scm:git:ssh://github.com:CitrineInformatics/lolo.git</developerConnection>
-        <url>http://github.com/CitrineInformatics/lolo/tree/master</url>
-    </scm>
-
-    <properties>
-        <project.build.sourceEncoding>UTF-8</project.build.sourceEncoding>
-        <project.reporting.outputEncoding>UTF-8</project.reporting.outputEncoding>
-    </properties>
-
-    <dependencies>
-        <dependency>
-            <groupId>com.github.fommil.netlib</groupId>
-            <artifactId>all</artifactId>
-            <version>1.1.2</version>
-            <type>pom</type>
-        </dependency>
-        <dependency>
-            <groupId>junit</groupId>
-            <artifactId>junit</artifactId>
-            <version>4.12</version>
-            <scope>test</scope>
-        </dependency>
-        <dependency>
-            <groupId>io.citrine</groupId>
-            <artifactId>theta</artifactId>
-            <version>1.0.0</version>
-        </dependency>
-        <dependency>
-            <groupId>org.scalanlp</groupId>
-            <artifactId>breeze_2.12</artifactId> <!-- or 2.11 -->
-            <version>0.13.2</version>
-        </dependency>
-        <dependency>
-            <groupId>org.knowm.xchart</groupId>
-            <artifactId>xchart</artifactId>
-            <version>3.5.2</version>
-        </dependency>
-        <dependency>
-            <groupId>org.scalactic</groupId>
-            <artifactId>scalactic_2.12</artifactId>
-            <version>3.0.5</version>
-        </dependency>
-        <dependency>
-            <groupId>org.scalatest</groupId>
-            <artifactId>scalatest_2.12</artifactId>
-            <version>3.0.5</version>
-            <scope>test</scope>
-        </dependency>
-        <dependency>
-            <groupId>org.scala-lang</groupId>
-            <artifactId>scala-library</artifactId>
-            <version>2.12.4</version>
-        </dependency>
-    </dependencies>
-
-    <build>
-        <plugins>
-            <plugin>
-                <groupId>org.apache.maven.plugins</groupId>
-                <artifactId>maven-surefire-plugin</artifactId>
-                <version>2.19.1</version>
-                <configuration>
-                    <argLine>-Xmx1024m</argLine>
-                    <parallel>methods</parallel>
-                    <threadCount>1</threadCount>
-                    <perCoreThreadCount>true</perCoreThreadCount>
-                    <excludedGroups>io.citrine.lolo.SlowTest</excludedGroups>
-                </configuration>
-            </plugin>
-            <plugin>
-                <groupId>org.apache.maven.plugins</groupId>
-                <artifactId>maven-compiler-plugin</artifactId>
-                <version>3.3</version>
-                <configuration>
-                    <source>1.8</source>
-                    <target>1.8</target>
-                </configuration>
-            </plugin>
-            <plugin>
-                <groupId>org.sonatype.plugins</groupId>
-                <artifactId>nexus-staging-maven-plugin</artifactId>
-                <version>1.6.7</version>
-                <extensions>true</extensions>
-                <configuration>
-                    <serverId>ossrh</serverId>
-                    <nexusUrl>https://oss.sonatype.org/</nexusUrl>
-                    <autoReleaseAfterClose>true</autoReleaseAfterClose>
-                </configuration>
-            </plugin>
-            <plugin>
-                <groupId>org.scalastyle</groupId>
-                <artifactId>scalastyle-maven-plugin</artifactId>
-                <version>0.8.0</version>
-                <executions>
-                    <execution>
-                        <id>validate</id>
-                        <phase>validate</phase>
-                        <configuration>
-                            <verbose>true</verbose>
-                            <failOnViolation>true</failOnViolation>
-                            <includeTestSourceDirectory>true</includeTestSourceDirectory>
-
-                            <failOnWarning>false</failOnWarning>
-
-                            <sourceDirectory>${basedir}/src/main/scala</sourceDirectory>
-                            <testSourceDirectory>${basedir}/src/test/scala</testSourceDirectory>
-                            <configLocation>
-                                https://raw.githubusercontent.com/CitrineInformatics/citrine-style/master/citrine-scala-style.xml
-                            </configLocation>
-                            <outputFile>${project.basedir}/scalastyle-output.xml</outputFile>
-                            <outputEncoding>UTF-8</outputEncoding>
-                        </configuration>
-                        <goals>
-                            <goal>check</goal>
-                        </goals>
-                    </execution>
-                </executions>
-            </plugin>
-            <plugin>
-                <!-- see http://davidb.github.com/scala-maven-plugin -->
-                <groupId>net.alchim31.maven</groupId>
-                <artifactId>scala-maven-plugin</artifactId>
-                <version>3.2.1</version>
-                <executions>
-                    <execution>
-                        <id>scala-compile-first</id>
-                        <phase>process-resources</phase>
-                        <goals>
-                            <goal>add-source</goal>
-                            <goal>compile</goal>
-                        </goals>
-                    </execution>
-                    <execution>
-                        <id>Compile</id>
-                        <goals>
-                            <goal>compile</goal>
-                            <goal>testCompile</goal>
-                        </goals>
-                    </execution>
-                    <execution>
-                        <id>attach-javadocs</id>
-                        <goals>
-                            <goal>doc-jar</goal>
-                        </goals>
-                    </execution>
-                </executions>
-            </plugin>
-        </plugins>
-    </build>
-
-    <profiles>
-        <profile>
-            <id>lolopy</id>
-            <build>
-                <plugins>
-                    <plugin>
-                        <groupId>org.apache.maven.plugins</groupId>
-                        <artifactId>maven-shade-plugin</artifactId>
-                        <version>2.4.3</version>
-                        <!-- http://maven.apache.org/plugins/maven-shade-plugin/examples/attached-artifact.html -->
-                        <configuration>
-                            <shadedArtifactAttached>true</shadedArtifactAttached>
-                            <shadedClassifierName>jar-with-dependencies</shadedClassifierName>
-                        </configuration>
-                        <executions>
-                            <execution>
-                                <phase>package</phase>
-                                <goals>
-                                    <goal>shade</goal>
-                                </goals>
-                            </execution>
-                        </executions>
-                    </plugin>
-                </plugins>
-            </build>
-        </profile>
-        <profile>
-            <id>sign</id>
-            <build>
-                <plugins>
-                    <plugin>
-                        <groupId>org.apache.maven.plugins</groupId>
-                        <artifactId>maven-gpg-plugin</artifactId>
-                        <version>1.6</version>
-                        <executions>
-                            <execution>
-                                <id>sign-artifacts</id>
-                                <phase>verify</phase>
-                                <goals>
-                                    <goal>sign</goal>
-                                </goals>
-                            </execution>
-                        </executions>
-                    </plugin>
-                </plugins>
-            </build>
-        </profile>
-        <profile>
-            <id>build-extras</id>
-            <activation>
-                <activeByDefault>true</activeByDefault>
-            </activation>
-            <build>
-                <plugins>
-                    <plugin>
-                        <groupId>org.apache.maven.plugins</groupId>
-                        <artifactId>maven-source-plugin</artifactId>
-                        <version>2.4</version>
-                        <executions>
-                            <execution>
-                                <id>attach-sources</id>
-                                <goals>
-                                    <goal>jar-no-fork</goal>
-                                </goals>
-                            </execution>
-                        </executions>
-                    </plugin>
-                    <plugin>
-                        <groupId>org.apache.maven.plugins</groupId>
-                        <artifactId>maven-javadoc-plugin</artifactId>
-                        <version>2.10.3</version>
-                        <executions>
-                            <execution>
-                                <id>attach-javadocs</id>
-                                <goals>
-                                    <goal>jar</goal>
-                                </goals>
-                            </execution>
-                        </executions>
-                    </plugin>
-
-                </plugins>
-            </build>
-        </profile>
-    </profiles>
-
-    <distributionManagement>
-        <snapshotRepository>
-            <id>ossrh</id>
-            <url>https://oss.sonatype.org/content/repositories/snapshots</url>
-        </snapshotRepository>
-        <repository>
-            <id>ossrh</id>
-            <url>https://oss.sonatype.org/service/local/staging/deploy/maven2/</url>
-        </repository>
-    </distributionManagement>
-
-    <repositories>
-        <repository>
-            <id>scala</id>
-            <name>Scala Tools</name>
-            <url>http://scala-tools.org/repo-releases/</url>
-            <releases>
-                <enabled>true</enabled>
-            </releases>
-            <snapshots>
-                <enabled>false</enabled>
-            </snapshots>
-        </repository>
-    </repositories>
-
-    <pluginRepositories>
-        <pluginRepository>
-            <id>scala</id>
-            <name>Scala Tools</name>
-            <url>http://scala-tools.org/repo-releases/</url>
-            <releases>
-                <enabled>true</enabled>
-            </releases>
-            <snapshots>
-                <enabled>false</enabled>
-            </snapshots>
-        </pluginRepository>
-    </pluginRepositories>
-
-</project>
diff --git a/project/Dependencies.scala b/project/Dependencies.scala
new file mode 100644
index 00000000..e2164cd7
--- /dev/null
+++ b/project/Dependencies.scala
@@ -0,0 +1,22 @@
+import sbt._
+
+object Dependencies {
+  lazy val netlibVersion = "1.1.2"
+  lazy val junitVersion = "4.13.1"
+  lazy val thetaVersion = "1.1.2"
+  lazy val breezeVersion = "1.1"
+  lazy val scalaTestVersion = "3.2.2"
+
+  // this is used in build.sbt so that it works with both 2.12 and 2.13
+  lazy val scalaParallelCollectionsVersion = "1.0.0"
+
+  val loloDeps = Seq(
+    "com.github.fommil.netlib" % "all" % netlibVersion,
+    "junit" % "junit" % junitVersion % "test",
+    "org.scalanlp" %% "breeze" % breezeVersion,
+    "io.citrine" %% "theta" % thetaVersion,
+    "org.scalatest" %% "scalatest" % scalaTestVersion % "test",
+    "com.novocode" % "junit-interface" % "0.11" % "test",
+    "org.knowm.xchart" % "xchart" % "3.5.2"
+  )
+}
diff --git a/project/build.properties b/project/build.properties
new file mode 100644
index 00000000..a919a9b5
--- /dev/null
+++ b/project/build.properties
@@ -0,0 +1,1 @@
+sbt.version=1.3.8
diff --git a/project/ghpages.sbt b/project/ghpages.sbt
new file mode 100644
index 00000000..350a9557
--- /dev/null
+++ b/project/ghpages.sbt
@@ -0,0 +1,2 @@
+addSbtPlugin("com.typesafe.sbt" % "sbt-ghpages" % "0.6.2")
+
diff --git a/project/plugins.sbt b/project/plugins.sbt
new file mode 100644
index 00000000..9ae90e80
--- /dev/null
+++ b/project/plugins.sbt
@@ -0,0 +1,3 @@
+addSbtPlugin("com.geirsson" % "sbt-ci-release" % "1.5.3")
+addSbtPlugin("com.typesafe.sbt" % "sbt-site" % "1.3.2")
+addSbtPlugin("com.eed3si9n" % "sbt-assembly" % "0.14.10")
diff --git a/python/Makefile b/python/Makefile
index 1e315233..03f911ca 100644
--- a/python/Makefile
+++ b/python/Makefile
@@ -1,3 +1,3 @@
 all:
-	cd ..; mvn -B -P lolopy clean package -DskipTests=true -Dmaven.javadoc.skip=True
+	cd ..; sbt assembly
 	pip install -e .
diff --git a/python/README.md b/python/README.md
index 78f89108..2bf31e3f 100644
--- a/python/README.md
+++ b/python/README.md
@@ -19,7 +19,7 @@ Lolo is a Scala library that contains a variety of machine learning algorithms,
 
 ### Development 
 
-Lolopy requires Python >= 3.6, Java JDK >= 1.8, and Maven to be installed on your system when developing lolopy.  
+Lolopy requires Python >= 3.6, Java JDK >= 1.8, and sbt to be installed on your system when developing lolopy.  
 
 Before developing `lolopy`, compile `lolo` on your system using Maven.
 We have provided a `Makefile` that contains the needed operations.
diff --git a/python/lolopy/learners.py b/python/lolopy/learners.py
index 9e724a47..1882c68b 100644
--- a/python/lolopy/learners.py
+++ b/python/lolopy/learners.py
@@ -3,10 +3,11 @@
 import numpy as np
 from lolopy.loloserver import get_java_gateway
 from lolopy.utils import send_feature_array, send_1D_array
+import random
 from sklearn.base import BaseEstimator, RegressorMixin, ClassifierMixin, is_regressor
 from sklearn.exceptions import NotFittedError
 
-__all__ = ['RandomForestRegressor', 'RandomForestClassifier']
+__all__ = ['RandomForestRegressor', 'RandomForestClassifier', 'ExtraRandomTreesRegressor', 'ExtraRandomTreesClassifier']
 
 
 class BaseLoloLearner(BaseEstimator, metaclass=ABCMeta):
@@ -253,11 +254,11 @@ class in Lolo"""
     def __init__(self, num_trees=-1, use_jackknife=True, bias_learner=None,
                  leaf_learner=None, subset_strategy="auto", min_leaf_instances=1,
                  max_depth=2**30, uncertainty_calibration=False, randomize_pivot_location=False,
-                 randomly_rotate_features=False):
+                 randomly_rotate_features=False, random_seed=None):
         """Initialize the RandomForest
 
         Args:
-            num_trees (int): Number of trees to use in the forest
+            num_trees (int): Number of trees to use in the forest (default of -1 sets the number of trees to the number of training rows)
             use_jackknife (bool): Whether to use jackknife based variance estimates
             bias_learner (BaseLoloLearner): Algorithm used to model bias (default: no model)
             leaf_learner (BaseLoloLearner): Learner used at each leaf of the random forest (default: GuessTheMean)
@@ -273,6 +274,7 @@ def __init__(self, num_trees=-1, use_jackknife=True, bias_learner=None,
             uncertainty_calibration (bool): whether to re-calibrate the predicted uncertainty based on out-of-bag residuals
             randomize_pivot_location (bool): whether to draw pivots randomly or always select the midpoint
             randomly_rotate_features (bool): whether to randomly rotate real features for each tree in the forest
+            random_seed (int): random number generator seed used for nondeterministic functionality
         """
         super(RandomForestMixin, self).__init__()
 
@@ -287,8 +289,11 @@ def __init__(self, num_trees=-1, use_jackknife=True, bias_learner=None,
         self.uncertainty_calibration = uncertainty_calibration
         self.randomize_pivot_location = randomize_pivot_location
         self.randomly_rotate_features = randomly_rotate_features
+        self.random_seed = random_seed
 
     def _make_learner(self):
+        rng = self.gateway.jvm.scala.util.Random(self.random_seed) if self.random_seed else self.gateway.jvm.scala.util.Random()
+
         #  TODO: Figure our a more succinct way of dealing with optional arguments/Option values
         #  TODO: that ^^, please
         learner = self.gateway.jvm.io.citrine.lolo.learners.RandomForest(
@@ -304,7 +309,8 @@ def _make_learner(self):
             self.max_depth,
             self.uncertainty_calibration,
             self.randomize_pivot_location,
-            self.randomly_rotate_features
+            self.randomly_rotate_features,
+            rng
         )
         return learner
 
@@ -317,10 +323,89 @@ class RandomForestClassifier(BaseLoloClassifier, RandomForestMixin):
     """Random Forest model used for classification"""
 
 
+class ExtraRandomTreesMixIn(BaseLoloLearner):
+    """Extra Random Trees base class
+
+    Implements the _make_learner operation and the __init__ function with options specific to the ExtraRandomTrees
+    class in Lolo"""
+
+    def __init__(self, num_trees=-1, use_jackknife=False, bias_learner=None,
+                 leaf_learner=None, subset_strategy="auto", min_leaf_instances=1,
+                 max_depth=2**30, uncertainty_calibration=False, disable_bootstrap=True,
+                 randomly_rotate_features=False, random_seed=None):
+        """Initialize the ExtraRandomTrees ensemble
+
+        Args:
+            num_trees (int): Number of trees to use in the forest (default of -1 sets the number of trees to the number of training rows)
+            use_jackknife (bool): Whether to use jackknife based variance estimates (default: False)
+            bias_learner (BaseLoloLearner): Algorithm used to model bias (default: no model)
+            leaf_learner (BaseLoloLearner): Learner used at each leaf of the random forest (default: GuessTheMean)
+            subset_strategy (Union[string,int,float]): Strategy used to determine number of features used at each split
+                Available options:
+                    "auto": Use the default for lolo (all features for regression; classification not supported)
+                    "log2": Use the base 2 log of the number of features
+                    "sqrt": Use the square root of the number of features
+                    integer: Set the number of features explicitly
+                    float: Use a certain fraction of the features
+            min_leaf_instances (int): Minimum number of features used at each leaf
+            max_depth (int): Maximum depth to which to allow the decision trees to grow
+            uncertainty_calibration (bool): whether to re-calibrate the predicted uncertainty based on out-of-bag residuals
+            randomize_pivot_location (bool): whether to draw pivots randomly or always select the midpoint
+            disable_bootstrap (bool): whether to disable bootstrapping (default: true)
+            randomly_rotate_features (bool): whether to randomly rotate real features for each tree in the forest
+        """
+        super(ExtraRandomTreesMixIn, self).__init__()
+
+        # Store the variables
+        self.num_trees = num_trees
+
+        self.use_jackknife = use_jackknife
+        self.bias_learner = bias_learner
+        self.leaf_learner = leaf_learner
+        self.subset_strategy = subset_strategy
+        self.min_leaf_instances = min_leaf_instances
+        self.max_depth = max_depth
+        self.uncertainty_calibration = uncertainty_calibration
+        self.randomly_rotate_features = randomly_rotate_features
+        self.random_seed = random_seed
+
+    def _make_learner(self):
+        #  TODO: Figure our a more succinct way of dealing with optional arguments/Option values
+        #  TODO: that ^^, please
+
+        rng = self.gateway.jvm.scala.util.Random(self.random_seed) if self.random_seed else self.gateway.jvm.scala.util.Random()
+
+        learner = self.gateway.jvm.io.citrine.lolo.learners.ExtraRandomTrees(
+            self.num_trees, self.use_jackknife,
+            getattr(self.gateway.jvm.io.citrine.lolo.learners.ExtraRandomTrees,
+                    "$lessinit$greater$default$3")() if self.bias_learner is None
+            else self.gateway.jvm.scala.Some(self.bias_learner._make_learner()),
+            getattr(self.gateway.jvm.io.citrine.lolo.learners.ExtraRandomTrees,
+                    "$lessinit$greater$default$4")() if self.leaf_learner is None
+            else self.gateway.jvm.scala.Some(self.leaf_learner._make_learner()),
+            self.subset_strategy,
+            self.min_leaf_instances,
+            self.max_depth,
+            self.uncertainty_calibration,
+            self.disable_bootstrap,
+            self.randomly_rotate_features,
+            rng
+        )
+        return learner
+
+
+class ExtraRandomTreesRegressor(BaseLoloRegressor, RandomForestMixin):
+    """Random Forest model used for regression"""
+
+
+class ExtraRandomTreesClassifier(BaseLoloClassifier, RandomForestMixin):
+    """Random Forest model used for regression"""
+
+
 class RegressionTreeLearner(BaseLoloRegressor):
     """Regression tree learner, based on the RandomTree algorithm."""
 
-    def __init__(self, num_features=-1, max_depth=30, min_leaf_instances=1, leaf_learner=None):
+    def __init__(self, num_features=-1, max_depth=30, min_leaf_instances=1, leaf_learner=None, random_seed=None):
         """Initialize the learner
 
         Args:
@@ -334,6 +419,7 @@ def __init__(self, num_features=-1, max_depth=30, min_leaf_instances=1, leaf_lea
         self.max_depth = max_depth
         self.min_leaf_instances = min_leaf_instances
         self.leaf_learner = leaf_learner
+        self.random_seed = random_seed
 
     def _make_learner(self):
         if self.leaf_learner is None:
@@ -351,9 +437,11 @@ def _make_learner(self):
             "$lessinit$greater$default$5"
         )()
 
+        rng = self.gateway.jvm.scala.util.Random() if self.random_seed is None else self.gateway.jvm.scala.util.Random(self.random_seed)
+
         return self.gateway.jvm.io.citrine.lolo.trees.regression.RegressionTreeLearner(
             self.num_features, self.max_depth, self.min_leaf_instances,
-            leaf_learner, splitter
+            leaf_learner, splitter, rng
         )
 
 class LinearRegression(BaseLoloRegressor):
diff --git a/python/lolopy/loloserver.py b/python/lolopy/loloserver.py
index eebe057c..acd9dec7 100644
--- a/python/lolopy/loloserver.py
+++ b/python/lolopy/loloserver.py
@@ -1,6 +1,5 @@
 """Methods related to starting and stopping the Java Gateway"""
 from py4j.java_gateway import JavaGateway
-from xml.etree import ElementTree
 import sys
 import os
 
@@ -28,12 +27,8 @@ def find_lolo_jar(skip_devel_version=False):
     """
 
     if not skip_devel_version and _is_development_installation():
-        # Get the version of lolo from the pom.xml file
-        pom = ElementTree.parse(os.path.join(_lolo_root, 'pom.xml'))
-        lolo_version = pom.find('{http://maven.apache.org/POM/4.0.0}version').text
-
         # Get the appropriate Jar
-        jar_path = os.path.join(_lolo_root, 'target', 'lolo-{}-jar-with-dependencies.jar'.format(lolo_version))
+        jar_path = os.path.join(_lolo_root, 'target', 'scala-2.13', 'lolo-jar-with-dependencies.jar')
         if not os.path.isfile(jar_path):
             raise RuntimeError('Current version of lolo jar not found. Try re-building project with make')
 
diff --git a/python/lolopy/metrics.py b/python/lolopy/metrics.py
index 2f85ac18..77bc5c87 100644
--- a/python/lolopy/metrics.py
+++ b/python/lolopy/metrics.py
@@ -8,7 +8,7 @@
 import numpy as np
 
 
-def _call_lolo_merit(metric_name, y_true, y_pred, y_std=None, *args):
+def _call_lolo_merit(metric_name, y_true, y_pred, rng, y_std=None, *args):
     """Call a metric from lolopy
     
     Args:
@@ -16,7 +16,7 @@ def _call_lolo_merit(metric_name, y_true, y_pred, y_std=None, *args):
         y_true ([double]): True value
         y_pred ([double]): Predicted values
         y_uncert ([double]): Prediction uncertainties
-        *args (list): Any paramters to the constructor of the Metric
+        *args (list): Any parameters to the constructor of the Metric
     Returns:
         (double): Metric score
     """
@@ -25,8 +25,10 @@ def _call_lolo_merit(metric_name, y_true, y_pred, y_std=None, *args):
     if y_std is None:
         y_std = np.ones(len(y_true))
 
-    # Get the metric object
     gateway = get_java_gateway()
+    # Get default rng
+    rng = rng if rng else gateway.jvm.scala.util.Random()
+    # Get the metric object
     metric = getattr(gateway.jvm.io.citrine.lolo.validation, metric_name)
     if len(args) > 0:
         metric = metric(*args)
@@ -40,10 +42,10 @@ def _call_lolo_merit(metric_name, y_true, y_pred, y_std=None, *args):
     pred_result = gateway.jvm.io.citrine.lolo.util.LoloPyDataLoader.makeRegressionPredictionResult(y_pred_java, y_std_java)
 
     # Run the prediction result through the metric
-    return metric.evaluate(pred_result, y_true_java)
+    return metric.evaluate(pred_result, y_true_java, rng)
 
 
-def root_mean_squared_error(y_true, y_pred):
+def root_mean_squared_error(y_true, y_pred, rng=None):
     """Compute the root mean squared error
     
     Args:
@@ -53,10 +55,10 @@ def root_mean_squared_error(y_true, y_pred):
         (double): RMSE
     """
 
-    return _call_lolo_merit('RootMeanSquareError', y_true, y_pred)
+    return _call_lolo_merit('RootMeanSquareError', y_true, y_pred, rng)
 
 
-def standard_confidence(y_true, y_pred, y_std):
+def standard_confidence(y_true, y_pred, y_std, rng=None):
     """Fraction of entries that have errors within the predicted confidence interval. 
     
     Args:
@@ -67,10 +69,10 @@ def standard_confidence(y_true, y_pred, y_std):
         (double): standard confidence
     """
 
-    return _call_lolo_merit('StandardConfidence', y_true, y_pred, y_std)
+    return _call_lolo_merit('StandardConfidence', y_true, y_pred, rng, y_std)
 
 
-def standard_error(y_true, y_pred, y_std, rescale=1.0):
+def standard_error(y_true, y_pred, y_std, rng=None, rescale=1.0):
     """Root mean square of the error divided by the predicted uncertainty 
     
     Args:
@@ -82,10 +84,10 @@ def standard_error(y_true, y_pred, y_std, rescale=1.0):
         (double): standard error
     """
 
-    return _call_lolo_merit('StandardError', y_true, y_pred, y_std, float(rescale))
+    return _call_lolo_merit('StandardError', y_true, y_pred, rng, y_std, float(rescale))
 
 
-def uncertainty_correlation(y_true, y_pred, y_std):
+def uncertainty_correlation(y_true, y_pred, y_std, rng=None):
     """Measure of the correlation between the predicted uncertainty and error magnitude
     
     Args:
@@ -95,7 +97,4 @@ def uncertainty_correlation(y_true, y_pred, y_std):
     Returns:
         (double):
     """
-
-    return _call_lolo_merit('UncertaintyCorrelation', y_true, y_pred, y_std)
-
-
+    return _call_lolo_merit('UncertaintyCorrelation', y_true, y_pred, rng, y_std)
diff --git a/python/lolopy/version.py b/python/lolopy/version.py
index 013fa1d4..b736fa8e 100644
--- a/python/lolopy/version.py
+++ b/python/lolopy/version.py
@@ -1,3 +1,3 @@
 # single source of truth for package version,
 # see https://packaging.python.org/en/latest/single_source_version/
-__version__ = "1.1.1"
+__version__ = "1.2.0"
diff --git a/python/setup.py b/python/setup.py
index ff13257a..f097797d 100644
--- a/python/setup.py
+++ b/python/setup.py
@@ -11,11 +11,11 @@
 version = version_ns['__version__']
 
 # Find the lolo jar
-JAR_FILE = glob(os.path.join('..', 'target', 'lolo-*-jar-with-dependencies.jar'))
+JAR_FILE = glob(os.path.join('..', 'target', 'scala-2.13', 'lolo-jar-with-dependencies.jar'))
 if len(JAR_FILE) == 0:
-    raise Exception('No Jar files found. Build lolo first by calling "make" or "cd ..; mvn -P lolopy clean package"')
+    raise Exception('No Jar files found. Build lolo first by calling "make" or "cd ..; sbt assembly"')
 elif len(JAR_FILE) > 1:
-    raise Exception('Found >1 Jar file. Clean and rebuild lolopy: cd ..; mvn -P lolopy clean package')
+    raise Exception('Found >1 Jar file. Clean and rebuild lolopy: cd ..; sbt assembly')
 
 # Copy the jar file to a directory at the same level as the package
 jar_path = os.path.join('lolopy', 'jar')
diff --git a/src/main/scala/io/citrine/lolo/Model.scala b/src/main/scala/io/citrine/lolo/Model.scala
index dea41dad..794e3915 100644
--- a/src/main/scala/io/citrine/lolo/Model.scala
+++ b/src/main/scala/io/citrine/lolo/Model.scala
@@ -1,5 +1,7 @@
 package io.citrine.lolo
 
+import breeze.linalg.DenseMatrix
+
 /**
   * Created by maxhutch on 11/14/16.
   */
@@ -9,7 +11,22 @@ trait Model[+T <: PredictionResult[Any]] extends Serializable {
     * Apply the model to a seq of inputs
     *
     * @param inputs to apply the model to
-    * @return a predictionresult which includes, at least, the expected outputs
+    * @return a PredictionResult which includes, at least, the expected outputs
     */
   def transform(inputs: Seq[Vector[Any]]): T
+
+  /**
+    * Compute Shapley feature attributions for a given input in this node's subtree
+    *
+    * @param input for which to compute feature attributions.
+    * @param omitFeatures feature indices to omit in computing Shapley values
+    * @return matrix of attributions for each feature and output
+    *         One row per feature, each of length equal to the output dimension.
+    *         The output dimension is 1 for single-task regression, or equal to the number of classification categories.
+    */
+  def shapley(
+               input: Vector[Any],
+               omitFeatures: Set[Int] = Set()
+             ): Option[DenseMatrix[Double]] = None
+
 }
diff --git a/src/main/scala/io/citrine/lolo/TrainingResult.scala b/src/main/scala/io/citrine/lolo/TrainingResult.scala
index b9819ed9..07b525c0 100644
--- a/src/main/scala/io/citrine/lolo/TrainingResult.scala
+++ b/src/main/scala/io/citrine/lolo/TrainingResult.scala
@@ -1,5 +1,7 @@
 package io.citrine.lolo
 
+import breeze.linalg.DenseVector
+
 /**
   * Created by maxhutch on 12/4/16.
   */
diff --git a/src/main/scala/io/citrine/lolo/bags/BaggedResult.scala b/src/main/scala/io/citrine/lolo/bags/BaggedResult.scala
index f16fee7b..4821565b 100644
--- a/src/main/scala/io/citrine/lolo/bags/BaggedResult.scala
+++ b/src/main/scala/io/citrine/lolo/bags/BaggedResult.scala
@@ -45,7 +45,8 @@ case class BaggedSingleResult(
                                predictions: Seq[PredictionResult[Double]],
                                NibIn: Vector[Vector[Int]],
                                bias: Option[Double] = None,
-                               rescale: Double = 1.0
+                               rescale: Double = 1.0,
+                               disableBootstrap: Boolean = false
                              ) extends BaggedResult[Double] with RegressionResult {
   private lazy val treePredictions: Array[Double] = predictions.map(_.getExpected().head).toArray
 
@@ -62,13 +63,28 @@ case class BaggedSingleResult(
     treePredictions.map(x => Math.pow(x - expected, 2.0)).sum / (treePredictions.length - 1)
   }
 
-  override def getStdDevMean(): Option[Seq[Double]] = Some(Seq(stdDevMean))
+  override def getStdDevMean(): Option[Seq[Double]] = {
+    if (disableBootstrap) {
+      // If bootstrap is disabled, rescale is unity and treeVariance is our only option for UQ.
+      // Since it's not recalibrated, it's best considered to be a confidence interval of the underlying weak learner.
+      assert(rescale == 1.0)
+      Some(Seq(Math.sqrt(treeVariance)))
+    } else {
+      Some(Seq(stdDevMean))
+    }
+  }
 
-  override def getStdDevObs(): Option[Seq[Double]] = Some(Seq(stdDevObs))
+  override def getStdDevObs(): Option[Seq[Double]] = {
+    if (disableBootstrap) {
+      None
+    } else {
+      Some(Seq(stdDevObs))
+    }
+  }
 
   /**
-   * For the sake of parity, we were using this method
-   */
+    * For the sake of parity, we were using this method
+    */
   override def getUncertainty(observational: Boolean): Option[Seq[Any]] = {
     if (observational) {
       getStdDevObs()
@@ -119,7 +135,7 @@ case class BaggedSingleResult(
           tNotCount = tNotCount + 1
         }
       }
-      // Compute the infinitessimal jackknife estimate
+      // Compute the infinitesimal jackknife estimate
       val varIJ = Math.pow(cov / vecN.size, 2.0)
 
       if (tNotCount > 0) {
@@ -142,12 +158,13 @@ case class BaggedSingleResult(
 }
 
 case class BaggedClassificationResult(
-                                  predictions: Seq[PredictionResult[Any]]
-                                ) extends BaggedResult[Any] {
+                                       predictions: Seq[PredictionResult[Any]]
+                                     ) extends BaggedResult[Any] {
+  val predictionEnsemble = predictions.map{ p => p.getExpected() }
   lazy val expectedMatrix: Seq[Seq[Any]] = predictions.map(p => p.getExpected()).transpose
 
   lazy val expected: Seq[Any] = expectedMatrix.map(ps => ps.groupBy(identity).maxBy(_._2.size)._1).seq
-  lazy val uncertainty: Seq[Map[Any, Double]] = expectedMatrix.map(ps => ps.groupBy(identity).mapValues(_.size.toDouble / ps.size))
+  lazy val uncertainty: Seq[Map[Any, Double]] = expectedMatrix.map(ps => ps.groupBy(identity).mapValues(_.size.toDouble / ps.size).toMap)
 
   /**
    * Return the majority vote vote
@@ -174,7 +191,8 @@ case class BaggedMultiResult(
                          predictions: Seq[PredictionResult[Double]],
                          NibIn: Vector[Vector[Int]],
                          bias: Option[Seq[Double]] = None,
-                         rescale: Double = 1.0
+                         rescale: Double = 1.0,
+                         disableBootstrap: Boolean = false
                        ) extends BaggedResult[Double] with RegressionResult {
 
   /**
@@ -184,9 +202,24 @@ case class BaggedMultiResult(
     */
   override def getExpected(): Seq[Double] = expected
 
-  override def getStdDevObs(): Option[Seq[Double]] = Some(varObs.map{v => Math.sqrt(v)})
+  override def getStdDevObs(): Option[Seq[Double]] = {
+    if (disableBootstrap) {
+      None
+    } else {
+      Some(varObs.map { v => Math.sqrt(v) })
+    }
+  }
 
-  override def getStdDevMean(): Option[Seq[Double]] = Some(stdDevMean)
+  override def getStdDevMean(): Option[Seq[Double]] = {
+    if (disableBootstrap) {
+      // If bootstrap is disabled, rescale is unity and treeVariance is our only option for UQ.
+      // Since it's not recalibrated, it's best considered to be a confidence interval of the underlying weak learner.
+      assert(rescale == 1.0)
+      Some(varObs.map{v => Math.sqrt(v)})
+    } else {
+      Some(stdDevMean)
+    }
+  }
 
   /**
    * For the sake of parity, we were using this method
@@ -297,11 +330,11 @@ case class BaggedMultiResult(
     /* These operations are pulled out of the loop and extra-verbose for performance */
     val JMat = NibJ.t * predMat
     Async.canStop()
-    val JMat2 = JMat :* JMat * ((Nib.size - 1.0) / Nib.size)
+    val JMat2 = JMat *:* JMat * ((Nib.size - 1.0) / Nib.size)
     Async.canStop()
     val IJMat = NibIJ.t * predMat
     Async.canStop()
-    val IJMat2 = IJMat :* IJMat
+    val IJMat2 = IJMat *:* IJMat
     Async.canStop()
     val arg = IJMat2 + JMat2
     Async.canStop()
diff --git a/src/main/scala/io/citrine/lolo/bags/Bagger.scala b/src/main/scala/io/citrine/lolo/bags/Bagger.scala
index 67f81c36..717152a0 100644
--- a/src/main/scala/io/citrine/lolo/bags/Bagger.scala
+++ b/src/main/scala/io/citrine/lolo/bags/Bagger.scala
@@ -1,12 +1,13 @@
 package io.citrine.lolo.bags
 
-import breeze.stats.distributions.Poisson
+import breeze.linalg.DenseMatrix
+import breeze.stats.distributions.{Poisson, Rand, RandBasis}
 import io.citrine.lolo.stats.metrics.ClassificationMetrics
 import io.citrine.lolo.util.{Async, InterruptibleExecutionContext}
 import io.citrine.lolo.{Learner, Model, PredictionResult, RegressionResult, TrainingResult}
 
 import scala.collection.parallel.ExecutionContextTaskSupport
-import scala.collection.parallel.immutable.ParSeq
+import scala.collection.parallel.immutable.{ParRange, ParSeq}
 import scala.reflect._
 
 /**
@@ -15,15 +16,25 @@ import scala.reflect._
   * Created by maxhutch on 11/14/16.
   *
   * @param method  learner to train each model in the ensemble
-  * @param numBags number of models in the ensemble
-  */
+  * @param numBags number of base models to aggregate (default of -1 sets the number of models to the number of training rows)
+  * @param useJackknife whether to enable jackknife uncertainty estimate
+  * @param biasLearner learner to use for estimating bias
+  * @param uncertaintyCalibration whether to enable empirical uncertainty calibration
+  * @param disableBootstrap whether to disable bootstrap (useful when `method` implements its own randomization)
+  * @param randBasis breeze RandBasis to use for generating breeze random numbers
+  * */
 case class Bagger(
                    method: Learner,
                    numBags: Int = -1,
                    useJackknife: Boolean = true,
                    biasLearner: Option[Learner] = None,
-                   uncertaintyCalibration: Boolean = false
+                   uncertaintyCalibration: Boolean = false,
+                   disableBootstrap: Boolean = false,
+                   randBasis: RandBasis = Rand
                  ) extends Learner {
+  require(
+    !(uncertaintyCalibration && disableBootstrap),
+    "Options uncertaintyCalibration and disableBootstrap are incompatible. At most one may be set true.")
 
   /**
     * Draw with replacement from the training data for each model
@@ -66,22 +77,26 @@ case class Bagger(
     )
 
     /* Compute the number of instances of each training row in each training sample */
-    val dist = new Poisson(1.0)
-    val Nib: Vector[Vector[Int]] = Iterator.continually{
-      // Generate Poisson distributed weights, filtering out any that don't have the minimum required number
-      // of non-zero training weights
-      Iterator.continually {
-        Vector.fill(trainingData.size)(dist.draw())
-      }.filter(_.count(_ > 0) >= Bagger.minimumNonzeroWeightSize).take(actualBags).toVector
-    }.filter{nMat =>
-      lazy val noAlwaysPresentTrainingData = nMat.transpose.forall{vec => vec.contains(0)}
-      // Make sure that at least one learner is missing each training point
-      // This prevents a divide-by-zero error in the jackknife-after-bootstrap calcualtion
-      !useJackknife || noAlwaysPresentTrainingData
-    }.next()
+    val dist = new Poisson(1.0)(randBasis)
+    val Nib: Vector[Vector[Int]] = if (disableBootstrap) {
+      Vector.fill[Vector[Int]](actualBags)(Vector.fill[Int](trainingData.size)(1))
+    } else {
+      Iterator.continually{
+        // Generate Poisson distributed weights, filtering out any that don't have the minimum required number
+        // of non-zero training weights
+        Iterator.continually {
+          Vector.fill(trainingData.size)(dist.draw())
+        }.filter(_.count(_ > 0) >= Bagger.minimumNonzeroWeightSize).take(actualBags).toVector
+      }.filter{nMat =>
+        lazy val noAlwaysPresentTrainingData = nMat.transpose.forall{vec => vec.contains(0)}
+        // Make sure that at least one learner is missing each training point
+        // This prevents a divide-by-zero error in the jackknife-after-bootstrap calculation
+        !useJackknife || noAlwaysPresentTrainingData
+      }.next()
+    }
 
     /* Learn the actual models in parallel */
-    val parIterator = (0 until actualBags).par
+    val parIterator = new ParRange(0 until actualBags)
     parIterator.tasksupport = new ExecutionContextTaskSupport(InterruptibleExecutionContext())
     val (models, importances: ParSeq[Option[Vector[Double]]]) = parIterator.map { i =>
       // get weights
@@ -109,7 +124,7 @@ case class Bagger(
         Async.canStop()
         val model = new BaggedModel(oobModels, Nib.filter {
           _ (idx) == 0
-        }, useJackknife)
+        }, useJackknife, disableBootstrap = disableBootstrap)
         val predicted = model.transform(Seq(trainingData(idx)._1))
         val error = predicted.getExpected().head - trainingData(idx)._2.asInstanceOf[Double]
         val uncertainty = predicted match {
@@ -138,9 +153,9 @@ case class Bagger(
     if (biasLearner.isEmpty || oobErrors.isEmpty) {
       Async.canStop()
       if (isRegression) {
-        new BaggedTrainingResult(models.asInstanceOf[ParSeq[Model[PredictionResult[Double]]]], averageImportance, Nib, trainingData, useJackknife, None, ratio)
+        new BaggedTrainingResult(models.asInstanceOf[ParSeq[Model[PredictionResult[Double]]]], averageImportance, Nib, trainingData, useJackknife, None, ratio, disableBootstrap)
       } else {
-        new BaggedTrainingResult(models, averageImportance, Nib, trainingData, useJackknife, None, ratio)
+        new BaggedTrainingResult(models, averageImportance, Nib, trainingData, useJackknife, None, ratio, disableBootstrap)
       }
     } else {
       val biasTraining = oobErrors.map { case (f, e, u) =>
@@ -154,9 +169,9 @@ case class Bagger(
       Async.canStop()
 
       if (isRegression) {
-        new BaggedTrainingResult[Double](models.asInstanceOf[ParSeq[Model[PredictionResult[Double]]]], averageImportance, Nib, trainingData, useJackknife, Some(biasModel), ratio)
+        new BaggedTrainingResult[Double](models.asInstanceOf[ParSeq[Model[PredictionResult[Double]]]], averageImportance, Nib, trainingData, useJackknife, Some(biasModel), ratio, disableBootstrap)
       } else {
-        new BaggedTrainingResult[Any](models, averageImportance, Nib, trainingData, useJackknife, None, ratio)
+        new BaggedTrainingResult[Any](models, averageImportance, Nib, trainingData, useJackknife, None, ratio, disableBootstrap)
       }
     }
   }
@@ -169,15 +184,21 @@ class BaggedTrainingResult[+T : ClassTag](
                             trainingData: Seq[(Vector[Any], Any)],
                             useJackknife: Boolean,
                             biasModel: Option[Model[PredictionResult[T]]] = None,
-                            rescale: Double = 1.0
+                            rescale: Double = 1.0,
+                            disableBootstrap: Boolean = false
                           )
   extends TrainingResult {
 
   lazy val NibT = Nib.transpose
-  lazy val model = new BaggedModel[T](models, Nib, useJackknife, biasModel, rescale)
+  lazy val model = new BaggedModel[T](models, Nib, useJackknife, biasModel, rescale, disableBootstrap)
   lazy val rep = trainingData.find(_._2 != null).get._2
   lazy val predictedVsActual = trainingData.zip(NibT).flatMap { case ((f, l), nb) =>
-    val oob = models.zip(nb).filter(_._2 == 0)
+    val oob = if (disableBootstrap) {
+      models.zip(nb)
+    } else {
+      models.zip(nb).filter(_._2 == 0)
+    }
+
     if (oob.isEmpty || l == null || (l.isInstanceOf[Double] && l.asInstanceOf[Double].isNaN)) {
       Seq()
     } else {
@@ -207,7 +228,13 @@ class BaggedTrainingResult[+T : ClassTag](
 
   override def getPredictedVsActual(): Option[Seq[(Vector[Any], Any, Any)]] = Some(predictedVsActual)
 
-  override def getLoss(): Option[Double] = Some(loss)
+  override def getLoss(): Option[Double] = {
+    if (predictedVsActual.nonEmpty) {
+      Some(loss)
+    } else {
+      None
+    }
+  }
 }
 
 /**
@@ -221,10 +248,10 @@ class BaggedModel[+T: ClassTag](
                    Nib: Vector[Vector[Int]],
                    useJackknife: Boolean,
                    biasModel: Option[Model[PredictionResult[T]]] = None,
-                   rescale: Double = 1.0
+                   rescale: Double = 1.0,
+                   disableBootstrap: Boolean = false
                  ) extends Model[BaggedResult[T]] {
 
-
   /**
     * Apply each model to the outputs and wrap them up
     *
@@ -244,15 +271,43 @@ class BaggedModel[+T: ClassTag](
 
     val res = if (inputs.size == 1 && isRegression) {
       // In the special case of a single prediction on a real value, emit an optimized BaggedSingleResult
-      BaggedSingleResult(ensemblePredictions.map(_.asInstanceOf[PredictionResult[Double]]), Nib, bias.map(_.head), rescale)
+      BaggedSingleResult(ensemblePredictions.map(_.asInstanceOf[PredictionResult[Double]]), Nib, bias.map(_.head), rescale, disableBootstrap)
     } else if (isRegression) {
-      new BaggedMultiResult(ensemblePredictions.map(_.asInstanceOf[PredictionResult[Double]]), Nib, bias, rescale)
+      BaggedMultiResult(ensemblePredictions.map(_.asInstanceOf[PredictionResult[Double]]), Nib, bias, rescale, disableBootstrap)
     } else {
-      new BaggedClassificationResult(ensemblePredictions)
+      BaggedClassificationResult(ensemblePredictions)
     }
 
     res.asInstanceOf[BaggedResult[T]]
   }
+
+  /**
+    * Compute Shapley feature attributions for a given input
+    *
+    * @param input for which to compute feature attributions.
+    * @param omitFeatures feature indices to omit in computing Shapley values
+    * @return matrix of attributions for each feature and output
+    *         One row per feature, each of length equal to the output dimension.
+    *         The output dimension is 1 for single-task regression, or equal to the number of classification categories.
+    */
+  override def shapley(input: Vector[Any], omitFeatures: Set[Int] = Set()): Option[DenseMatrix[Double]] = {
+    val ensembleShapley = models.map(model => model.shapley(input, omitFeatures))
+    if (!ensembleShapley.head.isDefined) {
+      None
+    }
+    assert(ensembleShapley.forall(x => x.isDefined))
+
+    def sumReducer(a: Option[DenseMatrix[Double]],
+                   b: Option[DenseMatrix[Double]]): Option[DenseMatrix[Double]] = {
+      Some(a.get + b.get)
+    }
+    val scale = 1.0 / ensembleShapley.length
+
+    Some(scale * ensembleShapley.reduce(sumReducer).get)
+  }
+
+  // Accessor useful for testing.
+  private[bags] def getModels(): ParSeq[Model[PredictionResult[T]]] = models
 }
 
 object Bagger {
diff --git a/src/main/scala/io/citrine/lolo/bags/MultiTaskBagger.scala b/src/main/scala/io/citrine/lolo/bags/MultiTaskBagger.scala
index 7c20f0ff..a708a323 100644
--- a/src/main/scala/io/citrine/lolo/bags/MultiTaskBagger.scala
+++ b/src/main/scala/io/citrine/lolo/bags/MultiTaskBagger.scala
@@ -1,11 +1,11 @@
 package io.citrine.lolo.bags
 
-import breeze.stats.distributions.Poisson
+import breeze.stats.distributions.{Poisson, Rand, RandBasis}
 import io.citrine.lolo._
 import io.citrine.lolo.util.{Async, InterruptibleExecutionContext}
 
 import scala.collection.parallel.ExecutionContextTaskSupport
-import scala.collection.parallel.immutable.ParSeq
+import scala.collection.parallel.immutable.{ParRange, ParSeq}
 
 /**
   * Create an ensemble of multi-task models
@@ -18,7 +18,8 @@ case class MultiTaskBagger(
                             numBags: Int = -1,
                             useJackknife: Boolean = true,
                             biasLearner: Option[Learner] = None,
-                            uncertaintyCalibration: Boolean = false
+                            uncertaintyCalibration: Boolean = false,
+                            randBasis: RandBasis = Rand
                           ) extends MultiTaskLearner {
 
   private def combineImportance(v1: Option[Vector[Double]], v2: Option[Vector[Double]]): Option[Vector[Double]] = {
@@ -52,7 +53,7 @@ case class MultiTaskBagger(
     }
 
     /* Compute the number of instances of each training row in each training sample */
-    val dist = new Poisson(1.0)
+    val dist = new Poisson(1.0)(randBasis)
     val Nib = Vector.tabulate(actualBags) { _ =>
       Vector.tabulate(inputs.size) { _ =>
         dist.draw()
@@ -60,7 +61,7 @@ case class MultiTaskBagger(
     }
 
     /* Learn the actual models in parallel */
-    val parIterator = (0 until actualBags).par
+    val parIterator = new ParRange(0 until actualBags)
     parIterator.tasksupport = new ExecutionContextTaskSupport(InterruptibleExecutionContext())
     val (models: ParSeq[Seq[Model[PredictionResult[Any]]]], importances: ParSeq[Seq[Option[Vector[Double]]]]) = parIterator.map { i =>
       // Train the model
@@ -74,45 +75,50 @@ case class MultiTaskBagger(
      * Transpose the models and importances so the bags are the inner index and the labels are the outer index.
      * Foreach label, emit a BaggedTrainingResult
      */
-    models.transpose.zip(importances.seq.transpose).zipWithIndex.map { case ((m, i: Seq[Option[Vector[Double]]]), k) =>
-      val averageImportance: Option[Vector[Double]] = i.reduce {
-        combineImportance
-      }.map(_.map(_ / importances.size))
-      val trainingData = inputs.zip(labels(k))
-      val helper = BaggerHelper(m, trainingData, Nib, useJackknife, uncertaintyCalibration)
+    type arg = ((ParSeq[Model[PredictionResult[Any]]], Seq[Option[Vector[Double]]]), Int)
+    models.transpose.zip(importances.seq.transpose).zipWithIndex.map { xx: arg =>
+      xx match {
+        case ((m: ParSeq[Model[PredictionResult[Any]]], i: Seq[Option[Vector[Double]]]), k: Int) =>
+          val averageImportance: Option[Vector[Double]] = i.reduce {
+            combineImportance
+          }.map(_.map(_ / importances.size))
+          val trainingData = inputs.zip(labels(k))
+          val helper = BaggerHelper(m, trainingData, Nib, useJackknife, uncertaintyCalibration)
 
-      Async.canStop()
-      if (!helper.isRegression) {
-        new BaggedTrainingResult[Any](m, averageImportance, Nib, inputs.zip(labels(k)), useJackknife)
-      } else {
-        if (biasLearner.isEmpty) {
-          new BaggedTrainingResult[Double](m.asInstanceOf[ParSeq[Model[PredictionResult[Double]]]], averageImportance, Nib, trainingData, useJackknife, None, helper.ratio)
-        } else {
           Async.canStop()
-          val baggedModel = new BaggedModel[Double](m.asInstanceOf[ParSeq[Model[PredictionResult[Double]]]], Nib, useJackknife, None, helper.ratio)
-          Async.canStop()
-          val baggedRes = baggedModel.transform(trainingData.map(_._1))
-          Async.canStop()
-          val foo = baggedRes.getUncertainty()
-          val biasTraining = trainingData.zip(
-            baggedRes.getExpected().zip(baggedRes.getUncertainty().get)
-          ).flatMap { case ((f, a), (p, u)) =>
-            if (a == null || (a.isInstanceOf[Double] && a.asInstanceOf[Double].isNaN)) {
-              None
+          if (!helper.isRegression) {
+            new BaggedTrainingResult[Any](m, averageImportance, Nib, inputs.zip(labels(k)), useJackknife)
+          } else {
+            if (biasLearner.isEmpty) {
+              new BaggedTrainingResult[Double](m.asInstanceOf[ParSeq[Model[PredictionResult[Double]]]], averageImportance, Nib, trainingData, useJackknife, None, helper.ratio)
             } else {
-              // Math.E is only statistically correct.  It should be actualBags / Nib.transpose(i).count(_ == 0)
-              // Or, better yet, filter the bags that don't include the training example
-              val bias = Math.E * Math.max(Math.abs(p.asInstanceOf[Double] - a.asInstanceOf[Double]) - u.asInstanceOf[Double], 0.0)
-              Some((f, bias))
+              Async.canStop()
+              val baggedModel = new BaggedModel[Double](m.asInstanceOf[ParSeq[Model[PredictionResult[Double]]]], Nib, useJackknife, None, helper.ratio)
+              Async.canStop()
+              val baggedRes = baggedModel.transform(trainingData.map(_._1))
+              Async.canStop()
+              val foo = baggedRes.getUncertainty()
+              val biasTraining = trainingData.zip(
+                baggedRes.getExpected().zip(baggedRes.getUncertainty().get)
+              ).flatMap { case ((f, a), (p, u)) =>
+                if (a == null || (a.isInstanceOf[Double] && a.asInstanceOf[Double].isNaN)) {
+                  None
+                } else {
+                  // Math.E is only statistically correct.  It should be actualBags / Nib.transpose(i).count(_ == 0)
+                  // Or, better yet, filter the bags that don't include the training example
+                  val bias = Math.E * Math.max(Math.abs(p.asInstanceOf[Double] - a.asInstanceOf[Double]) - u.asInstanceOf[Double], 0.0)
+                  Some((f, bias))
+                }
+              }
+              Async.canStop()
+              val biasModel = biasLearner.get.train(helper.biasTraining).getModel().asInstanceOf[Model[PredictionResult[Double]]]
+              Async.canStop()
+
+              new BaggedTrainingResult[Double](m.asInstanceOf[ParSeq[Model[PredictionResult[Double]]]], averageImportance, Nib, trainingData, useJackknife, Some(biasModel), helper.ratio)
             }
           }
-          Async.canStop()
-          val biasModel = biasLearner.get.train(helper.biasTraining).getModel().asInstanceOf[Model[PredictionResult[Double]]]
-          Async.canStop()
-
-          new BaggedTrainingResult[Double](m.asInstanceOf[ParSeq[Model[PredictionResult[Double]]]], averageImportance, Nib, trainingData, useJackknife, Some(biasModel), helper.ratio)
-        }
       }
-    }.seq
+      }
+      .seq
   }
 }
diff --git a/src/main/scala/io/citrine/lolo/encoders/CategoricalEncoder.scala b/src/main/scala/io/citrine/lolo/encoders/CategoricalEncoder.scala
index 405788a7..a0bda2ed 100644
--- a/src/main/scala/io/citrine/lolo/encoders/CategoricalEncoder.scala
+++ b/src/main/scala/io/citrine/lolo/encoders/CategoricalEncoder.scala
@@ -10,7 +10,7 @@ package io.citrine.lolo.encoders
 class CategoricalEncoder[T](encoding: Map[T, Char]) extends Serializable {
 
   /** Inverse of the encoding */
-  lazy val decoding: Map[Char, T] = encoding.groupBy(_._2).mapValues(_.keys.head).map(identity)
+  lazy val decoding: Map[Char, T] = encoding.groupBy(_._2).mapValues(_.keys.head).toMap
 
   /**
     * Just call the encoding.  Use 0 for unknown inputs
diff --git a/src/main/scala/io/citrine/lolo/hypers/RandomHyperOptimizer.scala b/src/main/scala/io/citrine/lolo/hypers/RandomHyperOptimizer.scala
index 3f99cb4b..37dfce92 100644
--- a/src/main/scala/io/citrine/lolo/hypers/RandomHyperOptimizer.scala
+++ b/src/main/scala/io/citrine/lolo/hypers/RandomHyperOptimizer.scala
@@ -11,7 +11,7 @@ import scala.util.Random
   * calls.
   * Created by maxhutch on 12/7/16.
   */
-class RandomHyperOptimizer() extends HyperOptimizer {
+class RandomHyperOptimizer(rng: Random = Random) extends HyperOptimizer {
 
   /** Keep track of the best hypers outside of the optimize call so it persists across calls */
   var best: Map[String, Any] = Map()
@@ -29,7 +29,7 @@ class RandomHyperOptimizer() extends HyperOptimizer {
     /* Just draw numIteration times */
     (0 until numIterations).foreach { i =>
       val testHypers = hyperGrids.map { case (n, v) =>
-        n -> Random.shuffle(v).head
+        n -> rng.shuffle(v).head
       }
       val testLearner = builder(testHypers)
       val res = testLearner.train(trainingData)
diff --git a/src/main/scala/io/citrine/lolo/learners/ExtraRandomTrees.scala b/src/main/scala/io/citrine/lolo/learners/ExtraRandomTrees.scala
new file mode 100644
index 00000000..b7f3424b
--- /dev/null
+++ b/src/main/scala/io/citrine/lolo/learners/ExtraRandomTrees.scala
@@ -0,0 +1,125 @@
+package io.citrine.lolo.learners
+
+import breeze.stats.distributions.RandBasis
+import io.citrine.lolo.bags.Bagger
+import io.citrine.lolo.transformers.FeatureRotator
+import io.citrine.lolo.trees.classification.ClassificationTreeLearner
+import io.citrine.lolo.trees.regression.RegressionTreeLearner
+import io.citrine.lolo.trees.splits.{ExtraRandomClassificationSplitter, ExtraRandomRegressionSplitter}
+import io.citrine.lolo.{Learner, TrainingResult}
+import org.apache.commons.math3.random.MersenneTwister
+
+import scala.util.Random
+
+/**
+  * Extremely randomized tree ensemble
+  *
+  * This is based on Geurts, P., Ernst, D. & Wehenkel, L. Extremely randomized trees. Mach Learn 63, 3–42 (2006).
+  * https://doi.org/10.1007/s10994-006-6226-1.
+  *
+  * @param numTrees       number of trees to use (-1 => number of training instances)
+  * @param useJackknife   whether to use jackknife based variance estimates
+  * @param biasLearner    learner to model bias (absolute residual)
+  * @param leafLearner    learner to use at the leaves of the trees
+  * @param subsetStrategy for random feature selection at each split
+  *                       (auto => all features for regression)
+  * @param minLeafInstances minimum number of instances per leave in each tree
+  * @param maxDepth       maximum depth of each tree in the forest (default: unlimited)
+  * @param uncertaintyCalibration whether to empirically recalibrate the predicted uncertainties (default: false)
+  * @param disableBootstrap whether to disable bootstrap (default: true)
+  * @param randomlyRotateFeatures whether to randomly rotate real features for each tree in the forest (default: false)
+  * @param rng            random number generator to use
+  */
+case class ExtraRandomTrees(
+                             numTrees: Int = -1,
+                             useJackknife: Boolean = false,
+                             biasLearner: Option[Learner] = None,
+                             leafLearner: Option[Learner] = None,
+                             subsetStrategy: Any = "auto",
+                             minLeafInstances: Int = 1,
+                             maxDepth: Int = Integer.MAX_VALUE,
+                             uncertaintyCalibration: Boolean = false,
+                             disableBootstrap: Boolean = true,
+                             randomlyRotateFeatures: Boolean = false,
+                             rng: Random = Random
+                           ) extends Learner {
+    /**
+    * Train an extremely randomized tree ensemble model
+    *
+    * @param trainingData to train on
+    * @param weights      for the training rows, if applicable
+    * @return training result containing a model
+    */
+  override def train(trainingData: Seq[(Vector[Any], Any)], weights: Option[Seq[Double]]): TrainingResult = {
+    val breezeRandBasis: RandBasis = new RandBasis(new MersenneTwister(rng.nextLong()))
+
+    trainingData.head._2 match {
+      case _: Double => {
+        val numFeatures: Int = subsetStrategy match {
+          case x: String => {
+            x match {
+              case "auto" => trainingData.head._1.size
+              case "sqrt" => Math.ceil(Math.sqrt(trainingData.head._1.size)).toInt
+              case "log2" => Math.ceil(Math.log(trainingData.head._1.size) / Math.log(2)).toInt
+              case x: String =>
+                println(s"Unrecognized subsetStrategy ${x}; using auto")
+                trainingData.head._1.size
+            }
+          }
+          case x: Int => x
+          case x: Double => (trainingData.head._1.size * x).toInt
+        }
+        val DTLearner = RegressionTreeLearner(
+          leafLearner = leafLearner,
+          numFeatures = numFeatures,
+          minLeafInstances = minLeafInstances,
+          maxDepth = maxDepth,
+          splitter = ExtraRandomRegressionSplitter(rng)
+        )
+        val bagger = Bagger(
+          if (randomlyRotateFeatures) FeatureRotator(DTLearner) else DTLearner,
+          numBags = numTrees,
+          useJackknife = useJackknife,
+          biasLearner = biasLearner,
+          uncertaintyCalibration = uncertaintyCalibration,
+          disableBootstrap = disableBootstrap,
+          randBasis = breezeRandBasis
+        )
+        bagger.train(trainingData, weights)
+      }
+      case _: Any => {
+        val numFeatures: Int = subsetStrategy match {
+          case x: String => {
+            x match {
+              case "auto" => Math.ceil(Math.sqrt(trainingData.head._1.size)).toInt
+              case "sqrt" => Math.ceil(Math.sqrt(trainingData.head._1.size)).toInt
+              case "log2" => Math.ceil(Math.log(trainingData.head._1.size) / Math.log(2)).toInt
+              case x: String =>
+                println(s"Unrecognized subsetStrategy ${x}; using auto")
+                trainingData.head._1.size
+            }
+          }
+          case x: Int => x
+          case x: Double => (trainingData.head._1.size * x).toInt
+        }
+        val DTLearner = ClassificationTreeLearner(
+          leafLearner = leafLearner,
+          numFeatures = numFeatures,
+          minLeafInstances = minLeafInstances,
+          maxDepth = maxDepth,
+          splitter = ExtraRandomClassificationSplitter(rng)
+        )
+        val bagger = Bagger(
+          if (randomlyRotateFeatures) FeatureRotator(DTLearner) else DTLearner,
+          numBags = numTrees,
+          useJackknife = useJackknife,
+          biasLearner = biasLearner,
+          uncertaintyCalibration = uncertaintyCalibration,
+          disableBootstrap = disableBootstrap,
+          randBasis = breezeRandBasis
+        )
+        bagger.train(trainingData, weights)
+      }
+    }
+  }
+}
diff --git a/src/main/scala/io/citrine/lolo/learners/RandomForest.scala b/src/main/scala/io/citrine/lolo/learners/RandomForest.scala
index 6e129da7..68faeecb 100644
--- a/src/main/scala/io/citrine/lolo/learners/RandomForest.scala
+++ b/src/main/scala/io/citrine/lolo/learners/RandomForest.scala
@@ -1,11 +1,15 @@
 package io.citrine.lolo.learners
 
+import breeze.stats.distributions.{RandBasis, ThreadLocalRandomGenerator}
 import io.citrine.lolo.bags.Bagger
 import io.citrine.lolo.transformers.FeatureRotator
 import io.citrine.lolo.trees.classification.ClassificationTreeLearner
 import io.citrine.lolo.trees.regression.RegressionTreeLearner
 import io.citrine.lolo.trees.splits.{ClassificationSplitter, RegressionSplitter}
 import io.citrine.lolo.{Learner, TrainingResult}
+import org.apache.commons.math3.random.MersenneTwister
+
+import scala.util.Random
 
 /**
   * Standard random forest as a wrapper around bagged decision trees
@@ -22,6 +26,7 @@ import io.citrine.lolo.{Learner, TrainingResult}
   * @param uncertaintyCalibration whether to empirically recalibrate the predicted uncertainties (default: false)
   * @param randomizePivotLocation whether generate splits randomly between the data points (default: false)
   * @param randomlyRotateFeatures whether to randomly rotate real features for each tree in the forest (default: false)
+  * @param rng            random number generator to use for stochastic functionality
   */
 case class RandomForest(
                          numTrees: Int = -1,
@@ -33,8 +38,10 @@ case class RandomForest(
                          maxDepth: Int = Integer.MAX_VALUE,
                          uncertaintyCalibration: Boolean = false,
                          randomizePivotLocation: Boolean = false,
-                         randomlyRotateFeatures: Boolean = false
+                         randomlyRotateFeatures: Boolean = false,
+                         rng: Random = Random
                        ) extends Learner {
+
   /**
     * Train a random forest model
     *
@@ -46,6 +53,8 @@ case class RandomForest(
     * @return training result containing a model
     */
   override def train(trainingData: Seq[(Vector[Any], Any)], weights: Option[Seq[Double]]): TrainingResult = {
+    val breezeRandBasis: RandBasis = new RandBasis(new ThreadLocalRandomGenerator(new MersenneTwister(rng.nextLong())))
+
     trainingData.head._2 match {
       case _: Double =>
         val numFeatures: Int = subsetStrategy match {
@@ -68,13 +77,16 @@ case class RandomForest(
           numFeatures = numFeatures,
           minLeafInstances = minLeafInstances,
           maxDepth = maxDepth,
-          splitter = RegressionSplitter(randomizePivotLocation)
+          splitter = RegressionSplitter(randomizePivotLocation),
+          rng = rng
         )
+
         val bagger = Bagger(if (randomlyRotateFeatures) FeatureRotator(DTLearner) else DTLearner,
           numBags = numTrees,
           useJackknife = useJackknife,
           biasLearner = biasLearner,
-          uncertaintyCalibration = uncertaintyCalibration
+          uncertaintyCalibration = uncertaintyCalibration,
+          randBasis = breezeRandBasis
         )
         bagger.train(trainingData, weights)
       case _: Any =>
@@ -98,10 +110,13 @@ case class RandomForest(
           numFeatures = numFeatures,
           minLeafInstances = minLeafInstances,
           maxDepth = maxDepth,
-          splitter = ClassificationSplitter(randomizePivotLocation)
+          splitter = ClassificationSplitter(randomizePivotLocation),
+          rng = rng
         )
-        val bagger = Bagger(if (randomlyRotateFeatures) FeatureRotator(DTLearner) else DTLearner,
-          numBags = numTrees
+        val bagger = Bagger(
+          if (randomlyRotateFeatures) FeatureRotator(DTLearner) else DTLearner,
+          numBags = numTrees,
+          randBasis = breezeRandBasis
         )
         bagger.train(trainingData, weights)
     }
diff --git a/src/main/scala/io/citrine/lolo/linear/GuessTheMean.scala b/src/main/scala/io/citrine/lolo/linear/GuessTheMean.scala
index 67f02e44..57562d6e 100644
--- a/src/main/scala/io/citrine/lolo/linear/GuessTheMean.scala
+++ b/src/main/scala/io/citrine/lolo/linear/GuessTheMean.scala
@@ -7,7 +7,7 @@ import scala.util.Random
 /**
   * Created by maxhutch on 11/15/16.
   */
-case class GuessTheMeanLearner() extends Learner {
+case class GuessTheMeanLearner(rng: Random = Random) extends Learner {
 
   /**
     * Train a model
@@ -20,7 +20,7 @@ case class GuessTheMeanLearner() extends Learner {
     val data = trainingData.map(_._2).zip(weights.getOrElse(Seq.fill(trainingData.size)(1.0)))
     val mean = data.head._1 match {
       case x: Double => data.asInstanceOf[Seq[(Double, Double)]].map(p => p._1 * p._2).sum / data.map(_._2).sum
-      case x: Any => Random.shuffle(data.groupBy(_._1).mapValues(_.map(_._2).sum).toSeq).maxBy(_._2)._1
+      case x: Any => rng.shuffle(data.groupBy(_._1).mapValues(_.map(_._2).sum).toSeq).maxBy(_._2)._1
     }
 
     new GuessTheMeanTrainingResult(new GuessTheMeanModel(mean))
diff --git a/src/main/scala/io/citrine/lolo/stats/correlations/DistanceCorrelation.scala b/src/main/scala/io/citrine/lolo/stats/correlations/DistanceCorrelation.scala
index 829fa12a..415b2e2e 100644
--- a/src/main/scala/io/citrine/lolo/stats/correlations/DistanceCorrelation.scala
+++ b/src/main/scala/io/citrine/lolo/stats/correlations/DistanceCorrelation.scala
@@ -24,7 +24,7 @@ object DistanceCorrelation {
     }
     val grandMean = sum(pairDistances) / (x.size * x.size)
     val colSum: DenseVector[Double] = sum(pairDistances(*, ::))
-    val colMean: DenseVector[Double] = colSum :* (1.0 / x.size)
+    val colMean: DenseVector[Double] = colSum *:* (1.0 / x.size)
     DenseMatrix.tabulate(x.size, x.size) { case (i, j) =>
       pairDistances(i, j) - colMean(i) - colMean(j) + grandMean
     }
@@ -42,7 +42,7 @@ object DistanceCorrelation {
   def distanceCovariance[T](x: Seq[T], y: Seq[T], distance: (T, T) => Double): Double = {
     val A: DenseMatrix[Double] = doubleCenter(x, distance)
     val B: DenseMatrix[Double] = doubleCenter(y, distance)
-    Math.sqrt(sum(A :* B) / (x.size * x.size))
+    Math.sqrt(sum(A *:* B) / (x.size * x.size))
   }
 
   /**
diff --git a/src/main/scala/io/citrine/lolo/stats/functions/Linear.scala b/src/main/scala/io/citrine/lolo/stats/functions/Linear.scala
index 7d73bcc8..7a142af7 100644
--- a/src/main/scala/io/citrine/lolo/stats/functions/Linear.scala
+++ b/src/main/scala/io/citrine/lolo/stats/functions/Linear.scala
@@ -21,10 +21,10 @@ object Linear {
     * @param magnitude of the gradient (default: 1)
     * @return linear function with specified gradient magnitude in a random direction
     */
-  def randomDirection(nDim: Int, magnitude: Double = 1.0): Linear = {
+  def randomDirection(nDim: Int, magnitude: Double = 1.0, rng: Random = Random): Linear = {
     val gradient = {
       // Draw guassian so the direction is random rather than preferring corners
-      val unnormalized = Seq.fill(nDim){Random.nextGaussian()}
+      val unnormalized = Seq.fill(nDim){rng.nextGaussian()}
       val norm = Math.sqrt(unnormalized.map(Math.pow(_, 2)).sum)
       unnormalized.map(x => x * magnitude / norm)
     }
diff --git a/src/main/scala/io/citrine/lolo/stats/metrics/ClassificationMetrics.scala b/src/main/scala/io/citrine/lolo/stats/metrics/ClassificationMetrics.scala
index 8d162969..0d8c469b 100644
--- a/src/main/scala/io/citrine/lolo/stats/metrics/ClassificationMetrics.scala
+++ b/src/main/scala/io/citrine/lolo/stats/metrics/ClassificationMetrics.scala
@@ -14,7 +14,9 @@ object ClassificationMetrics {
     * @return the weighted average f1 score
     */
   def f1scores(predictedVsActual: Seq[(Vector[Any], Any, Any)]): Double = {
-    val labels = predictedVsActual.map(_._3).distinct
+    // the test set and training set might have different labels, so we might
+    // have predictions that are never actual and vice versa
+    val labels = (predictedVsActual.map(_._2) ++ predictedVsActual.map(_._3)).distinct
     val index = labels.zipWithIndex.toMap
     val numLabels = labels.size
     val confusionMatrix = DenseMatrix.zeros[Int](numLabels, numLabels)
diff --git a/src/main/scala/io/citrine/lolo/trees/FeaturePowerSetTerms.scala b/src/main/scala/io/citrine/lolo/trees/FeaturePowerSetTerms.scala
new file mode 100644
index 00000000..e93358e4
--- /dev/null
+++ b/src/main/scala/io/citrine/lolo/trees/FeaturePowerSetTerms.scala
@@ -0,0 +1,150 @@
+package io.citrine.lolo.trees
+
+/**
+  * Description of a feature's effect on the weight ascribed to the node
+  *
+  * Decision trees can be thought of as weighted sums over their leaves.  When there is knowledge of every feature,
+  * all of the weight goes into a single leaf.  However, when feature are excluded from the prediction, then both
+  * paths of decisions that depend on those excluded features are assigned non-zero weights that are proportional to
+  * the share of the training data that followed that path.
+  *
+  * This class captures those weights for a single feature and a single node. The `weightWhenIncluded` gives the factor
+  * of the weight of the node due to this feature when the feature is included (i.e. known).  It is always 0 or 1.
+  * The `weightWhenExcluded` gives the factor of the weight of the node due to this feature when the feature is
+  * excluded (i.e. unknown).  It is always strictly greater than 0 and strictly less than 1.
+  *
+  * @param weightWhenExcluded factor of the weight of the leaf when the feature is excluded
+  * @param weightWhenIncluded factor of the weight of the leaf when the feature is included
+  */
+case class FeatureWeightFactor(
+                                weightWhenExcluded: Double,
+                                weightWhenIncluded: Double
+                              ) {
+  require(weightWhenIncluded == 0.0 || weightWhenIncluded == 1.0,
+    s"Got weightWhenIncluded=$weightWhenIncluded, but should only ever be 0.0 or 1.0")
+  require(weightWhenExcluded > 0 && weightWhenExcluded < 1,
+    s"Got weightWhenExcluded=$weightWhenExcluded, but should be > 0 and < 1")
+}
+
+/**
+  * Weighted combinatorial factors corresponding to the power-set of features that are included
+  *
+  * The set is extended as new features are encountered.  Each element in `weightBySubsetSize` corresponds to all of the
+  * subsets of the features that are included vs excluded of a given size:
+  *  - the 0th element corresponds to the term where all of the encountered features are excluded,
+  *  - the 1st element sums the terms that have exactly one of the encountered feature included,
+  *  - the 2nd element sums the terms that have exactly two of the encountered features included,
+  * and so on and so forth.  The final element represents the single term where all of the features are "turned on".
+  *
+  * The calculation is performed efficiently using a dynamic programming technique that runs in quadratic time with
+  * the number of features.  It is described in the `extend` and `unwind` methods.
+  *
+  * @param maxFeatures number of features to allocate space for
+  */
+class FeaturePowerSetTerms(maxFeatures: Int) {
+  // pre-allocation of this whole array is an attempted performance optimization.
+  val weightBySubsetSize: Array[Double] = Array.fill[Double](maxFeatures + 1)(0.0)
+  weightBySubsetSize(0) = 1.0
+  var size: Int = 0
+
+  /**
+    * Extend the power set by adding a new feature (in-place)
+    *
+    * Each subset weight is updated with a contribution from the subset of the same size multiplied by the weight when
+    * the feature is excluded and the subset of size one smaller multiplied by the weight when the feature is included.
+    * The additional terms in the expression track the combinatorial factors |S|! and  (M - |S| - 1)! from eq. 2 of
+    * https://arxiv.org/pdf/1802.03888.pdf
+    *
+    * @param weightWhenExcluded weight factor for this node when the feature is excluded (unknown)
+    * @param weightWhenIncluded weight factor for this node when the feature is included (known)
+    * @return this (in-place)
+    */
+  def extend(
+              weightWhenExcluded: Double,
+              weightWhenIncluded: Double
+            ): FeaturePowerSetTerms = {
+    size += 1
+
+    // use a single scratch variable so we can do a forward in-place update
+    var register: Double = weightBySubsetSize(0)
+
+    // there is no subset with size -1, so special case i = 0
+    weightBySubsetSize(0) = weightWhenExcluded * weightBySubsetSize(0) * ((size - 0).toDouble / (size + 1))
+
+    (1 until size + 1).foreach { i =>
+      val contrib = weightWhenIncluded * register * (i.toDouble / (size + 1))
+      // stash the weight for use in the next iteration
+      register = weightBySubsetSize(i)
+      weightBySubsetSize(i) = weightWhenExcluded * weightBySubsetSize(i) * ((size - i).toDouble / (size + 1)) + contrib
+    }
+
+    this
+  }
+
+  /**
+    * Remove a feature from the power set, reversing the action of the `extend` method
+    *
+    * This method is probably better called "remove", but it is called unwind in the paper.  It's not pretty, but it
+    * simply inverts `extend`.  If you think about extend as multiplying by a matrix with weightWhenExcluded... on the
+    * diagonal and weightWhenIncluded... on the lower diagonal.  This method is performing a fast linear solve on that
+    * banded system, with the two logical branches corresponding to the banded vs purely diagonal case.
+    * The ... above denotes some combinatorial terms also present.
+    *
+    * @return feature power set with a feature removed (out-of-place)
+    */
+  def unwind(weightWhenExcluded: Double, weightWhenIncluded: Double): FeaturePowerSetTerms = {
+    // make a copy so this is out of place
+    val out = this.copy()
+
+    // reverse the procedure in extend
+    var n = out.weightBySubsetSize(size)
+    (size - 1 to 0 by -1).foreach { j =>
+      if (weightWhenIncluded != 0.0) {
+        val t = out.weightBySubsetSize(j)
+        out.weightBySubsetSize(j) = n * (size + 1) / ((j + 1) * weightWhenIncluded)
+        n = t - out.weightBySubsetSize(j) * weightWhenExcluded * ((size - j).toDouble / (size + 1))
+      } else {
+        out.weightBySubsetSize(j) = out.weightBySubsetSize(j) * (size + 1).toDouble / (weightWhenExcluded * (size - j))
+      }
+    }
+
+    // bookkeeping
+    out.size -= 1
+    out
+  }
+
+  /**
+    * Get the total weight of the feature power set, as a sum over the weights by size of the sets in the power set
+    */
+  def totalWeight: Double = weightBySubsetSize.take(size + 1).sum
+
+  /**
+    * This is an optimized version of unwind().totalWeight.  Because the unwound feature power set isn't used,
+    * it is not even allocated.  That avoided allocation is a significant performance improvement.
+    */
+  def unwoundTotalWeight(weightWhenExcluded: Double, weightWhenIncluded: Double): Double = {
+    var res = 0.0
+    // reverse the procedure in extend
+
+    if (weightWhenIncluded != 0.0) {
+      var n = weightBySubsetSize(size)
+      (size - 1 to 0 by -1).foreach { j =>
+        val x = n * (size + 1) / ((j + 1) * weightWhenIncluded)
+        res += x
+        n = weightBySubsetSize(j) - x * weightWhenExcluded * (size - j) / (size + 1)
+      }
+    } else {
+      (0 until size).foreach { j =>
+        res += weightBySubsetSize(j) * (size + 1) / (weightWhenExcluded * (size - j))
+      }
+    }
+    res
+  }
+
+  def copy(): FeaturePowerSetTerms = {
+    val newSet = new FeaturePowerSetTerms(this.maxFeatures)
+    newSet.size = this.size
+    this.weightBySubsetSize.zipWithIndex.foreach { case (x, i) => newSet.weightBySubsetSize(i) = x }
+    newSet
+  }
+}
diff --git a/src/main/scala/io/citrine/lolo/trees/Nodes.scala b/src/main/scala/io/citrine/lolo/trees/Nodes.scala
index 19653502..a3dd5582 100644
--- a/src/main/scala/io/citrine/lolo/trees/Nodes.scala
+++ b/src/main/scala/io/citrine/lolo/trees/Nodes.scala
@@ -1,5 +1,6 @@
 package io.citrine.lolo.trees
 
+import breeze.linalg.DenseMatrix
 import io.citrine.lolo.trees.splits.Split
 import io.citrine.lolo.{Learner, Model, PredictionResult}
 
@@ -35,20 +36,65 @@ abstract class TrainingNode[T <: AnyVal, S](
 
 trait ModelNode[T <: PredictionResult[Any]] extends Serializable {
   def transform(input: Vector[AnyVal]): (T, TreeMeta)
+
+  /**
+    * Compute Shapley feature attributions for a given input in this node's subtree
+    *
+    * @param input for which to compute feature attributions.
+    * @return array of vector-valued attributions for each feature
+    *         One Vector[Double] per feature, each of length equal to the output dimension.
+    */
+  def shapley(input: Vector[AnyVal], omitFeatures: Set[Int] = Set()): Option[DenseMatrix[Double]]
+
+  /**
+    * Get the Shapley feature attribution from a subtree
+    *
+    * This uses an algorithm based on TreeShap from https://arxiv.org/abs/1802.03888
+    *
+    * The algorithm recursively tracks the feature that is used for each split and the share of the training data
+    * that went to the left and the right of the split.  In the leaves, these features and their weights are used
+    * to construct the leaf's contribution to the SHAP values, which are then summed.
+    *
+    * @param input          for which to compute feature attributions.
+    * @param featureWeights Map from feature index to [[FeatureWeightFactor]], which stores the weight of the child
+    *                       of the split when the feature is known vs unknown
+    * @return matrix of attributions for each feature and output
+    *         One row per feature, each of length equal to the output dimension.
+    *         The output dimension is 1 for single-task regression, or equal to the number of classification categories.
+    */
+  private[lolo] def shapleyRecurse(
+                                    input: Vector[AnyVal],
+                                    omitFeatures: Set[Int] = Set(),
+                                    featureWeights: Map[Int, FeatureWeightFactor]
+                                  ): DenseMatrix[Double]
+
+  /**
+    * Weight of training data in subtree, specifically the number of data for unweighted training sets
+    *
+    * @return total weight of training weight in subtree
+    */
+  private[lolo] def getTrainingWeight(): Double
 }
 
 /**
   * Internal node in the decision tree
   *
-  * @param split to decide which branch to take
-  * @param left  branch node
-  * @param right branch node
+  * @param split           to decide which branch to take
+  * @param left            branch node
+  * @param right           branch node
+  * @param numFeatures     number of features, used for Shapley computation
+  * @param outputDimension dimension of model output, used for Shapley computation
+  *                        1 for single-task regression, or equal to the number of classification categories.
+  * @param trainingWeight  weight of training data in subtree (i.e. size of unweighted training set)
   * @tparam T type of the output
   */
 class InternalModelNode[T <: PredictionResult[Any]](
                                                      split: Split,
                                                      left: ModelNode[T],
-                                                     right: ModelNode[T]
+                                                     right: ModelNode[T],
+                                                     numFeatures: Int,
+                                                     outputDimension: Int,
+                                                     trainingWeight: Double
                                                    ) extends ModelNode[T] {
   /**
     * Just propagate the prediction call through the appropriate child
@@ -63,6 +109,76 @@ class InternalModelNode[T <: PredictionResult[Any]](
       right.transform(input)
     }
   }
+
+  override def shapley(input: Vector[AnyVal], omitFeatures: Set[Int] = Set()): Option[DenseMatrix[Double]] = {
+    // Kick off a recursive procedure, the base case of which is an empty map
+    Some(shapleyRecurse(input, omitFeatures, Map()))
+  }
+
+  /**
+    * On the way down: append this node to the set of features that have been encountered in this path
+    * through the decision tree unless this feature is in the omitted features list.
+    *
+    * On the way up: sum the contributions from the two children of this node.  If the feature that this node splits
+    * on is in the omitted features list, then multiply the contributions by the share of the training data that
+    * went to each child in lieu of including the feature on the way down.
+    */
+  def shapleyRecurse(
+                      input: Vector[AnyVal],
+                      omitFeatures: Set[Int],
+                      featureWeights: Map[Int, FeatureWeightFactor]
+                    ): DenseMatrix[Double] = {
+    val featureIndex = split.getIndex()
+
+    // The hot node is the one that is selected when the feature is present
+    val (hot, cold) = if (this.split.turnLeft(input)) {
+      (left, right)
+    } else {
+      (right, left)
+    }
+
+    // If the feature is omitted or unknown, then the weight assigned to each child is taken
+    // as the share of the training data that got split into each child
+    val hotPortion = hot.getTrainingWeight() / trainingWeight
+    val coldPortion = cold.getTrainingWeight() / trainingWeight
+
+    if (omitFeatures.contains(featureIndex)) {
+      // Don't add the feature to the featureWeights
+      // Multiply the child contributions by the hotPortion and coldPortion instead
+      val hotContrib = hotPortion * hot.shapleyRecurse(input, omitFeatures, featureWeights)
+      val coldContrib = coldPortion * cold.shapleyRecurse(input, omitFeatures, featureWeights)
+
+      // Finally sum
+      hotContrib + coldContrib
+    } else {
+      // Check if the feature in this node has been used already (on the way down)
+      val (hotFactor, coldFactor) = featureWeights.get(featureIndex) match {
+        case None =>
+          // If this is the first time the feature is being split on,
+          // then when the feature is excluded the weights are hotPortion and coldPortion
+          // and when the feature is included they are 1.0 and 0.0
+          (
+            FeatureWeightFactor(hotPortion, 1.0),
+            FeatureWeightFactor(coldPortion, 0.0)
+          )
+        case Some(previousWeights) =>
+          // If the feature has been split on before, then multiply the previous weights by the new ones
+          (
+            FeatureWeightFactor(hotPortion * previousWeights.weightWhenExcluded, previousWeights.weightWhenIncluded),
+            FeatureWeightFactor(coldPortion * previousWeights.weightWhenExcluded, 0.0)
+          )
+      }
+
+      // Add the features to the featureWeights map and recurse down
+      val hotContrib = hot.shapleyRecurse(input, omitFeatures, featureWeights.updated(featureIndex, hotFactor))
+      val coldContrib = cold.shapleyRecurse(input, omitFeatures, featureWeights.updated(featureIndex, coldFactor))
+
+      // Simple sum on the way up
+      coldContrib + hotContrib
+    }
+  }
+
+  override def getTrainingWeight(): Double = trainingWeight
 }
 
 /**
@@ -84,14 +200,53 @@ class TrainingLeaf[T](
     * @return lightweight prediction node
     */
   def getNode(): ModelNode[PredictionResult[T]] = {
-    new ModelLeaf(leafLearner.train(trainingData).getModel().asInstanceOf[Model[PredictionResult[T]]], depth)
+    new ModelLeaf(leafLearner.train(trainingData).getModel().asInstanceOf[Model[PredictionResult[T]]], depth, trainingData.head._1.length, trainingData.size.toDouble)
   }
 
   override def getFeatureImportance(): mutable.ArraySeq[Double] = mutable.ArraySeq.fill(trainingData.head._1.size)(0.0)
 }
 
-class ModelLeaf[T](model: Model[PredictionResult[T]], depth: Int) extends ModelNode[PredictionResult[T]] {
+class ModelLeaf[T](model: Model[PredictionResult[T]], depth: Int, numFeatures: Int, trainingWeight: Double) extends ModelNode[PredictionResult[T]] {
   override def transform(input: Vector[AnyVal]): (PredictionResult[T], TreeMeta) = {
     (model.transform(Seq(input)), TreeMeta(depth))
   }
+
+  /**
+    * Compute the contribution to SHAP in the leaf based on the features that were encountered between the root node
+    * and this leaf.  Note that the order of these features does *not* matter.
+    *
+    * The contributions are based on a tricky combinatorial factor that can be computed using dynamic programming.
+    * For details of this procedure, see [[FeaturePowerSetTerms]].
+    */
+  def shapleyRecurse(
+                      input: Vector[AnyVal],
+                      omitFeatures: Set[Int],
+                      featureWeights: Map[Int, FeatureWeightFactor]
+                    ): DenseMatrix[Double] = {
+    // Start with an empty matrix, into which we'll set the non-zero contributions
+    val shapValues = DenseMatrix.zeros[Double](1, input.length)
+
+    // Load all of the features into the set, which performs a dynamic programming calculation
+    val set = new FeaturePowerSetTerms(featureWeights.size)
+    featureWeights.values.foreach { case FeatureWeightFactor(exclude, include) => set.extend(exclude, include) }
+
+    // The contribution is proportional to the leaf's prediction, so grab that
+    this.model.transform(Seq(input)).getExpected().head match {
+      case v: Double =>
+        // For each feature, compute the contribution and store it in shapValues
+        featureWeights.foreach { case (featureIndex, node) =>
+          // Compute the weight of the contribution by removing this feature from the set and then computing the weights
+          // These two steps are fused to avoid an extra memory allocation
+          // This is equivalent to `set.unwind(node.weightWhenExcluded, node.weightWhenIncluded).totalWeight`
+          val w = set.unwoundTotalWeight(node.weightWhenExcluded, node.weightWhenIncluded)
+          shapValues(0, featureIndex) = w * (node.weightWhenIncluded - node.weightWhenExcluded) * v
+        }
+      case _ => throw new NotImplementedError()
+    }
+    shapValues
+  }
+
+  override def getTrainingWeight(): Double = trainingWeight
+
+  override def shapley(input: Vector[AnyVal], omitFeatures: Set[Int] = Set()): Option[DenseMatrix[Double]] = None
 }
diff --git a/src/main/scala/io/citrine/lolo/trees/classification/ClassificationTrainingNode.scala b/src/main/scala/io/citrine/lolo/trees/classification/ClassificationTrainingNode.scala
index a38f76f9..224c2010 100644
--- a/src/main/scala/io/citrine/lolo/trees/classification/ClassificationTrainingNode.scala
+++ b/src/main/scala/io/citrine/lolo/trees/classification/ClassificationTrainingNode.scala
@@ -1,11 +1,23 @@
 package io.citrine.lolo.trees.classification
 
+import breeze.linalg.DenseVector
 import io.citrine.lolo.trees.splits.{ClassificationSplitter, NoSplit, Split, Splitter}
 import io.citrine.lolo.trees.{InternalModelNode, ModelNode, TrainingLeaf, TrainingNode}
 import io.citrine.lolo.{Learner, PredictionResult}
 
 /**
   * Created by maxhutch on 1/12/17.
+  *
+  * @param trainingData     that this node sees
+  * @param leafLearner      to train on leaves
+  * @param split            used to split this node's direct descendants
+  * @param deltaImpurity    change in training set impurity between this node vs direct descendants
+  * @param numFeatures      number of features used in selecting splits
+  * @param minLeafInstances minimum training instances per node
+  * @param remainingDepth   to stop growing the node
+  * @param maxDepth         maximum depth of tree
+  * @param numClasses       total number of classes, used as an output dimensionality
+  * @param splitter         used to select splits
   */
 class ClassificationTrainingNode(
                                   trainingData: Seq[(Vector[AnyVal], Char, Double)],
@@ -16,6 +28,7 @@ class ClassificationTrainingNode(
                                   minLeafInstances: Int,
                                   remainingDepth: Int,
                                   maxDepth: Int,
+                                  numClasses: Int,
                                   splitter: Splitter[Char]
                                 ) extends TrainingNode(trainingData, remainingDepth) {
 
@@ -25,9 +38,9 @@ class ClassificationTrainingNode(
   lazy val (leftTrain, rightTrain) = trainingData.partition(r => split.turnLeft(r._1))
   assert(leftTrain.size > 0 && rightTrain.size > 0, s"Split ${split} resulted in zero size: ${trainingData.map(_._1(split.getIndex()))}")
 
-  lazy val leftChild = ClassificationTrainingNode.buildChild(leftTrain, leafLearner, minLeafInstances, remainingDepth, maxDepth, numFeatures, splitter)
+  lazy val leftChild = ClassificationTrainingNode.buildChild(leftTrain, leafLearner, minLeafInstances, remainingDepth, maxDepth, numFeatures, numClasses, splitter)
 
-  lazy val rightChild = ClassificationTrainingNode.buildChild(rightTrain, leafLearner, minLeafInstances, remainingDepth, maxDepth, numFeatures, splitter)
+  lazy val rightChild = ClassificationTrainingNode.buildChild(rightTrain, leafLearner, minLeafInstances, remainingDepth, maxDepth, numFeatures, numClasses, splitter)
 
 
   /**
@@ -36,7 +49,7 @@ class ClassificationTrainingNode(
     * @return lightweight prediction node
     */
   override def getNode(): ModelNode[PredictionResult[Char]] = new InternalModelNode(
-    split, leftChild.getNode(), rightChild.getNode()
+    split, leftChild.getNode(), rightChild.getNode(), numClasses, numFeatures, trainingData.size.toDouble
   )
 
   override def getFeatureImportance(): scala.collection.mutable.ArraySeq[Double] = {
@@ -60,6 +73,8 @@ object ClassificationTrainingNode {
     * @param remainingDepth   the number of splits left
     * @param maxDepth         to compute depth via remainingDepth
     * @param numFeatures      to consider in the split
+    * @param numClasses       total number of classes, used as an output dimensionality
+    * @param splitter         used to select splits
     * @return the child node, either a RegressionTrainingNode or TrainingLeaf
     */
   def buildChild(
@@ -69,12 +84,13 @@ object ClassificationTrainingNode {
                   remainingDepth: Int,
                   maxDepth: Int,
                   numFeatures: Int,
+                  numClasses: Int,
                   splitter: Splitter[Char]
                 ): TrainingNode[AnyVal, Char] = {
     if (trainingData.size >= 2 * minLeafInstances && remainingDepth > 0 && trainingData.exists(_._2 != trainingData.head._2)) {
       val (leftSplit, leftDelta) = splitter.getBestSplit(trainingData, numFeatures, minLeafInstances)
       if (!leftSplit.isInstanceOf[NoSplit]) {
-        new ClassificationTrainingNode(trainingData, leafLearner, leftSplit, leftDelta, numFeatures, minLeafInstances, remainingDepth - 1, maxDepth, splitter)
+        new ClassificationTrainingNode(trainingData, leafLearner, leftSplit, leftDelta, numFeatures, minLeafInstances, remainingDepth - 1, maxDepth, numClasses, splitter)
       } else {
         new TrainingLeaf(trainingData, leafLearner, maxDepth - remainingDepth)
       }
diff --git a/src/main/scala/io/citrine/lolo/trees/classification/ClassificationTree.scala b/src/main/scala/io/citrine/lolo/trees/classification/ClassificationTree.scala
index 138b441d..9f196ab1 100644
--- a/src/main/scala/io/citrine/lolo/trees/classification/ClassificationTree.scala
+++ b/src/main/scala/io/citrine/lolo/trees/classification/ClassificationTree.scala
@@ -6,21 +6,28 @@ import io.citrine.lolo.trees.splits.{ClassificationSplitter, NoSplit, Splitter}
 import io.citrine.lolo.trees.{ModelNode, TrainingLeaf, TrainingNode, TreeMeta}
 import io.citrine.lolo.{Learner, Model, PredictionResult, TrainingResult}
 
+import scala.util.Random
+
 
 /**
   * Created by maxhutch on 12/2/16.
   *
   * @param numFeatures subset of features to select splits from
+  * @param maxDepth maximum depth of tree
+  * @param minLeafInstances minimum training instances per node
+  * @param leafLearner to train on leaves
+  * @param splitter used to select splits
   */
 case class ClassificationTreeLearner(
                                       numFeatures: Int = -1,
                                       maxDepth: Int = 30,
                                       minLeafInstances: Int = 1,
                                       leafLearner: Option[Learner] = None,
-                                      splitter: Splitter[Char] = ClassificationSplitter()
+                                      splitter: Splitter[Char] = ClassificationSplitter(),
+                                      rng: Random = Random
                                     ) extends Learner {
 
-  @transient private lazy val myLeafLearner: Learner = leafLearner.getOrElse(new GuessTheMeanLearner)
+  @transient private lazy val myLeafLearner: Learner = leafLearner.getOrElse(GuessTheMeanLearner(rng = rng))
 
   /**
     * Train classification tree
@@ -75,6 +82,7 @@ case class ClassificationTreeLearner(
         remainingDepth = maxDepth - 1,
         maxDepth = maxDepth,
         minLeafInstances = minLeafInstances,
+        numClasses = trainingData.map{_._2}.distinct.length,
         splitter = splitter
       )
     }
@@ -154,7 +162,3 @@ class ClassificationResult(
     predictions.map(_._2.depth)
   }
 }
-
-
-
-
diff --git a/src/main/scala/io/citrine/lolo/trees/multitask/MultiTaskTrainingNode.scala b/src/main/scala/io/citrine/lolo/trees/multitask/MultiTaskTrainingNode.scala
index a404cec8..e080ea47 100644
--- a/src/main/scala/io/citrine/lolo/trees/multitask/MultiTaskTrainingNode.scala
+++ b/src/main/scala/io/citrine/lolo/trees/multitask/MultiTaskTrainingNode.scala
@@ -5,6 +5,8 @@ import io.citrine.lolo.linear.GuessTheMeanLearner
 import io.citrine.lolo.trees.splits.{MultiTaskSplitter, NoSplit}
 import io.citrine.lolo.trees.{InternalModelNode, ModelNode, TrainingLeaf}
 
+import scala.util.Random
+
 /**
   * Node in a multi-task training tree, which can produce nodes for its model trees
   *
@@ -12,18 +14,19 @@ import io.citrine.lolo.trees.{InternalModelNode, ModelNode, TrainingLeaf}
   */
 class MultiTaskTrainingNode(
                              inputs: Seq[(Vector[AnyVal], Array[AnyVal], Double)],
-                             randomizePivotLocation: Boolean = false
+                             randomizePivotLocation: Boolean = false,
+                             rng: Random = Random
                            ) {
 
   // Compute a split
-  val (split, _) = MultiTaskSplitter.getBestSplit(inputs, inputs.head._1.size, 1, randomizePivotLocation)
+  val (split, _) = MultiTaskSplitter(rng = rng).getBestSplit(inputs, inputs.head._1.size, 1, randomizePivotLocation)
 
   // Try to construct left and right children
   val (leftChild: Option[MultiTaskTrainingNode], rightChild: Option[MultiTaskTrainingNode]) = split match {
     case _: NoSplit => (None, None)
     case _: Any =>
       val (leftData, rightData) = inputs.partition(row => split.turnLeft(row._1))
-      (Some(new MultiTaskTrainingNode(leftData, randomizePivotLocation)), Some(new MultiTaskTrainingNode(rightData, randomizePivotLocation)))
+      (Some(new MultiTaskTrainingNode(leftData, randomizePivotLocation, rng = rng)), Some(new MultiTaskTrainingNode(rightData, randomizePivotLocation, rng = rng)))
   }
 
   // Construct the model node for the `index`th label
@@ -46,14 +49,20 @@ class MultiTaskTrainingNode(
         new InternalModelNode[PredictionResult[Double]](
           split,
           leftChild.get.getNode(index).asInstanceOf[ModelNode[PredictionResult[Double]]],
-          rightChild.get.getNode(index).asInstanceOf[ModelNode[PredictionResult[Double]]]
+          rightChild.get.getNode(index).asInstanceOf[ModelNode[PredictionResult[Double]]],
+          numFeatures=inputs.head._1.length,
+          outputDimension = 0,  // Don't support multitask SHAP at this time.
+          trainingWeight = reducedData.length.toDouble
         )
       } else {
         if (!label.isInstanceOf[Char]) throw new IllegalArgumentException("Training data wasn't double or char")
         new InternalModelNode[PredictionResult[Char]](
           split,
           leftChild.get.getNode(index).asInstanceOf[ModelNode[PredictionResult[Char]]],
-          rightChild.get.getNode(index).asInstanceOf[ModelNode[PredictionResult[Char]]]
+          rightChild.get.getNode(index).asInstanceOf[ModelNode[PredictionResult[Char]]],
+          numFeatures=inputs.head._1.length,
+          outputDimension = 0,  // Don't support multitask SHAP at this time.
+          trainingWeight = reducedData.length.toDouble
         )
       }
     } else if (leftChild.isDefined && left.nonEmpty) {
diff --git a/src/main/scala/io/citrine/lolo/trees/multitask/MultiTaskTree.scala b/src/main/scala/io/citrine/lolo/trees/multitask/MultiTaskTree.scala
index ba08b849..ef713f9f 100644
--- a/src/main/scala/io/citrine/lolo/trees/multitask/MultiTaskTree.scala
+++ b/src/main/scala/io/citrine/lolo/trees/multitask/MultiTaskTree.scala
@@ -6,12 +6,15 @@ import io.citrine.lolo.trees.classification.ClassificationTree
 import io.citrine.lolo.trees.regression.RegressionTree
 import io.citrine.lolo.{Model, MultiTaskLearner, PredictionResult, TrainingResult}
 
+import scala.util.Random
+
 /**
   * Multi-task tree learner, which produces multiple decision trees with the same split structure
   *
   */
 case class MultiTaskTreeLearner(
-                                 randomizePivotLocation: Boolean = false
+                                 randomizePivotLocation: Boolean = false,
+                                 rng: Random = Random
                                ) extends MultiTaskLearner {
 
   /**
diff --git a/src/main/scala/io/citrine/lolo/trees/regression/RegressionTrainingLeaf.scala b/src/main/scala/io/citrine/lolo/trees/regression/RegressionTrainingLeaf.scala
index 3e3f21ed..9e7c6d40 100644
--- a/src/main/scala/io/citrine/lolo/trees/regression/RegressionTrainingLeaf.scala
+++ b/src/main/scala/io/citrine/lolo/trees/regression/RegressionTrainingLeaf.scala
@@ -21,7 +21,7 @@ class RegressionTrainingLeaf(
     * @return lightweight prediction node
     */
   def getNode(): ModelNode[PredictionResult[Double]] = {
-    new ModelLeaf(model.asInstanceOf[Model[PredictionResult[Double]]], depth)
+    new ModelLeaf(model.asInstanceOf[Model[PredictionResult[Double]]], depth, 1, trainingData.size.toDouble)
   }
 
   /**
diff --git a/src/main/scala/io/citrine/lolo/trees/regression/RegressionTrainingNode.scala b/src/main/scala/io/citrine/lolo/trees/regression/RegressionTrainingNode.scala
index 8395fd27..e7c81dd8 100644
--- a/src/main/scala/io/citrine/lolo/trees/regression/RegressionTrainingNode.scala
+++ b/src/main/scala/io/citrine/lolo/trees/regression/RegressionTrainingNode.scala
@@ -1,8 +1,7 @@
 package io.citrine.lolo.trees.regression
 
-import io.citrine.lolo.trees.splits.{NoSplit, RegressionSplitter, Split, Splitter}
+import io.citrine.lolo.trees.splits.{NoSplit, Split, Splitter}
 import io.citrine.lolo.trees.{InternalModelNode, ModelNode, TrainingNode}
-import io.citrine.lolo.trees.{InternalModelNode, ModelNode, TrainingLeaf, TrainingNode}
 import io.citrine.lolo.{Learner, PredictionResult}
 
 /**
@@ -40,7 +39,7 @@ class RegressionTrainingNode(
     * @return lightweight prediction node
     */
   override def getNode(): ModelNode[PredictionResult[Double]] = {
-    new InternalModelNode[PredictionResult[Double]](split, leftChild.getNode(), rightChild.getNode())
+    new InternalModelNode[PredictionResult[Double]](split, leftChild.getNode(), rightChild.getNode(), numFeatures, 1, trainingData.size.toDouble)
   }
 
   /**
@@ -53,14 +52,14 @@ class RegressionTrainingNode(
     */
   override def getFeatureImportance(): scala.collection.mutable.ArraySeq[Double] = {
     val improvement = deltaImpurity
-    var ans = leftChild.getFeatureImportance().zip(rightChild.getFeatureImportance()).map(p => p._1 + p._2)
-    ans(split.getIndex) = ans(split.getIndex) + improvement
+    val ans = leftChild.getFeatureImportance().zip(rightChild.getFeatureImportance()).map(p => p._1 + p._2)
+    ans(split.getIndex()) = ans(split.getIndex()) + improvement
     ans
   }
 }
 
 /**
-  * Companion objet to hold helper functions
+  * Companion object to hold helper functions
   */
 object RegressionTrainingNode {
   /**
diff --git a/src/main/scala/io/citrine/lolo/trees/regression/RegressionTree.scala b/src/main/scala/io/citrine/lolo/trees/regression/RegressionTree.scala
index c23d3428..97b79be9 100644
--- a/src/main/scala/io/citrine/lolo/trees/regression/RegressionTree.scala
+++ b/src/main/scala/io/citrine/lolo/trees/regression/RegressionTree.scala
@@ -1,12 +1,14 @@
 package io.citrine.lolo.trees.regression
 
+import breeze.linalg.DenseMatrix
 import io.citrine.lolo.encoders.CategoricalEncoder
 import io.citrine.lolo.linear.GuessTheMeanLearner
 import io.citrine.lolo.trees.splits.{NoSplit, RegressionSplitter, Splitter}
-import io.citrine.lolo.trees.splits.{NoSplit, RegressionSplitter}
 import io.citrine.lolo.trees.{ModelNode, TrainingNode, TreeMeta}
 import io.citrine.lolo.{Learner, Model, PredictionResult, TrainingResult}
 
+import scala.util.Random
+
 /**
   * Learner for regression trees
   *
@@ -22,10 +24,11 @@ case class RegressionTreeLearner(
                                   maxDepth: Int = 30,
                                   minLeafInstances: Int = 1,
                                   leafLearner: Option[Learner] = None,
-                                  splitter: Splitter[Double] = RegressionSplitter()
+                                  splitter: Splitter[Double] = RegressionSplitter(),
+                                  rng: Random = Random
                                 ) extends Learner {
   /** Learner to use for training the leaves */
-  @transient private lazy val myLeafLearner = leafLearner.getOrElse(GuessTheMeanLearner())
+  @transient private lazy val myLeafLearner = leafLearner.getOrElse(GuessTheMeanLearner(rng = rng))
 
   /**
     * Train the tree by recursively partitioning (splitting) the training data on a single feature
@@ -139,6 +142,20 @@ class RegressionTree(
       inputs.map(inp => root.transform(CategoricalEncoder.encodeInput(inp, encoders)))
     )
   }
+
+  /**
+    * Compute Shapley feature attributions for a given input
+    *
+    * @param input for which to compute feature attributions.
+    * @param omitFeatures feature indices to omit in computing Shapley values
+    * @return array of Shapley feature attributions, one per input feature, each a vector of
+    *         One Vector[Double] per feature, each of length equal to the output dimension.
+    *         The output dimension is 1 for single-task regression, or equal to the number of classification categories.
+    */
+  override def shapley(input: Vector[Any], omitFeatures: Set[Int] = Set()): Option[DenseMatrix[Double]] = {
+    root.shapley(CategoricalEncoder.encodeInput(input, encoders), omitFeatures)
+  }
+
 }
 
 /**
diff --git a/src/main/scala/io/citrine/lolo/trees/splits/BoltzmannSplitter.scala b/src/main/scala/io/citrine/lolo/trees/splits/BoltzmannSplitter.scala
index af7921b5..08866c57 100644
--- a/src/main/scala/io/citrine/lolo/trees/splits/BoltzmannSplitter.scala
+++ b/src/main/scala/io/citrine/lolo/trees/splits/BoltzmannSplitter.scala
@@ -30,7 +30,7 @@ import scala.util.Random
   * @param temperature used to control how sensitive the probability of a split is to its change in variance.
   *                    The temperature can be thought of as a hyperparameter.
   */
-case class BoltzmannSplitter(temperature: Double) extends Splitter[Double] {
+case class BoltzmannSplitter(temperature: Double, rng: Random = Random) extends Splitter[Double] {
   require(temperature >= Float.MinPositiveValue, s"Temperature must be >= ${Float.MinPositiveValue} to avoid numerical underflows")
 
   /**
@@ -58,11 +58,11 @@ case class BoltzmannSplitter(temperature: Double) extends Splitter[Double] {
     /* Try every feature index */
     val featureIndices: Seq[Int] = rep._1.indices
 
-    val possibleSplits: Seq[SplitterResult] = Random.shuffle(featureIndices).take(numFeatures).flatMap { index =>
+    val possibleSplits: Seq[SplitterResult] = rng.shuffle(featureIndices).take(numFeatures).flatMap { index =>
       /* Use different spliters for each type */
       rep._1(index) match {
-        case _: Double => BoltzmannSplitter.getBestRealSplit(data, calculator, index, minInstances, beta)
-        case _: Char => BoltzmannSplitter.getBestCategoricalSplit(data, calculator, index, minInstances, beta)
+        case _: Double => BoltzmannSplitter.getBestRealSplit(data, calculator, index, minInstances, beta, rng)
+        case _: Char => BoltzmannSplitter.getBestCategoricalSplit(data, calculator, index, minInstances, beta, rng)
         case _: Any => throw new IllegalArgumentException("Trying to split unknown feature type")
       }
     }
@@ -79,7 +79,7 @@ case class BoltzmannSplitter(temperature: Double) extends Splitter[Double] {
 
     // select from a discrete probability distribution by drawing a random number and then computing the CDF
     // where the "draw" is the bin for which the CDF crosses the drawn number
-    val draw = Random.nextDouble() * totalProbability
+    val draw = rng.nextDouble() * totalProbability
     // could be a scanLeft + find, but this is more readable
     var cumSum: Double = 0.0
     possibleSplits.foreach { case SplitterResult(split, variance, score, base) =>
@@ -125,7 +125,8 @@ object BoltzmannSplitter {
                         calculator: VarianceCalculator,
                         index: Int,
                         minCount: Int,
-                        beta: Double
+                        beta: Double,
+                        rng: Random
                       ): Option[SplitterResult] = {
     /* Pull out the feature that's considered here and sort by it */
     val thinData = data.map(dat => (dat._1(index).asInstanceOf[Double], dat._2, dat._3)).sortBy(_._1)
@@ -144,7 +145,7 @@ object BoltzmannSplitter {
       val right = thinData(j)._1
       if (j + 1 >= minCount && Splitter.isDifferent(left, right)) {
         val score = -totalVariance * beta
-        val pivot = (left - right) * Random.nextDouble() + right
+        val pivot = (left - right) * rng.nextDouble() + right
         Some(score, pivot, totalVariance)
       } else {
         None
@@ -157,7 +158,7 @@ object BoltzmannSplitter {
 
     val base: Double = possibleSplits.map(_._1).max
     val totalScore = possibleSplits.map { case (s, _, _) => Math.exp(s - base) }.sum
-    val draw = Random.nextDouble() * totalScore
+    val draw = rng.nextDouble() * totalScore
     var cumSum: Double = 0.0
     possibleSplits.foreach { case (score, pivot, variance) =>
       cumSum = cumSum + Math.exp(score - base)
@@ -184,14 +185,15 @@ object BoltzmannSplitter {
                                calculator: VarianceCalculator,
                                index: Int,
                                minCount: Int,
-                               beta: Double
+                               beta: Double,
+                               rng: Random
                              ): Option[SplitterResult] = {
     /* Extract the features at the index */
     val thinData = data.map(dat => (dat._1(index).asInstanceOf[Char], dat._2, dat._3))
     val totalWeight = thinData.map(_._3).sum
 
     /* Group the data by categorical feature and compute the weighted sum and sum of the weights for each */
-    val groupedData = thinData.groupBy(_._1).mapValues(g => (g.map(v => v._2 * v._3).sum, g.map(_._3).sum, g.size))
+    val groupedData = thinData.groupBy(_._1).mapValues(g => (g.map(v => v._2 * v._3).sum, g.map(_._3).sum, g.size)).toMap
 
     /* Make sure there is more than one member for most of the classes */
     val nonTrivial: Double = groupedData.filter(_._2._3 > 1).map(_._2._2).sum
@@ -200,7 +202,7 @@ object BoltzmannSplitter {
     }
 
     /* Compute the average label for each categorical value */
-    val categoryAverages: Map[Char, Double] = groupedData.mapValues(p => p._1 / p._2)
+    val categoryAverages: Map[Char, Double] = groupedData.mapValues(p => p._1 / p._2).toMap
 
     /* Create an orderd list of the categories by average label */
     val orderedNames = categoryAverages.toSeq.sortBy(_._2).map(_._1)
@@ -230,7 +232,7 @@ object BoltzmannSplitter {
 
     val base: Double = possibleSplits.map(_._1).max
     val totalScore = possibleSplits.map { case (s, _, _) => Math.exp(s - base) }.sum
-    val draw = Random.nextDouble() * totalScore
+    val draw = rng.nextDouble() * totalScore
     var cumSum: Double = 0.0
     possibleSplits.foreach { case (score, includeSet, variance) =>
       cumSum = cumSum + Math.exp(score - base)
diff --git a/src/main/scala/io/citrine/lolo/trees/splits/ClassificationSplitter.scala b/src/main/scala/io/citrine/lolo/trees/splits/ClassificationSplitter.scala
index 9a24be8b..095384a9 100644
--- a/src/main/scala/io/citrine/lolo/trees/splits/ClassificationSplitter.scala
+++ b/src/main/scala/io/citrine/lolo/trees/splits/ClassificationSplitter.scala
@@ -9,7 +9,7 @@ import scala.util.Random
   *
   * Created by maxhutch on 12/2/16.
   */
-case class ClassificationSplitter(randomizedPivotLocation: Boolean = false) extends Splitter[Char]{
+case class ClassificationSplitter(randomizedPivotLocation: Boolean = false, rng: Random = Random) extends Splitter[Char]{
 
   /**
     * Get the best split, considering numFeature random features (w/o replacement)
@@ -33,7 +33,7 @@ case class ClassificationSplitter(randomizedPivotLocation: Boolean = false) exte
     val rep = data.head
     /* Try every feature index */
     val featureIndices: Seq[Int] = rep._1.indices
-    Random.shuffle(featureIndices).take(numFeatures).foreach { index =>
+    rng.shuffle(featureIndices).take(numFeatures).foreach { index =>
 
       /* Use different spliters for each type */
       val (possibleSplit, possibleImpurity) = rep._1(index) match {
@@ -98,7 +98,7 @@ case class ClassificationSplitter(randomizedPivotLocation: Boolean = false) exte
         val left = thinData(j + 1)._1
         val right = thinData(j)._1
         bestPivot = if (randomizePivotLocation) {
-          (left - right) * Random.nextDouble() + right
+          (left - right) * rng.nextDouble() + right
         } else {
           (left + right) / 2.0
         }
diff --git a/src/main/scala/io/citrine/lolo/trees/splits/ExtraRandomClassificationSplitter.scala b/src/main/scala/io/citrine/lolo/trees/splits/ExtraRandomClassificationSplitter.scala
new file mode 100644
index 00000000..cf131f86
--- /dev/null
+++ b/src/main/scala/io/citrine/lolo/trees/splits/ExtraRandomClassificationSplitter.scala
@@ -0,0 +1,140 @@
+package io.citrine.lolo.trees.splits
+
+import io.citrine.lolo.trees.impurity.{GiniCalculator, ImpurityCalculator}
+
+import scala.util.Random
+
+/**
+  * A splitter that defines Extremely Randomized Trees
+  *
+  * This is based on Geurts, P., Ernst, D. & Wehenkel, L. Extremely randomized trees. Mach Learn 63, 3–42 (2006).
+  * https://doi.org/10.1007/s10994-006-6226-1.
+  *
+  * @param rng random number generator
+  */
+case class ExtraRandomClassificationSplitter(
+                                              rng: Random = Random
+                                            ) extends Splitter[Char] {
+
+  /**
+    * Get the best split, considering numFeature random features (w/o replacement)
+    *
+    * @param data        to split
+    * @param numFeatures to consider, randomly
+    * @return a split object that optimally divides data
+    */
+  def getBestSplit(
+                    data: Seq[(Vector[AnyVal], Char, Double)],
+                    numFeatures: Int,
+                    minInstances: Int
+                  ): (Split, Double) = {
+
+    val calculator = GiniCalculator.build(data.map{ p => (p._2, p._3) })
+    val initialVariance = calculator.getImpurity
+    var bestSplit: Split = new NoSplit()
+    var bestVariance = Double.MaxValue
+
+    val rep = data.head
+
+    /* Try every feature index */
+    val featureIndices: Seq[Int] = rep._1.indices
+    rng.shuffle(featureIndices.toVector).take(numFeatures).foreach { index =>
+      /* Use different spliters for each type */
+      val (possibleSplit, possibleVariance) = rep._1(index) match {
+        case _: Double => getBestRealSplit(data, calculator.copy(), index, minInstances)
+        case _: Char => getBestCategoricalSplit(data, calculator.copy(), index, minInstances)
+        case _: Any => throw new IllegalArgumentException("Trying to split unknown feature type")
+      }
+
+      /* Keep track of the best split */
+      if (possibleVariance < bestVariance) {
+        bestVariance = possibleVariance
+        bestSplit = possibleSplit
+      }
+    }
+    if (bestVariance >=  initialVariance) {
+      (new NoSplit(), 0.0)
+    } else {
+      val deltaImpurity = initialVariance - bestVariance
+      (bestSplit, deltaImpurity)
+    }
+  }
+
+  /**
+    * Find the best split on a continuous variable.
+    *
+    * @param data  to split
+    * @param index of the feature to split on
+    * @param minCount minimum number of data points to allow in each of the resulting splits
+    * @return the best split of this feature
+    */
+  def getBestRealSplit(
+                        data: Seq[(Vector[AnyVal], Char, Double)],
+                        calculator: ImpurityCalculator[Char],
+                        index: Int,
+                        minCount: Int
+                      ): (RealSplit, Double) = {
+    /* Pull out the feature that's considered here and sort by it */
+    val axis: Seq[Double] = data.map(_._1(index).asInstanceOf[Double])
+    val lowerBound = axis.min
+    val upperBound = axis.max
+    val pivot = lowerBound + (upperBound - lowerBound) * rng.nextDouble()
+
+    val split = new RealSplit(index, pivot)
+
+    /* Move the data from the right to the left partition one value at a time */
+    calculator.reset()
+    data.foreach{ dat =>
+      if (split.turnLeft(dat._1)) {
+        calculator.add(dat._2, dat._3)
+      }
+    }
+
+    (split, calculator.getImpurity)
+  }
+
+  /**
+    * Get find the best categorical splitter.
+    *
+    * @param data  to split
+    * @param index of the feature to split on
+    * @return the best split of this feature
+    */
+  def getBestCategoricalSplit(
+                               data: Seq[(Vector[AnyVal], Char, Double)],
+                               calculator: ImpurityCalculator[Char],
+                               index: Int,
+                               minCount: Int
+                             ): (Split, Double) = {
+    /* Extract the features at the index */
+    val thinData = data.map(dat => (dat._1(index).asInstanceOf[Char], dat._2, dat._3))
+
+    /* Group the data by categorical feature and compute the weighted sum and sum of the weights for each */
+    val groupedData = thinData.groupBy(_._1).mapValues(g => (g.map(v => v._2 * v._3).sum, g.map(_._3).sum, g.size))
+
+    /* Compute the average label for each categorical value */
+    val categoryAverages: Map[Char, Double] = groupedData.mapValues(p => p._1 / p._2).toMap
+
+    /* Create an ordered list of the categories by average label */
+    val orderedNames: Seq[Char] = categoryAverages.toSeq.sortBy(_._2).map(_._1)
+    if (orderedNames.size == 1) {
+      calculator.reset()
+      return (new NoSplit(), calculator.getImpurity)
+    }
+
+    val pivot = rng.nextInt(orderedNames.size - 1)
+    val bestSet = orderedNames.slice(0, pivot + 1).toSet
+    val split = new CategoricalSplit(index, new scala.collection.mutable.BitSet() ++ bestSet.map(_.toInt))
+
+    /* Move the data from the right to the left partition one value at a time */
+    calculator.reset()
+    data.foreach{ dat =>
+      if (split.turnLeft(dat._1)) {
+        calculator.add(dat._2, dat._3)
+      }
+    }
+
+    (split, calculator.getImpurity)
+  }
+
+}
diff --git a/src/main/scala/io/citrine/lolo/trees/splits/ExtraRandomRegressionSplitter.scala b/src/main/scala/io/citrine/lolo/trees/splits/ExtraRandomRegressionSplitter.scala
new file mode 100644
index 00000000..9c2a35d8
--- /dev/null
+++ b/src/main/scala/io/citrine/lolo/trees/splits/ExtraRandomRegressionSplitter.scala
@@ -0,0 +1,140 @@
+package io.citrine.lolo.trees.splits
+
+import io.citrine.lolo.trees.impurity.{ImpurityCalculator, VarianceCalculator}
+
+import scala.util.Random
+
+/**
+  * A splitter that defines Extremely Randomized Trees
+  *
+  * This is based on Geurts, P., Ernst, D. & Wehenkel, L. Extremely randomized trees. Mach Learn 63, 3–42 (2006).
+  * https://doi.org/10.1007/s10994-006-6226-1.
+  *
+  * @param rng random number generator
+  */
+case class ExtraRandomRegressionSplitter(
+                                          rng: Random = Random
+                                        ) extends Splitter[Double] {
+
+  /**
+    * Get the best split, considering numFeature random features (w/o replacement)
+    *
+    * @param data        to split
+    * @param numFeatures to consider, randomly
+    * @return a split object that optimally divides data
+    */
+  def getBestSplit(
+                    data: Seq[(Vector[AnyVal], Double, Double)],
+                    numFeatures: Int,
+                    minInstances: Int
+                  ): (Split, Double) = {
+
+    val calculator = VarianceCalculator.build(data.map(_._2), data.map(_._3))
+    val initialVariance = calculator.getImpurity
+    var bestSplit: Split = new NoSplit()
+    var bestVariance = Double.MaxValue
+
+    val rep = data.head
+
+    /* Try every feature index */
+    val featureIndices: Seq[Int] = rep._1.indices
+    rng.shuffle(featureIndices.toVector).take(numFeatures).foreach { index =>
+      /* Use different spliters for each type */
+      val (possibleSplit, possibleVariance) = rep._1(index) match {
+        case _: Double => getBestRealSplit(data, calculator.copy(), index, minInstances)
+        case _: Char => getBestCategoricalSplit(data, calculator.copy(), index, minInstances)
+        case _: Any => throw new IllegalArgumentException("Trying to split unknown feature type")
+      }
+
+      /* Keep track of the best split */
+      if (possibleVariance < bestVariance) {
+        bestVariance = possibleVariance
+        bestSplit = possibleSplit
+      }
+    }
+    if (bestVariance >=  initialVariance) {
+      (new NoSplit(), 0.0)
+    } else {
+      val deltaImpurity = initialVariance - bestVariance
+      (bestSplit, deltaImpurity)
+    }
+  }
+
+  /**
+    * Find the best split on a continuous variable.
+    *
+    * @param data  to split
+    * @param index of the feature to split on
+    * @param minCount minimum number of data points to allow in each of the resulting splits
+    * @return the best split of this feature
+    */
+  def getBestRealSplit(
+                        data: Seq[(Vector[AnyVal], Double, Double)],
+                        calculator: ImpurityCalculator[Double],
+                        index: Int,
+                        minCount: Int
+                      ): (RealSplit, Double) = {
+    /* Pull out the feature that's considered here and sort by it */
+    val axis: Seq[Double] = data.map(_._1(index).asInstanceOf[Double])
+    val lowerBound = axis.min
+    val upperBound = axis.max
+    val pivot = lowerBound + (upperBound - lowerBound) * rng.nextDouble()
+
+    val split = new RealSplit(index, pivot)
+
+    /* Move the data from the right to the left partition one value at a time */
+    calculator.reset()
+    data.foreach{ dat =>
+      if (split.turnLeft(dat._1)) {
+        calculator.add(dat._2, dat._3)
+      }
+    }
+
+    (split, calculator.getImpurity)
+  }
+
+  /**
+    * Get find the best categorical splitter.
+    *
+    * @param data  to split
+    * @param index of the feature to split on
+    * @return the best split of this feature
+    */
+  def getBestCategoricalSplit(
+                               data: Seq[(Vector[AnyVal], Double, Double)],
+                               calculator: ImpurityCalculator[Double],
+                               index: Int,
+                               minCount: Int
+                             ): (Split, Double) = {
+    /* Extract the features at the index */
+    val thinData = data.map(dat => (dat._1(index).asInstanceOf[Char], dat._2, dat._3))
+
+    /* Group the data by categorical feature and compute the weighted sum and sum of the weights for each */
+    val groupedData = thinData.groupBy(_._1).mapValues(g => (g.map(v => v._2 * v._3).sum, g.map(_._3).sum, g.size))
+
+    /* Compute the average label for each categorical value */
+    val categoryAverages: Map[Char, Double] = groupedData.mapValues(p => p._1 / p._2).toMap
+
+    /* Create an ordered list of the categories by average label */
+    val orderedNames: Seq[Char] = categoryAverages.toSeq.sortBy(_._2).map(_._1)
+    if (orderedNames.size == 1) {
+      calculator.reset()
+      return (new NoSplit(), calculator.getImpurity)
+    }
+
+    val pivot = rng.nextInt(orderedNames.size - 1)
+    val bestSet = orderedNames.slice(0, pivot + 1).toSet
+    val split = new CategoricalSplit(index, new scala.collection.mutable.BitSet() ++ bestSet.map(_.toInt))
+
+    /* Move the data from the right to the left partition one value at a time */
+    calculator.reset()
+    data.foreach{ dat =>
+      if (split.turnLeft(dat._1)) {
+        calculator.add(dat._2, dat._3)
+      }
+    }
+
+    (split, calculator.getImpurity)
+  }
+
+}
diff --git a/src/main/scala/io/citrine/lolo/trees/splits/MultiTaskSplitter.scala b/src/main/scala/io/citrine/lolo/trees/splits/MultiTaskSplitter.scala
index 006d333d..11608d13 100644
--- a/src/main/scala/io/citrine/lolo/trees/splits/MultiTaskSplitter.scala
+++ b/src/main/scala/io/citrine/lolo/trees/splits/MultiTaskSplitter.scala
@@ -8,7 +8,7 @@ import scala.util.Random
   *
   * Created by maxhutch on 11/29/16.
   */
-object MultiTaskSplitter {
+case class MultiTaskSplitter(rng: Random = Random) {
 
   /**
     * Get the best split, considering numFeature random features (w/o replacement)
@@ -33,7 +33,7 @@ object MultiTaskSplitter {
 
     /* Try every feature index */
     val featureIndices: Seq[Int] = rep._1.indices
-    Random.shuffle(featureIndices).take(numFeatures).foreach { index =>
+    rng.shuffle(featureIndices).take(numFeatures).foreach { index =>
 
       /* Use different spliters for each type */
       val (possibleSplit, possibleImpurity) = rep._1(index) match {
@@ -82,7 +82,7 @@ object MultiTaskSplitter {
         val left = features(j + 1)
         val right = features(j)
         val pivot = if (randomizePivotLocation) {
-          (left - right) * Random.nextDouble() + right
+          (left - right) * rng.nextDouble() + right
         } else {
           (left + right) / 2.0
         }
@@ -119,8 +119,8 @@ object MultiTaskSplitter {
 
     /* Group the data by categorical feature and compute the weighted sum and sum of the weights for each */
     val groupedData: Map[Char, (Double, Double, Double)] = thinData.groupBy(_._1).mapValues(g =>
-      (computeImpurity(g.map(x => (x._2, x._3))), g.map(_._3).sum, g.size)
-    )
+      (computeImpurity(g.map(x => (x._2, x._3))), g.map(_._3).sum, g.size.toDouble)
+    ).toMap
 
     /* Make sure there is more than one member for most of the classes */
     val nonTrivial: Double = groupedData.filter(_._2._3 > 1).map(_._2._2).sum
diff --git a/src/main/scala/io/citrine/lolo/trees/splits/RegressionSplitter.scala b/src/main/scala/io/citrine/lolo/trees/splits/RegressionSplitter.scala
index 987aba6f..5a6b7881 100644
--- a/src/main/scala/io/citrine/lolo/trees/splits/RegressionSplitter.scala
+++ b/src/main/scala/io/citrine/lolo/trees/splits/RegressionSplitter.scala
@@ -20,7 +20,7 @@ import scala.util.Random
   *
   * Created by maxhutch on 11/29/16.
   */
-case class RegressionSplitter(randomizePivotLocation: Boolean = false) extends Splitter[Double] {
+case class RegressionSplitter(randomizePivotLocation: Boolean = false, rng: Random = Random) extends Splitter[Double] {
 
   /**
     * Get the best split, considering numFeature random features (w/o replacement)
@@ -44,7 +44,7 @@ case class RegressionSplitter(randomizePivotLocation: Boolean = false) extends S
 
     /* Try every feature index */
     val featureIndices: Seq[Int] = rep._1.indices
-    Random.shuffle(featureIndices).take(numFeatures).foreach { index =>
+    rng.shuffle(featureIndices).take(numFeatures).foreach { index =>
 
       /* Use different spliters for each type */
       val (possibleSplit, possibleVariance) = rep._1(index) match {
@@ -111,7 +111,7 @@ case class RegressionSplitter(randomizePivotLocation: Boolean = false) extends S
         bestVariance = totalVariance
         /* Try pivots at the midpoints between consecutive member values */
         bestPivot = if (randomizePivotLocation) {
-          (left - right) * Random.nextDouble() + right
+          (left - right) * rng.nextDouble() + right
         } else {
           (left + right) / 2.0
         }
@@ -147,7 +147,7 @@ case class RegressionSplitter(randomizePivotLocation: Boolean = false) extends S
     }
 
     /* Compute the average label for each categorical value */
-    val categoryAverages: Map[Char, Double] = groupedData.mapValues(p => p._1 / p._2)
+    val categoryAverages: Map[Char, Double] = groupedData.mapValues(p => p._1 / p._2).toMap
 
     /* Create an orderd list of the categories by average label */
     val orderedNames = categoryAverages.toSeq.sortBy(_._2).map(_._1)
diff --git a/src/main/scala/io/citrine/lolo/validation/CrossValidation.scala b/src/main/scala/io/citrine/lolo/validation/CrossValidation.scala
index 6c1163bf..6ad81b40 100644
--- a/src/main/scala/io/citrine/lolo/validation/CrossValidation.scala
+++ b/src/main/scala/io/citrine/lolo/validation/CrossValidation.scala
@@ -7,7 +7,7 @@ import scala.util.Random
 /**
   * Methods tha use cross-validation to calculate predicted-vs-actual data and metric estimates
   */
-object CrossValidation {
+case object CrossValidation {
 
   /**
     * Driver to apply named metrics to k-fold cross-validated predicted-vs-actual
@@ -17,6 +17,7 @@ object CrossValidation {
     * @param metrics      apply to the predicted-vs-actual data
     * @param k            number of folds
     * @param nTrial       number of times to refold the data to improve statistics
+    * @param rng          random number generator to use in choosing folds
     * @tparam T type of the response, e.g. Double for scalar regression
     * @return a Map from the name of the metric to its mean value and the error in that mean
     */
@@ -25,11 +26,13 @@ object CrossValidation {
                                learner: Learner,
                                metrics: Map[String, Merit[T]],
                                k: Int = 8,
-                               nTrial: Int = 1
+                               nTrial: Int = 1,
+                               rng: Random = Random
                              ): Map[String, (Double, Double)] = {
     Merit.estimateMerits(
-      kFoldPvA(trainingData, learner, k, nTrial).iterator,
-      metrics
+      kFoldPvA(trainingData, learner, k, nTrial, rng).iterator,
+      metrics,
+      rng
     )
   }
 
@@ -47,11 +50,12 @@ object CrossValidation {
                    trainingData: Seq[(Vector[Any], T)],
                    learner: Learner,
                    k: Int = 8,
-                   nTrial: Int = 1
+                   nTrial: Int = 1,
+                   rng: Random = Random
                  ): Iterable[(PredictionResult[T], Seq[T])] = {
     val nTest: Int = Math.ceil(trainingData.size.toDouble / k).toInt
     (0 until nTrial).flatMap { _ =>
-      val folds: Seq[Seq[(Vector[Any], T)]] = Random.shuffle(trainingData).grouped(nTest).toSeq
+      val folds: Seq[Seq[(Vector[Any], T)]] = rng.shuffle(trainingData).grouped(nTest).toSeq
 
       folds.indices.map { idx =>
         val (testFolds, trainFolds) = folds.zipWithIndex.partition(_._2 == idx)
diff --git a/src/main/scala/io/citrine/lolo/validation/Merit.scala b/src/main/scala/io/citrine/lolo/validation/Merit.scala
index 2399290e..a28e9f68 100644
--- a/src/main/scala/io/citrine/lolo/validation/Merit.scala
+++ b/src/main/scala/io/citrine/lolo/validation/Merit.scala
@@ -18,7 +18,7 @@ trait Merit[T] {
     *
     * @return the value of the figure of merit
     */
-  def evaluate(predictionResult: PredictionResult[T], actual: Seq[T]): Double
+  def evaluate(predictionResult: PredictionResult[T], actual: Seq[T], rng: Random = Random): Double
 
   /**
     * Estimate the merit and the uncertainty in the merit over batches of predicted and ground-truth values
@@ -26,8 +26,8 @@ trait Merit[T] {
     * @param pva predicted-vs-actual data as an iterable over [[PredictionResult]] and ground-truth tuples
     * @return the estimate of the merit value and the uncertainty in that estimate
     */
-  def estimate(pva: Iterable[(PredictionResult[T], Seq[T])]): (Double, Double) = {
-    val samples = pva.map { case (prediction, actual) => evaluate(prediction, actual) }
+  def estimate(pva: Iterable[(PredictionResult[T], Seq[T])], rng: Random = Random): (Double, Double) = {
+    val samples = pva.map { case (prediction, actual) => evaluate(prediction, actual, rng) }
     val mean: Double = samples.sum / samples.size
     val variance: Double = (samples.size / (samples.size - 1)) * samples.map(x => Math.pow(x - mean, 2)).sum / samples.size
     (mean, Math.sqrt(variance / samples.size))
@@ -38,7 +38,7 @@ trait Merit[T] {
   * Square root of the mean square error. For an unbiased estimator, this is equal to the standard deviation of the difference between predicted and actual values.
   */
 case object RootMeanSquareError extends Merit[Double] {
-  override def evaluate(predictionResult: PredictionResult[Double], actual: Seq[Double]): Double = {
+  override def evaluate(predictionResult: PredictionResult[Double], actual: Seq[Double], rng: Random = Random): Double = {
     Math.sqrt(
       (predictionResult.getExpected(), actual).zipped.map {
         case (x, y) => Math.pow(x - y, 2)
@@ -51,7 +51,7 @@ case object RootMeanSquareError extends Merit[Double] {
   * R2 = 1 - MSE(y) / Var(y), where y is the predicted variable
   */
 case object CoefficientOfDetermination extends Merit[Double] {
-  override def evaluate(predictionResult: PredictionResult[Double], actual: Seq[Double]): Double = {
+  override def evaluate(predictionResult: PredictionResult[Double], actual: Seq[Double], rng: Random = Random): Double = {
     val averageActual = actual.sum / actual.size
     val sumOfSquares = actual.map(x => Math.pow(x - averageActual, 2)).sum
     val sumOfResiduals = predictionResult.getExpected().zip(actual).map { case (x, y) => Math.pow(x - y, 2.0) }.sum
@@ -63,7 +63,7 @@ case object CoefficientOfDetermination extends Merit[Double] {
   * The fraction of predictions that fall within the predicted uncertainty
   */
 case object StandardConfidence extends Merit[Double] {
-  override def evaluate(predictionResult: PredictionResult[Double], actual: Seq[Double]): Double = {
+  override def evaluate(predictionResult: PredictionResult[Double], actual: Seq[Double], rng: Random = Random): Double = {
     if (predictionResult.getUncertainty().isEmpty) return 0.0
 
     (predictionResult.getExpected(), predictionResult.getUncertainty().get, actual).zipped.count {
@@ -76,7 +76,7 @@ case object StandardConfidence extends Merit[Double] {
   * Root mean square of (the error divided by the predicted uncertainty)
   */
 case class StandardError(rescale: Double = 1.0) extends Merit[Double] {
-  override def evaluate(predictionResult: PredictionResult[Double], actual: Seq[Double]): Double = {
+  override def evaluate(predictionResult: PredictionResult[Double], actual: Seq[Double], rng: Random = Random): Double = {
     if (predictionResult.getUncertainty().isEmpty) return Double.PositiveInfinity
     val standardized = (predictionResult.getExpected(), predictionResult.getUncertainty().get, actual).zipped.map {
       case (x, sigma: Double, y) => (x - y) / sigma
@@ -97,7 +97,7 @@ case class StandardError(rescale: Double = 1.0) extends Merit[Double] {
   * an "ideal" error series from which the correlation coefficient can be estimated.
   */
 case object UncertaintyCorrelation extends Merit[Double] {
-  override def evaluate(predictionResult: PredictionResult[Double], actual: Seq[Double]): Double = {
+  override def evaluate(predictionResult: PredictionResult[Double], actual: Seq[Double], rng: Random = Random): Double = {
     val predictedUncertaintyActual: Seq[(Double, Double, Double)] = (
       predictionResult.getExpected(),
       predictionResult.getUncertainty().get.asInstanceOf[Seq[Double]],
@@ -105,7 +105,7 @@ case object UncertaintyCorrelation extends Merit[Double] {
     ).zipped.toSeq
 
     val ideal = predictedUncertaintyActual.map { case (_, uncertainty, actual) =>
-      val error = Random.nextGaussian() * uncertainty
+      val error = rng.nextGaussian() * uncertainty
       (actual + error, uncertainty, actual)
     }
 
@@ -145,18 +145,19 @@ object Merit {
     */
   def estimateMerits[T](
                           pva: Iterator[(PredictionResult[T], Seq[T])],
-                          merits: Map[String, Merit[T]]
+                          merits: Map[String, Merit[T]],
+                          rng: Random = Random
                         ): Map[String, (Double, Double)] = {
 
     pva.flatMap { case (predictions, actual) =>
       // apply all the merits to the batch at the same time so the batch can fall out of memory
-      merits.mapValues(f => f.evaluate(predictions, actual)).toSeq
+      merits.mapValues(f => f.evaluate(predictions, actual, rng)).toSeq
     }.toIterable.groupBy(_._1).mapValues { x =>
       val meritResults = x.map(_._2)
       val mean = meritResults.sum / meritResults.size
       val variance = meritResults.map(y => Math.pow(y - mean, 2)).sum / meritResults.size
       (mean, Math.sqrt(variance / meritResults.size))
-    }
+    }.toMap
   }
 
   /**
@@ -167,6 +168,7 @@ object Merit {
     * @param merits          to apply at each parameter value
     * @param logScale        whether the parameters should be plotted on a log scale
     * @param pvaBuilder      function that takes the parameter to predicted-vs-actual data
+    * @param rng             random number generator to use
     * @return an [[XYChart]] that plots the merits vs the parameter value
     */
   def plotMeritScan[T](
@@ -175,7 +177,8 @@ object Merit {
                          merits: Map[String, Merit[T]],
                          logScale: Boolean = false,
                          yMin: Option[Double] = None,
-                         yMax: Option[Double] = None
+                         yMax: Option[Double] = None,
+                         rng: Random = Random
                        )(
                          pvaBuilder: Double => Iterator[(PredictionResult[T], Seq[T])]
                        ): XYChart = {
@@ -189,7 +192,7 @@ object Merit {
 
     parameterValues.foreach { param =>
       val pva = pvaBuilder(param)
-      val meritResults = Merit.estimateMerits(pva, merits)
+      val meritResults = Merit.estimateMerits(pva, merits, rng)
       meritResults.foreach { case (name, (mean, err)) =>
         seriesData(name).add(mean)
         seriesData(s"${name}_err").add(err)
diff --git a/src/main/scala/io/citrine/lolo/validation/StatisticalValidation.scala b/src/main/scala/io/citrine/lolo/validation/StatisticalValidation.scala
index c6d83f8d..69f51f7f 100644
--- a/src/main/scala/io/citrine/lolo/validation/StatisticalValidation.scala
+++ b/src/main/scala/io/citrine/lolo/validation/StatisticalValidation.scala
@@ -7,7 +7,7 @@ import scala.util.Random
 /**
   * Methods that draw data from a distribution and compute predicted-vs-actual data
   */
-object StatisticalValidation {
+case class StatisticalValidation(rng: Random = Random) {
 
   /**
     * Generate predicted-vs-actual data given a source of ground truth data and a learner
@@ -69,7 +69,7 @@ object StatisticalValidation {
                                nRound: Int
                              ): Iterator[(PredictionResult[T], Seq[T])] = {
     Iterator.tabulate(nRound) { _ =>
-      val subset = Random.shuffle(source).take(nTrain + nTest)
+      val subset = rng.shuffle(source).take(nTrain + nTest)
       val (trainingData: Seq[(Vector[Any], T)], testData: Seq[(Vector[Any], T)]) = subset.splitAt(nTrain)
       val model = learner.train(trainingData).getModel()
       val predictions: PredictionResult[T] = model.transform(testData.map(_._1)).asInstanceOf[PredictionResult[T]]
diff --git a/src/site/index.html b/src/site/index.html
new file mode 100644
index 00000000..80b876b8
--- /dev/null
+++ b/src/site/index.html
@@ -0,0 +1,21 @@
+<!DOCTYPE html>
+<html lang="en">
+<head>
+    <meta charset="UTF-8">
+    <title>Project Documentation</title>
+    <script language="JavaScript">
+    <!--
+    function doRedirect()
+    {
+        window.location.replace("latest/api");
+    }
+
+    doRedirect();
+    //-->
+    </script>
+</head>
+<body>
+<a href="latest/api">Go to the project documentation
+</a>
+</body>
+</html>


Test patch:
diff --git a/python/lolopy/tests/test_learners.py b/python/lolopy/tests/test_learners.py
index 37db6940..c82c9aae 100644
--- a/python/lolopy/tests/test_learners.py
+++ b/python/lolopy/tests/test_learners.py
@@ -1,4 +1,11 @@
-from lolopy.learners import RandomForestRegressor, RandomForestClassifier, RegressionTreeLearner, LinearRegression
+from lolopy.learners import (
+    RandomForestRegressor,
+    RandomForestClassifier,
+    RegressionTreeLearner,
+    LinearRegression,
+    ExtraRandomTreesRegressor,
+    ExtraRandomTreesClassifier
+)
 from sklearn.exceptions import NotFittedError
 from sklearn.metrics import r2_score, accuracy_score, log_loss
 from sklearn.datasets import load_boston, load_iris
@@ -25,7 +32,7 @@ def _make_linear_data():
 class TestRF(TestCase):
 
     def test_rf_regressor(self):
-        rf = RandomForestRegressor()
+        rf = RandomForestRegressor(random_seed = 31247895)
 
         # Train the model
         X, y = load_boston(True)
@@ -51,7 +58,7 @@ def test_rf_regressor(self):
         self.assertGreater(score, 0.98)
 
         # Test with weights (make sure it doesn't crash)
-        rf.fit(X, y, [1.0]*len(y))
+        rf.fit(X, y, [2.0]*len(y))
 
         # Make sure feature importances are stored
         self.assertEqual(np.shape(rf.feature_importances_), (X.shape[1],))
@@ -67,9 +74,8 @@ def test_rf_regressor(self):
         rf.clear_model()
         self.assertIsNone(rf.model_)
 
-
     def test_classifier(self):
-        rf = RandomForestClassifier()
+        rf = RandomForestClassifier(random_seed = 34789)
 
         # Load in the iris dataset
         X, y = load_iris(True)
@@ -89,10 +95,10 @@ def test_classifier(self):
         self.assertEqual(len(X), len(y_pred))
         acc = accuracy_score(y, y_pred)
         print('Accuracy:', acc)
-        self.assertAlmostEqual(acc, 1)  # Given default settings, we should get perfect fittness to training data
+        self.assertAlmostEqual(acc, 1)  # Given default settings, we should get perfect fitness to training data
 
     def test_serialization(self):
-        rf = RandomForestClassifier()
+        rf = RandomForestClassifier(random_seed = 234785)
 
         # Make sure this doesn't error without a model training
         data = pkl.dumps(rf)
@@ -173,14 +179,103 @@ def test_adjust_rtree_learners(self):
         self.assertAlmostEqual(1.0, r2_score(y, tree.predict(X)))  # Linear leaves means perfect fit
 
         # Test whether changing leaf learner does something
-        rf = RandomForestRegressor(leaf_learner=LinearRegression(), min_leaf_instances=16)
+        rf = RandomForestRegressor(leaf_learner=LinearRegression(), min_leaf_instances=16, random_seed = 23478)
         rf.fit(X[:16, :], y[:16])  # Train only on a subset
         self.assertAlmostEqual(1.0, r2_score(y, rf.predict(X)))  # Should fit perfectly on whole dataset
 
-        rf = RandomForestRegressor()
+        rf = RandomForestRegressor(random_seed = 7834)
         rf.fit(X[:16, :], y[:16])
         self.assertLess(r2_score(y, rf.predict(X)), 1.0)  # Should not fit the whole dataset perfectly
 
 
+class TestExtraRandomTrees(TestCase):
+
+    def test_extra_random_trees_regressor(self):
+        rf = ExtraRandomTreesRegressor(random_seed = 378456)
+
+        # Train the model
+        X, y = load_boston(True)
+
+        # Make sure we get a NotFittedError
+        with self.assertRaises(NotFittedError):
+            rf.predict(X)
+
+        # Fit the model
+        rf.fit(X, y)
+
+        # Run some predictions
+        y_pred = rf.predict(X)
+        self.assertEqual(len(y_pred), len(y))
+
+        # Test the ability to get importance scores
+        y_import = rf.get_importance_scores(X[:100, :])
+        self.assertEqual((100, len(X)), y_import.shape)
+
+        # Basic test for functionality. R^2 above 0.98 was measured on 2020-04-06
+        score = r2_score(y_pred, y)
+        print("R2: ", score)
+        self.assertGreater(score, 0.98)
+
+        # Test with weights (make sure it doesn't crash)
+        rf.fit(X, y, [2.0]*len(y))
+
+        # Make sure feature importances are stored
+        self.assertEqual(np.shape(rf.feature_importances_), (X.shape[1],))
+        self.assertAlmostEqual(1.0, np.sum(rf.feature_importances_))
+
+        # Run predictions with std dev
+        y_pred, y_std = rf.predict(X, return_std=True)
+        self.assertEqual(len(y_pred), len(y_std))
+        self.assertTrue((y_std >= 0).all())  # They must be positive
+        self.assertGreater(np.std(y_std), 0)  # Must have a variety of values
+
+        # Make sure the detach operation functions
+        rf.clear_model()
+        self.assertIsNone(rf.model_)
+
+    def test_extra_random_trees_classifier(self):
+        rf = ExtraRandomTreesClassifier(random_seed = 378456)
+
+        # Load in the iris dataset
+        X, y = load_iris(True)
+        rf.fit(X, y)
+
+        # Predict the probability of membership in each class
+        y_prob = rf.predict_proba(X)
+        self.assertEqual((len(X), 3), np.shape(y_prob))
+        self.assertAlmostEqual(len(X), np.sum(y_prob))
+        ll_score = log_loss(y, y_prob)
+        print('Log loss:', ll_score)
+        self.assertLess(ll_score, 0.03)  # Measured at 0.026 on 2020-04-06
+
+        # Test just getting the predicted class
+        y_pred = rf.predict(X)
+        self.assertTrue(np.isclose(np.argmax(y_prob, 1), y_pred).all())
+        self.assertEqual(len(X), len(y_pred))
+        acc = accuracy_score(y, y_pred)
+        print('Accuracy:', acc)
+        self.assertAlmostEqual(acc, 1)  # Given default settings, we should get perfect fitness to training data
+
+    def test_serialization(self):
+        rf = ExtraRandomTreesClassifier(random_seed = 378945)
+
+        # Make sure this doesn't error without a model training
+        data = pkl.dumps(rf)
+        rf2 = pkl.loads(data)
+
+        # Load in the iris dataset and train model
+        X, y = load_iris(True)
+        rf.fit(X, y)
+
+        # Try saving and loading the model
+        data = pkl.dumps(rf)
+        rf2 = pkl.loads(data)
+
+        # Make sure it yields the same predictions as the first model
+        probs1 = rf.predict_proba(X)
+        probs2 = rf2.predict_proba(X)
+        self.assertTrue(np.isclose(probs1, probs2).all())
+
+
 if __name__ == "__main__":
     main()
diff --git a/python/lolopy/tests/test_metrics.py b/python/lolopy/tests/test_metrics.py
index 08b49402..b0a669d7 100644
--- a/python/lolopy/tests/test_metrics.py
+++ b/python/lolopy/tests/test_metrics.py
@@ -1,3 +1,4 @@
+from lolopy.loloserver import get_java_gateway
 from lolopy.metrics import (root_mean_squared_error, standard_confidence, standard_error, uncertainty_correlation)
 from numpy.random import multivariate_normal, uniform, normal, seed
 from unittest import TestCase
@@ -10,16 +11,20 @@ def test_rmse(self):
         self.assertAlmostEqual(root_mean_squared_error([4, 5], [1, 2]), 3)
 
     def test_standard_confidene(self):
-        self.assertAlmostEqual(standard_confidence([1, 2], [2, 3], [1.5, 0.9]), 0.5)
-        self.assertAlmostEqual(standard_confidence([1, 2], [2, 3], [1.5, 1.1]), 1)
+        gateway = get_java_gateway()
+        rng = gateway.jvm.scala.util.Random(367894)
+        self.assertAlmostEqual(standard_confidence([1, 2], [2, 3], [1.5, 0.9]), 0.5, rng)
+        self.assertAlmostEqual(standard_confidence([1, 2], [2, 3], [1.5, 1.1]), 1, rng)
 
     def test_standard_error(self):
         self.assertAlmostEqual(standard_error([1, 2], [1, 2], [1, 1]), 0)
         self.assertAlmostEqual(standard_error([4, 5], [1, 2], [3, 3]), 1)
 
     def test_uncertainty_correlation(self):
-        seed(1)
+        seed(3893789455)
         sample_size = 2 ** 15
+        gateway = get_java_gateway()
+        rng = gateway.jvm.scala.util.Random(783245)
         for expected in [0, 0.75]:
             # Make the error distribution
             y_true = uniform(0, 1, sample_size)
@@ -32,6 +37,6 @@ def test_uncertainty_correlation(self):
             y_std = [abs(d[1]) for d in draw]
 
             # Test with a very large tolerance for now
-            measured_corr = uncertainty_correlation(y_true, y_pred, y_std)
+            measured_corr = uncertainty_correlation(y_true, y_pred, y_std, rng = rng)
             corr_error = abs(measured_corr - expected)
             self.assertLess(corr_error, 0.25, 'Error for {:.2f}: {:.2f}'.format(expected, corr_error))
diff --git a/python/test_requirements.txt b/python/test_requirements.txt
new file mode 100644
index 00000000..f3c7e8e6
--- /dev/null
+++ b/python/test_requirements.txt
@@ -0,0 +1,1 @@
+nose
diff --git a/src/test/scala/io/citrine/lolo/AccuracyTest.scala b/src/test/scala/io/citrine/lolo/AccuracyTest.scala
index 71529510..01cba66e 100644
--- a/src/test/scala/io/citrine/lolo/AccuracyTest.scala
+++ b/src/test/scala/io/citrine/lolo/AccuracyTest.scala
@@ -6,10 +6,13 @@ import io.citrine.lolo.trees.splits.{BoltzmannSplitter, RegressionSplitter}
 import io.citrine.theta.Stopwatch
 import org.junit.Test
 
+import scala.util.Random
+
 /**
   * Created by maxhutch on 7/10/17.
   */
 class AccuracyTest {
+  val rng = new Random(234785L)
 
   val noiseLevel: Double = 0.000
   val nFeat: Int = 10
@@ -27,12 +30,62 @@ class AccuracyTest {
     learner.train(trainingData).getLoss().get
   }
 
+  /**
+    * Test that setting rng yields repeatable results.
+    */
+  /**
+    * DISABLED until implementing perfect reproducibility.
+  @Test
+  def testRepeatable(): Unit = {
+    val (errorsStandardTree, errorsAnnealingTree, errorsUnrandomizedTree) = (1 to 4).map { _ =>
+      val seed = 354128L
+      val rng = new Random(seed)
+      val errorStandardTree = {
+        val baseLearner = RegressionTreeLearner(
+          numFeatures = nFeat,
+          splitter = RegressionSplitter(randomizePivotLocation = true, rng = rng),
+          rng = rng
+        )
+        val learner = new Bagger(baseLearner, numBags = nRow * 8, randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
+        computeMetrics(learner)
+      }
+      rng.setSeed(seed)
+      val errorAnnealingTree = {
+        val baseLearner = RegressionTreeLearner(
+          numFeatures = nFeat,
+          splitter = BoltzmannSplitter(temperature = Float.MinPositiveValue, rng = rng),
+          rng = rng
+        )
+        val learner = new Bagger(baseLearner, numBags = nRow * 8, randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
+        computeMetrics(learner)
+      }
+      rng.setSeed(seed)
+      val randBasis = TestUtils.getBreezeRandBasis(rng.nextLong())
+      val errorUnrandomizedTree = {
+        val baseLearner = RegressionTreeLearner(
+          numFeatures = nFeat,
+          splitter = RegressionSplitter(randomizePivotLocation = false, rng = rng),
+          rng = rng
+        )
+        val learner = new Bagger(baseLearner, numBags = nRow * 8, randBasis = randBasis)
+        computeMetrics(learner)
+      }
+      (errorStandardTree, errorAnnealingTree, errorUnrandomizedTree)
+    }.unzip3
+
+    val atol = 1e-12
+    assert(errorsStandardTree.forall{ e => Math.abs(e - errorsStandardTree.head) < atol })
+    assert(errorsAnnealingTree.forall{ e => Math.abs(e - errorsStandardTree.head) < atol })
+    assert(errorsUnrandomizedTree.forall{ e => Math.abs(e - errorsUnrandomizedTree.head) < atol })
+  }
+  */
+
   /**
     * Quick sanity check of the test setup
     */
   @Test
   def testRandomForest(): Unit = {
-    val baseLearner = RegressionTreeLearner(numFeatures = nFeat / 3, minLeafInstances = minInstances)
+    val baseLearner = RegressionTreeLearner(numFeatures = nFeat / 3, minLeafInstances = minInstances, rng = rng)
     val learner = new Bagger(baseLearner, numBags = nRow * nScal)
     val error = computeMetrics(learner)
     assert(error > noiseLevel, s"Can't do better than noise")
@@ -47,18 +100,20 @@ class AccuracyTest {
     val errorStandardTree = {
       val baseLearner = RegressionTreeLearner(
         numFeatures = nFeat,
-        splitter = RegressionSplitter(randomizePivotLocation = true)
+        splitter = RegressionSplitter(randomizePivotLocation = true, rng = rng),
+        rng = rng
       )
-      val learner = new Bagger(baseLearner, numBags = nRow * 16)
+      val learner = new Bagger(baseLearner, numBags = nRow * 16, randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
       // println(s"Normal train time: ${Stopwatch.time(computeMetrics(learner))}")
       computeMetrics(learner)
     }
     val errorAnnealingTree = {
       val baseLearner = RegressionTreeLearner(
         numFeatures = nFeat,
-        splitter = BoltzmannSplitter(temperature = Float.MinPositiveValue)
+        splitter = BoltzmannSplitter(temperature = Float.MinPositiveValue, rng = rng),
+        rng = rng
       )
-      val learner = new Bagger(baseLearner, numBags = nRow * 16)
+      val learner = new Bagger(baseLearner, numBags = nRow * 16, randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
       // println(s"Annealing train time: ${Stopwatch.time(computeMetrics(learner))}")
       computeMetrics(learner)
     }
@@ -78,7 +133,7 @@ class AccuracyTest {
 object AccuracyTest {
 
   val trainingDataFull: Seq[(Vector[Any], Double)] = TestUtils.binTrainingData(
-    TestUtils.generateTrainingData(2048, 48),
+    TestUtils.generateTrainingData(2048, 48, seed = 283467L),
     inputBins = Seq((2, 32)) // bin the 3rd feature into a categorical
   ).asInstanceOf[Seq[(Vector[Any], Double)]]
 
@@ -102,7 +157,8 @@ object AccuracyTest {
     val baseLearner = RegressionTreeLearner(
       numFeatures = nFeatSub,
       splitter = splitter,
-      minLeafInstances = minInstances
+      minLeafInstances = minInstances,
+      rng = new Random(247895L)
     )
     val learner = new Bagger(baseLearner, numBags = nRow * nScal, biasLearner = None)
     val model = learner.train(trainingData).getModel()
diff --git a/src/test/scala/io/citrine/lolo/TestUtils.scala b/src/test/scala/io/citrine/lolo/TestUtils.scala
index 00de7efc..d6e3b6a6 100644
--- a/src/test/scala/io/citrine/lolo/TestUtils.scala
+++ b/src/test/scala/io/citrine/lolo/TestUtils.scala
@@ -1,7 +1,11 @@
 package io.citrine.lolo
 
+import breeze.stats.distributions.{RandBasis, ThreadLocalRandomGenerator}
 import io.citrine.lolo.stats.functions.Friedman
+import org.apache.commons.math3.random.MersenneTwister
 
+import scala.collection.generic.CanBuildFrom
+import scala.collection.mutable
 import scala.util.{Random, Try}
 
 /**
@@ -77,4 +81,24 @@ object TestUtils {
     outputData
   }
 
+  /**
+    * Enumerate the cartesian product of items in baseGrids.
+    *
+    * @param baseGrids a sequence of 1-d mesh specifications, one for each dimension of the output vectors
+    * @return a sequence of vectors enumerating the cartesian product of items in baseGrids
+    */
+  def enumerateGrid(baseGrids: Seq[Seq[Double]]): Seq[Vector[Double]] = {
+    if (baseGrids.length == 1) {
+      baseGrids.head.map { x => Vector(x) }
+    } else {
+      baseGrids.head.flatMap { x =>
+        enumerateGrid(baseGrids.takeRight(baseGrids.length - 1)).map { n =>
+          x +: n
+        }
+      }
+    }
+  }
+
+  def getBreezeRandBasis(seed: Long): RandBasis = new RandBasis(new ThreadLocalRandomGenerator(new MersenneTwister(seed)))
 }
+
diff --git a/src/test/scala/io/citrine/lolo/bags/BaggedResultTest.scala b/src/test/scala/io/citrine/lolo/bags/BaggedResultTest.scala
index f031375c..66e01d9d 100644
--- a/src/test/scala/io/citrine/lolo/bags/BaggedResultTest.scala
+++ b/src/test/scala/io/citrine/lolo/bags/BaggedResultTest.scala
@@ -11,6 +11,7 @@ import org.scalatest.Assertions._
 import scala.util.Random
 
 class BaggedResultTest {
+  val rng = new Random(278945L)
 
   @Test
   def testSingleMultiConsistency(): Unit = {
@@ -19,8 +20,8 @@ class BaggedResultTest {
       inputBins = Seq((0, 8))
     )
 
-    val DTLearner = RegressionTreeLearner(numFeatures = 12)
-    val biasLearner = RegressionTreeLearner(maxDepth = 5, leafLearner = Some(GuessTheMeanLearner()))
+    val DTLearner = RegressionTreeLearner(numFeatures = 12, rng = rng)
+    val biasLearner = RegressionTreeLearner(maxDepth = 5, leafLearner = Some(GuessTheMeanLearner(rng = rng)))
 
     Array(
       Bagger(DTLearner, numBags = 64, biasLearner = None, uncertaintyCalibration = false, useJackknife = true),
@@ -38,9 +39,8 @@ class BaggedResultTest {
   @Test
   def testBaggedSingleResultGetUncertainty(): Unit = {
     val noiseLevel = 100.0
-    val rng = new Random(237485L)
 
-    Seq(RegressionTreeLearner(), GuessTheMeanLearner()).foreach { baseLearner =>
+    Seq(RegressionTreeLearner(), GuessTheMeanLearner(rng = rng)).foreach { baseLearner =>
       // These are in Seqs as a convenience for repurposing this test as a diagnostic tool.
       Seq(128).foreach { nRows =>
         Seq(16).foreach { nCols =>
@@ -126,7 +126,7 @@ class BaggedResultTest {
     * @param model        The trained model
     */
   private def testConsistency(trainingData: Seq[(Vector[Any], Any)], model: BaggedModel[Any]): Unit = {
-    val testSubset = Random.shuffle(trainingData).take(16)
+    val testSubset = rng.shuffle(trainingData).take(16)
     val (singleValues, singleUncertainties) = testSubset.map { case (x, _) =>
       val res = model.transform(Seq(x))
       (res.getExpected().head.asInstanceOf[Double], res.getUncertainty().get.head.asInstanceOf[Double])
diff --git a/src/test/scala/io/citrine/lolo/bags/BaggerTest.scala b/src/test/scala/io/citrine/lolo/bags/BaggerTest.scala
index a13f1e27..8b9764f3 100644
--- a/src/test/scala/io/citrine/lolo/bags/BaggerTest.scala
+++ b/src/test/scala/io/citrine/lolo/bags/BaggerTest.scala
@@ -1,7 +1,9 @@
 package io.citrine.lolo.bags
 
+import java.io.{File, PrintWriter}
 import java.util.concurrent.{Callable, CancellationException, Executors, Future, TimeUnit}
 
+import breeze.linalg.DenseMatrix
 import io.citrine.lolo.TestUtils
 import io.citrine.lolo.linear.GuessTheMeanLearner
 import io.citrine.lolo.stats.functions.Friedman
@@ -19,18 +21,21 @@ import scala.util.Random
   */
 @Test
 class BaggerTest {
+  val rng = new Random(97793500L)
 
   /**
     * Test the fit performance of the regression bagger
     */
   @Test
   def testRegressionBagger(): Unit = {
+    rng.setSeed(24795L)
+
     val trainingData = TestUtils.binTrainingData(
       TestUtils.generateTrainingData(1024, 12, noise = 0.1, function = Friedman.friedmanSilverman),
       inputBins = Seq((0, 8))
     )
-    val DTLearner = RegressionTreeLearner(numFeatures = 3)
-    val baggedLearner = Bagger(DTLearner, numBags = trainingData.size)
+    val DTLearner = RegressionTreeLearner(numFeatures = 3, rng = rng)
+    val baggedLearner = Bagger(DTLearner, numBags = trainingData.size, randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
     val RFMeta = baggedLearner.train(trainingData)
     val RF = RFMeta.getModel()
 
@@ -53,12 +58,14 @@ class BaggerTest {
     */
   @Test
   def testClassificationBagger(): Unit = {
+    rng.setSeed(24795L)
+
     val trainingData = TestUtils.binTrainingData(
-      TestUtils.generateTrainingData(1024, 12, noise = 0.1, function = Friedman.friedmanSilverman),
+      TestUtils.generateTrainingData(1024, 12, noise = 0.1, function = Friedman.friedmanSilverman, seed = rng.nextLong()),
       inputBins = Seq((0, 8)), responseBins = Some(8)
     )
     val DTLearner = ClassificationTreeLearner()
-    val baggedLearner = Bagger(DTLearner, numBags = trainingData.size / 2)
+    val baggedLearner = Bagger(DTLearner, numBags = trainingData.size / 2, randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
     val RFMeta = baggedLearner.train(trainingData)
     val RF = RFMeta.getModel()
 
@@ -90,18 +97,20 @@ class BaggerTest {
     * set independently
     */
   def testUncertaintyCalibration(): Unit = {
+    rng.setSeed(24795L)
+
     val width = 0.10 // make the function more linear
     val nFeatures = 5
     val bagsPerRow = 4 // picked to be large enough that bias correction is small but model isn't too expensive
-    val trainingData = TestUtils.generateTrainingData(128, nFeatures, xscale = width, seed = Random.nextLong())
-    val DTLearner = RegressionTreeLearner(numFeatures = nFeatures)
-    val bias = RegressionTreeLearner(maxDepth = 4)
-    val baggedLearner = Bagger(DTLearner, numBags = bagsPerRow * trainingData.size, biasLearner = Some(bias))
+    val trainingData = TestUtils.generateTrainingData(128, nFeatures, xscale = width, seed = rng.nextLong())
+    val DTLearner = RegressionTreeLearner(numFeatures = nFeatures, rng = rng)
+    val bias = RegressionTreeLearner(maxDepth = 4, rng = rng)
+    val baggedLearner = Bagger(DTLearner, numBags = bagsPerRow * trainingData.size, biasLearner = Some(bias), randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
     val RFMeta = baggedLearner.train(trainingData)
     val RF = RFMeta.getModel()
 
-    val interiorTestSet = TestUtils.generateTrainingData(128, nFeatures, xscale = width / 2.0, xoff = width / 4.0, seed = Random.nextLong())
-    val fullTestSet = TestUtils.generateTrainingData(128, nFeatures, xscale = width, seed = Random.nextLong())
+    val interiorTestSet = TestUtils.generateTrainingData(128, nFeatures, xscale = width / 2.0, xoff = width / 4.0, seed = rng.nextLong())
+    val fullTestSet = TestUtils.generateTrainingData(128, nFeatures, xscale = width, seed = rng.nextLong())
 
     val interiorStandardRMSE = BaggerTest.getStandardRMSE(interiorTestSet, RF)
     val fullStandardRMSE = BaggerTest.getStandardRMSE(fullTestSet, RF)
@@ -117,17 +126,20 @@ class BaggerTest {
     */
   @Test
   def testUncertaintyCalibrationWithConstantResponse(): Unit = {
+    rng.setSeed(24795L)
+
     // setup some training data with constant labels
     val nFeatures = 5
-    val X: Vector[Vector[Any]] = TestUtils.generateTrainingData(128, nFeatures, xscale = 0.5, seed = Random.nextLong()).map(_._1)
+    val X: Vector[Vector[Any]] = TestUtils.generateTrainingData(128, nFeatures, xscale = 0.5, seed = rng.nextLong()).map(_._1)
     val y: Vector[Any] = X.map(_ => 0.0)
 
     // setup a relatively complicated random forest (turn a bunch of stuff on)
     val DTLearner = RegressionTreeLearner(
       numFeatures = nFeatures,
-      leafLearner = Some(GuessTheMeanLearner()),
+      leafLearner = Some(GuessTheMeanLearner(rng = rng)),
       maxDepth = 30,
-      splitter = RegressionSplitter(randomizePivotLocation = true)
+      splitter = RegressionSplitter(randomizePivotLocation = true, rng = rng),
+      rng = rng
     )
 
     val bagger = new Bagger(
@@ -136,17 +148,19 @@ class BaggerTest {
       useJackknife = true,
       biasLearner = Some(RegressionTreeLearner(
         maxDepth = 3,
-        leafLearner = Some(GuessTheMeanLearner()),
-        splitter = RegressionSplitter(randomizePivotLocation = true))
+        leafLearner = Some(GuessTheMeanLearner(rng = rng)),
+        splitter = RegressionSplitter(randomizePivotLocation = true),
+        rng = rng)
       ),
-      uncertaintyCalibration = true
+      uncertaintyCalibration = true,
+      randBasis = TestUtils.getBreezeRandBasis(rng.nextLong())
     )
 
     // Make sure the model trains
     val model = bagger.train(X.zip(y)).getModel()
 
     // Generate a new test set and make sure the predictions are 0 +/- 0
-    val testX: Vector[Vector[Any]] = TestUtils.generateTrainingData(128, nFeatures, xscale = 0.5, seed = Random.nextLong()).map(_._1)
+    val testX: Vector[Vector[Any]] = TestUtils.generateTrainingData(128, nFeatures, xscale = 0.5, seed = rng.nextLong()).map(_._1)
     val predictions = model.transform(testX)
     assert(predictions.getExpected().forall(_ == 0.0))
     assert(predictions.getUncertainty().get.forall(_ == 0.0))
@@ -163,10 +177,12 @@ class BaggerTest {
     */
   @Test
   def testScores(): Unit = {
+    rng.setSeed(24795L)
+
     val csv = TestUtils.readCsv("double_example.csv")
     val trainingData = csv.map(vec => (vec.init, vec.last.asInstanceOf[Double]))
-    val DTLearner = RegressionTreeLearner()
-    val baggedLearner = Bagger(DTLearner, numBags = trainingData.size * 16) // use lots of trees to reduce noise
+    val DTLearner = RegressionTreeLearner(rng = rng)
+    val baggedLearner = Bagger(DTLearner, numBags = trainingData.size * 16, randBasis = TestUtils.getBreezeRandBasis(rng.nextLong())) // use lots of trees to reduce noise
     val RF = baggedLearner.train(trainingData).getModel()
 
     /* Call transform on the training data */
@@ -185,19 +201,21 @@ class BaggerTest {
     */
   @Test
   def calibrationTimeTest(): Unit = {
+    rng.setSeed(24795L)
+
     val trainingData = TestUtils.binTrainingData(
       TestUtils.generateTrainingData(1024, 12, noise = 0.1, function = Friedman.friedmanSilverman),
       inputBins = Seq((0, 8))
     )
-    val DTLearner = RegressionTreeLearner(numFeatures = 3)
+    val DTLearner = RegressionTreeLearner(numFeatures = 3, rng = rng)
     val start = System.nanoTime()
-    Bagger(DTLearner, numBags = trainingData.size, uncertaintyCalibration = false)
+    Bagger(DTLearner, numBags = trainingData.size, uncertaintyCalibration = false, randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
       .train(trainingData)
       .getModel()
     val unCalibratedTime = 1.0e-9 * (System.nanoTime() - start)
 
     val startAgain = System.nanoTime()
-    Bagger(DTLearner, numBags = 64, uncertaintyCalibration = true)
+    Bagger(DTLearner, numBags = 64, uncertaintyCalibration = true, randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
       .train(trainingData)
       .getModel()
     val calibratedTime = 1.0e-9 * (System.nanoTime() - startAgain)
@@ -210,9 +228,11 @@ class BaggerTest {
     */
   @Test
   def testInterrupt(): Unit = {
+    rng.setSeed(24795L)
+
     val trainingData = TestUtils.generateTrainingData(2048, 12, noise = 0.1, function = Friedman.friedmanSilverman)
-    val DTLearner = RegressionTreeLearner(numFeatures = 3)
-    val baggedLearner = Bagger(DTLearner, numBags = trainingData.size)
+    val DTLearner = RegressionTreeLearner(numFeatures = 3, rng = rng)
+    val baggedLearner = Bagger(DTLearner, numBags = trainingData.size, randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
 
     // Create a future to run train
     val tmpPool = Executors.newFixedThreadPool(1)
@@ -260,6 +280,8 @@ class BaggerTest {
     */
   @Test
   def testSmallDataRecalibration(): Unit = {
+    rng.setSeed(24795L)
+
     // Define a simple, binary function and create training data
     def stepFunction(x: Seq[Double]): Double = Math.floor(2 * x(0))
 
@@ -270,29 +292,41 @@ class BaggerTest {
      * then the model will fail to train
      */
     val DTLearner = RegressionTreeLearner(
-      leafLearner = Some(GuessTheMeanLearner()),
+      leafLearner = Some(GuessTheMeanLearner(rng = rng)),
       numFeatures = 2,
-      splitter = RegressionSplitter(randomizePivotLocation = true)
+      splitter = RegressionSplitter(randomizePivotLocation = true, rng = rng),
+      rng = rng
     )
-    val trainedModel: BaggedModel[Any] = Bagger(DTLearner,
-      numBags = 16, useJackknife = true, uncertaintyCalibration = true)
+    val trainedModel: BaggedModel[Any] = Bagger(DTLearner, numBags = 16, useJackknife = true,
+      uncertaintyCalibration = true, randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
       .train(trainingData)
       .getModel()
+
+    try {
+      val _: BaggedModel[Any] = Bagger(DTLearner, numBags = 16, useJackknife = true,
+        uncertaintyCalibration = true, disableBootstrap = true, randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
+        .train(trainingData)
+        .getModel()
+      fail("Setting both uncertaintyCalibration and disableBootstrap should throw an exception.")
+    } catch {
+      case _: Throwable =>
+    }
   }
 
   /**
     * Test that the uncertainty is always positive (and non-zero)
     *
     * This happens randomly, so let's repeat a test many times to make sure we catch it.  On my machine, this fails
-    * in the first couple thousand times and takes runs for 13 seconds once its resolved, so I don't think
+    * in the first couple thousand times and takes runs for 13 seconds once it's resolved, so I don't think
     * that's too much overhead.
     */
   @Test
   def testUncertaintyFloor(): Unit = {
+    rng.setSeed(24795L)
     (0 until 16384).foreach { idx =>
-      val trainingData = TestUtils.generateTrainingData(16, 5, noise = 0.0, function = Friedman.friedmanSilverman, seed = Random.nextLong())
-      val DTLearner = RegressionTreeLearner(numFeatures = 2)
-      val sigma = Bagger(DTLearner, numBags = 7)
+      val trainingData = TestUtils.generateTrainingData(16, 5, noise = 0.0, function = Friedman.friedmanSilverman, seed = rng.nextLong())
+      val DTLearner = RegressionTreeLearner(numFeatures = 2, rng = rng)
+      val sigma = Bagger(DTLearner, numBags = 7, randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
         .train(trainingData)
         .getModel()
         .transform(trainingData.map(_._1))
@@ -308,10 +342,11 @@ class BaggerTest {
     */
   @Test
   def testUncertaintyFloorWithBias(): Unit = {
+    rng.setSeed(24795L)
     (0 until 1024).foreach { idx =>
-      val trainingData = TestUtils.generateTrainingData(16, 5, noise = 0.0, function = Friedman.friedmanSilverman, seed = Random.nextLong())
-      val DTLearner = RegressionTreeLearner(numFeatures = 2)
-      val sigma = Bagger(DTLearner, numBags = 7, biasLearner = Some(GuessTheMeanLearner()))
+      val trainingData = TestUtils.generateTrainingData(16, 5, noise = 0.0, function = Friedman.friedmanSilverman, seed = rng.nextLong())
+      val DTLearner = RegressionTreeLearner(numFeatures = 2, rng = rng)
+      val sigma = Bagger(DTLearner, numBags = 7, biasLearner = Some(GuessTheMeanLearner(rng = rng)), randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
         .train(trainingData)
         .getModel()
         .transform(trainingData.map(_._1))
@@ -319,6 +354,39 @@ class BaggerTest {
       assert(sigma.forall(_ > 0.0), s"Found an predicted uncertainty of ${sigma.min} during trial $idx")
     }
   }
+
+  /**
+    * Test Shapley values are correctly averaged over trees.
+    */
+  @Test
+  def testShapley(): Unit = {
+    rng.setSeed(24795L)
+    val nCols = 5
+    val trainingData = TestUtils.generateTrainingData(64, nCols, noise = 0.0, function = Friedman.friedmanSilverman, seed = rng.nextLong())
+    val DTLearner = RegressionTreeLearner(numFeatures = nCols, rng = rng)
+    val model = Bagger(DTLearner, randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
+      .train(trainingData)
+      .getModel()
+    val trees = model.getModels()
+    trainingData.foreach { case (x, _) =>
+      val shapley = model.shapley(x).get
+
+      // Do a quick sanity check on the output format.
+      assert(shapley.cols == nCols)
+      assert(shapley.rows == 1)
+
+      // Compute the mean shap value over trees and ensure the bagged model gives the same result.
+      val treeMean = (1.0 / trees.length) * trees.map { t =>
+        t.shapley(x).get
+      }.reduce[DenseMatrix[Double]] { case (a: DenseMatrix[Double], b: DenseMatrix[Double]) =>
+        a +:+ b
+      }
+      val atol = 1e-8
+      assert(
+        (treeMean - shapley).toDenseVector.toScalaVector.forall { x => Math.abs(x) < atol }
+      )
+    }
+  }
 }
 
 /**
@@ -331,10 +399,40 @@ object BaggerTest {
     * @param argv args
     */
   def main(argv: Array[String]): Unit = {
-    new BaggerTest()
-      .testUncertaintyCalibration()
+    measureShapleyPerf()
   }
 
+  /**
+    * Simple driver for running a performance test of BaggedModel.shapley().
+    */
+  def measureShapleyPerf(): Unit = {
+    val rng = new Random(278345L)
+    val pw = new PrintWriter(new File(s"/tmp/shapley-perf_${new Random().nextInt()}.tsv"))
+    pw.write("nCols\tnRows\trepNum\trowIdx\tns\n")
+    (256 to 1024 by 256).foreach { nCols =>
+      (4 to 8 by 2).foreach { nRowsLog: Int =>
+        val nRows = 1 << nRowsLog
+        (1 to 3).foreach { repNum =>
+          val trainingData = TestUtils.generateTrainingData(nRows, nCols, noise = 0.0, function = Friedman.friedmanSilverman, seed = rng.nextLong())
+          val DTLearner = RegressionTreeLearner(numFeatures = nCols, rng = rng)
+          println(s"Training model nCols=$nCols\tnRows=$nRows\trepNum=$repNum")
+          val model = Bagger(DTLearner, randBasis = TestUtils.getBreezeRandBasis(rng.nextLong()))
+            .train(trainingData)
+            .getModel()
+          println(s"Trained")
+
+          rng.shuffle(trainingData).take(16).zipWithIndex.foreach { case (x, i) =>
+            val t0 = System.nanoTime()
+            val shapley = model.shapley(x._1).get
+            val t1 = System.nanoTime()
+            pw.write(s"$nCols\t$nRows\t$repNum\t$i\t${t1-t0}\n")
+            pw.flush()
+          }
+        }
+      }
+    }
+    pw.close()
+  }
 
   def getStandardRMSE(testSet: Seq[(Vector[Any], Double)], model: BaggedModel[Any]): Double = {
     val predictions = model.transform(testSet.map(_._1))
diff --git a/src/test/scala/io/citrine/lolo/bags/MultiTaskBaggerTest.scala b/src/test/scala/io/citrine/lolo/bags/MultiTaskBaggerTest.scala
index 990f5f8c..1f86d998 100644
--- a/src/test/scala/io/citrine/lolo/bags/MultiTaskBaggerTest.scala
+++ b/src/test/scala/io/citrine/lolo/bags/MultiTaskBaggerTest.scala
@@ -18,6 +18,7 @@ import scala.util.Random
   */
 @Test
 class MultiTaskBaggerTest {
+  val rng = new Random(37895L)
 
   /**
     * Test that we get a reasonable output on a single regression problem
@@ -31,7 +32,7 @@ class MultiTaskBaggerTest {
     val inputs = trainingData.map(_._1)
     val labels = trainingData.map(_._2)
     val DTLearner = MultiTaskTreeLearner()
-    val baggedLearner = MultiTaskBagger(DTLearner, numBags = trainingData.size)
+    val baggedLearner = MultiTaskBagger(DTLearner, numBags = trainingData.size, randBasis = TestUtils.getBreezeRandBasis(10478L))
     val RFMeta = baggedLearner.train(inputs, Seq(labels)).head
     val RF = RFMeta.getModel()
 
@@ -51,11 +52,10 @@ class MultiTaskBaggerTest {
   @Test
   def testBaggedMultiTaskGetUncertainty(): Unit = {
     val noiseLevel = 100.0
-    val rng = new Random(237485L)
 
     Seq(MultiTaskTreeLearner()).foreach { baseLearner =>
       // These are in Seqs as a convenience for repurposing this test as a diagnostic tool.
-      Seq(128).foreach { nRows =>
+      Seq(64).foreach { nRows =>
         Seq(16).foreach { nCols =>
           Seq(2).map { n => n * nRows }.foreach { nBags =>
             // Used for error output.
@@ -66,7 +66,7 @@ class MultiTaskBaggerTest {
               val trainingData = trainingDataTmp.map { x => (x._1, x._2 + noiseLevel * rng.nextDouble()) }
               val inputs = trainingData.map(_._1)
               val labels = trainingData.map(_._2)
-              val baggedLearner = MultiTaskBagger(baseLearner, numBags = nBags, uncertaintyCalibration = true)
+              val baggedLearner = MultiTaskBagger(baseLearner, numBags = nBags, uncertaintyCalibration = true, randBasis = TestUtils.getBreezeRandBasis(7835178L))
               val RFMeta = baggedLearner.train(inputs, Seq(labels)).head
               val RF = RFMeta.getModel()
               val results = RF.transform(trainingData.take(4).map(_._1))
@@ -144,7 +144,7 @@ class MultiTaskBaggerTest {
     val inputs = trainingData.map(_._1)
     val labels = trainingData.map(_._2)
     val DTLearner = MultiTaskTreeLearner()
-    val baggedLearner = MultiTaskBagger(DTLearner, numBags = trainingData.size)
+    val baggedLearner = MultiTaskBagger(DTLearner, numBags = trainingData.size, randBasis = TestUtils.getBreezeRandBasis(478L))
     val RFMeta = baggedLearner.train(inputs, Seq(labels)).head
     val RF = RFMeta.getModel()
 
@@ -176,7 +176,7 @@ class MultiTaskBaggerTest {
     val realLabel: Seq[Double] = raw.map(_._2)
     val catLabel: Seq[Boolean] = raw.map(_._2 > realLabel.max / 2.0)
     val DTLearner = MultiTaskTreeLearner()
-    val baggedLearner = MultiTaskBagger(DTLearner, numBags = inputs.size, biasLearner = Some(new RegressionTreeLearner(maxDepth = 2)))
+    val baggedLearner = MultiTaskBagger(DTLearner, numBags = inputs.size, biasLearner = Some(new RegressionTreeLearner(maxDepth = 2)), randBasis = TestUtils.getBreezeRandBasis(78495L))
     val RFMeta = baggedLearner.train(inputs, Seq(realLabel, catLabel)).last
     val RF = RFMeta.getModel()
 
@@ -196,14 +196,14 @@ class MultiTaskBaggerTest {
     val realLabel: Seq[Double] = raw.map(_._2)
     val catLabel: Seq[Boolean] = raw.map(_._2 > realLabel.max / 2.0)
     val sparseCat = catLabel.map(x =>
-      if (Random.nextDouble() > 0.125) {
+      if (rng.nextDouble() > 0.125) {
         null
       } else {
         x
       }
     )
     val sparseReal = realLabel.map(x =>
-      if (Random.nextDouble() > 0.5) {
+      if (rng.nextDouble() > 0.5) {
         Double.NaN
       } else {
         x
@@ -211,7 +211,7 @@ class MultiTaskBaggerTest {
     )
 
     val DTLearner = MultiTaskTreeLearner()
-    val baggedLearner = MultiTaskBagger(DTLearner, numBags = inputs.size, biasLearner = Some(new GuessTheMeanLearner))
+    val baggedLearner = MultiTaskBagger(DTLearner, numBags = inputs.size, biasLearner = Some(GuessTheMeanLearner(rng = rng)), randBasis = TestUtils.getBreezeRandBasis(7839L))
     val trainingResult = baggedLearner.train(inputs, Seq(sparseReal, sparseCat))
     val RFMeta = trainingResult.last
     val RF = RFMeta.getModel()
diff --git a/src/test/scala/io/citrine/lolo/learners/ExtraRandomTreesTest.scala b/src/test/scala/io/citrine/lolo/learners/ExtraRandomTreesTest.scala
new file mode 100644
index 00000000..c5401c9f
--- /dev/null
+++ b/src/test/scala/io/citrine/lolo/learners/ExtraRandomTreesTest.scala
@@ -0,0 +1,182 @@
+package io.citrine.lolo.learners
+
+import breeze.stats.distributions.Beta
+import io.citrine.lolo.TestUtils
+import io.citrine.lolo.stats.functions.Friedman
+import org.junit.Test
+
+import scala.util.Random
+
+@Test
+class ExtraRandomTreesTest {
+  val rng = new Random(92345L)
+
+  /**
+    * Test that the regression does the same thing as the regression bagger
+    */
+  @Test
+  def testRegression(): Unit = {
+    rng.setSeed(982534L)
+    val trainingData = TestUtils.binTrainingData(
+      TestUtils.generateTrainingData(256, 5, noise = 0.1, function = Friedman.friedmanSilverman, seed = rng.nextLong()),
+      inputBins = Seq((0, 8))
+    )
+
+    Seq(true, false).foreach { randomlyRotateFeatures =>
+      val RFMeta = ExtraRandomTrees(randomlyRotateFeatures = randomlyRotateFeatures, rng = rng)
+        .train(trainingData)
+      val RF = RFMeta.getModel()
+
+      val loss = RFMeta.getLoss().get
+      assert(0 <= loss && loss < 1e-8, "Expected zero loss.")
+
+      val results = RF.transform(trainingData.map(_._1))
+      // val means = results.getExpected()
+      val sigma: Seq[Double] = results.getUncertainty(observational = false).get.asInstanceOf[Seq[Double]]
+      assert(sigma.forall(_ >= 0.0))
+
+      assert(results.getUncertainty(observational = true).isEmpty,
+        "Observational uncertainty should be empty, since jackknife is disabled.")
+
+      assert(results.getGradient().isEmpty, "Returned a gradient when there shouldn't be one")
+
+      if (!randomlyRotateFeatures) {
+        /* The first feature should be the most important */
+        val importances = RFMeta.getFeatureImportance().get
+        assert(importances(1) == importances.max)
+      }
+    }
+  }
+
+  /**
+    * Test that the regression does the same thing as the classification bagger
+    */
+  @Test
+  def testClassification(): Unit = {
+    rng.setSeed(25723478834L)
+
+    val nTrain = 64
+    val nTest = nTrain
+    val nBins = 8
+    val trainingData = TestUtils.binTrainingData(
+      TestUtils.generateTrainingData(nTrain, 5, noise = 0.0, function = Friedman.friedmanSilverman, seed = rng.nextLong()),
+      responseBins = Some(nBins)
+    )
+
+    /* Generate small perturbations from the training data */
+    val testData = trainingData.map { r =>
+      (r._1.map { _.asInstanceOf[Double] + 0.05 * rng.nextGaussian() }, r._2)
+    }.take(nTest)
+
+    Seq(true, false).foreach { randomlyRotateFeatures =>
+      Seq(true, false).foreach { disableBootstrap =>
+        Seq(1, 2).foreach { minLeafInstances =>
+          rng.setSeed(238834L)
+
+          val RFMeta = ExtraRandomTrees(numTrees = trainingData.size * 4, randomlyRotateFeatures = randomlyRotateFeatures,
+            disableBootstrap = disableBootstrap, rng = rng, minLeafInstances = minLeafInstances
+          ).train(trainingData)
+          val RF = RFMeta.getModel()
+
+          /* Inspect the results on the training set */
+          val trueRateTrainingSet = trainingData.map(_._2).zip(RF.transform(trainingData.map(_._1)).getExpected()).count { case (a, p) => a == p }.toDouble / nTrain
+          if (minLeafInstances == 1 && disableBootstrap) {
+            assert(trueRateTrainingSet == 1.0)
+          } else {
+            assert(trueRateTrainingSet > 0.9)
+          }
+
+          /* Inspect the results on the perturbed test set */
+          val results = RF.transform(testData.map(_._1))
+          val means = results.getExpected()
+          val trueRateTestSet = testData.map(_._2).zip(means).count { case (a, p) => a == p }.toDouble / nTest
+          assert(trueRateTestSet > 0.5)
+
+          /* Check that class probabilities are reasonable */
+          val uncertainty = results.getUncertainty()
+          assert(uncertainty.isDefined)
+          assert(testData.map(_._2).zip(uncertainty.get).forall { case (a, probs) =>
+            val classProbabilities = probs.asInstanceOf[Map[Any, Double]]
+            val maxProb = classProbabilities(a)
+            maxProb > 1.0/(2*nBins) && maxProb < 1.0 && Math.abs(classProbabilities.values.sum - 1.0) < 1.0e-6
+          })
+        }
+      }
+    }
+  }
+
+  /**
+    * Ensure classification isn't biased toward one or the other ordering of class labels
+    * when there are duplicate inputs.
+    */
+  @Test
+  def testClassificationUnbiased(): Unit = {
+    rng.setSeed(257834L)
+    val numTrials = 20
+    val (winsSuffixed, winsPrefixed): (Int,Int) = (0 until numTrials).map { _ =>
+      val nTrain = 64
+      val nTest = 16
+      val (mainTrainingData, testData) = {
+        val allData = TestUtils.binTrainingData(
+          TestUtils.generateTrainingData(nTrain + nTest, 5, noise = 0.1, function = Friedman.friedmanSilverman, seed = rng.nextLong()),
+          responseBins = Some(2)
+        )
+        (allData.take(nTrain), allData.takeRight(nTest))
+      }
+      val dupeLabel = "DUPE"
+      val trainingDataSuffixed = mainTrainingData ++ Seq(
+        (mainTrainingData.head._1, dupeLabel)
+      )
+      val trainingDataPrefixed = Seq(
+        (mainTrainingData.head._1, dupeLabel)
+      ) ++ mainTrainingData
+
+      val RFSuffixed = ExtraRandomTrees(numTrees = trainingDataSuffixed.size * 2, rng = rng).train(trainingDataSuffixed)
+      val RFPrefixed = ExtraRandomTrees(numTrees = trainingDataPrefixed.size * 2, rng = rng).train(trainingDataPrefixed)
+      val predictedSuffixed = RFSuffixed.getModel().transform(testData.map(_._1))
+      val predictedPrefixed = RFPrefixed.getModel().transform(testData.map(_._1))
+      val extraLabelCountSuffixed = predictedSuffixed.getExpected().count { case p: String => p == dupeLabel }
+      val extraLabelCountPrefixed = predictedPrefixed.getExpected().count { case p: String => p == dupeLabel }
+
+      if (extraLabelCountSuffixed > extraLabelCountPrefixed) {
+        (1,0)
+      } else if (extraLabelCountSuffixed < extraLabelCountPrefixed) {
+        (0,1)
+      } else {
+        (0,0)
+      }
+    }.asInstanceOf[Seq[(Int,Int)]].reduce{(a: (Int,Int),b: (Int,Int))=>(a._1 + b._1, a._2 + b._2)}
+
+    // Posterior beta distribution with Jeffreys prior.
+    val d = new Beta(winsSuffixed + 0.5, winsPrefixed + 0.5)
+    val l = d.inverseCdf(2e-6)
+    val r = d.inverseCdf(1-2e-6)
+    val tol = 1e-2
+    assert(l < 0.5 - tol, f"Bias detected toward prefixed duplicate rows: rate ${d.mean}%.3f (1e-6 CI ${l}%.3f - ${r}%.3f) should be close to 0.5")
+    assert(r > 0.5 + tol, f"Bias detected toward suffixed duplicate rows: rate ${d.mean}%.3f (1e-6 CI ${l}%.3f - ${r}%.3f) should be close to 0.5")
+  }
+
+  /**
+    * Make sure that we can draw training weights consistently even when the training size is small
+    */
+  @Test
+  def testWeightsWithSmallData(): Unit = {
+    rng.setSeed(7834L)
+    val trainingData = TestUtils.generateTrainingData(8, 1, seed = rng.nextLong())
+    // the number of trees is the number of times we generate weights
+    // so this has the effect of creating lots of different sets of weights
+    val learner = ExtraRandomTrees(numTrees = 16384, rng = rng)
+    // the test is that this training doesn't throw an exception
+    learner.train(trainingData).getModel()
+  }
+
+}
+
+object ExtraRandomTreesTest {
+  def main(argv: Array[String]): Unit = {
+    new ExtraRandomTreesTest()
+      .testClassification()
+    new ExtraRandomTreesTest()
+      .testClassificationUnbiased()
+  }
+}
diff --git a/src/test/scala/io/citrine/lolo/learners/RandomForestTest.scala b/src/test/scala/io/citrine/lolo/learners/RandomForestTest.scala
index 108be621..86abcf25 100644
--- a/src/test/scala/io/citrine/lolo/learners/RandomForestTest.scala
+++ b/src/test/scala/io/citrine/lolo/learners/RandomForestTest.scala
@@ -1,15 +1,20 @@
 package io.citrine.lolo.learners
 
+import breeze.linalg.DenseMatrix
 import breeze.stats.distributions.Beta
 import io.citrine.lolo.TestUtils
 import io.citrine.lolo.stats.functions.Friedman
 import org.junit.Test
+import org.scalatest.Assertions._
+
+import scala.util.Random
 
 /**
   * Created by maxhutch on 1/9/17.
   */
 @Test
 class RandomForestTest {
+  val rng = new Random(92345L)
 
   /**
     * Test that the regression forest does the same thing as the regression bagger
@@ -17,19 +22,19 @@ class RandomForestTest {
   @Test
   def testRegressionForest(): Unit = {
     val trainingData = TestUtils.binTrainingData(
-      TestUtils.generateTrainingData(1024, 12, noise = 0.1, function = Friedman.friedmanSilverman),
+      TestUtils.generateTrainingData(1024, 12, noise = 0.1, function = Friedman.friedmanSilverman, seed = rng.nextLong()),
       inputBins = Seq((0, 8))
     )
 
     Seq(true, false).foreach { randomlyRotateFeatures =>
-      val RFMeta = RandomForest(randomlyRotateFeatures = randomlyRotateFeatures)
+      val RFMeta = RandomForest(randomlyRotateFeatures = randomlyRotateFeatures, rng = rng)
         .train(trainingData)
       val RF = RFMeta.getModel()
 
       assert(RFMeta.getLoss().get < 1.0, "Loss of bagger is larger than expected")
 
       val results = RF.transform(trainingData.map(_._1))
-      val means = results.getExpected()
+      // val means = results.getExpected()
       val sigma: Seq[Double] = results.getUncertainty().get.asInstanceOf[Seq[Double]]
       assert(sigma.forall(_ >= 0.0))
 
@@ -49,12 +54,12 @@ class RandomForestTest {
   @Test
   def testClassificationForest(): Unit = {
     val trainingData = TestUtils.binTrainingData(
-      TestUtils.generateTrainingData(1024, 12, noise = 0.1, function = Friedman.friedmanSilverman),
+      TestUtils.generateTrainingData(128, 12, noise = 0.1, function = Friedman.friedmanSilverman, seed = rng.nextLong()),
       inputBins = Seq((0, 8)), responseBins = Some(8)
     )
 
     Seq(true, false).foreach { randomlyRotateFeatures =>
-      val RFMeta = new RandomForest(numTrees = trainingData.size * 2, randomlyRotateFeatures = randomlyRotateFeatures)
+      val RFMeta = RandomForest(numTrees = trainingData.size * 2, randomlyRotateFeatures = randomlyRotateFeatures, rng = rng)
         .train(trainingData)
       val RF = RFMeta.getModel()
 
@@ -79,10 +84,10 @@ class RandomForestTest {
     */
   @Test
   def testClassificationForestUnbiased(): Unit = {
-    val numTrials = 30
-    val (winsSuffixed, winsPrefixed): (Int,Int) = (0 until numTrials).map {case i: Int =>
+    val numTrials = 20
+    val (winsSuffixed, winsPrefixed): (Int,Int) = (0 until numTrials).map { _ =>
       val mainTrainingData = TestUtils.binTrainingData(
-        TestUtils.generateTrainingData(128, 5, noise = 0.1, function = Friedman.friedmanSilverman, seed = i),
+        TestUtils.generateTrainingData(64, 5, noise = 0.1, function = Friedman.friedmanSilverman, seed = rng.nextLong()),
         responseBins = Some(2)
       )
       val dupeLabel = "DUPE"
@@ -93,8 +98,8 @@ class RandomForestTest {
         (mainTrainingData.head._1, dupeLabel)
       ) ++ mainTrainingData
 
-      val RFSuffixed = RandomForest(numTrees = trainingDataSuffixed.size * 2).train(trainingDataSuffixed)
-      val RFPrefixed = RandomForest(numTrees = trainingDataPrefixed.size * 2).train(trainingDataPrefixed)
+      val RFSuffixed = RandomForest(numTrees = trainingDataSuffixed.size * 2, rng = rng).train(trainingDataSuffixed)
+      val RFPrefixed = RandomForest(numTrees = trainingDataPrefixed.size * 2, rng = rng).train(trainingDataPrefixed)
       val predictedSuffixed = RFSuffixed.getModel().transform(mainTrainingData.map(_._1))
       val predictedPrefixed = RFPrefixed.getModel().transform(mainTrainingData.map(_._1))
       val extraLabelCountSuffixed = predictedSuffixed.getExpected().count { case p: String => p == dupeLabel }
@@ -127,10 +132,10 @@ class RandomForestTest {
     // Generate a linear signal in one dimension: 2 * x
     val trainingData: Seq[(Vector[Double], Double)] = TestUtils.generateTrainingData(32, 1, function = {x =>
       x.head * 2.0
-    })
+    }, seed = rng.nextLong())
 
     // Create a consistent set of parameters
-    val baseForest = RandomForest(numTrees = 16384, useJackknife = false)
+    val baseForest = RandomForest(numTrees = 16384, useJackknife = false, rng = rng)
 
     // Turn off split randomization and compute the loss (out-of-bag error)
     val lossWithoutRandomization: Double = baseForest.copy(randomizePivotLocation = false)
@@ -153,11 +158,85 @@ class RandomForestTest {
     val trainingData = TestUtils.generateTrainingData(8, 1)
     // the number of trees is the number of times we generate weights
     // so this has the effect of creating lots of different sets of weights
-    val learner = RandomForest(numTrees = 16384)
+    val learner = RandomForest(numTrees = 16384, rng = rng)
     // the test is that this training doesn't throw an exception
     learner.train(trainingData).getModel()
   }
 
+  def shapleyCompare(
+                      trainingData: Seq[(Vector[Any],Double)],
+                      evalLocation: Vector[Any],
+                      expected: Vector[Double],
+                      rtol: Double = 5e-2
+                    ): Unit = {
+    val actual = RandomForest(rng = rng).train(trainingData).getModel().shapley(evalLocation) match {
+      case None => fail("Unexpected None returned by shapley.")
+      case x: Option[DenseMatrix[Double]] => {
+        val a = x.get
+        assert(a.cols == trainingData.head._1.length, "Expected one Shapley value per feature.")
+        assert(a.rows == 1, "Expected a single output dimension.")
+        a.toDenseVector.toScalaVector
+      }
+      case _ => fail("Unexpected return type.")
+    }
+    expected.zip(actual).foreach {
+      case (e: Double, a: Double) => assert(Math.abs(e - a)/a < rtol)
+    }
+  }
+
+  /**
+    * Test Shapley value for  a simple tree.
+    */
+  @Test
+  def testShapley(): Unit = {
+    rng.setSeed(3751L)
+
+    // Example from Lundberg paper (https://arxiv.org/pdf/1802.03888.pdf)
+    val trainingData1 = Seq(
+      (Vector(1.0, 1.0), 80.0),
+      (Vector(1.0, 0.0), 0.0),
+      (Vector(0.0, 1.0), 0.0),
+      (Vector(0.0, 0.0), 0.0)
+    )
+    val expected1 = Vector(30.0, 30.0)
+    shapleyCompare((1 to 8).map(_=>trainingData1).flatten, Vector[Any](1.0, 1.0), expected1)
+
+    // Second example from Lundberg paper (https://arxiv.org/pdf/1802.03888.pdf)
+    val trainingData2 = Seq(
+      (Vector(1.0, 1.0), 90.0),
+      (Vector(1.0, 0.0), 10.0),
+      (Vector(0.0, 1.0), 0.0),
+      (Vector(0.0, 0.0), 0.0)
+    )
+    val expected2 = Vector(35.0, 30.0)
+    shapleyCompare((1 to 8).map(_=>trainingData2).flatten, Vector[Any](1.0, 1.0), expected2)
+
+    // Example with two splits on one feature
+    // Worked out with pen-and-paper from Lundberg Equation 2.
+    val trainingData3 = Seq(
+      (Vector(1.0, 1.0), 100.0),
+      (Vector(1.0, 0.0), 80.0),
+      (Vector(1.0, 0.2), 70.0),
+      (Vector(0.0, 1.0), 0.0),
+      (Vector(0.0, 0.2), 0.0),
+      (Vector(0.0, 0.0), 0.0)
+    )
+    val expected3 = Vector(45.8333333333333, 12.5)
+    shapleyCompare((1 to 8).map(_=>trainingData3).flatten, Vector[Any](1.0, 1.0), expected3)
+
+    // Example with 5 features, to exercise all the factorials in Lundberg equation 2.
+    // Referenced against the shap package on a sklearn decision tree.
+    val trainingData4 = Seq(
+      (Vector(0.0, 0.0, 0.0, 0.0, 0.0), 1.0),
+      (Vector(1.0, 0.0, 0.0, 0.0, 0.0), 2.0),
+      (Vector(0.0, 1.0, 0.0, 0.0, 0.0), 4.0),
+      (Vector(0.0, 0.0, 1.0, 0.0, 0.0), 8.0),
+      (Vector(0.0, 0.0, 0.0, 1.0, 0.0), 16.0),
+      (Vector(0.0, 0.0, 0.0, 0.0, 1.0), 32.0)
+    )
+    val expected4 = Vector(0.0333333333333333, 0.2, 0.8666666666666667, 3.533333333333333, 16.866666666666667)
+    shapleyCompare((1 to 8).map{_=>trainingData4}.flatten, Vector.fill[Any](5)(1.0), expected4)
+  }
 }
 
 object RandomForestTest {
@@ -166,5 +245,7 @@ object RandomForestTest {
       .testClassificationForest()
     new RandomForestTest()
       .testClassificationForestUnbiased()
+    new RandomForestTest()
+      .testShapley()
   }
 }
diff --git a/src/test/scala/io/citrine/lolo/linear/LinearRegressionTest.scala b/src/test/scala/io/citrine/lolo/linear/LinearRegressionTest.scala
index 9b30a1ad..bdb3a79e 100644
--- a/src/test/scala/io/citrine/lolo/linear/LinearRegressionTest.scala
+++ b/src/test/scala/io/citrine/lolo/linear/LinearRegressionTest.scala
@@ -10,6 +10,7 @@ import scala.util.Random
   */
 @Test
 class LinearRegressionTest {
+  val rng = new Random(783578L)
   /** Number of training rows */
   val n = 6
   /** Number of features in each row */
@@ -46,7 +47,7 @@ class LinearRegressionTest {
     */
   @Test
   def testRegression(): Unit = {
-    val int0 = Random.nextDouble()
+    val int0 = rng.nextDouble()
     val result = data * beta0 + int0
 
     val trainingData = (0 until n).map { i =>
@@ -69,7 +70,7 @@ class LinearRegressionTest {
     */
   @Test
   def testWeightedRegression(): Unit = {
-    val int0 = Random.nextDouble()
+    val int0 = rng.nextDouble()
     val result = data * beta0 + int0
 
     val trainingData = (0 until n).map { i =>
@@ -77,7 +78,7 @@ class LinearRegressionTest {
     }
 
     val lr = new LinearRegressionLearner()
-    val lrm = lr.train(trainingData, weights = Some(Seq.tabulate(n)(i => Math.abs(Random.nextDouble()))))
+    val lrm = lr.train(trainingData, weights = Some(Seq.tabulate(n)(i => Math.abs(rng.nextDouble()))))
     val model = lrm.getModel()
     val output = model.transform(trainingData.map(_._1))
     val predicted = output.getExpected()
@@ -121,7 +122,7 @@ class LinearRegressionTest {
     */
   @Test
   def testCategoricalValue(): Unit = {
-    val int0 = Random.nextDouble()
+    val int0 = rng.nextDouble()
     val result = data * beta0 + int0
 
     val trainingData = (0 until n).map { i =>
diff --git a/src/test/scala/io/citrine/lolo/stats/correlations/DistanceCorrelationTest.scala b/src/test/scala/io/citrine/lolo/stats/correlations/DistanceCorrelationTest.scala
index 798ef1ab..d59e741c 100644
--- a/src/test/scala/io/citrine/lolo/stats/correlations/DistanceCorrelationTest.scala
+++ b/src/test/scala/io/citrine/lolo/stats/correlations/DistanceCorrelationTest.scala
@@ -9,6 +9,7 @@ import scala.util.Random
   */
 @Test
 class DistanceCorrelationTest {
+  val rng = new Random(36783677836L)
 
   val N: Int = 512
 
@@ -16,8 +17,8 @@ class DistanceCorrelationTest {
 
   @Test
   def testIndependent(): Unit = {
-    val X: Seq[Double] = Seq.tabulate(N)(i => Random.nextGaussian())
-    val Y: Seq[Double] = Seq.tabulate(N)(i => Random.nextGaussian())
+    val X: Seq[Double] = Seq.tabulate(N)(i => rng.nextGaussian())
+    val Y: Seq[Double] = Seq.tabulate(N)(i => rng.nextGaussian())
     val dcorr = DistanceCorrelation.distanceCorrelation(X, Y, dist)
 
     assert(dcorr < 0.15, s"dCorr for independent vars is ${dcorr}")
@@ -26,7 +27,7 @@ class DistanceCorrelationTest {
 
   @Test
   def testLinear(): Unit = {
-    val X: Seq[Double] = Seq.tabulate(N)(i => Random.nextDouble())
+    val X: Seq[Double] = Seq.tabulate(N)(i => rng.nextDouble())
     val Y: Seq[Double] = X.map(v => v * 10)
     val dcorr = DistanceCorrelation.distanceCorrelation(X, Y, dist)
 
@@ -35,7 +36,7 @@ class DistanceCorrelationTest {
 
   @Test
   def testQuadratic(): Unit = {
-    val X: Seq[Double] = Seq.tabulate(N)(i => Random.nextDouble() - .5)
+    val X: Seq[Double] = Seq.tabulate(N)(i => rng.nextDouble() - .5)
     val Y: Seq[Double] = X.map(v => v * v * 2)
     val dcorr = DistanceCorrelation.distanceCorrelation(X, Y, dist)
 
@@ -44,7 +45,7 @@ class DistanceCorrelationTest {
 
   @Test
   def testSin(): Unit = {
-    val X: Seq[Double] = Seq.tabulate(N)(i => Random.nextDouble() - .5)
+    val X: Seq[Double] = Seq.tabulate(N)(i => rng.nextDouble() - .5)
     val Y: Seq[Double] = X.map(v => Math.sin(2 * Math.PI * v))
     val dcorr = DistanceCorrelation.distanceCorrelation(X, Y, dist)
 
diff --git a/src/test/scala/io/citrine/lolo/stats/metrics/ClassificationMetricsTest.scala b/src/test/scala/io/citrine/lolo/stats/metrics/ClassificationMetricsTest.scala
index 6bcc4db4..04d5aa08 100644
--- a/src/test/scala/io/citrine/lolo/stats/metrics/ClassificationMetricsTest.scala
+++ b/src/test/scala/io/citrine/lolo/stats/metrics/ClassificationMetricsTest.scala
@@ -7,23 +7,23 @@ import scala.util.Random
 /**
   * Created by maxhutch on 12/28/16.
   */
-@Test
 class ClassificationMetricsTest {
+  val rng = new Random(9230L)
 
   /**
     * Test that the metric works on a sparse confusion matrix
     */
+  @Test
   def testSparse(): Unit = {
     val N = 512
     /* Make random predictions */
     val pva = Seq.tabulate(N) { i =>
-      (Vector(0.0), Random.nextInt(N).toString, Random.nextInt(N).toString)
+      (Vector(0.0), rng.nextInt(N).toString, rng.nextInt(N).toString)
     }
     val loss = ClassificationMetrics.f1scores(pva)
     /* The loss should be in (0.0, 1.0) */
-    assert(!loss.isNaN)
-    assert(loss > 0.0)
-    assert(loss < 1.0)
+    assert(!loss.isNaN, s"Loss was NaN")
+    assert(loss >= 0.0, s"Loss was less than 0: ${loss}")
+    assert(loss <= 1.0, s"Loss was greater than 1: ${loss}")
   }
-
 }
diff --git a/src/test/scala/io/citrine/lolo/transformers/FeatureRotatorTest.scala b/src/test/scala/io/citrine/lolo/transformers/FeatureRotatorTest.scala
index 269f4232..c677d072 100644
--- a/src/test/scala/io/citrine/lolo/transformers/FeatureRotatorTest.scala
+++ b/src/test/scala/io/citrine/lolo/transformers/FeatureRotatorTest.scala
@@ -18,11 +18,13 @@ import scala.util.Random
 @Test
 class FeatureRotatorTest {
 
+  val rng = new Random(1297843L)
+
   val data: Seq[(Vector[Any], Any)] = TestUtils.binTrainingData(
     TestUtils.generateTrainingData(1024, 12, noise = 0.1, function = Friedman.friedmanSilverman),
     inputBins = Seq((0, 8))
   )
-  val weights: Vector[Double] = Vector.fill(data.size)(if (Random.nextBoolean()) Random.nextDouble() else 0.0)
+  val weights: Vector[Double] = Vector.fill(data.size)(if (rng.nextBoolean()) rng.nextDouble() else 0.0)
 
   @Test
   def testRandomRotation(): Unit = {
@@ -106,11 +108,11 @@ class FeatureRotatorTest {
     */
   @Test
   def testRotatedGTM(): Unit = {
-    val learner = GuessTheMeanLearner()
+    val learner = GuessTheMeanLearner(rng = rng)
     val model = learner.train(data).getModel()
     val result = model.transform(data.map(_._1)).getExpected()
 
-    val rotatedLearner = FeatureRotator(GuessTheMeanLearner())
+    val rotatedLearner = FeatureRotator(GuessTheMeanLearner(rng = rng))
     val rotatedModel = rotatedLearner.train(data).getModel()
     val rotatedResult = rotatedModel.transform(data.map(_._1)).getExpected()
 
@@ -245,7 +247,7 @@ class FeatureRotatorTest {
     val catLabel = inputs.map(x => Friedman.friedmanGrosseSilverman(x) > 15.0)
 
     // Sparsify the categorical labels
-    val sparseCatLabel = catLabel.map(x => if (Random.nextBoolean()) x else null)
+    val sparseCatLabel = catLabel.map(x => if (rng.nextBoolean()) x else null)
 
     // Train and evaluate rotated models on original and rotated features
     val baseLearner = new MultiTaskTreeLearner()
diff --git a/src/test/scala/io/citrine/lolo/transformers/StandardizerTest.scala b/src/test/scala/io/citrine/lolo/transformers/StandardizerTest.scala
index b32246ca..045b4178 100644
--- a/src/test/scala/io/citrine/lolo/transformers/StandardizerTest.scala
+++ b/src/test/scala/io/citrine/lolo/transformers/StandardizerTest.scala
@@ -15,9 +15,10 @@ import scala.util.Random
   */
 @Test
 class StandardizerTest {
+  val rng = new Random(823478686L)
 
   val data: Vector[(Vector[Double], Double)] = TestUtils.generateTrainingData(1024, 12, noise = 0.1, function = Friedman.friedmanSilverman)
-  val weights: Vector[Double] = Vector.fill(data.size)(if (Random.nextBoolean()) Random.nextDouble() else 0.0)
+  val weights: Vector[Double] = Vector.fill(data.size)(if (rng.nextBoolean()) rng.nextDouble() else 0.0)
 
   // Creating another dataset which has 1 feature that has 0 variance.
   val dataWithConstant: Vector[(Vector[Double], Double)] = data.map(d => (0.0 +: d._1, d._2))
@@ -53,11 +54,11 @@ class StandardizerTest {
     */
   @Test
   def testStandardGTM(): Unit = {
-    val learner = GuessTheMeanLearner()
+    val learner = GuessTheMeanLearner(rng = rng)
     val model = learner.train(data).getModel()
     val result = model.transform(data.map(_._1)).getExpected()
 
-    val standardLearner = Standardizer(GuessTheMeanLearner())
+    val standardLearner = Standardizer(GuessTheMeanLearner(rng = rng))
     val standardModel = standardLearner.train(data).getModel()
     val standardResult = standardModel.transform(data.map(_._1)).getExpected()
 
@@ -177,7 +178,7 @@ class StandardizerTest {
     val catLabel = inputs.map(x => Friedman.friedmanGrosseSilverman(x) > 15.0)
 
     // Sparsify the categorical labels
-    val sparseCatLabel = catLabel.map(x => if (Random.nextBoolean()) x else null)
+    val sparseCatLabel = catLabel.map(x => if (rng.nextBoolean()) x else null)
 
     // Screw up the scale of the double labels
     val scale = 10000.0
diff --git a/src/test/scala/io/citrine/lolo/trees/FeaturePowerSetTermsTest.scala b/src/test/scala/io/citrine/lolo/trees/FeaturePowerSetTermsTest.scala
new file mode 100644
index 00000000..c6c13503
--- /dev/null
+++ b/src/test/scala/io/citrine/lolo/trees/FeaturePowerSetTermsTest.scala
@@ -0,0 +1,102 @@
+package io.citrine.lolo.trees
+
+import breeze.linalg.{DenseVector, norm}
+import org.junit.Test
+import org.scalatest.Assertions._
+
+import scala.util.Random
+
+@Test
+class FeaturePowerSetTermsTest {
+
+  /**
+    * Test a few simple extend cases explicitly
+    */
+  @Test
+  def testExtend() {
+    val path = new FeaturePowerSetTerms(2)
+    assert(path.size == 0)
+
+    path.extend(0.5, 1.0)
+    assert(path.size == 1)
+    assert(Math.abs(path.weightBySubsetSize(0) - 0.25) < 1e-9)
+    assert(Math.abs(path.weightBySubsetSize(1) - 0.5) < 1e-9)
+
+    path.extend(0.1, 0.0)
+    assert(path.size == 2)
+    assert(Math.abs(path.weightBySubsetSize(0) - 1.0 / 60.0) < 1e-9)
+    assert(Math.abs(path.weightBySubsetSize(1) - 1.0 / 60.0) < 1e-9)
+    assert(Math.abs(path.weightBySubsetSize(2) - 0.0) < 1e-9)
+  }
+
+  /**
+    * Test that the order that we extend by doesn't matter
+    */
+  @Test
+  def testExtendOrderInvariant(): Unit = {
+    val N = 8
+    val rng = new Random(378654L)
+    val factors = Seq.fill(N)((rng.nextDouble(), if (rng.nextBoolean()) 1.0 else 0.0))
+
+    val results = Seq.fill(N) {
+      val set = new FeaturePowerSetTerms(N)
+      rng.shuffle(factors).foreach(x => set.extend(x._1, x._2))
+      // make dense vectors so we can take norms
+      new DenseVector(set.weightBySubsetSize)
+    }
+
+    assert(results.forall { x => norm(x - results.head) < 1.0e-9 })
+  }
+
+  @Test
+  def testCopy(): Unit = {
+    val path = new FeaturePowerSetTerms(3)
+    var newPath = path.copy()
+    assert(path.size == newPath.size)
+
+    path.extend(0.99, 1.0)
+    assert(path.size == newPath.size + 1)
+
+    path.extend(0.5, 0.0)
+    assert(path.size == newPath.size + 2)
+
+    path.extend(0.1, 1.0)
+    assert(path.size == newPath.size + 3)
+
+    newPath = path.copy()
+    path.weightBySubsetSize.zip(newPath.weightBySubsetSize).foreach { case (p, np) =>
+      assert(p == np)
+    }
+  }
+
+  /**
+    * Test that unwind will reverse an extend
+    */
+  @Test
+  def testUnwind(): Unit = {
+    val N = 8
+    val rng = new Random(83647L)
+    val firstFactors = Seq.fill(N)((rng.nextDouble(), if (rng.nextBoolean()) 1.0 else 0.0))
+    val secondFactors = Seq.fill(N)((rng.nextDouble(), if (rng.nextBoolean()) 1.0 else 0.0))
+
+    // add both sets of factors to the set (note that we test that extend order doesn't matter already
+    val set = new FeaturePowerSetTerms(firstFactors.size + secondFactors.size)
+    firstFactors.foreach(x => set.extend(x._1, x._2))
+    secondFactors.foreach(x => set.extend(x._1, x._2))
+
+    // Now only add the second set
+    val anotherSet = new FeaturePowerSetTerms(secondFactors.size)
+    secondFactors.foreach(x => anotherSet.extend(x._1, x._2))
+
+    val thinedSet = firstFactors.foldLeft(set) { case (tmp, x) =>
+      tmp.unwind(x._1, x._2)
+    }
+
+    val anotherWeights = new DenseVector(anotherSet.weightBySubsetSize)
+    val weights = new DenseVector(thinedSet.weightBySubsetSize.take(anotherWeights.length))
+
+    assert(norm(weights - anotherWeights) < 1.0e-9)
+    assert(thinedSet.weightBySubsetSize.drop(anotherWeights.size).forall(x => x == 0.0))
+    assert(thinedSet.size == anotherSet.size)
+  }
+}
diff --git a/src/test/scala/io/citrine/lolo/trees/classification/ClassificationTreeTest.scala b/src/test/scala/io/citrine/lolo/trees/classification/ClassificationTreeTest.scala
index 0db00743..383b1926 100644
--- a/src/test/scala/io/citrine/lolo/trees/classification/ClassificationTreeTest.scala
+++ b/src/test/scala/io/citrine/lolo/trees/classification/ClassificationTreeTest.scala
@@ -89,13 +89,7 @@ class ClassificationTreeTest {
 
     /* The first feature should be the most important */
     val importances = DTMeta.getFeatureImportance().get
-    println(importances)
     assert(importances.slice(0, 5).min > importances.slice(5, importances.size).max)
-
-    /* Test serialization */
-    val tmpFile: File = File.createTempFile("tmp", ".csv")
-    val oos = new ObjectOutputStream(new FileOutputStream(tmpFile))
-    oos.writeObject(DT)
   }
 
   /**
diff --git a/src/test/scala/io/citrine/lolo/trees/multitask/MultiTaskTreeTest.scala b/src/test/scala/io/citrine/lolo/trees/multitask/MultiTaskTreeTest.scala
index f13f52c2..cc07b733 100644
--- a/src/test/scala/io/citrine/lolo/trees/multitask/MultiTaskTreeTest.scala
+++ b/src/test/scala/io/citrine/lolo/trees/multitask/MultiTaskTreeTest.scala
@@ -15,9 +15,10 @@ import scala.util.Random
   */
 @Test
 class MultiTaskTreeTest {
+  val rng = new Random(1012795L)
 
   /* Setup some data */
-  val raw: Seq[(Vector[Double], Double)] = TestUtils.generateTrainingData(1024, 12, noise = 0.1, function = Friedman.friedmanSilverman)
+  val raw: Seq[(Vector[Double], Double)] = TestUtils.generateTrainingData(1024, 12, noise = 0.1, function = Friedman.friedmanSilverman, seed = rng.nextLong())
   val inputs: Seq[Vector[Double]] = raw.map(_._1)
   val realLabel: Seq[Double] = raw.map(_._2)
   val catLabel: Seq[Boolean] = raw.map(_._2 > realLabel.max / 2.0)
@@ -27,7 +28,7 @@ class MultiTaskTreeTest {
     */
   @Test
   def testTwoLabels(): Unit = {
-    val learner = new MultiTaskTreeLearner()
+    val learner = MultiTaskTreeLearner()
     val models = learner.train(inputs, Seq(realLabel, catLabel)).map(_.getModel())
     assert(models.size == 2)
     assert(models.head.isInstanceOf[RegressionTree])
@@ -45,14 +46,14 @@ class MultiTaskTreeTest {
   @Test
   def testSparseExact(): Unit = {
     val sparseCat = catLabel.map(x =>
-      if (Random.nextBoolean()) {
+      if (rng.nextBoolean()) {
         null
       } else {
         x
       }
     )
 
-    val learner = new MultiTaskTreeLearner()
+    val learner = MultiTaskTreeLearner()
     val models = learner.train(inputs, Seq(realLabel, sparseCat)).map(_.getModel())
     val realResults = models.head.transform(inputs).getExpected().asInstanceOf[Seq[Double]]
     val catResults = models.last.transform(inputs).getExpected().asInstanceOf[Seq[Boolean]]
@@ -67,7 +68,7 @@ class MultiTaskTreeTest {
   @Test
   def testSparseQuantitative(): Unit = {
     val sparseCat = catLabel.map(x =>
-      if (Random.nextDouble() > 0.125) {
+      if (rng.nextDouble() > 0.125) {
         null
       } else {
         x
diff --git a/src/test/scala/io/citrine/lolo/trees/regression/RegressionTreeTest.scala b/src/test/scala/io/citrine/lolo/trees/regression/RegressionTreeTest.scala
index 07c35a97..ea3cb2ea 100644
--- a/src/test/scala/io/citrine/lolo/trees/regression/RegressionTreeTest.scala
+++ b/src/test/scala/io/citrine/lolo/trees/regression/RegressionTreeTest.scala
@@ -2,6 +2,7 @@ package io.citrine.lolo.trees.regression
 
 import java.io.{File, FileOutputStream, ObjectOutputStream}
 
+import breeze.linalg.DenseMatrix
 import io.citrine.lolo.TestUtils
 import io.citrine.lolo.linear.LinearRegressionLearner
 import io.citrine.lolo.stats.functions.Friedman
@@ -16,6 +17,7 @@ import scala.util.Random
   */
 @Test
 class RegressionTreeTest {
+  val rng = new Random(90476L)
 
   /**
     * Trivial models with no splits should have finite feature importance.
@@ -215,7 +217,7 @@ class RegressionTreeTest {
     val linearLearner = LinearRegressionLearner(regParam = Some(1.0))
     val DTLearner = RegressionTreeLearner(leafLearner = Some(linearLearner), maxDepth = 1)
     val DTMeta = DTLearner.train(trainingData, weights = Some(Seq.fill(trainingData.size) {
-      Random.nextInt(8)
+      rng.nextInt(8)
     }))
     val DT = DTMeta.getModel()
 
@@ -245,6 +247,94 @@ class RegressionTreeTest {
     }
   }
 
+  /**
+    * Convenience method to compare the Shapley value at a test location for a tree trained on trainingData to a known reference.
+    *
+    * @param trainingData on which to train a regression tree
+    * @param evalLocation at which to evaluate the Shapley value
+    * @param expected     known reference Shapley value to which to compare (fails upon mismatch)
+    */
+  def shapleyCompare(
+                     trainingData: Seq[(Vector[Double],Double)],
+                     evalLocation: Vector[Any],
+                     expected: Vector[Double],
+                     omitFeatures: Set[Int] = Set()
+                    ): Unit = {
+    val actual = RegressionTreeLearner().train(trainingData).getModel().shapley(evalLocation, omitFeatures) match {
+      case None => fail("Unexpected None returned by shapley.")
+      case x: Option[DenseMatrix[Double]] => {
+        val a = x.get
+        assert(a.cols == trainingData.head._1.length, "Expected one Shapley value per feature.")
+        assert(a.rows == 1, "Expected a single output dimension.")
+        a.toDenseVector.toScalaVector
+      }
+      case _ => fail("Unexpected return type.")
+    }
+    assert(expected.zip(actual).forall{
+      case (e: Double, a: Double) => Math.abs(e - a) < 1e-12
+    }, s"Shapley value ${actual} does not match reference ${expected}.")
+  }
+
+  /**
+    * Test Shapley value for  a simple tree.
+    */
+  @Test
+  def testShapley(): Unit = {
+    // Example from Lundberg paper (https://arxiv.org/pdf/1802.03888.pdf)
+    val trainingData1 = Seq(
+      (Vector(1.0, 1.0), 80.0),
+      (Vector(1.0, 0.0), 0.0),
+      (Vector(0.0, 1.0), 0.0),
+      (Vector(0.0, 0.0), 0.0)
+    )
+    val expected1 = Vector(30.0, 30.0)
+    shapleyCompare(trainingData1, Vector[Any](1.0, 1.0), expected1)
+
+    // Second example from Lundberg paper (https://arxiv.org/pdf/1802.03888.pdf)
+    val trainingData2 = Seq(
+      (Vector(1.0, 1.0), 90.0),
+      (Vector(1.0, 0.0), 10.0),
+      (Vector(0.0, 1.0), 0.0),
+      (Vector(0.0, 0.0), 0.0)
+    )
+    val expected2 = Vector(35.0, 30.0)
+    shapleyCompare(trainingData2, Vector[Any](1.0, 1.0), expected2)
+
+    // Example with two splits on one feature
+    // Worked out with pen-and-paper from Lundberg Equation 2.
+    val trainingData3 = Seq(
+      (Vector(1.0, 1.0), 100.0),
+      (Vector(1.0, 0.0), 80.0),
+      (Vector(1.0, 0.2), 70.0),
+      (Vector(0.0, 1.0), 0.0),
+      (Vector(0.0, 0.2), 0.0),
+      (Vector(0.0, 0.0), 0.0)
+    )
+    val expected3 = Vector(45.8333333333333, 12.5)
+    shapleyCompare(trainingData3, Vector[Any](1.0, 1.0), expected3)
+
+    // Example with 5 features, to exercise all the factorials in Lundberg equation 2.
+    // Referenced against the shap package on a sklearn decision tree.
+    val trainingData4 = Seq(
+      (Vector(0.0, 0.0, 0.0, 0.0, 0.0), 1.0),
+      (Vector(1.0, 0.0, 0.0, 0.0, 0.0), 2.0),
+      (Vector(0.0, 1.0, 0.0, 0.0, 0.0), 4.0),
+      (Vector(0.0, 0.0, 1.0, 0.0, 0.0), 8.0),
+      (Vector(0.0, 0.0, 0.0, 1.0, 0.0), 16.0),
+      (Vector(0.0, 0.0, 0.0, 0.0, 1.0), 32.0)
+    )
+    val expected4 = Vector(0.0333333333333333, 0.2, 0.8666666666666667, 3.533333333333333, 16.866666666666667)
+    shapleyCompare(trainingData4, Vector.fill[Any](5)(1.0), expected4)
+
+    // Test omitted features
+    val expected2a = Vector(0.0, 20.0)
+    shapleyCompare(trainingData2, Vector[Any](1.0, 1.0), expected2a, omitFeatures = Set(0))
+    val expected2b = Vector(25.0, 0.0)
+    shapleyCompare(trainingData2, Vector[Any](1.0, 1.0), expected2b, omitFeatures = Set(1))
+
+    // Ensure we don't crash when restricting number of features.
+    RegressionTreeLearner(numFeatures = 1).train(trainingData4).getModel().shapley(Vector.fill[Any](5)(0.0), Set())
+  }
 }
 
 /** Companion driver */
diff --git a/src/test/scala/io/citrine/lolo/trees/splits/BoltzmannSplitterTest.scala b/src/test/scala/io/citrine/lolo/trees/splits/BoltzmannSplitterTest.scala
index 167f62e2..4d83a43c 100644
--- a/src/test/scala/io/citrine/lolo/trees/splits/BoltzmannSplitterTest.scala
+++ b/src/test/scala/io/citrine/lolo/trees/splits/BoltzmannSplitterTest.scala
@@ -5,15 +5,16 @@ import org.junit.Test
 import scala.util.Random
 
 class BoltzmannSplitterTest {
+  val rng = new Random(45321L)
 
   /**
     * Test that uniform labels result in "NoSplit" with zero reduced impurity
     */
   @Test
   def testZeroVariance(): Unit = {
-    val splitter = BoltzmannSplitter(1.0e-9)
+    val splitter = BoltzmannSplitter(1.0e-9, rng = rng)
     val testData = Seq.fill(64){
-      val x = Random.nextDouble()
+      val x = rng.nextDouble()
       val y = 1.0
       val weight = 1.0
       (Vector(x), y, weight)
@@ -32,10 +33,10 @@ class BoltzmannSplitterTest {
     */
   @Test
   def testLowVariance(): Unit = {
-    val splitter = BoltzmannSplitter(1.0e-18)
+    val splitter = BoltzmannSplitter(1.0e-18, rng = rng)
     val testData = Seq.fill(256){
-      val x = Random.nextDouble()
-      val y = Random.nextGaussian() * 1.0e-9 + 1.0
+      val x = rng.nextDouble()
+      val y = rng.nextGaussian() * 1.0e-9 + 1.0
       val weight = 1.0
       (Vector(x), y, weight)
     }
diff --git a/src/test/scala/io/citrine/lolo/trees/splits/ExtraRandomRegressionSplitterTest.scala b/src/test/scala/io/citrine/lolo/trees/splits/ExtraRandomRegressionSplitterTest.scala
new file mode 100644
index 00000000..0f2e08bd
--- /dev/null
+++ b/src/test/scala/io/citrine/lolo/trees/splits/ExtraRandomRegressionSplitterTest.scala
@@ -0,0 +1,134 @@
+package io.citrine.lolo.trees.splits
+
+import io.citrine.lolo.TestUtils
+import org.junit.Test
+
+import scala.util.Random
+
+class ExtraRandomRegressionSplitterTest {
+
+  /**
+    * Test that uniform labels result in "NoSplit" with zero reduced impurity
+    */
+  @Test
+  def testZeroVariance(): Unit = {
+    val rng = new Random(278345L)
+    val splitter = ExtraRandomRegressionSplitter(rng)
+    val testData = Seq.fill(64){
+      val x = rng.nextDouble()
+      val y = 1.0
+      val weight = 1.0
+      (Vector(x), y, weight)
+    }
+
+    val (bestSplit, improvement) = splitter.getBestSplit(testData, 1, 1)
+    assert(bestSplit.isInstanceOf[NoSplit], "Can't split on zero variance")
+    assert(improvement == 0, "Can't improve on no variance")
+  }
+
+  /**
+    * Test that very small random signals don't lead to exceptions.
+    *
+    * Note that this test will return NoSplit occasionally, when the starting impurity is numerically zero.  That's OK
+    * though.
+    */
+  @Test
+  def testLowVariance(): Unit = {
+    val rng = new Random(278345L)
+    val splitter = ExtraRandomRegressionSplitter(rng)
+    val testData = Seq.fill(256){
+      val x = rng.nextDouble()
+      val y = rng.nextGaussian() * 1.0e-9 + 1.0
+      val weight = 1.0
+      (Vector(x), y, weight)
+    }
+
+    splitter.getBestSplit(testData, 1, 1)
+  }
+
+  /**
+    * Test correctness of the split.
+    *
+    * Specifically, this test ensures that ExtraRandomSplitter.getBestSplit selects the split
+    */
+  @Test
+  def testCorrectSplit(): Unit = {
+    val rng = new Random(34682L)
+    val baseGrid = Seq(0.0, 1.2)
+    Seq(1, 5, 9).foreach { numFeatures =>
+      // Arrange training input data as a regular grid, and then as a Gaussian-distributed sample.
+      Seq(TestUtils.enumerateGrid(Seq.fill(numFeatures)(baseGrid)), Seq.fill(64)(Vector.fill(numFeatures)(rng.nextGaussian()))).foreach { xTrain =>
+        // Have ExtraRandomSplitter.getBestSplit choose a split from this number of randomly-chosen features.
+        (1 to numFeatures by Math.max(numFeatures/2,1)).foreach { numFeaturesToConsider =>
+          // Repeat with the same training inputs a few times with different random number streams.
+          (1 to 4).foreach { repetitionNumber =>
+            // Feature indices, which will be used to determine which features will be considered by the splitter (and in which order).
+            // This is stored as a Vector because rng.shuffle was having trouble with a Seq; possibly a type erasure issue.
+            val featureIndices = (0 until numFeatures).toVector
+
+            // Generate linear training data, using xTrain from above as inputs, starting by selecting coefficients for each feature.
+            val featureCoefficients = rng.shuffle(featureIndices).map { i => Math.pow((i + 1).toDouble, 2.0) }
+            val trainingData = xTrain.map { x =>
+              val y = x.zip(featureCoefficients).map { case (a, b) => a * b }.sum
+              val weight = 1.0
+              (x, y, weight)
+            }
+
+            // Determine exactly what shuffle() and nextDouble() will return within getBestSplit by resetting rng to a common rng seed.
+            val sharedSeed = 238745L + repetitionNumber
+            rng.setSeed(sharedSeed)
+            // This is the set of shuffled indices to return when getBestSplit calls the RNG's shuffle method.
+            val shuffledFeatureIndices = rng.shuffle(featureIndices)
+            // Precompute the sequence of numbers returned when getBestSplit calls the RNG's nextDouble method to select cut points.
+            val randomUniforms = Seq.fill(numFeatures)(rng.nextDouble())
+
+            // Compute where the cut points should be placed, based on the the sequence of randomUniforms.
+            val cutPoints = shuffledFeatureIndices.zip(randomUniforms).map { case (i, u) =>
+              val x = trainingData.map(_._1(i))
+              val xmin = x.min
+              val xmax = x.max
+              val cutPoint = xmin + u * (xmax - xmin)
+              (i, cutPoint)
+            }.sortBy(_._1).map(_._2)
+
+            // Compute the sum of variances across partitions of each possible cut.
+            val varianceSums = featureIndices.map { k =>
+              trainingData.groupBy { v => v._1(k) < cutPoints(k) }.flatMap { case (_, subset) =>
+                val mean = subset.map(_._2).sum / subset.length
+                subset.map { case (_, yi, _) =>
+                  Math.pow(yi - mean, 2)
+                }
+              }.sum
+            }
+
+            // Choose the index on which a split results in the lowest sum of variances.
+            val indexOfBest = shuffledFeatureIndices.take(numFeaturesToConsider).map { i =>
+              (i, varianceSums(i))
+            }.minBy(_._2)._1
+
+            // Instantiate the splitter to test, passing in the random number generator that is reset to its former state used above.
+            rng.setSeed(sharedSeed)
+            val splitter = ExtraRandomRegressionSplitter(rng)
+            // Ask the splitter what it chose as a split.
+            val bestSplit = splitter.getBestSplit(trainingData, numFeaturesToConsider, 1)
+
+            // Finally, we can ensure that the index on which we split is correct...
+            val testCaveatMessage = "NOTE: this test may inaccurately fail due to changes in the sequence of rng calls."
+            assert(bestSplit._1.getIndex() == indexOfBest, s"Incorrect index of best split. $testCaveatMessage")
+            // Do a sanity check about the directionality of turnLeft to ensure this test is valid.
+            assert(bestSplit._1.turnLeft(Vector.fill(numFeatures)(xTrain.flatten.min)),
+              s"Unexpected directionality of turnLeft, which is probably a bug with the test itself.")
+            assert(!bestSplit._1.turnLeft(Vector.fill(numFeatures)(xTrain.flatten.max)),
+              s"Unexpected directionality of turnLeft, which is probably a bug with the test itself.")
+            // And ensure that we're turning at the right place.
+            assert(bestSplit._1.turnLeft(Vector.fill(numFeatures)(cutPoints(indexOfBest) - 1e-8)),
+              s"Split not placed at correct location. $testCaveatMessage")
+            assert(!bestSplit._1.turnLeft(Vector.fill(numFeatures)(cutPoints(indexOfBest) + 1e-8)),
+              s"Split not placed at correct location. $testCaveatMessage")
+          }
+        }
+      }
+    }
+  }
+
+}
diff --git a/src/test/scala/io/citrine/lolo/trees/splits/MultiTaskSplitterTest.scala b/src/test/scala/io/citrine/lolo/trees/splits/MultiTaskSplitterTest.scala
index 73ae73e1..4505cc20 100644
--- a/src/test/scala/io/citrine/lolo/trees/splits/MultiTaskSplitterTest.scala
+++ b/src/test/scala/io/citrine/lolo/trees/splits/MultiTaskSplitterTest.scala
@@ -3,12 +3,15 @@ package io.citrine.lolo.trees.splits
 import io.citrine.lolo.trees.impurity.MultiImpurityCalculator
 import org.junit.Test
 
+import scala.util.Random
+
 
 /**
   * Created by maxhutch on 12/1/16.
   */
 @Test
 class MultiTaskSplitterTest {
+  val rng = new Random(38945L)
 
   /**
     * Test that the real split goes in the right place
@@ -23,7 +26,8 @@ class MultiTaskSplitterTest {
     )
     val calculator = MultiImpurityCalculator.build(data.map(_._2), data.map(_._3))
 
-    val (pivot, impurity) = MultiTaskSplitter.getBestRealSplit(data, calculator, 0, 1)
+    val splitter = MultiTaskSplitter(rng = rng)
+    val (pivot, impurity) = splitter.getBestRealSplit(data, calculator, 0, 1)
     val turns = data.map(x => pivot.turnLeft(x._1))
     assert(turns == Seq(true, false, false, false))
     assert(impurity == 0.5)
@@ -42,7 +46,8 @@ class MultiTaskSplitterTest {
     )
     val calculator = MultiImpurityCalculator.build(data.map(_._2), data.map(_._3))
 
-    val (pivot, impurity) = MultiTaskSplitter.getBestCategoricalSplit(data, calculator, 0, 1)
+    val splitter = MultiTaskSplitter(rng = rng)
+    val (pivot, impurity) = splitter.getBestCategoricalSplit(data, calculator, 0, 1)
     val turns = data.map(x => pivot.turnLeft(x._1)).take(3)
     println(turns)
     assert(turns == Seq(true, false, false))
@@ -63,7 +68,8 @@ class MultiTaskSplitterTest {
     val data = inputs.indices.map { i =>
       (inputs(i), labels(i).asInstanceOf[Array[AnyVal]], weights(i))
     }
-    val (pivot, impurity) = MultiTaskSplitter.getBestSplit(data, data.head._1.size, 1)
+    val splitter = MultiTaskSplitter(rng = rng)
+    val (pivot, impurity) = splitter.getBestSplit(data, data.head._1.size, 1)
     assert(!pivot.isInstanceOf[NoSplit])
   }
 }
diff --git a/src/test/scala/io/citrine/lolo/trees/splits/SplitterTest.scala b/src/test/scala/io/citrine/lolo/trees/splits/SplitterTest.scala
index abc8e436..54cc99aa 100644
--- a/src/test/scala/io/citrine/lolo/trees/splits/SplitterTest.scala
+++ b/src/test/scala/io/citrine/lolo/trees/splits/SplitterTest.scala
@@ -10,14 +10,15 @@ import scala.util.Random
   */
 @Test
 class SplitterTest {
+  val rng = new Random(39578L)
 
   /**
     * Test different variance formulations
     */
   @Test
   def testTotalVarianceCalculation(): Unit = {
-    val values = (0 until 16).map(i => Random.nextDouble())
-    val (left, right) = values.partition(f => Random.nextBoolean())
+    val values = (0 until 16).map(i => rng.nextDouble())
+    val (left, right) = values.partition(f => rng.nextBoolean())
 
     val totalSquareSum = values.map(v => v * v).sum
     val totalNum = values.size
@@ -54,7 +55,7 @@ class SplitterTest {
   def testLargeDuplicates(): Unit = {
     val base: Double = 3.0e9
     val trainingData = Seq.fill(8) {
-      (Vector(base + Random.nextDouble()), Random.nextDouble(), 1.0)
+      (Vector(base + rng.nextDouble()), rng.nextDouble(), 1.0)
     }
 
     val calculator = new VarianceCalculator(0.0, 0.0, 8.0)
diff --git a/src/test/scala/io/citrine/lolo/validation/CalibrationStudy.scala b/src/test/scala/io/citrine/lolo/validation/CalibrationStudy.scala
index 50b9778e..491f7bcb 100644
--- a/src/test/scala/io/citrine/lolo/validation/CalibrationStudy.scala
+++ b/src/test/scala/io/citrine/lolo/validation/CalibrationStudy.scala
@@ -12,6 +12,7 @@ import scala.collection.JavaConverters._
 import scala.util.Random
 
 object CalibrationStudy {
+  val rng = new Random(372845L)
 
   def main(args: Array[String]): Unit = {
 
@@ -48,7 +49,7 @@ object CalibrationStudy {
                                  ratios: Seq[Int] = Seq(1, 2, 4)
                                ): Unit = {
     val nFeature = 8
-    val data = TestUtils.iterateTrainingData(nFeature, Friedman.friedmanSilverman, seed = Random.nextLong())
+    val data = TestUtils.iterateTrainingData(nFeature, Friedman.friedmanSilverman, seed = rng.nextLong())
 
     ratios.foreach { ratio =>
       val chart = Merit.plotMeritScan(
@@ -57,7 +58,8 @@ object CalibrationStudy {
         Map("R2" -> CoefficientOfDetermination, "confidence" -> StandardConfidence, "standard error / 4" -> StandardError(0.25), "sigmaCorr" -> UncertaintyCorrelation),
         logScale = true,
         yMin = Some(0.0),
-        yMax = Some(1.0)
+        yMax = Some(1.0),
+        rng = rng
       ) { nTrain: Double =>
         val nTree = ratio * nTrain
         val learner = Bagger(
@@ -68,7 +70,7 @@ object CalibrationStudy {
           useJackknife = true,
           uncertaintyCalibration = false
         )
-        StatisticalValidation.generativeValidation[Double](
+        StatisticalValidation(rng = rng).generativeValidation[Double](
           data,
           learner,
           nTrain = nTrain.toInt,
@@ -84,7 +86,7 @@ object CalibrationStudy {
                                         calibrated: Boolean = false
                                       ): Unit = {
     val nFeature = 8
-    val data = TestUtils.iterateTrainingData(nFeature, Friedman.friedmanSilverman, seed = Random.nextLong())
+    val data = TestUtils.iterateTrainingData(nFeature, Friedman.friedmanSilverman, seed = rng.nextLong())
     sizes.foreach { nTrain =>
       val chart = Merit.plotMeritScan(
         "Number of trees",
@@ -92,7 +94,8 @@ object CalibrationStudy {
         Map("R2" -> CoefficientOfDetermination, "confidence" -> StandardConfidence, "standard error / 4" -> StandardError(0.25), "sigmaCorr" -> UncertaintyCorrelation),
         logScale = true,
         yMin = Some(0.0),
-        yMax = Some(1.0)
+        yMax = Some(1.0),
+        rng = rng
       ) { nTree: Double =>
         val learner = Bagger(
           RegressionTreeLearner(
@@ -102,7 +105,7 @@ object CalibrationStudy {
           useJackknife = true,
           uncertaintyCalibration = calibrated
         )
-        StatisticalValidation.generativeValidation[Double](
+        StatisticalValidation(rng = rng).generativeValidation[Double](
           data,
           learner,
           nTrain = nTrain.toInt,
@@ -125,7 +128,7 @@ object CalibrationStudy {
                                      funcName: String = "fs"
                                    ): Unit = {
     val nFeature = 8
-    val data = TestUtils.iterateTrainingData(nFeature, func, seed = Random.nextLong())
+    val data = TestUtils.iterateTrainingData(nFeature, func, seed = rng.nextLong())
     sizes.foreach { nTree =>
       val chart = Merit.plotMeritScan(
         "Number of training points",
@@ -133,7 +136,8 @@ object CalibrationStudy {
         Map("R2" -> CoefficientOfDetermination, "confidence" -> StandardConfidence, "standard error / 4" -> StandardError(0.25), "sigmaCorr" -> UncertaintyCorrelation),
         logScale = true,
         yMin = Some(0.0),
-        yMax = Some(1.0)
+        yMax = Some(1.0),
+        rng = rng
       ) { nTrain: Double =>
         val learner = Bagger(
           RegressionTreeLearner(
@@ -143,7 +147,7 @@ object CalibrationStudy {
           useJackknife = true,
           uncertaintyCalibration = calibrated
         )
-        StatisticalValidation.generativeValidation[Double](
+        StatisticalValidation(rng = rng).generativeValidation[Double](
           data,
           learner,
           nTrain = nTrain.toInt,
@@ -170,7 +174,7 @@ object CalibrationStudy {
                                ignoreDims: Int = 0
                              ): Map[String, (Double, Double)] = {
 
-    val data = TestUtils.iterateTrainingData(nFeature, func, seed = Random.nextLong())
+    val data = TestUtils.iterateTrainingData(nFeature, func, seed = rng.nextLong())
       .map{case (x, y) => (x.drop(ignoreDims), y)}
     val learner = Bagger(
       RegressionTreeLearner(
@@ -182,7 +186,7 @@ object CalibrationStudy {
       biasLearner = None // Some(RegressionTreeLearner(maxDepth = (Math.log(nTrain) / Math.log(2) / 2).toInt))
     )
 
-    val fullStream = StatisticalValidation.generativeValidation[Double](
+    val fullStream = StatisticalValidation(rng = rng).generativeValidation[Double](
       data,
       learner,
       nTrain,
@@ -198,7 +202,7 @@ object CalibrationStudy {
     )
     BitmapEncoder.saveBitmap(chart, s"./stdres_${funcName}-nTrain.${nTrain}-nTree.${nTree}", BitmapFormat.PNG)
 
-    val dataStream = StatisticalValidation.generativeValidation[Double](
+    val dataStream = StatisticalValidation(rng = rng).generativeValidation[Double](
       data,
       learner,
       nTrain,
@@ -208,7 +212,8 @@ object CalibrationStudy {
     BitmapEncoder.saveBitmap(pva, s"./pva_${funcName}-nTrain.${nTrain}-nTree.${nTree}", BitmapFormat.PNG)
     Merit.estimateMerits(
       dataStream.iterator,
-      Map("R2" -> CoefficientOfDetermination, "confidence" -> StandardConfidence, "error / 4" -> StandardError(0.25), "sigmaCorr" -> UncertaintyCorrelation)
+      Map("R2" -> CoefficientOfDetermination, "confidence" -> StandardConfidence, "error / 4" -> StandardError(0.25), "sigmaCorr" -> UncertaintyCorrelation),
+      rng = rng
     )
   }
 
@@ -282,7 +287,8 @@ object CalibrationStudy {
         Map("R2" -> CoefficientOfDetermination, "confidence" -> StandardConfidence, "standard error / 4" -> StandardError(0.25), "sigmaCorr" -> UncertaintyCorrelation),
         logScale = true,
         yMin = Some(0.0),
-        yMax = Some(1.0)
+        yMax = Some(1.0),
+        rng = rng
       ) { nTrain: Double =>
         val learner = Bagger(
           RegressionTreeLearner(
@@ -292,7 +298,7 @@ object CalibrationStudy {
           useJackknife = true,
           uncertaintyCalibration = calibrated
         )
-        StatisticalValidation.generativeValidation[Double](
+        StatisticalValidation(rng = rng).generativeValidation[Double](
           trainingData,
           learner,
           nTrain = nTrain.toInt,
@@ -327,7 +333,7 @@ object CalibrationStudy {
       uncertaintyCalibration = false
     )
 
-    val fullStream = StatisticalValidation.generativeValidation[Double](
+    val fullStream = StatisticalValidation(rng = rng).generativeValidation[Double](
       trainingData,
       learner,
       nTrain,
@@ -342,7 +348,7 @@ object CalibrationStudy {
     )
     BitmapEncoder.saveBitmap(chart, s"./stdres_hcep-nTrain.${nTrain}-nTree.${nTree}", BitmapFormat.PNG)
 
-    val dataStream = StatisticalValidation.generativeValidation[Double](
+    val dataStream = StatisticalValidation(rng = rng).generativeValidation[Double](
       trainingData,
       learner,
       nTrain,
@@ -352,7 +358,8 @@ object CalibrationStudy {
     BitmapEncoder.saveBitmap(pva, s"./pva_hcep-nTrain.${nTrain}-nTree.${nTree}", BitmapFormat.PNG)
     Merit.estimateMerits(
       dataStream.iterator,
-      Map("R2" -> CoefficientOfDetermination, "confidence" -> StandardConfidence, "error / 4" -> StandardError(0.25), "sigmaCorr" -> UncertaintyCorrelation)
+      Map("R2" -> CoefficientOfDetermination, "confidence" -> StandardConfidence, "error / 4" -> StandardError(0.25), "sigmaCorr" -> UncertaintyCorrelation),
+      rng = rng
     )
   }
 
diff --git a/src/test/scala/io/citrine/lolo/validation/CrossValidationTest.scala b/src/test/scala/io/citrine/lolo/validation/CrossValidationTest.scala
index cf0bbdae..e194066f 100644
--- a/src/test/scala/io/citrine/lolo/validation/CrossValidationTest.scala
+++ b/src/test/scala/io/citrine/lolo/validation/CrossValidationTest.scala
@@ -5,7 +5,10 @@ import io.citrine.lolo.learners.RandomForest
 import io.citrine.lolo.stats.functions.Friedman
 import org.junit.Test
 
+import scala.util.Random
+
 class CrossValidationTest {
+  val rng = new Random(92486L)
 
   /**
     * Test that CV results are consistent with out-of-bag estimates from the bagged learner
@@ -15,7 +18,7 @@ class CrossValidationTest {
   @Test
   def testCompareToOutOfBag(): Unit = {
     val learner = RandomForest()
-    val data = TestUtils.generateTrainingData(128, 8, Friedman.friedmanSilverman)
+    val data = TestUtils.generateTrainingData(128, 8, Friedman.friedmanSilverman, seed = rng.nextLong())
 
     val metric = RootMeanSquareError
     val (rmseFromCV, uncertainty) = CrossValidation.kFoldCrossvalidation(data, learner, Map("rmse" -> metric), k = 3)("rmse")
diff --git a/src/test/scala/io/citrine/lolo/validation/MeritTest.scala b/src/test/scala/io/citrine/lolo/validation/MeritTest.scala
index 8d77caa9..f7fd07d6 100644
--- a/src/test/scala/io/citrine/lolo/validation/MeritTest.scala
+++ b/src/test/scala/io/citrine/lolo/validation/MeritTest.scala
@@ -8,6 +8,7 @@ import org.scalatest.Assertions._
 import scala.util.{Random, Try}
 
 class MeritTest {
+  val rng = new Random(34578L)
 
   /**
     * Generate test data by adding Gaussian noise to a uniformly distributed response
@@ -38,9 +39,9 @@ class MeritTest {
 
     Seq.fill(numBatch) {
       val pua = Seq.fill(batchSize) {
-        val y = Random.nextDouble()
+        val y = rng.nextDouble()
         val draw = errorDistribution.sample().toSeq
-        val error: Double = draw(0) * Random.nextGaussian()
+        val error: Double = draw(0) * rng.nextGaussian()
         val uncertainty = if (uncertaintyCorrelation >= maximumCorrelation) {
           Math.abs(draw(0))
         } else {
diff --git a/src/test/scala/io/citrine/lolo/validation/StatisticalValidationTest.scala b/src/test/scala/io/citrine/lolo/validation/StatisticalValidationTest.scala
index ac3f09fb..31bcfb9a 100644
--- a/src/test/scala/io/citrine/lolo/validation/StatisticalValidationTest.scala
+++ b/src/test/scala/io/citrine/lolo/validation/StatisticalValidationTest.scala
@@ -5,19 +5,22 @@ import io.citrine.lolo.learners.RandomForest
 import io.citrine.lolo.stats.functions.Friedman
 import org.junit.Test
 
+import scala.util.Random
+
 class StatisticalValidationTest {
+  val rng = new Random(783443769L)
 
   @Test
   def testCompareToKFolds(): Unit = {
     val learner = RandomForest()
-    val dataSet = TestUtils.generateTrainingData(128, 8, Friedman.friedmanSilverman)
+    val dataSet = TestUtils.generateTrainingData(128, 8, Friedman.friedmanSilverman, seed = rng.nextLong())
     val dataGenerator = TestUtils.iterateTrainingData(8, Friedman.friedmanSilverman)
 
     val metrics = Map("rmse" -> RootMeanSquareError)
     val (rmseFromCV, uncertaintyFromCV) = CrossValidation.kFoldCrossvalidation(dataSet, learner, metrics, k = 4, nTrial = 4)("rmse")
 
     val (rmseFromStats, uncertaintyFromStats) = Merit.estimateMerits(
-      StatisticalValidation.generativeValidation(dataGenerator, learner, 96, 32, 16),
+      StatisticalValidation(rng = rng).generativeValidation(dataGenerator, learner, 96, 32, 16),
       metrics
     )("rmse")