function NetParams = rbfOLS(X,T,Goal,MaxNeurons)
% Modification of Matlab built-in newrb.
%  Simplifications:  Must use Gaussian transfer function
%                    Spread is a default value.
%                    The output will be the weights and centers.
%
%  X, T are dim x numpts

[r1,c1]=size(X);
[r2,c2]=size(T);
if c1~=c2
   error("Dimension mismatch in rbfOLS.  Data should be dim x numpts");
end
     
    [C,b1,W,b,tr] = OLS(X,T,Goal,1,MaxNeurons);
    spr=b1(1);
    W=W';
    b=b';
    NetParams.C=C;  %Centers for the RBF:  dim x numcenters
    NetParams.W=W;  %Weights:  out dim x in dim
    NetParams.b=b;  % bias:    out dim
    NetParams.spr=spr;  %Spread (actually, the reciprocal of spread)
    NetParams.tr=tr;   % Training error
    
end

%======================================================
function [w1,b1,w2,b2,tr] = OLS(p,t,eg,sp,mn)

  [r,q] = size(p);
  [s2,q] = size(t);
  b = sqrt(-log(.5))/sp;

  % RADIAL BASIS LAYER OUTPUTS
  P = radbas1(dist(p',p)*b);
  PP = sum(P.*P)';
  d = t';
  dd = sum(d.*d)';

  % CALCULATE "ERRORS" ASSOCIATED WITH VECTORS
  e = ((P' * d)' .^ 2) ./ (dd * PP');

  % PICK VECTOR WITH MOST "ERROR"
  pick = findLargeColumn(e);
  used = [];
  left = 1:q;
  W = P(:,pick);
  P(:,pick) = []; PP(pick,:) = [];
  e(:,pick) = [];
  used = [used left(pick)];
  left(pick) = [];

  % CALCULATE ACTUAL ERROR
  w1 = p(:,used)';
  a1 = radbas1(dist(w1,p)*b);
  [w2,b2] = solvelin2(a1,t);
  a2 = w2*a1 + b2*ones(1,q);
  MSE = mse(t,a2);

  tr(1) = mse(t,repmat(mean(t,2),1,q));  %Tracking error
  tr(2) = MSE;

  iterations = min(mn,q);
  for k = 2:iterations

    % CALCULATE "ERRORS" ASSOCIATED WITH VECTORS
    wj = W(:,k-1);
    a = wj' * P / (wj'*wj);
    P = P - wj * a;
    PP = sum(P.*P)';
    e = ((P' * d)' .^ 2) ./ (dd * PP');

    % PICK VECTOR WITH MOST "ERROR"
    pick = findLargeColumn(e);
    W = [W, P(:,pick)];
    P(:,pick) = []; PP(pick,:) = [];
    e(:,pick) = [];
    used = [used left(pick)];
    left(pick) = [];

    % CALCULATE ACTUAL ERROR
    w1 = p(:,used)';
    a1 = radbas1(dist(w1,p)*b);
    [w2,b2] = solvelin2(a1,t);
    a2 = w2*a1 + b2*ones(1,q);
    MSE = mse(t,a2);

    % PROGRESS
    tr(k+1) = MSE;

    % CHECK ERROR
    if (MSE < eg), break, end

  end  %End of k loop

  [S1,R] = size(w1);
  b1 = ones(S1,1)*b;

end  %End of function

%======================================================

function i = findLargeColumn(m)
  replace = find(isnan(m));
  m(replace) = zeros(size(replace));
  m = sum(m .^ 2,1);
  i = find(m == max(m));
  i = i(1);
end

%======================================================

function [w,b] = solvelin2(p,t)
  if nargout <= 1
    w= t/p;
  else
    [pr,pc] = size(p);
    x = t/[p; ones(1,pc)];
    w = x(:,1:pr);
    b = x(:,pr+1);
  end
end

%======================================================
function a = radbas1(n)
    a = exp(-(n.*n));
end

%======================================================
function out=mse(t,y)
%  function mse(t,y)
%    These should be dim x numpts
%    Output is then scalar
  temp=sum((t-y).^2);
  out=mean(temp);

end

%======================================================
function z=dist(w,p)
% To be consistent with Matlab's dist function, the first
% matrix is transposed with what it ought to be.
% Incoming is numpts x dim and dim x numpts
p=p'; %This should convert "dist" to "edm" input

% EDM starts here; w and p should both be numpts x dim at this point
[S,R] = size(w);
[Q,R2] = size(p);
p=p';
if (R ~= R2), error('Inner matrix dimensions do not match.\n'), end
z = zeros(S,Q);
if (Q<S)
p = p';
copies = zeros(1,S);
for q=1:Q
z(:,q) = sum((w-p(q+copies,:)).^2,2);
end
else
w = w';
copies = zeros(1,Q);
for i=1:S
z(i,:) = sum((w(:,i+copies)-p).^2,1);
end
end
z = z.^0.5;
end