classdef Volterra_model < handle
    %Class structure with different methods and properties
    %The idea is that all data that is processed stays as information in
    %the class. 
    %To get the properties the code is . based:  Volterra.Order = 2
    %An overview of the properties, see under %Properties below
    %
    %All methods can be run by running: 
    % example: Volterra_model.FindLambda(Input, Output, ...)
    % For explenation of the methods, see the documentation on the
    % different methods
    %
    
    
    
    properties
        %%
        %Declaring properties of the class, do it is easily accessible
        Memory=NaN;         %Memory
        Order = NaN;        %Order
                
        %Structure of the regularization for training data
        XX=NaN;  

        %Kernels
        HH=NaN;             %Full kernel structure(Flat, Nx1)
        H1=NaN;             %1st-order kernel(M+1,1)
        HH2=NaN;            %2nd-order kernel(M+1,M+1)
        HH3=NaN;            %3rd-order kernel(M+1,M+1,M+1)
        
        %Structure with different predictions of validation data
        YPred={};        %YPred dumps the last predicted value at the end of YPred cell structure. 
        
        
        
        %Vectors for plotting the single and double L-curves
        
        %For single L-curve, 
        Fit_err_norm = NaN;                 %||r*H-F1|||     (1xlength(LambdaInterval))
        LH_norm = NaN;                      %||L*H||            (1xlength(LambdaInterval))
        LambdaInterval = NaN;              %LambdaInterval vector that fits to the different points, Lambda^2   (1xlength(LambdaInterval))
       
               
       %For double L-curve 
        Fit_error_normFull = NaN;                %||r*H-F1|||   ((length(LambdaDecayInt21)xlength(LambdaDecayInt22))          
        LambdaDecayInt21 = NaN;                   %Lambda Decay Interval i direction 1
        LambdaDecayInt22 = NaN;                   %Lambda Decay Interval i direction 2
        LambdaNormFullDecay1 = NaN;              %||L1*H||, ((length(LambdaDecayInt21)xlength(LambdaDecayInt22))     
        LambdaNormFullDecay2 = NaN;              %||L2*H||, ((length(LambdaDecayInt21)xlength(LambdaDecayInt22))     
        
        
        
    end
    %%
    methods  
        function FindLambda(obj,Inp,Out,Memory,Order,LambdaIntervalReg1,LambdaDecayInterval1,LambdaDecayInterval2,Sym,RegType)
            %FindLambda(Input,Output,Memory,Order,LambdaIntReg,LambdeDecayInt1,LambdeDecayInt2,Sym,type)
            %Inputs:
            %
            %Input = Training input(motion)
            %Output = Training output(force)
            %Order = model order(1,2,3) ps. max 2 for decay type
            %LambdaIntReg = The Interval Lambda is looped over, for Tikhonov regularization(1st and 0th)
            %LambdeDecayInt1 = The Interval Lambda decay is looped over, axis 1
            %LambdeDecayInt2 = The Interval Lambda decay is looped over, axis 2
            %Sym = true/false. If one uses symmetry reduction(true false). 
            %RegType = 'Reg1', 'Reg0' or 'decay' for 1st order Tikhonov, 0th order Tikhonov or decay type regularization    
            
            
            %Outputs:
            % 
            %For Single L-curves: In other words Tikhonov regularization
            %Fit_err_norm            %||r*H-F1|||     (1xlength(LambdaInterval))
            %LH_norm                 %||L*H||            (1xlength(LambdaInterval))
            %LambdaInterval        %%LambdaInterval vector that fits to the different points, Lambda^2   (1xlength(LambdaInterval)) = LambdaIntervalReg1
            %
            %For double L-curves(RegType = 'decay' and order =2):
            %FittError_normFull         %||r*H-F1|||   ((length(LambdaDecayInt21)xlength(LambdaDecayInt22))   
            %LambdaNormFullDecay1       %Lambda Decay Interval i direction 1
            %LambdaNormFullDecay2       %Lambda Decay Interval i direction 2
            %LambdaDecayInt21            %||L1*H||, ((length(LambdaDecayInt21)xlength(LambdaDecayInt22)) 
            %LambdaDecayInt22            %||L2*H||, ((length(LambdaDecayInt21)xlength(LambdaDecayInt22))        
            [obj.XX]=InpVolterraSeries(Inp,Memory,Order);            
            
            [obj.Fit_err_norm,obj.LH_norm,obj.LambdaInterval,obj.Fit_error_normFull,obj.LambdaNormFullDecay1,obj.LambdaNormFullDecay2,obj.LambdaDecayInt21,obj.LambdaDecayInt22...
                ]=FindLambda(obj.XX,Out,Memory,Order,LambdaIntervalReg1,LambdaDecayInterval1,LambdaDecayInterval2,Sym,RegType);      
        end       
   
 
        function TrainVolterra(obj,Inp,Out,Memory,Order,Sym,LambdaDecay1,LambdaDecay2,LambdaReg1,LambdaReg0)
            %TrainVolterra(Input,Output,Memory,Order,Symmetry,LambdaDecay1,LambdaDecay2,LambdaReg1,LambdaReg0)
            %Inputs:
            %Input = Training input(motion)
            %Output = Training output(force)
            %Order = model order(1,2,3) ps. max 2 for decay type
            %Sym = true/false. If one uses symmetry reduction(true false). 
            %LambdaDecay1 = Lambda decay 1st value
            %LambdaDecay2 = Lambda decay 2nd value
            %LambdaReg1 = Lambda Tikhonov 1st
            %LambdaReg0 = Lambda Tikhonov 0th
            
            %Outputs:
            %HH = Full structure of output kernel matrix, se Eq XXX. 
            %H1 = First kernel(1x(Memory + 1))
            %HH2 = Matrix of second kernel((Memory++1)x(Memory+1))
            %HH3 = 3dim Matrix of third kernel((Memory+1)x(Memory+1)x(Memory+1))
            
            [obj.XX]=InpVolterraSeries(Inp,Memory,Order);           
            [obj.HH, obj.H1 ,obj.HH2, obj.HH3] = TrainVolterra(obj.XX,Out,Memory,Order,Sym,LambdaDecay1,LambdaDecay2,LambdaReg1,LambdaReg0);
            
            obj.Memory = Memory;
            obj.Order = Order;
        end
        
        
        function ValidateVolterra(obj,InpValid)
            %Method that makes time-realization of independent motion series
            %ValidateVolterra(InputValid)
            %InputValid = Input validation series
            %Dumps results in YPred{end}, Note that training a new model does not
            %delete the old, but places it at the end. 
            
            %PS. The model needs to be trained first!
            
            %Input: 
            %InputValid = Validation input series(1xN)
            
            %Output: 
            %YPredL{end} = Predicted output from the validation series(1xN)
           [XXL]=InpVolterraSeries(InpValid,obj.Memory,obj.Order);          
            
           if isnan(obj.HH)
               error('Model needs to be trained first, use the .TrainVolterra() - function')
           end
           YPredL= XXL*obj.HH;
           obj.YPred{end+1}=YPredL;
        end
    end
    
end