classdef symmetric_chf
   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 = symmetric_chf(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  = 0.5*(log(mean(cos(v))^2 + mean(sin(v))^2) + ...
                   u'*obj.data_cov*u);
      end
      %% Function to evaluate gradient
      function r = get_gradient(obj, u)
        % Input : obj - object
        %         u - (k) dimensional vector
        
        v = obj.data*u;
        denom = mean(cos(v))^2 + mean(sin(v))^2;
        numer = (mean(cos(v))*mean(-obj.data.*sin(v), 1) + ...
                 mean(sin(v))*mean(obj.data.*cos(v), 1))';
        r = numer/denom + obj.data_cov*u;
%         keyboard;
      end
      %% Function to evaluate hessian
      function r = get_hessian(obj, u)
        % Input : obj - object
        %         u - (k) dimensional vector
        % keyboard;
        v     = obj.data*u;
        cos_v = cos(v);
        sin_v = sin(v);
        denom = mean(cos_v)^2 + mean(sin_v)^2;
        numer = (mean(cos(v))*mean(-obj.data.*sin(v), 1) + ...
                 mean(sin(v))*mean(obj.data.*cos(v), 1))';
        denom_grad = 2*numer';
        temp1 = (obj.data'*(-cos_v.*obj.data))/size(obj.data,1);
        temp2 = (obj.data'*(-sin_v.*obj.data))/size(obj.data,1);
        temp3 = mean(obj.data.*sin_v, 1)';
        temp4 = mean(obj.data.*cos_v, 1)';
        numer_grad = mean(cos_v)*temp1 + mean(sin_v)*temp2 + ...
                     temp3*temp3' + temp4*temp4';

        r = (denom*numer_grad - numer*denom_grad)/denom^2 + obj.data_cov;
      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
        if 0
            k = size(obj.data,2);
            C = zeros(k,k);
            for i = 1:m
                u = randn(k,1); u = u/norm(u);
                C = C + obj.get_hessian(u);
            end
            r = (C/m);
        else
            k = size(obj.data,2);
            C = zeros(k,k);
            for i = 1:m
                % u = randn(k,1); u = u/k;
                u = randn(k,1); u = u/norm(u);
                % C = C + obj.get_hessian(u);
                % u = zeros(k,1); u(i) = 1;
                C = C + obj.get_hessian(u);
            end
            r = (C/m);
            % r = (C/m);
            % r = r/norm(r);
            % [V,E]=eig(r);e=diag(E);
            % e(abs(e)<0.1*max(abs(e)))=0;
            % r=V*diag(e)*V';
        end
         %r = (C/m);
%          r = r/max(r, [], 'all');
%         for i = 1:k
%             u = zeros(k,1); u(i) = 1;
%             C = C + obj.get_hessian(u);
%         end
%         r = C/12;
      end
   end
end