#ifndef EM_SQDIST_H
#define EM_SQDIST_H

#include <stdint.h>
#include <vector>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_vector.h>

#include "em_circular.h"

class EMSqdist : public EMCircular
{
  public:
    EMSqdist();
    EMSqdist(int cluster_count_);

    void train();
    void add_time(double phase);

    void print();

    double get_density_at(double phase) const;

  private:
    void expectation();
    double maximisation();

    class Cluster {
      public:
        Cluster();
        Cluster(double xx_, double yy_, double weight_);
        double xx;
        double yy;
        double weight;

        void print();

        double density_at(double phase) const;
    };

    class MinParams {
      public:
        MinParams(std::vector<Timestamp>* timestamps_, int k_);
        std::vector<Timestamp>* timestamps;
        int k;
    };

    std::vector<Cluster> clusters;

    static double my_f (const gsl_vector *v, void *params);
    static void my_df (const gsl_vector *v, void *params, gsl_vector *df);
    static void my_fdf(const gsl_vector *v, void *params, double *f, gsl_vector *df);
    static void tanh_ext(double& x, double& y);
    static void inv_tanh_ext(double& x, double& y);

    void U(double zr, double zi, double fr, double fi, double w, double& vr, double& vi);
    void mle(int c, double s, double& xhat, double& yhat);
};

#endif // EM_SQDIST_H
