function [ width ] = NicForest_CalculateWidth(xtrain,ytrain,opt)
%NICFOREST_CALCULATEWIDTH	Calculate width using a mini-forest or
%regression
%	[ width ] = NicForest_CalculateWidth(xtrain,ytrain,opt) calculates a 
%   reasonable starting value for the width parameter for the ensemble 
%   forest development.
%
%	Inputs:
%		xtrain      - Training features
%		ytrain      - Training targets
%       opt         - Number of trees to be used in development
%
%	Outputs:
%		opt.Width   - Scalar initial value for the intercept prior's width
%		
%	Example
%		[ width ] = NicForest_CalculateWidth(xtrain)
%	
%	See also NICFOREST NICFOREST_TRAIN

%	$LastChangedBy: alistair $
%	$LastChangedDate: 2012-05-30 12:21:30 +0100 (Wed, 30 May 2012) $
%	$Revision: 21 $
%	Originally written on GLNXA64 by Alistair Johnson, 09-May-2012 16:26:13
%	Contact: alistairewj@gmail.com

Ntrees = opt.Trees;
num_tar = numel(unique(ytrain));
if num_tar==1
    error('Only one class provided');
elseif num_tar==2 % binary classification
    %=== Do a quick MCMC to find a reasonable width
    opt = forest_opt_set(opt,...
        'Iterations',20000,...
        'Save',2000,...
        'Resets', 1,...
        'UpdatedTrees', 2,...
        'BurnIn', 20, ...
        'Width', []);
    
    %=== Split into 2 folds + train 2 models
    idxSplit = false(size(xtrain,1),1);
    idxSplit(1:2:end) = true;
    
%     group=opt.Group;
%     group_uniq = unique(group,'first');
    
    %=== if there is more than 1 group, ignore it for width estimation
    % this is probably temporary
    opt1 = opt; opt2 = opt;
    N1 = sum(idxSplit); N2 = numel(idxSplit)-N1;
    opt1.Group = ones(N1,1);
    opt2.Group = ones(N2,1);
    
    [ forests1  ] = NicForest_train(xtrain(idxSplit,:),ytrain(idxSplit,:),opt1);
    [ forests2  ] = NicForest_train(xtrain(~idxSplit,:),ytrain(~idxSplit,:),opt2);
    
    ypred = zeros(size(ytrain,1),1);
    [ ypred1 ] = NicForest_apply_quick( forests1 , xtrain(~idxSplit,:) );
    ypred(~idxSplit) = ypred1;
    [ ypred2 ] = NicForest_apply_quick( forests2 , xtrain(idxSplit,:) );
    ypred(idxSplit) = ypred2;
    
    % width^2*Ntrees/4 ~ var(logit(Pi)) where Pi is the pred from reasonable model.
    width = 2*sqrt(var(logit(ypred))/Ntrees);
else
    
    width = 2*Ntrees^-0.5; % ytrain variance is normalized to be 1 in fcn
    % width = sqrt(4*var(ytrain)/Ntrees);
    
end

end

function [p] = logit(p) % logit
p = log(p) - log(1-p);
end