classdef cgf
   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
      imag_data %whether data is imaginary or not.
   end
   methods
      %% Constructor to set input parameters
      function obj = cgf(val)
        if nargin == 1
            obj.data = val;
            obj.data_cov = cov(val);
            obj.imag_data = 1*(~isreal(val));
        end
      end
      %% Function to evaluate objective function value
      function r = get_value(obj, u)
         % Input : obj - object
         % u - (k) dimensional vector
         u=u*3;
         v  = obj.data*u;
         pv = exp(v);
         t  = mean(pv);
         r  = log(t) - (1-2*obj.imag_data)*u'*obj.data_cov*u/2;
        %r  = log(t)+log(mean(exp(-v))) - 2*(1-2*obj.imag_data)*u'*obj.data_cov*u/2;
      end
      %% Function to evaluate gradient
      function r = get_gradient(obj, u)
        % Input : obj - object
        % u - (k) dimensional vector
        % v     = obj.data*u;
        % exp_v = exp(v-max(v));
        % pv    = exp_v/sum(exp_v);
        % t     = sum(obj.data'*diag(pv),2);
        % r     = t - obj.data_cov*u;
        % if obj.imag_data==1
        %     r=-imag(r);
        % end
        %u=u;
        v = obj.data*u;v=v-max(v);
        
        r = transpose(mean(obj.data.*exp(v), 1)/mean(exp(v))) ...
            - obj.data_cov*u;
        
        % if(isnan(r))
        %     [~,k] = size(obj.data);
        %     r = 0.01*randn(k,1);
        % end
      end
      %% Function to evaluate hessian
      function r = get_hessian(obj, u)
        % Input : obj - object
        %         u - (k) dimensional vector
        % v     = obj.data*u;
        % exp_v = exp(v-max(v));
        % pv    = exp_v/sum(exp_v);
        % t     = diag(pv)*obj.data;
        % C1    = t'*obj.data;
        % tmean = sum(t,1)';
        % C2    = tmean*tmean';
        % r     = C1 - C2 - obj.data_cov;
        v    = obj.data*u;
        %C1   = (transpose(obj.data)*diag(exp(v))*obj.data)/(size(obj.data,1)*mean(exp(v)));
        C1   = (transpose(obj.data)*(exp(v).*obj.data))/(size(obj.data,1)*mean(exp(v)));
        
        temp = transpose(mean(obj.data.*exp(v), 1)/mean(exp(v)));
        C2   = temp*transpose(temp);
        r    = C1 - C2 - obj.data_cov;
%         fprintf("============\n");
%         display(C1);
%         display(tmean);
%         display(C2);
%         display(obj.data_cov);
%         display(r);
%         fprintf("============\n");
        if obj.imag_data==1
            r=-real(r);
        end
      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:m
            % u = randn(k,1); u = u/k;
            u = randn(k,1); u = u/norm(u);
            % u = randn(k,1); u = u/sqrt(k);
            % u = zeros(k,1); u(i) = 1;
            C = C + obj.get_hessian(u);
        end
         % r = (C/k);
         r = C/max(abs(C(:)));
         % 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';

%          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