classdef kurtosis
   properties
      data % data X on which all computation takes place
           % The expected shape of the data is (n,k)
           % where n denotes the number of data points
           % and k denotes the dimensions of each data-point
      data_cov % covariance matrix of the data
   end
   methods
      %% Constructor to set input parameters
      function obj = kurtosis(val)
        if nargin == 1
            obj.data = val;
            obj.data_cov = cov(val);
        end
      end
      %% Function to evaluate objective function value
      function r = get_value(obj, u)
         % Input : obj - object
         %         u - (k) dimensional vector
         v = obj.data * u;
         r = mean(v.^4) - 3*(mean(v.^2)^2);
      end
      %% Function to evaluate gradient
      function r = get_gradient(obj, u)
        % Input : obj - object
        %         u - (k) dimensional vector
        % v = obj.data * u;
        % r = 4*mean(diag(v.^3)*obj.data,1) ...
        %     - 12*mean(diag(v)*obj.data,1)*mean(v.^2);
        % r = r';
        v = obj.data * u;
        v_3 = v.^3;
        r = 4*mean(v_3.*obj.data,1) ...
            - 12*mean(v.*obj.data,1)*mean(v.^2);
        r = r';
      end
      %% Function to evaluate hessian
      function r = get_hessian(obj, u)
        % Input : obj - object
        %         u - (k) dimensional vector
        n = size(obj.data,1);
        v = obj.data * u;
        % D = diag(v);
        X1 = v .* obj.data;
        L1 = 12*(X1')*X1/n;
        L2 = 12*mean(v.^2)*obj.data_cov;
        L3 = 24*mean(X1,1)'*mean(X1,1);
%         r = 12*(X1')*X1/n - 12*mean(v.^2)*obj.data_cov  ...
%             - 24*mean(X1,1)'*mean(X1,1);
        r = L1 - L2 - L3;
      end
      %% Function to estimate whitening matrix C
      function r = estimate_C(obj, m)
        % Input : obj - object
        %         m - number of samples over which to calculate C
        k = size(obj.data,2);
        C = zeros(k,k);
        for i = 1:k
            u = zeros(k,1); u(i) = 1;
            C = C + obj.get_hessian(u);
        end
        r = C/12;
      end
   end
end