function [TrainingTime, TestingTime, TrainingAccuracy, TestingAccuracy] = LCIELM(TrainingInput, TrainingTarget, TestingInput, TestingTarget, k, lambda, MaxHiddenNeurons, epsilon, ActivationFunType)

[TrainingDataNum, InputDimension] = size(TrainingInput);
start_time_training = cputime;
f = zeros(TrainingDataNum,1);
W = [];
b = [];
TrainingAccuracy = sqrt(mse(TrainingTarget));
OutputWeights = [];
num = 0;
count = 0;
%Training phase
while num <= MaxHiddenNeurons & TrainingAccuracy > epsilon
    len = floor(k*exp(-lambda*count))+1; %calculate the number of newly added hidden nodes for this loop
    wi = rand(InputDimension,len)*2-1;
    switch ActivationFunType
        case 0
            bi = rand(1,len)*2-1;
            H = activeadd(TrainingInput,wi,bi);
        case 1
            bi = rand(1,len)*0.5;
            H = activerbf(TrainingInput,wi,bi);
        otherwise
    end
    beta = pinv([f H])*TrainingTarget; %calculate the output weights for newly added hidden nodes
    f = f*beta(1) + H*beta(2:len+1); %update the output function
    OutputWeights = [OutputWeights*beta(1);beta(2:len+1)]; %update the output weights
    TrainingAccuracy = sqrt(mse(TrainingTarget-f));
    W = [W wi];
    b = [b bi];
    count = count+1;
    num = num+len;
end
end_time_training = cputime;
TrainingTime = end_time_training - start_time_training;
clear wi bi H beta f;
%Testing phase
start_time_testing = cputime;
switch ActivationFunType
    case 0
        O = activeadd(TestingInput,W,b)*OutputWeights;
    case 1
        O = activerbf(TestingInput,W,b)*OutputWeights;
    otherwise
end
end_time_testing = cputime;
TestingTime = end_time_testing - start_time_testing;
TestingAccuracy = sqrt(mse(TestingTarget - O));
%calculate the feature matrix for sigmoidal additive networks
function value = activeadd(x,wi,bi)
index = ones(1,size(x,1));
x = x*wi + bi(index,:);
value = 1./(1+exp(-x));
%calculate the feature matrix for Gaussian RBF networks
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);