Improve data loading and sampling for SGD (#296)

* Ensures source and target IDs are paired up in SGD for map-like sources
* Can validate map-like sources
* Can use a simple sampling scheme for non-map-like sources

Author: weefuzzy

include/flucoma/algorithms/public/MLP.hpp

@@ -30,8 +30,6 @@ class MLP
   using ArrayXXd = Eigen::ArrayXXd;
 
 public:
-  explicit MLP() = default;
-  ~MLP() = default;
 
   void init(index inputSize, index outputSize,
             FluidTensor<index, 1> hiddenSizes, index hiddenAct, index outputAct)

include/flucoma/algorithms/public/SGD.hpp

@@ -15,10 +15,10 @@ under the European Union’s Horizon 2020 research and innovation programme
 #include "../../data/FluidDataSet.hpp"
 #include "../../data/FluidIndex.hpp"
 #include "../../data/FluidTensor.hpp"
+#include "../../data/SimpleDataSampler.hpp"
 #include "../../data/TensorTypes.hpp"
 #include <Eigen/Core>
 #include <limits>
-#include <random>
 
 namespace fluid {
 namespace algorithm {
@@ -27,72 +27,63 @@ class SGD
 {
   using ArrayXd = Eigen::ArrayXd;
   using ArrayXXd = Eigen::ArrayXXd;
-  using Permutation = Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>;
 
 public:
-  explicit SGD() = default;
-  ~SGD() = default;
 
-  double train(MLP& model, const RealMatrixView in, RealMatrixView out,
+  double train(MLP& model, InputRealMatrixView in, RealMatrixView out,
                index nIter, index batchSize, double learningRate,
                double momentum, double valFrac)
+  {
+    return train(model, in, out,
+                 SimpleDataSampler(in.rows(), batchSize, valFrac, true), nIter,
+                 learningRate, momentum);
+  }
+
+  template <typename Sampler>
+  double train(MLP& model, InputRealMatrixView in, RealMatrixView out,
+               Sampler&& loader, index nIter, double learningRate,
+               double momentum)
   {
     using namespace _impl;
     using namespace std;
     using namespace Eigen;
-    index       nExamples = in.rows();
-    index       inputSize = in.cols();
-    index       outputSize = out.cols();
-    ArrayXXd    input = asEigen<Eigen::Array>(in);
-    ArrayXXd    output = asEigen<Eigen::Array>(out);
-    Permutation valPerm(nExamples);
-    valPerm.setIdentity();
-    shuffle(valPerm.indices().data(),
-            valPerm.indices().data() + valPerm.indices().size(),
-            mt19937{random_device{}()});
-    input = valPerm * input.matrix();
-    output = valPerm * output.matrix();
-    index nVal = std::lround(nExamples * valFrac);
-    index nTrain = nExamples - nVal;
+    index nExamples = in.rows();
+    // index inputSize = in.cols();
+    index outputSize = out.cols();
 
-    ArrayXXd trainInput = input.block(0, 0, nTrain, inputSize);
-    ArrayXXd trainOutput = output.block(0, 0, nTrain, outputSize);
-    ArrayXXd valInput = input.block(nTrain, 0, nVal, inputSize);
-    ArrayXXd valOutput = output.block(nTrain, 0, nVal, outputSize);
+    auto                    valIdx = loader.validationSet();
+    std::optional<ArrayXXd> valInput;
+    std::optional<ArrayXXd> valOutput;
+    if (valIdx)
+    {
+      valInput =
+          ArrayXXd(asEigen<Eigen::Array>(in)(valIdx->col(0), Eigen::all));
+      valOutput =
+          ArrayXXd(asEigen<Eigen::Array>(out)(valIdx->col(1), Eigen::all));
+    }
 
-    Permutation iterPerm(nTrain);
-    iterPerm.setIdentity();
     double error = 0;
     index  patience = mInitialPatience;
     double prevValLoss = std::numeric_limits<double>::max();
     while (nIter-- > 0)
     {
-      shuffle(iterPerm.indices().data(),
-              iterPerm.indices().data() + iterPerm.indices().size(),
-              mt19937{random_device{}()});
-      ArrayXXd inPerm = iterPerm * trainInput.matrix();
-      ArrayXXd outPerm = iterPerm * trainOutput.matrix();
-      for (index batchStart = 0; batchStart < inPerm.rows();
-           batchStart += batchSize)
+      for (auto batch : loader)
       {
-        index thisBatchSize = (batchStart + batchSize) <= nTrain
-                                  ? batchSize
-                                  : nTrain - batchStart;
-        ArrayXXd batchIn =
-            inPerm.block(batchStart, 0, thisBatchSize, inPerm.cols());
+        index    thisBatchSize = batch->rows();
+        ArrayXXd batchIn = asEigen<Eigen::Array>(in)(batch->col(0), Eigen::all);
         ArrayXXd batchOut =
-            outPerm.block(batchStart, 0, thisBatchSize, outPerm.cols());
+            asEigen<Eigen::Array>(out)(batch->col(1), Eigen::all);
         ArrayXXd batchPred = ArrayXXd::Zero(thisBatchSize, outputSize);
         model.forward(batchIn, batchPred);
         ArrayXXd diff = batchPred - batchOut;
         model.backward(diff);
         model.update(learningRate, momentum);
       }
-      if (nVal > 0)
+      if (valIdx)
       {
-        ArrayXXd valPred = ArrayXXd::Zero(nVal, outputSize);
-        model.forward(valInput, valPred);
-        double valLoss = model.loss(valPred, valOutput);
+        ArrayXXd valPred = ArrayXXd::Zero(valInput->rows(), outputSize);
+        model.forward(*valInput, valPred);
+        double valLoss = model.loss(valPred, *valOutput);
         if (valLoss < prevValLoss)
           patience = mInitialPatience;
         else
@@ -101,6 +92,13 @@ class SGD
         prevValLoss = valLoss;
       }
     }
+
+    auto trainingIdx = loader.trainingSet();
+    nExamples = trainingIdx->rows();
+    ArrayXXd input = asEigen<Eigen::Array>(in)(trainingIdx->col(0), Eigen::all);
+    ArrayXXd output =
+        asEigen<Eigen::Array>(out)(trainingIdx->col(1), Eigen::all);
+
     ArrayXXd finalPred = ArrayXXd::Zero(nExamples, outputSize);
     model.forward(input, finalPred);
     bool isNan = !((finalPred == finalPred)).all();

include/flucoma/clients/nrt/CommonResults.hpp

@@ -40,6 +40,7 @@ static const std::string FileRead{"Couldn't read file"};
 static const std::string FileWrite{"Couldn't write file"};
 static const std::string NotImplemented{"Not implemented"};
 static const std::string SizesDontMatch{"Sizes do not match"};
+static const std::string TooFewOutputPoints{"Not enough output points"};
 static const std::string DimensionsDontMatch{"Dimensions do not match"};
 
 template <typename T>

include/flucoma/clients/nrt/MLPClassifierClient.hpp

@@ -16,6 +16,7 @@ under the European Union’s Horizon 2020 research and innovation programme
 #include "../../algorithms/public/LabelSetEncoder.hpp"
 #include "../../algorithms/public/MLP.hpp"
 #include "../../algorithms/public/SGD.hpp"
+#include "../../data/FluidDataSetSampler.hpp"
 #include <string>
 
 namespace fluid {
@@ -149,8 +150,8 @@ class MLPClassifierClient : public FluidBaseClient,
     if (!targetClientPtr) return Error<double>(NoLabelSet);
     auto targetDataSet = targetClientPtr->getLabelSet();
     if (targetDataSet.size() == 0) return Error<double>(EmptyLabelSet);
-    if (sourceDataSet.size() != targetDataSet.size())
-      return Error<double>(SizesDontMatch);
+    if (sourceDataSet.size() > targetDataSet.size())
+      return Error<double>(TooFewOutputPoints);
 
     mAlgorithm.encoder.fit(targetDataSet);
 
@@ -162,24 +163,40 @@ class MLPClassifierClient : public FluidBaseClient,
                           mAlgorithm.encoder.numLabels(), get<kHidden>(),
                           get<kActivation>(), 1); // sigmoid output
     }
-    mAlgorithm.mlp.setTrained(false);
-    DataSet result(1);
-    auto    data = sourceDataSet.getData();
-    auto    tgt = targetDataSet.getData();
-
-    RealMatrix oneHot(targetDataSet.size(), mAlgorithm.encoder.numLabels());
-    oneHot.fill(0);
-    for (index i = 0; i < targetDataSet.size(); i++)
+    
+    if (auto missingIDs = sourceDataSet.checkIDs(targetDataSet);
+        missingIDs.size() == 0)
     {
-      mAlgorithm.encoder.encodeOneHot(tgt.row(i)(0), oneHot.row(i));
-    }
+      mAlgorithm.mlp.setTrained(false);
+      DataSet result(1);
+      auto    data = sourceDataSet.getData();
+      auto    tgt = targetDataSet.getData();
+  
+      RealMatrix oneHot(targetDataSet.size(), mAlgorithm.encoder.numLabels());
+      oneHot.fill(0);
+      for (index i = 0; i < targetDataSet.size(); i++)
+      {
+        mAlgorithm.encoder.encodeOneHot(tgt.row(i)(0), oneHot.row(i));
+      }
+  
+      FluidDataSetSampler sampler(sourceDataSet, targetDataSet,
+        get<kBatchSize>(), get<kVal>(), true);
 
-    algorithm::SGD sgd;
-    double         error =
-        sgd.train(mAlgorithm.mlp, data, oneHot, get<kIter>(), get<kBatchSize>(),
-                  get<kRate>(), get<kMomentum>(), get<kVal>());
+      algorithm::SGD sgd;
+      double error = sgd.train(mAlgorithm.mlp, data, oneHot, sampler, get<kIter>(),
+                               get<kRate>(), get<kMomentum>());
 
-    return error;
+      return error;
+    }
+    else
+    {
+      std::ostringstream oss;
+      oss << "Can't train because these IDs are missing from the output data: ";
+      std::copy(missingIDs.begin(), missingIDs.end() - 1,
+                std::ostream_iterator<std::string>(oss, ","));
+      oss << missingIDs.back();
+      return {Result::Status::kError, oss.str()};
+    }
   }
 
   MessageResult<void> predict(InputDataSetClientRef  srcClient,

include/flucoma/clients/nrt/MLPRegressorClient.hpp

@@ -14,8 +14,10 @@ under the European Union’s Horizon 2020 research and innovation programme
 #include "NRTClient.hpp"
 #include "../../algorithms/public/MLP.hpp"
 #include "../../algorithms/public/SGD.hpp"
+#include "../../data/FluidDataSetSampler.hpp"
 #include <string>
 
+
 namespace fluid {
 namespace client {
 namespace mlpregressor {
@@ -108,8 +110,8 @@ class MLPRegressorClient : public FluidBaseClient,
     if (!targetClientPtr) return Error<double>(NoDataSet);
     auto targetDataSet = targetClientPtr->getDataSet();
     if (targetDataSet.size() == 0) return Error<double>(EmptyDataSet);
-    if (sourceDataSet.size() != targetDataSet.size())
-      return Error<double>(SizesDontMatch);
+    if (sourceDataSet.size() > targetDataSet.size())
+      return Error<double>(TooFewOutputPoints);
     index outputAct = get<kOutputActivation>() == -1 ? get<kActivation>()
                                                      : get<kOutputActivation>();
     if (!mAlgorithm.initialized() ||
@@ -121,15 +123,29 @@ class MLPRegressorClient : public FluidBaseClient,
                       get<kHidden>(), get<kActivation>(), outputAct);
     }
 
-    mAlgorithm.setTrained(false);
-    DataSet        result(1);
-    auto           data = sourceDataSet.getData();
-    auto           tgt = targetDataSet.getData();
-    algorithm::SGD sgd;
-    double         error =
-        sgd.train(mAlgorithm, data, tgt, get<kIter>(), get<kBatchSize>(),
-                  get<kRate>(), get<kMomentum>(), get<kVal>());
-    return error;
+    if (auto missingIDs = sourceDataSet.checkIDs(targetDataSet);
+        missingIDs.size() == 0)
+    {
+      mAlgorithm.setTrained(false);
+      DataSet             result(1);
+      auto                data = sourceDataSet.getData();
+      auto                tgt = targetDataSet.getData();
+      FluidDataSetSampler sampler(sourceDataSet, targetDataSet,
+                                  get<kBatchSize>(), get<kVal>(), true);
+      algorithm::SGD      sgd;
+      double error = sgd.train(mAlgorithm, data, tgt, sampler, get<kIter>(),
+                               get<kRate>(), get<kMomentum>());
+      return error;
+    }
+    else
+    {
+      std::ostringstream oss;
+      oss << "Can't train because these IDs are missing from the output data: ";
+      std::copy(missingIDs.begin(), missingIDs.end() - 1,
+                std::ostream_iterator<std::string>(oss, ","));
+      oss << missingIDs.back();
+      return {Result::Status::kError, oss.str()};
+    }
   }
 
   MessageResult<void> predict(InputDataSetClientRef srcClient,

include/flucoma/data/FluidDataSet.hpp

@@ -13,6 +13,8 @@ namespace fluid {
 template <typename idType, typename dataType, index N>
 class FluidDataSet
 {
+  template<typename, typename, index> 
+  friend class FluidDataSet; 
 
 public:
   explicit FluidDataSet() = default;
@@ -64,7 +66,7 @@ class FluidDataSet
     }
   }
 
-  bool add(idType const& id, FluidTensorView<dataType, N> point)
+  bool add(idType const& id, FluidTensorView<const dataType, N> point)
   {
     assert(sameExtents(mDim, point.descriptor()));
     index pos = mData.rows();
@@ -196,6 +198,40 @@ class FluidDataSet
     return result.str();
   }
 
+  template<typename T, index M> 
+  auto indexMap(FluidDataSet<idType, T, M> const& x) const
+    -> std::pair<std::vector<index>,std::vector<index>>
+  {
+    using std::pair, std::vector, std::begin, std::end; 
+
+    pair<vector<index>, vector<index>> result; 
+    result.first.reserve(asUnsigned(x.size()));
+    result.second.reserve(asUnsigned(x.size()));
+
+    auto firstID = begin(x.getIds());
+    auto lastID = end(x.getIds());
+
+    std::transform(firstID, lastID, std::back_inserter(result.first),
+                   [this](auto const& id) { return mIndex.at(id); });
+    std::transform(firstID, lastID, std::back_inserter(result.second),
+                   [&x](auto const& id) { return x.mIndex.at(id); });
+
+    return result;
+  }
+
+  template <class U, index M>
+  std::vector<idType> checkIDs(FluidDataSet<idType, U, M> const& other) const
+  {
+    std::vector<idType> result;
+
+    std::for_each(mIndex.begin(), mIndex.end(), [&result, &other](auto& item) {
+      if (other.mIndex.find(item.first) == other.mIndex.end())
+        result.push_back(item.first);
+    });
+
+    return result;
+  }
+
 private:
   void initFromData()
   {

include/flucoma/data/FluidDataSetSampler.hpp

@@ -0,0 +1,55 @@
+/*
+Part of the Fluid Corpus Manipulation Project (http://www.flucoma.org/)
+Copyright University of Huddersfield.
+Licensed under the BSD-3 License.
+See license.md file in the project root for full license information.
+This project has received funding from the European Research Council (ERC)
+under the European Union’s Horizon 2020 research and innovation programme
+(grant agreement No 725899).
+*/
+
+/*
+Pairing two unordered DataSets for supervised use
+*/
+#pragma once
+
+#include "FluidDataSet.hpp"
+#include "FluidIndex.hpp"
+#include "detail/DataSampler.hpp"
+#include <optional>
+#include <random>
+
+
+namespace fluid {
+
+class FluidDataSetSampler : public detail::DataSampler<FluidDataSetSampler>
+{
+  friend detail::DataSampler<FluidDataSetSampler>;
+
+  std::pair<std::vector<index>, std::vector<index>> mIdxMaps;
+
+  template <class InputIter>
+  FluidTensorView<index, 2> map(InputIter start, InputIter end,
+                                FluidTensorView<index, 2> dst)
+  {
+    using std::begin;
+    auto inputSamples = dst.col(0);
+    auto outputSamples = dst.col(1);
+
+    transform(start, end, begin(inputSamples),
+              [&idx = mIdxMaps](index i) { return idx.first[asUnsigned(i)]; });
+    transform(start, end, begin(outputSamples),
+              [&idx = mIdxMaps](index i) { return idx.second[asUnsigned(i)]; });
+    return dst;
+  }
+
+public:
+  template <typename DataSetA, typename DataSetB>
+  FluidDataSetSampler(DataSetA const& in, DataSetB const& out, index batchSize,
+                      double validationFraction, bool shuffle = true)
+      : detail::DataSampler<FluidDataSetSampler>(in.size(), batchSize,
+                                                 validationFraction, shuffle),
+        mIdxMaps{in.indexMap(out)}
+  {}
+};
+} // namespace fluid