function [out1,out2] = getweights(maxiter, numofhidden, input, y, centers, sigma)
% [out1,out2] = getweights(1000000, 6, input', y, centers, sigma)
% maxiter = maximal number of iterations
% numofhidden = number of hidden neurons in the hidden layer
% input = column vector of input patterns
% y = the true output I am aiming for
% centers = centers computed by C-means clustering
% sigma = used to control smoothness of the interp fcn, also computed before

N = 1;              % dimension of true_out neurons
M = 1;              % dimension of output neurons
num = 100;          % number of training true_out
k = numofhidden;    % number of hidden neurons
E = zeros(maxiter/1000,1); % energy
true_out = y;       % y - matrix(num,M)
train_out = zeros(num,M);  % training output
w = zeros(k,M);     % weight matrix
r = zeros(k,1);     % euclidean norm of input - centers
deltaw = zeros(k,M);% will be computed in update rules
deltac = zeros(k,M);% will be computed in update rules
e = 0;

% NORMALIZE THE INPUT DATA
input(:,1)=input(:,1)-mean(input(:,1));
input(:,1)=input(:,1)/std(input(:,1));

% INITIALIZE THE PARAMETERS
eta1 = 0.02; % learning rate
eta2 = 0.02; % learning rate
w = 0.5*rand(k,M) - 0.25; % randomly initialize the weight vector
phi = zeros(k,M); % Gaussian RBF fcn
multip = 1; % serves as a counter for plotting the energy

for iter=1:maxiter
    mu=ceil(rand*num); % randomly pick a pattern
    for i=1:M
        for j=1:k
            r(j) = norm(input(mu,1)-centers(j,1)); % r - matrix (6,1)
        end
        phi(1:k) = exp(-(r).^2/(2*sigma^2));       % phi - matrix (6,1)
        train_out(mu,i) = w' * phi;                % train_out - matrix(1,1)
    end
    e = true_out(mu,i) - train_out(mu,i);          % compute the error between the original output and network output
    for i=1:M
        for j=1:k
            deltaw(j,i) = 2*eta1*e*phi(j,1);       % deltaw - matrix(6,1)
            deltac(j,i) = -2*eta2*w(j,i)*e*phi(j,1)*((input(mu,i)-centers(j,i))/(norm(input(mu,i)-centers(j,i))));
        end
    end
    w = w + deltaw;                                % update the weights
    centers = centers + deltac;                    % update the centers
    if iter == multip*1000                         % save energy only every 1000th iteration
        E(multip,1) = 1/2 * (sum((true_out-train_out).^2));
        multip = multip + 1;
    end
end
out1 = w;
out2 = centers;
xx=1:(maxiter/1000);
xx = xx';
plot(xx,E)
end