GP_Integrate.m

function [T,Y]=GP_Integrate(betas,matrix,addinp,norms,phis,start,stop,y0,h,used_inputs)
%Integrates a set of GP models of derivatives for a set of ODEs generated by emulator
%using a 4th order Runge-Kutta method

%Function Inputs:
%betas is a cell array in which each cell contains a specific row of the betas matrix, 
%or the mean of the the betas matrix for each model being integrated

%matrix is a cell array containing the interaction matrix of each model

%addinp is a matrix of of the values of all the other inputs to the model(s) (including
%any forcing functions) over the time period we integrate over. The length of addinp 
%should be equal to the number of points in the final time series (end-start)/h
%All values in addinp need to be normalized  with respect to the min and max values
%of their respective values in the training dataset

%h is the step size with respect to time

%norms is a matrix of the min and max values of all the inputs being 
%integrated (in the same order as y0). min values are in the top row, max values in the bottom.

%phis are the spline coefficients for the basis functions (cell array)

%Start is the time at which integration begins. Stop is the time to
%end integration.

%y0 is a vector of the inital conditions for the models being integrated

%Used inputs is a cell array containing the information as to what inputs
%are used in what model. Each cell should contain a vector corresponding to a different model. 
%Inputs should be referred to as those being integrated first, followed by
%those contained in addinp (in the same order as they appear in y0 and
%addinp respectively)
%For example, if two models were being integrated, with 3 other inputs total
%and the 1st model used both models outputs as inputs and the 1st and 3rd additional
%inputs, while the 2nd model used its own output as an input and the 2nd
%and 3rd additional inputs, used_inputs would be equal to
%{[1,1,1,0,1],[0,1,0,1,0]}. 
%If the models created do not follow this ordering scheme for their inputs
%the inputs can be rearranged based upon an alternate 
%numbering scheme provided to used_inputs. E.g. if the inputs need to be
%reordered the the 1st input should have a '1' in its place in the
%used_inputs vector, the 2nd input should have a '2' and so on. Using the
%same example as before, if the 1st models inputs needed rearranged so that
%the 3rd additional input came first, followed by the two model outputs in
%the same order as they are in y0, and ends with the 1st additional input,
%then the 1st cell in used_inputs would have the form [2,3,4,0,1]

%Function Outputs:

%T a vector of the time steps the models are integrated at.

%Y is a matrix of the models that have been integrated, at the time steps
%contained in T.


    function f = prediction(inputs)
        f=[];
        for kk=1:length(inputs)
            f= [f;bss_eval(inputs{kk}, betas{kk}, phis, matrix{kk})];
        end
    end

    function rein = reorder(used,inputs)
        used(used==0) = [];
        rein = zeros(size(inputs));
        for m=1:length(inputs)
            rein(used(m)) = inputs(m);
        end
    end

    function norm = norma(v,minim,maxim)
        %A 0 to 1 normalization function. v is a vector of the data to be
        %normalized. minim and maxim are the minimum and maximum values to normalize over.
        %If v is the only input the function uses the minimum and
        %maximum values of v as the global minimum and maximum values to normalize over.
        norm = zeros(length(v),1);

        for loop=1:length(v)
            norm(loop) = (v(loop) - minim)/(maxim - minim);
            if norm(loop) > 1
                norm(loop) = 1;
            else
                if norm(loop) < 0
                    norm(loop) = 0;
                end
            end
        end

    end

