function [C, R, w, info, conf] = rbf_rt_1(X, y, conf)%% Hybrid radial basis function network and regression tree.%% Solves a regression problem with inputs X and outputs y using a% regression tree and an RBF network selected using tree-guided% forward and backward subset selection. Returns the hidden unit% centres C, their radii R, the hidden-to-output weights w, some% additional information info and a fully instantiated configuration% structure conf.%% For further details of the function see:%%  'Matlab Routines for RBF Networks', 1999.%% For a description of the algorithm see:%%  'Further Advances in RBF Networks', 1999.%% Program name (for error messages).prog = 'rbf_rt_1';% Configuration specification.spec(1) = struct( ...  'comment', 'Scale length trial values', ...  'name',    'scales', ...  'type',    {{'vector', 'positive'}}, ...  'options', [], ...  'default', [1 2]);spec(2) = struct( ...  'comment', 'Minimum membership trial values', ...  'name', 'minmem', ...  'type', {{'vector', 'positive', 'integer'}}, ...  'options', [], ...  'default', 5);spec(3) = struct( ...  'comment', 'Model selection criterion', ...  'name', 'msc', ...  'type', 'string', ...  'options', {{'bic', 'gcv', 'loo'}}, ...  'default', 'bic');spec(4) = struct( ...  'comment', 'Factor for adjusted GCV', ...  'name', 'agcv', ...  'type', {{'number', 'positive'}}, ...  'options', [], ...  'default', 1);spec(5) = struct( ...  'comment', 'Basis function type', ...  'name', 'type', ...  'type', 'string', ...  'options', [], ...  'default', 'g');spec(6) = struct( ...  'comment', 'Basis function sharpness', ...  'name', 'exp', ...  'type', 'number', ...  'options', [], ...  'default', 2);spec(7) = struct( ...  'comment', 'Regularisation parameter', ...  'name', 'lambda', ...  'type', {{'number', 'nonnegative'}}, ...  'options', [], ...  'default', 0);spec(8) = struct( ...  'comment', 'Pseudo-inverse instead of normal', ...  'name', 'pseudo', ...  'type', 'number', ...  'options', {{0, 1}}, ...  'default', 0);spec(9) = struct( ...  'comment', 'Extra speed', ...  'name', 'speed', ...  'type', 'number', ...  'options', {{0, 1}}, ...  'default', 1);spec(10) = struct( ...  'comment', 'Resolve draws randomly', ...  'name', 'rand', ...  'type', 'number', ...  'options', {{0, 1}}, ...  'default', 1);spec(11) = struct( ...  'comment', 'Placement of centres', ...  'name', 'edge', ...  'type', 'number', ...  'options', {{0, 1}}, ...  'default', 0);spec(12) = struct( ...  'comment', 'Display graphical timer', ...  'name', 'timer', ...  'type', 'string', ...  'options', [], ...  'default', '');spec(13) = struct( ...  'comment', 'Verbose output', ...  'name', 'verb', ...  'type', 'number', ...  'options', {{0 1}}, ...  'default', 0);spec(14) = struct( ...  'comment', 'Print final summary', ...  'name', 'rprt', ...  'type', 'number', ...  'options', {{0, 1}}, ...  'default', 0);% Check number of arguments. Take special action if only one.switch nargincase 1  if isstring(X)    switch X    case 'conf'      conf_print(prog, spec)      return    case 'demo'      [mydemo, myclean] = eval(['demo_' prog]);      run_demo(mydemo, myclean)      return    otherwise      error([prog ': ''' X ''' unrecognised for single string argument'])    end  else    error([prog ': unrecognised type for single argument'])  endcase 2  if isstring(X) & isstring(y)    switch X    case 'conf'      conf_print(prog, spec, y)      return    otherwise      error([prog ': ''' X ''' unrecognised for double string argument'])    end  else    conf = [];  endcase 3otherwise  error([prog ': illegal number of arguments'])end% Check type of input arguments.if ~isnumeric(X) | ndims(X) ~= 2  error([prog ': first argument (X) should be a matrix'])endif ~isnumeric(y) | ndims(y) ~= 2 | size(y,2) ~= 1  error([prog ': second argument (y) should be a column vector'])end% Check for consistent size between X and y.[d,p] = size(X);if size(y,1) ~= p  error([prog ': number of samples inconsistent between X and y'])end% Check the configuration is okay and set defaults (if required).conf = conf_check(conf, spec, prog);% Scales and minimum memberships.scales = conf.scales;nscale = length(scales);mins = conf.minmem;nmin = length(mins);% Initialisation for flops and time.stats.comps = flops;stats.ticks = clock;% Set up static part of tree growing configuration.tcfg.verb = conf.verb;tcfg.rprt = conf.rprt;if conf.edge  tcfg.place = 'edge';else  tcfg.place = 'centre';end% Set up static part of RBF growing configuration.rcfg.msc = conf.msc;rcfg.agcv = conf.agcv;rcfg.type = conf.type;rcfg.exp = conf.exp;rcfg.lambda = conf.lambda;rcfg.pseudo = conf.pseudo;rcfg.speed = conf.speed;rcfg.rand = conf.rand;rcfg.rprt = conf.rprt;% Set up a timer for minmem (if required).if ~strcmp(conf.timer,'')  tconf.name = [prog ' minmem ' conf.timer];  tconf.n = nmin;  minm_tmr = get_tmr(tconf);end% Loop over mins.for m = 1:nmin  % Increment the timer for minmem.  if ~strcmp(conf.timer,'')    inc_tmr(minm_tmr)  end  % Grow a tree for this minmem.  tcfg.minm = mins(m);  tree = grow_tree(X, y, tcfg);  % Set up a timer for scales (if required).  if ~strcmp(conf.timer,'')    tconf.name = [prog ' scales ' conf.timer];    tconf.n = nscale;    scale_tmr = get_tmr(tconf);  end      for s = 1:nscale      % Increment the timer for scale.    if ~strcmp(conf.timer,'')      inc_tmr(scale_tmr)    end    % Grow an RBF on the current tree.    rcfg.scale = scales(s);    [rbf, tree] = growRBF(tree, rcfg);      % Is this the best tree so far.    if m == 1 & s == 1      btree = tree;      brbf = rbf;      bscale = rcfg.scale;      bminm = tcfg.minm;    else      if rbf.err < brbf.err        btree = tree;        brbf = rbf;        bscale = rcfg.scale;        bminm = tcfg.minm;      end    end  end  % Close the timer for scale.  if ~strcmp(conf.timer,'')    close(scale_tmr)  endend% Close the timer for minmem.if ~strcmp(conf.timer,'')  close(minm_tmr)end% Record the computations and time taken.stats.comps = flops - stats.comps;stats.ticks = etime(clock,stats.ticks);% Get the centres, radii and weights of the best tree.C = brbf.C;R = brbf.R;w = brbf.w;% Other potentially useful results go in 'info'.info.scale = bscale;info.minm = bminm;info.tree = btree;info.rbf = brbf;info.err = brbf.err;info.stats = stats;% Configuration for rbf_dm also goes in info.info.dmc = struct('type', conf.type, 'exp', conf.exp, 'bias', 0);function [rbf, tree] = growRBF(tree, conf)%% Grows an RBF over a tree.%% Initialisation for flops and time.if conf.rprt  comps = flops;  ticks = clock;end% Initialise and grow.[rbf, tree] = initNet(tree, conf);okay = 1;while okay  [rbf, okay] = splitNet(rbf, tree, conf);end% Record the computations and time taken.if conf.rprt  comps = flops - comps;  ticks = etime(clock,ticks);end% print a report.if conf.rprt  fprintf('-------- growRBF report --------\n')  fprintf('rbfs:       %d\n', rbf.m)  fprintf('actions:    '), fprintf('%d ', rbf.stats.actions), fprintf('\n')  fprintf('splits:     '), fprintf('%d ', rbf.stats.splitds), fprintf('\n')  fprintf('flops:      %d\n', comps)  fprintf('time:       %.3f\n', ticks)  fprintf('--------------------------------\n')endfunction [rbf, tree] = initNet(tree, conf)%% Initialise an RBF network prior to growing one over a tree.%% Calculate RBF columns for every node.dmcfg.type = conf.type;dmcfg.exp = conf.exp;for j = 1:tree.numn  n = tree.node(j);  tree.node(j).h = rbf_dm(tree.X, n.c, n.r * conf.scale, dmcfg);end% Initialise the network with the RBF from the root node of the tree.n = tree.node(1);rbf.C = n.c;rbf.R = n.r * conf.scale;rbf.H = n.h;rbf.m = 1;% MSC, variance, weights etc.[rbf.err, rbf.A, rbf.w, rbf.e, rbf.gam, rbf.AH] = learn(rbf.H, tree.y, conf);% Initialise the array of indexes which relate columns in H to tree nodes.rbf.node = 1;% Initialise the list of teminal tree nodes and their number.rbf.term = 1;rbf.numt = 1;% Stats.rbf.stats.actions = zeros(1,8);rbf.stats.splitds = zeros(tree.d,1);function [rbf, okay] = splitNet(rbf, tree, conf)%% Try to grow RBF further down the tree.%% For each terminal node, calculate the change in MSC% which would occur if that split was actually taken.% 8 possibilites: add both, add left, add right, add none;% replace both, replace left, replace right, remove the parent.dmsc = zeros(rbf.numt,8);for k = 1:rbf.numt    % Which node?  ind = rbf.term(k);  n = tree.node(ind);    % Has this node got children?  if isempty(n.split)      % No. It can't be split.    dmsc(k,:) = -1 * ones(1,8);  else    % Yes it can be split. What are its potential RBF contributions?    lh = tree.node(n.split.l).h;    rh = tree.node(n.split.r).h;    % Try adding both, left only, right only, or no change.    if (strcmp(conf.msc,'gcv') | strcmp(conf.msc,'bic')) & conf.speed == 1      dmsc(k,1:4) = learnQuick(rbf, tree, lh, rh, conf);    else      dmsc(k,1) = rbf.err - learn([rbf.H lh rh], tree.y, conf);      dmsc(k,2) = rbf.err - learn([rbf.H lh], tree.y, conf);      dmsc(k,3) = rbf.err - learn([rbf.H rh], tree.y, conf);      dmsc(k,4) = 0;    end    % Try replacing both.    if ~ismember(ind, rbf.node)      dmsc(k,5:7) = dmsc(k,1:3);      dmsc(k,8) = -1;    else      % Try replacing both, left only, right only, or just removing parent.      if (strcmp(conf.msc,'gcv') | strcmp(conf.msc,'bic')) & conf.speed == 1        % I don't know how to do this quickly yet, so avoid the issue for now.        % Meanwhile, if you want this feature, don't use conf.speed = 1.        dmsc(k,5:8) = -1 * ones(1,4);      else        keep = find(rbf.node ~= ind);        dmsc(k,5) = rbf.err - learn([rbf.H(:,keep) lh rh], tree.y, conf);        dmsc(k,6) = rbf.err - learn([rbf.H(:,keep) lh], tree.y, conf);        dmsc(k,7) = rbf.err - learn([rbf.H(:,keep) rh], tree.y, conf);        dmsc(k,8) = rbf.err - learn(rbf.H(:,keep), tree.y, conf);      end    end  end  end% What's the best action/terminal combo?[mdmsc, k] = max(dmsc, [], 1);[mmdmsc, act] = max(mdmsc, [], 2);% If there are no favourable changes possible, signal that and quit.if mmdmsc < 0  okay = 0;  returnelse  okay = 1;end% Record action in info stats.rbf.stats.actions(act) = rbf.stats.actions(act) + 1;% Index of the node (in tree.node) assuming no tie.ind = rbf.term(k(act));% In the case of a tie amongst two or more terminals,% randomly select one rather than always taking the first (k(act)).% Only do this if indeterminancy is permitted (by conf.rand == 1).if conf.rand  winners = find(dmsc(:,act) == mmdmsc);  numwins = length(winners);  if numwins > 1    winner = winners(ceil(rand * numwins));    ind = rbf.term(winner);  endend% Tree node being split and its children.n = tree.node(ind);lind = n.split.l;rind = n.split.r;ln = tree.node(lind);rn = tree.node(rind);% Record the split dimension. Maybe this should depend on action.dim = n.split.dim;rbf.stats.splitds(dim) = rbf.stats.splitds(dim) + 1;% Remove this parent from the terminal list.rbf.term(find(rbf.term == ind)) = [];rbf.numt = rbf.numt - 1;% Make the children of this node terminals.rbf.term = [rbf.term lind rind];rbf.numt = rbf.numt + 2;% Less awkward names.p = tree.p;y = tree.y;% Update the RBF network. Are we replacing the parent?if act <= 4  % No. Keep all the old columns.  m = rbf.m;  keep = 1:m;else  % Yes. Discard the parent if it's there.  if ismember(ind, rbf.node)    m = rbf.m - 1;    keep = find(rbf.node ~= ind);  else    m = rbf.m;    keep = 1:m;  endend% Are we adding both, one or none?if act == 1 | act == 5  % Both.  m = m + 2;  nc = [ln.c rn.c];  nr = [ln.r*conf.scale rn.r*conf.scale];  nh = [ln.h rn.h];  nt = [lind rind];elseif act == 2 | act == 6  % Left only.  m = m + 1;  nc = ln.c;  nr = ln.r*conf.scale;  nh = ln.h;  nt = lind;elseif act == 3 | act == 7  % Right only.  m = m + 1;  nc = rn.c;  nr = rn.r*conf.scale;  nh = rn.h;  nt = rind;else  % None.  nc = [];  nr = [];  nh = [];  nt = [];end% Update the RBF parameters.rbf.m = m;rbf.C = [rbf.C(:,keep) nc];rbf.R = [rbf.R(:,keep) nr];rbf.H = [rbf.H(:,keep) nh];[rbf.err, rbf.A, rbf.w, rbf.e, rbf.gam, rbf.AH] = learn(rbf.H, tree.y, conf);% Update the array of indexes which relate columns in H to tree nodes.rbf.node = [rbf.node(keep) nt];function [msc, A, w, e, gam, AH] = learn(H, y, conf)%% Get MSC and other things from from H and y.%[p,m] = size(H);if conf.pseudo  A = pinv(H'*H + conf.lambda * eye(m));else  A = inv(H'*H + conf.lambda * eye(m));endw = A * (H' * y);f = H * w;e = y - f;if conf.lambda == 0  gam = m;else  gam = m - conf.lambda * trace(A);endswitch conf.msccase 'gcv'  msc = p * (e' * e) / (p - conf.agcv * gam)^2;case 'loo'  dP = ones(p,1) - sum((H*A).*H,2);  msc = sum((e ./ dP).^2) / p;case 'bic'  msc = (p + (log(p) - 1) * gam) * (e' * e) / (p * (p - gam));otherwise  error(['growRBF: learn: bad MSC (' conf.msc ')'])endif nargout > 5  AH = A * H';endfunction changes = learnQuick(rbf, tree, h1, h2, conf)%% Speedy GCV or BIC calculation.%% Initialise.msc1 = rbf.err;H1 = rbf.H;A1H1 = rbf.AH;P1y = rbf.e;gam1 = rbf.gam;y = tree.y;lam = conf.lambda;changes = zeros(1,4);[p,m1] = size(H1);% Add h1 (left).m2 = m1 + 1;P1h1 = h1 - H1 * (A1H1 * h1);h1P1h1 = h1'*P1h1;yP1h1 = y'*P1h1;e2 = P1y - (yP1h1/(lam+h1P1h1)) * P1h1;sse2 = e2' * e2;if lam == 0  gam2 = m2;else  h1P1P1h1 = P1h1'*P1h1;  gam2 = gam1 + h1P1P1h1/(lam + h1P1h1);endswitch conf.msccase 'gcv'  msc2 = p * sse2 / (p - conf.agcv * gam2)^2;case 'bic'  msc2 = (p + (log(p) - 1) * gam2) * sse2 / (p * (p - gam2));otherwise  error(['growRBF: learnQuick: bad MSC (' conf.msc ')'])endchanges(2) = msc1 - msc2;% Add h2 (right) as well.m3 = m2 + 1;P2y = P1y - (yP1h1/(lam + h1P1h1)) * P1h1;A1H1h2 = A1H1 * h2;P1h2 = h2 - H1 * A1H1h2;h2P1h1 = h2' * P1h1;P2h2 = P1h2 - (h2P1h1/(lam + h1P1h1)) * P1h1;h2P2h2 = h2'*P2h2;yP2h2 = y'*P2h2;e3 = P2y - (yP2h2/(lam+h2P2h2)) * P2h2;sse3 = e3' * e3;if lam == 0  gam3 = m3;else  h2P2P2h2 = P2h2'*P2h2;  gam3 = gam2 + h2P2P2h2/(lam + h2P2h2);endswitch conf.msccase 'gcv'  msc3 = p * sse3 / (p - conf.agcv * gam3)^2;case 'bic'  msc3 = (p + (log(p) - 1) * gam3) * sse3 / (p * (p - gam3));endchanges(1) = msc1 - msc3;% What about adding h2 (right) only.h1 = h2;m2 = m1 + 1;P1h1 = h1 - H1 * (A1H1 * h1);h1P1h1 = h1'*P1h1;yP1h1 = y'*P1h1;e2 = P1y - (yP1h1/(lam+h1P1h1)) * P1h1;sse2 = e2' * e2;if lam == 0  gam2 = m2;else  h1P1P1h1 = P1h1'*P1h1;  gam2 = gam1 + h1P1P1h1/(lam + h1P1h1);endswitch conf.msccase 'gcv'  msc2 = p * sse2 / (p - conf.agcv * gam2)^2;case 'bic'  msc2 = (p + (log(p) - 1) * gam2) * sse2 / (p * (p - gam2));endchanges(3) = msc1 - msc2;