std::vector<double> NeuralNet::update(const std::vector<double> &inputs)
{
    using namespace util;
    std::copy(inputs.begin(), inputs.end(), mInputActivation.begin());

    for (int hiddenIdx = 0; hiddenIdx < mNumHidden; ++hiddenIdx)
    {
        mHiddenActivation[hiddenIdx] =
            sigmoid(multiplyVect(mInputActivation, mWin->getColVector(hiddenIdx)));
    }

    for (int outputIdx = 0; outputIdx < mNumOutputs; ++outputIdx)
    {
        mOutputActivation[outputIdx] =
            sigmoid(multiplyVect(mHiddenActivation, mWout->getColVector(outputIdx)));
    }

    return mOutputActivation;
}

double NeuralNet::backprop(const std::vector<double> &target, double learningRate)
{
    std::vector<double> output_deltas(mNumOutputs, 0.0);
    std::vector<double> hidden_deltas(mNumHidden, 0.0);

    // Compute output deltas
    for (int k = 0; k < mNumOutputs; k++)
    {
        double error = target[k] - mOutputActivation[k];
        output_deltas[k] = util::d_sigmoid(mOutputActivation[k]) * error;
    }

    // Compute hidden deltas
    for (int j = 0; j < mNumHidden; j++)
    {
        double error = util::multiplyVect(output_deltas, mWout->getRowVector(j));
        hidden_deltas[j] = util::d_sigmoid(mHiddenActivation[j]) * error;
    }

    // Update output weights
    for (int j = 0; j < mNumHidden; j++)
    {
        for (int i = 0; i < mNumOutputs; i++)
        {
            double new_weight = mWout->get(j, i) +
                learningRate * (output_deltas[i] * mOutputActivation[j]);

            mWout->setElement(j, i, new_weight);
        }
    }

    // Update hidden weights
    for (int k = 0; k < mNumInputs; k++)
    {
        for (int j = 0; j < mNumHidden; j++)
        {
            double new_weight = mWin->get(k, j) +
                learningRate * (hidden_deltas[j] * mInputActivation[k]);

            mWin->setElement(k, j, new_weight);
        }
    }

    return util::error(target, mOutputActivation);
}