T=start:h:stop;
y=y0;
Y=[y0,zeros(length(y0),length(T)-1)];
ind = 2;
for t=T(1:end-1)
    inputs1 = cell(length(y0),1);
    othinputs = cell(length(y0),1);
    for i=1:length(y) %Include values of the integrated functions in the inputs of the bss model
        for j=1:length(y)
            if used_inputs{i}(j) ~= 0
                inputs1{i} = [inputs1{i},norma(y(j),norms(1,j),norms(2,j))];
            end
        end
    end
    if isempty(addinp) == 0
        for ii=1:length(y0) %Determine the values of all other inputs at this time step
            for jj = length(y)+1:(size(addinp,2)+length(y))
                if used_inputs{ii}(jj) ~= 0
                    othinputs{ii} = [othinputs{ii},addinp(ind-1,jj-length(y0))];
                end
            end
        end
        for k=1:length(y0)
            inputs1{k} = [inputs1{k},othinputs{k}];
        end
    end
    for l=1:length(y0)
        if max(used_inputs{l}) > 1
            inputs1{l} = reorder(used_inputs{l},inputs1{l});
        end
    end

    dy1=prediction(inputs1)*h; %euler estimate of y(t+h)-y(t)
    %In between each estimate of the derivative we check to see if the
    %value of the outputs exceeds the range of the training data, and if
    %the derivative would produce a new prediction further outside the
    %accepted range. If both conditions are true, the derivative is set
    %equal to 0.
    for p=1:length(y)
        if y(p)>norms(2,p) && dy1(p)>0
            dy1(p)=0;
        else
            if y(p)<norms(1,p) && dy1(p)<0
                dy1(p)=0;
            end
        end
    end

    inputs2 = cell(length(y0),1);
    for i=1:length(y) %Include values of the integrated functions in the inputs of the bss model
        for j=1:length(y)
            if used_inputs{i}(j) ~= 0
                inputs2{i} = [inputs2{i},norma((y(j)+dy1(j)/2),norms(1,j),norms(2,j))];
            end
        end
    end
    for k=1:length(y)
        inputs2{k} = [inputs2{k},othinputs{k}];
    end
    for l=1:length(y0)
        if max(used_inputs{l}) > 1
            inputs2{l} = reorder(used_inputs{l},inputs2{l});
        end
    end
    dy2=prediction(inputs2)*h; %midpoint estimate
    for p=1:length(y)
        if (y(p)+dy1(p)/2)>norms(2,p) && dy2(p)>0
            dy2(p)=0;
        else
            if (y(p)+dy1(p)/2)<norms(1,p) && dy2(p)<0
                dy2(p)=0;
            end
        end
    end

    inputs3 = cell(length(y0),1);
    for i=1:length(y) %Include values of the integrated functions in the inputs of the bss model
        for j=1:length(y)
            if used_inputs{i}(j) ~= 0
                inputs3{i} = [inputs3{i},norma((y(j)+dy2(j)/2),norms(1,j),norms(2,j))];
            end
        end
    end
    for k=1:length(y)
        inputs3{k} = [inputs3{k},othinputs{k}];
    end
    for l=1:length(y0)
        if max(used_inputs{l}) > 1
            inputs3{l} = reorder(used_inputs{l},inputs3{l});
        end
    end
    dy3=prediction(inputs3)*h; %refined midpoint method
    for p=1:length(y)
        if (y(p)+dy2(p)/2)>norms(2,p) && dy3(p)>0
            dy3(p)=0;
        else
            if (y(p)+dy2(p)/2)<norms(1,p) && dy3(p)<0
                dy3(p)=0;
            end
        end
    end

    inputs4 = cell(length(y0),1);
    for i=1:length(y) %Include values of the integrated functions in the inputs of the bss model
        for j=1:length(y)
            if used_inputs{i}(j) ~= 0
                inputs4{i} = [inputs4{i},norma((y(j)+dy3(j)),norms(1,j),norms(2,j))];
            end
        end
    end
    for k=1:length(y)
        inputs4{k} = [inputs4{k},othinputs{k}];
    end
    for l=1:length(y0)
        if max(used_inputs{l}) > 1
            inputs4{l} = reorder(used_inputs{l},inputs4{l});
        end
    end
    dy4=prediction(inputs4)*h; %forward estimate
    for p=1:length(y)
        if (y(p)+dy3(p))>norms(2,p) && dy4(p)>0
            dy4(p)=0;
        else
            if (y(p)+dy3(p))<norms(1,p) && dy4(p)<0
                dy4(p)=0;
            end
        end
    end
    y=y+(dy1+2*dy2+2*dy3+dy4)/6; %averaging
    Y(:,ind)=y; %augment with a new column
    ind = ind+1;
end
end