X = [(o(:,1))',(x(:,1))';
    (o(:,2))',(x(:,2))'];
y = [-1.*ones(1,length(o(:,1))), ones(1,length(x(:,1)))]';
theta = rand(5,1); % [m1, m2, q1, p, q2]
m = theta(1:2);
Q = [theta(3), theta(4);
    theta(4), theta(5)];
error_log = zeros(1,1);

alpha = 0.1;
beta = 0.1;

error = 100;
i = 1;
while error > 185
    derror_ddelta = zeros(length(X),1);
    error_per_squared = zeros(length(X),1);
    ddelta_dm1 = zeros(length(X),1);
    ddelta_dm2 = zeros(length(X),1);
    ddelta_dq1 = zeros(length(X),1);
    ddelta_dp = zeros(length(X),1);
    ddelta_dq2 = zeros(length(X),1);
    for j = 1:length(X)
        delta = (X(:,j) - m)' * Q * (X(:,j) - m);
        error_per = hyp_tan(beta, delta) - y(j);
        error_per_squared(j) = error_per ^ 2;
        r_dot = hyp_tan_dot(beta, delta);
        derror_ddelta(j) = 2 * error_per * r_dot;

        ddelta_dm1(j) = -2 * (Q(1,1) * (X(1,j) - m(1)) + Q(1,2) * (X(2,j) - m(2)));
        ddelta_dm2(j) = -2 * (Q(2,2) * (X(2,j) - m(2)) + Q(1,2) * (X(1,j) - m(1)));
        ddelta_dq1(j) = (X(1,j) - m(1))^2;
        ddelta_dp(j) = 2 * (X(1,j) - m(1)) * (X(2) - m(2));
        ddelta_dq2(j) = (X(2,j) - m(2))^2;
    end
    error = 0.5 * sum(error_per_squared);
    theta = theta - alpha .* sum(derror_ddelta) .* [sum(ddelta_dm1);
                                                  sum(ddelta_dm2);
                                                  sum(ddelta_dq1);
                                                  sum(ddelta_dp);
                                                  sum(ddelta_dq2)];
    error_log(i) = error;
    i = i+1;
end

m = theta(1:2);
Q = [theta(3), theta(4);
    theta(4), theta(5)];

% ans_x2 = zeros(1, length(X(1,:)));
% for k = 1:length(X(1,:))
%     syms x2
%     x1 = X(1,k);
%     eqn = Q(1,1)*(x1 - m(1))^2 + 2*Q(1,2)*(x1 - m(1))*(x2 - m(2)) + Q(2,2)*(x2 - m(2))^2;
%     ans_x2(1,1) = solve(eqn, x2);
% end

close all
figure()
% scatter(o(:,1),o(:,2))
% hold
% scatter(x(:,1),x(:,2),'x')
% hold
x_plot = X(1,:);
% y_plot = ans_x2;
% line(x_plot, y_plot)
% syms x1
% syms y
eqn = Q(1,1)*(x_plot - m(1))^2 + 2*Q(1,2)*(x_plot - m(1))*(y_plot - m(2)) + Q(2,2)*(y_plot - m(2))^2;
% hold
f = @(x_plot,y_plot) eqn;
fimplicit(f,[-10 10 -10 10])


figure()
plot(linspace(0,length(error_log),length(error_log)), error_log)