classdef radialBasisNetwork2 %% this one uses the estimated variances for each cluster instead of b
    
    properties
        labData %% the experimental output data
        modelData %% the simulated output data
        error %% the difference between the experimental output data and the simulated output data
        u  % the input data, stores all the input data feed to it
        N  %The ammount of inputdata points stored in the object
        Nd  %backwards prediction horizon
        Nn  %number of nodes/clusters 
        mu  % the centers of the clusters
        Sigma % the inverse variance matrix(diagonal elements of covariance matrix), for each cluster
        w  % the weights in neural network
    end
    methods 
        function obj = radialBasisNetwork2(lData,mData,input,ND,NN)
            obj.labData=lData;
            obj.modelData=mData;
            obj.u=input;
            obj.Nd=ND;
            obj.Nn=NN;
            obj=uppdateParameters(obj);
        end
        
        function obj=uppdateParameters(obj)
            [indx,obj.mu] = kmeans(obj.u((end-obj.Nd+1):end,:),obj.Nn);
            obj.N = length(obj.u(:,1));
            obj.error = obj.labData - obj.modelData;
            obj=updateSigma(obj,indx);
            
            for k = (obj.N-obj.Nd+1):obj.N
                rSum = 0;
                for i = 1:obj.Nn
                    rSum = rSum + exp(-(1/8)* ((obj.u(k,:)-obj.mu(i,:)))*obj.Sigma{i}*((obj.u(k,:) - obj.mu(i,:))')); %% b = 1/(2*sigma^2)
                end
                for i =1:obj.Nn
                    
                    r((k-obj.N+obj.Nd),i) = exp(-(1/8)*((obj.u(k,:)-obj.mu(i,:)))*obj.Sigma{i}*((obj.u(k,:)-obj.mu(i,:))'))/rSum; %% some normalized densities
                end
            end
            R = [ones(obj.Nd,1),r];
            R=r;
            if (det(R'*R)>1e-8)
                Y = obj.error((end-obj.Nd+1):end,:);
                obj.w = (R'*R)\R'*Y;
            else
                fprintf("the matrix R'*R is badly scaled\n")
            end
        end
        
        function bias = predictError(obj,u)
                rSum = 0;
                for i = 1:obj.Nn
                    rSum = rSum + exp(-(1/8)*((u-obj.mu(i,:)))*obj.Sigma{i}*((u - obj.mu(i,:))')); %% b = 1/(2*sigma^2)
                end
                if rSum == 0
                    rSum =1;
                end
                for i =1:obj.Nn
                    r(i) = exp(-(1/8)*((u-obj.mu(i,:)))*obj.Sigma{i}*((u-obj.mu(i,:))'))/rSum; %% some normalized densities
                end
            R=r;
            bias = R*obj.w;
            if(isnan(bias(1)))
            end
        end
        
        function obj=uppdateWithNewData(obj,u,ylab,ymodel)
            obj=appendInput(obj,u);
            obj=appendLabdata(obj,ylab);
            obj=appendmodelData(obj,ymodel);
            obj=uppdateParameters(obj);
        end
        %% append functions
        function obj=appendInput(obj,newInp)
            obj.u = [obj.u;newInp];
        end
        function obj=appendLabdata(obj,newLabdata)
            obj.labData = [obj.labData;newLabdata];
        end
        function obj=appendmodelData(obj,newmodelData)
            obj.modelData = [obj.modelData;newmodelData];
        end
        function obj=updateSigma(obj,indx)
            uSplitt = cell(1,obj.Nn);
            for i = 1:obj.Nd
                uSplitt{indx(i)} = [uSplitt{indx(i)};obj.u((obj.N-obj.Nd+i),:)];
            end
            
            obj.Sigma = cell(obj.Nn,1);
            nk=zeros(obj.Nn,1);
            for i = 1:obj.Nn
                nk(i)=length(uSplitt{i}(:,1));
                if(nk(i)==1)
                    obj.Sigma{i}=eye(4);
                else
                    temp=0;
                    for j = 1:nk(i)
                        temp = temp + (uSplitt{i}(j,:) - obj.mu(i,:))'*(uSplitt{i}(j,:) - obj.mu(i,:));
                    end
                    temp2 = 1/(nk(i)-1)*temp;
                    temp3=diag(temp2);
                    for k = 1:4
                        if temp3(k)==0
                            temp3(k)=1;
                        end
                    end
                    obj.Sigma{i} = diag(1./(temp3));
                end      
            end
        end
        
    end
    
end