#include <math.h>
#include <stdio.h>
#include <iostream>

#include "logreg.hpp"

Logreg::Logreg() :
  w(1, 0.3, 0),
  points()
{

}

void Logreg::train() {
  Vec3D w_new(1, 0, 0);
  Vec3D w_prev(1, 0.3, 0);
  double step_size = 1.0;
  double E = loss(w);
  Vec3D g = gradient(w);
  double dist = 100;
  while (dist > 0.001) {
    w_new = w_prev - g * step_size;
    double E_new = loss(w_new);
    if (E_new != E_new) {
      w_new = w_prev;
      std::cout << "unlimited" << std::endl;
      break;
    }
    Vec3D g_new = gradient(w_new);

    if (E_new <= E) {
      dist = (w_new - w_prev).norm();
      w_prev = w_new;
      E = E_new;
      g = g_new;
      step_size *= 2;
    } else {
      step_size /= 2;
    }
  }

  w = w_new;
  std::cout << "w = ["<< w.x << ", " << w.y << ", " << w.z << "]" << std::endl;
}

void Logreg::add_measurement(double phase, bool positive) {
  points.push_back(Measurement(Vec3D(cos(phase), sin(phase), 1), positive ? 1 : -1));
}

double Logreg::loss(Vec3D at) const {
  double result = 0;
  for (int i = points.size(); i >= 0; --i) {
    result += log(1 + exp(-(at*points[i].x)*points[i].y));
  }
  return result / points.size();
}

Logreg::Vec3D Logreg::gradient(Vec3D at) const {
  Vec3D result(0, 0, 0);
  for (int i = points.size(); i >= 0; --i) {
    result += points[i].x * (points[i].y/(1 + exp((at*points[i].x)*points[i].y)));
  }
  result.x /= -int(points.size());
  result.y /= -int(points.size());
  result.z /= -int(points.size());
  return result;
}
