#ifndef EUCLIDEAN_LOST_HPP
#define EUCLIDEAN_LOST_HPP

#include <dlib/dnn.h>

namespace dlib{

class loss_euclidean_dist_
{
public:

    typedef matrix<float> training_label_type;
    typedef matrix<float> output_label_type;

    template <
            typename SUB_TYPE,
            typename label_iterator
            >
    void to_label (
            const tensor& input_tensor,
            const SUB_TYPE& sub,
            label_iterator iter
            ) const
    {
        //static_assert(SUB_TYPE::num_layers-1 > 0, "SUB_TYPE::num_layers-1 > 0");

        std::cout<<__func__<<" : num of layers "<<sub.num_layers<<std::endl;
        DLIB_CASSERT(input_tensor.num_samples() != 0);

        //std::cout<<sub.num_layers<<std::endl;
        auto const &output_tensor = layer<sub.num_layers-1>(sub).get_output();

        DLIB_CASSERT(output_tensor.num_samples() != 0);
        //this for loop should(?) make sure the iter is same as the
        //input labels, how could I do that?
        for (long i = 0; i < output_tensor.num_samples(); ++i)
        {
            *iter++ = image_plane(output_tensor, i);
        }
    }

    template <
            typename const_label_iterator,
            typename SUBNET
            >
    double compute_loss_value_and_gradient (
            const tensor& input_tensor,
            const_label_iterator truth,
            SUBNET& sub
            ) const
    {
        DLIB_CASSERT(input_tensor.num_samples() != 0);

        auto const &output_tensor = sub.get_output();
        tensor& grad = sub.get_gradient_input();
        DLIB_CASSERT(input_tensor.num_samples() == output_tensor.num_samples());
        DLIB_CASSERT(grad.nc() == input_tensor.nc() && grad.nr() == input_tensor.nr());
        DLIB_CASSERT(input_tensor.nr() == output_tensor.nr() &&
                     input_tensor.nc() == output_tensor.nc() &&
                     input_tensor.k() == output_tensor.k());

        double loss = 0;
        output_label_type update_gradient(grad.nr(), grad.nc());
        update_gradient = 0;
        for(long i = 0; i != output_tensor.num_samples(); ++i){
            auto const predict_residual = image_plane(output_tensor, i);
            DLIB_CASSERT(predict_residual.nr() == truth->nr() &&
                         predict_residual.nc() == truth->nc());
            //I set the label(truth) as
            //[high_res_img - low_res_img] == [truth residual_img].
            auto const error_residual = *truth - predict_residual;
            auto const eu_dist = sum(pointwise_multiply(error_residual, error_residual));
            //after derivative, gradient of euclidean distance are
            //input - reference, since output_tensor is a mini batch
            //I sum all of them together
            update_gradient += error_residual;
            loss += eu_dist;
            ++truth;
        }
        //I take the average as the gradient value
        update_gradient /= output_tensor.num_samples();
        auto *grad_ptr = grad.host_write_only();
        std::copy(&update_gradient(0), &update_gradient(0) + update_gradient.size(),
                  grad_ptr);

        return loss / 2.0 / output_tensor.num_samples();
    }

    friend void serialize(const loss_euclidean_dist_& , std::ostream& out)
    {
        serialize("loss_euclidean_dist_", out);
    }

    friend void deserialize(loss_euclidean_dist_& , std::istream& in)
    {
        std::string version;
        deserialize(version, in);
        if (version != "loss_euclidean_dist_")
            throw serialization_error("Unexpected version found while deserializing dlib::loss_euclidean_dist_.");
    }

    friend std::ostream& operator<<(std::ostream& out, const loss_euclidean_dist_& )
    {
        out << "loss_euclidean_dist_";
        return out;
    }

    friend void to_xml(const loss_euclidean_dist_&, std::ostream& out)
    {
        out << "<loss_euclidean_dist_/>";
    }

};

template <typename SUBNET>
using loss_euclidean_dist = add_loss_layer<loss_euclidean_dist_, SUBNET>;//*/

}

#endif // EUCLIDEAN_LOST_HPP