% This function gives the testing RMSE curve of EIR-ELM.

% Input:

%     N: Number of training samples

%     n: Number of input attributes

%     X: Input of raining data set, a N*n matrix

%     T: Target of training data set, a N*1 vector

%     tX: Input of testing data set

%     tT: Target of testing data set

%     maxnum: maximum number of network hidden nodes - 1

%     fun: Type of activation function: 0 for sigmoidal additive function and 1 for Gaussian radial basis function

%     init: Initial number of hidden nodes

%     k and C: Parameters of EIR-ELM

% Output:

%     E: A vector which contains the testing RMSE of each time when adding a chunk of hidden nodes

%     gap: A vector which contains the hidden layer size of each time when adding a chunk of hidden nodes

function [gap,E] = eir_elm(N,n,X,T,tX,tT,maxnum,fun,init,k,C)

W = rand(n,init)*2-1;

switch fun

    case 0

        b = rand(1,init)*2-1;

        H = activeadd(X,W,b);

    case 1

        b = rand(1,init)*0.5;

        H = activerbf(X,W,b);

    otherwise

end

D = inv(C*eye(init) + H'*H)*H';

beta = D*T;

switch fun

    case 0

        O = activeadd(tX,W,b)*beta;

    case 1

        O = activerbf(tX,W,b)*beta;

    otherwise

end

error = sqrt(mse(tT - O));

E = error;

num = init;

gap = num;

 

while num <= maxnum

    wi = rand(n,k)*2-1;

    switch fun

        case 0

            bi = rand(1,k)*2-1;

            dH = activeadd(X,wi,bi);

        case 1

            bi = rand(1,k)*0.5;

            dH = activerbf(X,wi,bi);

        otherwise

    end

    low = inf;

    pos = 0;

    for i = 1:k

        h = dH(:,i);

        temp = h'*(eye(N)-H*D);

        M = temp/(temp*h);

        L = D*(eye(N)-h*M);

        tD = [L;M];

        tbeta = tD*T;

        e = [H h]*tbeta-T;

        del = e'*e + C*(tbeta'*tbeta);

        if del < low

            low = del;

            pos = i;

            ttD = tD;

            beta = tbeta;

        end

    end

    W = [W wi(:,pos)];

    b = [b bi(:,pos)];

    H = [H dH(:,pos)];

    D = ttD;

    switch fun

        case 0

            O = activeadd(tX,W,b)*beta;

        case 1

            O = activerbf(tX,W,b)*beta;

        otherwise

    end

    error = sqrt(mse(tT - O));

    E = [E error];

    num = num+1;

    gap = [gap num];

end

 

function value = activeadd(x,wi,bi)

index = ones(1,size(x,1));

x = x*wi + bi(index,:);

value = 1./(1+exp(-x));

 

function value = activerbf(x,wi,bi)

n = size(x,1);

m = size(wi,2);

value = zeros(n,m);

for i = 1:m

    w = wi(:,i)';

    index = ones(n,1);

    wm = w(index,:);

   value(:,i) = -bi(i)*sum((x-wm).^2,2);

end

value = exp(value);