function error = van_bnpm_analysis

SCCSid = '1.4';
SPMid  = spm('SFnBanner','bnpm toolbox','SCCSid');
warning('off')
H = spm('FnUIsetup','VAN BnPM Analysis',0);
error = 0;

%----- Result directory ---------%
if exist('spm_get')
    BnPMstr.result_dir  = spm_str_manip(spm_get(-1,'*','Results directory',pwd),'r');
else
    BnPMstr.result_dir  = spm_str_manip(spm_select(1, 'dir', 'Results directory', [], pwd), '.*');
end
wdir = BnPMstr.result_dir ;
disp(['Results directory:  ' wdir]);
cd(wdir);

%------ Only do regression, BPM.type = 'REGRESSION' ------%
swd = pwd;        
[flist{1}{1}, swd] = wfu_bpm_get_any_file(sprintf('Dependent modality file list for %s',  'REGRESSION'), swd);
title_group{1}{1}  = spm_input('Name of the main modality ','+2','s');

% ------ Reading the main modality files ------%
yfile_names = wfu_bpm_read_flist(flist{1}{1});
nsubj = size(yfile_names,1); % number of subjects
for k = 1:nsubj
    yfile_names_subjs{k} = yfile_names(k,:); % name of each main modality image file
    index = find(yfile_names_subjs{k}==' ');
    if ~isempty(index)
        yfile_names_subjs{k} = [yfile_names(k,1:index-1),yfile_names(k,index+1:size(yfile_names(k)))];
    end
end
BnPMstr.fname = yfile_names_subjs{1};

% ------- determining the size of the images ---------------------------%
Vtemp = spm_vol(yfile_names_subjs{1});
nRows      = Vtemp.dim(1);
nCols      = Vtemp.dim(2);
nSlices = Vtemp.dim(3);
BnPMstr.dim = Vtemp.dim;

% ------ imageing covariates ------ %        
nModalities = spm_input('How many imaging covariates? ','+1');
for cv = 1:nModalities
    title_img_cov{cv} = spm_input(sprintf('Name of imaging covariate#%d ',cv),'+1','s');
end
       
for m = 2:nModalities+1                
    [flist{m}{1}, swd] = wfu_bpm_get_any_file(sprintf('Modality %d file list ', m), swd);                                     
end
% ------ Reading the regressors files ------%
nr = nModalities; % number of regressors
xfile_names = cell(1,nr);
for m = 1:nr
    xfile_names{m} = wfu_bpm_read_flist(flist{m+1}{1});
    for k = 1:nsubj
        xfile_names_subjs{m}{k} = xfile_names{m}(k,:);
        index = find(xfile_names_subjs{m}{k}==' ');
        if ~isempty(index)
            xfile_names_subjs{m}{k} = [xfile_names{m}(k,1:index-1),xfile_names{m}(k,index+1:size(xfile_names{m}(k)))];
        end
    end
end


% ------ non-imaging covariates ------ %         
niCov = spm_input('Any non-imaging covariates? ','+1','y/n',[1,0],2);
        
% ----- No regressors error ---------------------%
if (nModalities == 0) & (niCov == 0)
    display('Regression requires at least one covariate')
    error = 1;
    return
end
        
% --- Selecting the file with non-imaging covariates information--%       
% conf is no-image covariates        
if niCov > 0
    [conf, swd] = wfu_bpm_get_any_file(sprintf('Non-imaging covariates file for %s', 'REGRESSION'), swd);
    col_cof = load(conf);
    nIC = size(col_cof,2);
    for cv = 1:nIC
        title_ni_cov{cv} = spm_input(sprintf('Name of non-imaging covariates#%d ',cv),'+1','s');
    end
else
    conf = [];
end

        
%----- brain mask -----% 
if spm_input('Apply a predefined brain mask? ','+1','y/n',[1 0],0)
     %   ---------- Enter the mask ----------------%
    [mask, swd]  = wfu_bpm_get_image_file('Select mask', swd);
else
    tt = spm_input('Threshold type', '+1', 'Proportional|Absolute',[],1);
    if strcmp(tt, 'Proportional')
        fthr = 1;
        thr = spm_input('threshold ','+1','e',0.1);   
    else
        fthr = 2;
        thr = spm_input('threshold ','+1','e',0.1);   
    end
    mask = [];
end

%-Ask about variance smoothing & volumetric computation
%-----------------------------------------------------------------------
vFWHM=[0,0,0];	% FWHM for variance smoothing
bVarSm=0;	% Flag for variance smoothing

vFWHM = spm_input('FWHM(mm) for Variance smooth','+1','e',0);
if length(vFWHM)==1
	vFWHM = vFWHM * ones(1,3);
elseif length(vFWHM)==2
	vFWHM = [vFWHM, 0];
else
	vFWHM = reshape(vFWHM(1:3),1,3);
end
if ~all(vFWHM==0), bVarSm=1; end

% ------ input permutation times ------ %
nPiCond = gamma(nsubj+1);
bAproxTst = spm_input(sprintf('%d Perms. Use approx. test?',nPiCond),...
							'+1','y/n')=='y';
if bAproxTst
      tmp = 0;
      while ((tmp>nPiCond) | (tmp==0) )
      tmp = spm_input(sprintf('# perms. to use? (Max %d)',nPiCond),'+0');
      tmp = floor(max([0,tmp]));
      end
      if (tmp==nPiCond), bAproxTst=0; else nPiCond=tmp; end
end

% ------ Contrast parameter ------ %
tt = spm_input('Contrast Type', '+1', 'F|T',[],1);
if strcmp(tt, 'F') 
    STAT = 'F';
    sr = spm_input(sprintf('Select the regressors (1/0)'),'+1','e'); 
    c = sr(:);
else
    ts = spm_input('T-sign', '+1', '+|-',[],1);
    STAT = 'T';
    c = spm_input(sprintf('Enter new contrast '),'+1','e');   
    c=c(:);
    if strcmp(ts,'-')
         c = -c;
    end
end

% ------ SPM, BnPM directory ------%
spm_dir  = spm_str_manip(spm_select(1, 'dir', 'SPM directory', [], pwd), '.*');
bnpm_dir  = spm_str_manip(spm_select(1, 'dir', 'BnPM directory', [], pwd), '.*');

% ------ Number of jobs ------%
BnPMstr.njob = spm_input(sprintf('Enter # of jobs '),'+1','e'); 

% input walltime
walltime = spm_input(sprintf('Enter walltime hour:minute:second'),'+1','s');

%=======================================================================
% - C O R R E C T   P E R M U T A T I O N
%=======================================================================
% Work out correct permuation completely. Separating the first
BPM.type = 'REGRESSION';
BPM.result_dir = BnPMstr.result_dir;
BPM.flist = wfu_bpm_write_flist_gen(wdir,flist,nModalities+1);
BPM.conf = conf;
BPM.mask = mask;
if fthr == 1
    BPM.mask_pthr = thr;
else
    BPM.mask_athr = thr;
end
BPM.DMS(1) = 1;
% Add mean as part of titles

BPM.titles{1} = title_group{1}{1};

BPM.robust = 0;
if niCov > 0
    BPM.DMS(2) = size(col_cof,2);
        
    for k = 1:BPM.DMS(2)
        BPM.titles{k+BPM.DMS(1)} = title_ni_cov{k};
    end
else
    BPM.DMS(2) = 0;
end
if nModalities > 0
    BPM.DMS(3) = nModalities;
    for k = 1:BPM.DMS(3)
        BPM.titles{1} = 'mean'; % add mean to contrast
        BPM.titles{k+BPM.DMS(1)+BPM.DMS(2)} = title_img_cov{k};
    end
else
    BPM.DMS(3) = 0;
end

% ----- Treatment of the mask ----------------%
if nargin < 1  
    master_name = fullfile(BPM.result_dir,'master_flist.txt');
end    
file_names_mod = wfu_bpm_read_flist(master_name);                     
[file_names_subjs,no_grp] = wfu_bpm_get_file_names( file_names_mod(1,:) ); 

if isempty(BPM.mask)
    %  ---- Generating the default mask ------------- %    
%     if nargin < 1  
%         master_name = fullfile(BPM.result_dir,'master_flist.txt');
%     end    
%     file_names_mod = wfu_bpm_read_flist(master_name) ;                     
%     [file_names_subjs,no_grp] = wfu_bpm_get_file_names( file_names_mod(1,:) ); 
    k = 1;        
    for m = 1:no_grp
        for n = 1:size(file_names_subjs{m},1)
            CP1{k} = file_names_subjs{m}(n,:);
            k = k + 1;
        end
    end    
    P1 = strvcat(CP1);
    [inter_mask1] = wfu_bpm_mask(P1, BPM);
    BPM.mask      = fullfile(BPM.result_dir,'mask.img')    ;
else
    % ----- mask provided by the user ----------%
    Vmask  = spm_vol(BPM.mask);
    Vsubj1 = spm_vol(file_names_subjs{1}(1,:));
    if ~isequal(Vmask.dim(1:3),Vsubj1.dim(1:3))   
        C{1}       = file_names_subjs{1}(1,:);
        C{2}       = BPM.mask;
        P          = strvcat(C);    
        flag.mean  = 0;
        flag.which = 1;
        spm_reslice(P,flag);  
        [filepath,fname,ext]  = fileparts(BPM.mask);
        rfname                = strcat('r',fname);
        BPM.mask              = fullfile(filepath,strcat(rfname,ext));
    end
end

% ------ Reslicing the ROI mask -----------------------------%

if isfield(BPM,'mask_ROI') | isfield(BPM,'mask_ancova_ROI')
    Vmask      = spm_vol(BPM.mask);
    C{1}       = BPM.mask;
    if isfield(BPM,'mask_ROI') 
        Vmask_ROI                = spm_vol(BPM.mask_ROI);
        C{2}                     = BPM.mask_ROI;   
    else
        Vmask_ROI                = spm_vol(BPM.mask_ancova_ROI);
        C{2}                     = BPM.mask_ancova_ROI;
    end
    
    if ~isequal(Vmask.dim(1:3),Vmask_ROI.dim(1:3))              
        P          = strvcat(C);    
        flag.mean  = 0;
        flag.which = 1;
        spm_reslice(P,flag);
        if isfield(BPM,'mask_ROI') 
            [filepath,fname,ext]     = fileparts(BPM.mask_ROI);
            rfname                   = strcat('r',fname);
            BPM.mask_ROI             = fullfile(filepath,strcat(rfname,ext));
        else
            [filepath,fname,ext]     = fileparts(BPM.mask_ancova_ROI);
            rfname                   = strcat('r',fname);
            BPM.mask_ancova_ROI      = fullfile(filepath,strcat(rfname,ext));
        end
        
    end
end

%  ------------ Reslicing the imaging modalities -------------%

[BPM] = wfu_bpm_reslice_newflists(BPM);   

% set file name for design matrix information file
BPM.XtX = fullfile(BPM.result_dir, 'XtX');
warning off MATLAB:divideByZero;

% ---- deleting the content of the result directory ----------%

ddir = pwd;
cd(BPM.result_dir);

delete('S*.img');delete('S*.hdr'); delete('S*.mat');
delete('R*.img');delete('R*.hdr'); delete('R*.mat');
delete('T*.img');delete('T*.hdr'); delete('T*.mat');
delete('b*.img');delete('b*.hdr'); delete('b*.mat');
delete('s*.img');delete('s*.hdr'); delete('s*.mat');
delete('X.mat') ; delete('*PM.mat'); delete('X*.*');

% -----execute after result dir files deleted----- %

[results] = wfu_bpm_execute(BPM);
results.V = wfu_bpm_hdr_struct(results.V);
% ---------- storing the beta coefficients file names ------------ %

%------ writing the results ------%        
for k = 1:results.nr
    BPM.beta{k} = fullfile(BPM.result_dir,sprintf('beta%03d.img', k));           
end
%         % ---------- storing the sig in an img file ---------------------- %
%         results.V.fname = fullfile(BPM.result_dir,'sig.img');   
%         BPM.sig    =  results.V.fname ;
%         spm_write_vol(results.V, results.sig) ;                         
        % --------- storing the residuals file names -------------------- % 
for k = 1:results.nsubj
    BPM.E{k} = fullfile(BPM.result_dir,sprintf('Res%03d.img', k));           
end 
P = strvcat(BPM.E);
% ---------- storing the sig2 in an img file ---------------------- %
results.V.fname = fullfile(BPM.result_dir,'sig2.img');   
BPM.sig2    =  results.V.fname ;
spm_write_vol(results.V, results.sig2) ;
% ----------- storing X file -------------------------------------- %
Xfname = fullfile(BPM.result_dir,'X');   
X = results.X;
save(Xfname, 'X')  ; 
% ----------- filling other BPM fields ----------------------------- %
BPM.X = Xfname     ;
BPM.dof = results.dof      ;
    
%%%%%%%%%%%%% t statistic
BPM.contrast{1} = c;
% ---------- Deleting the residual files -------------------------%
ddir = pwd ;
cd(BPM.result_dir);
delete('R*.img');
delete('R*.hdr');
delete('R*.mat');  
BPM.Type_STAT = STAT;
BPM.vFWHM = vFWHM;
BPM.bVarSm = bVarSm;
BPM = wfu_bpm_con_man(BPM);    
Vt   = spm_vol(BPM.Stat{1});
T0 = spm_read_vols(Vt);



% ------ Create jobs ------%

for part = 1:BnPMstr.njob
                % For results
            pwd;
            rfolder = ['Result_Part',num2str(part)];
            mkdir(rfolder);
            cd(rfolder); 
            % master file
            fid = fopen('master_flist.txt','w');
            numfmod = nr+1;
            fprintf(fid,'%d\n',numfmod);            
            for m = 1:nr+1
                fmod{m} = [...
                    pwd,'//file_mod',num2str(m),'.txt'];
                fprintf(fid,'%s\n',fmod{m});
            end
            fclose(fid);
            % fmod file for y
            yfile = [...
                pwd,'/ypart',num2str(part),'.flist'];
            fid = fopen('file_mod1.txt','w');
            fprintf(fid,'%s',yfile);
            fclose(fid);
            % y file
            fid = fopen(['ypart',num2str(part),'.flist'],'w');
            fprintf(fid,'%d\n',nsubj);
            for k = 1:nsubj
                yname = ['./../Y_Part',num2str(part),'/Y',num2str(k),'_part',num2str(part),'.nii'];
                fprintf(fid,'%s\n',yname);
            end
            fclose(fid);
            % fmod file for x & x file
            for m = 1:nr
                xfile = [...
                    pwd,'/Reg',num2str(m),'_xpart',num2str(part),'.flist'];
                fid = fopen(['file_mod',num2str(m+1),'.txt'],'w');
                fprintf(fid,'%s',xfile);
                fclose(fid);
                fid = fopen(['Reg',num2str(m),'_xpart',num2str(part),'.flist'],'w');
                fprintf(fid,'%d\n',nsubj);
                for k = 1:nsubj
                    xname = ['./../Regressor',num2str(m),'X_Part',num2str(part),'/Reg',num2str(m),'X',num2str(k),'_part',num2str(part),'.nii'];
                    fprintf(fid,'%s\n',xname);
                end
                fclose(fid);
            end

    cd ..
    
    % ------ write matlab command file ------%
    fid = fopen(['bnpm_cluster',num2str(part),'.m'],'w');
    path = ['addpath ',spm_dir,' ',bnpm_dir];
    fprintf(fid,'%s\n',path);
    fprintf(fid,'%s\n','bnpm_startup;');
        % Replace the default stream at MATLAB startup, using a stream whose
    % seed is based on clock, so that rand will return different values in
    % different MATLAB sessions. 
    randinit = 'RandStream.setDefaultStream(RandStream(''mt19937ar'',''seed'',sum(100*clock)))';
    fprintf(fid,'%s\n',randinit);
    
    if isempty(mask) % no mask
        if isempty(conf) % no nonimage
            commend = ['van_bnpm(''REGRESSION'',''master_flist.txt'',[',num2str(c(1)),' ',num2str(c(2)),'],''',STAT,...
                ''',[],[],[',num2str(fthr),' ',num2str(thr),'],''',BnPMstr.result_dir,'Result_Part',num2str(part),''',[],[],',num2str(nPiCond),')'];
        else % nonimage covariate
            commend = ['van_bnpm(''REGRESSION'',''master_flist.txt'',[',num2str(c(1)),' ',num2str(c(2)),'],''',STAT,...
                ''',[],''',conf,''',[',num2str(fthr),' ',num2str(thr),'],''',BnPMstr.result_dir,'Result_Part',num2str(part),''',[],[],',num2str(nPiCond),')'];
        end
    else % mask
        if isempty(conf) % no nonimage
            commend = ['van_bnpm(''REGRESSION'',''master_flist.txt'',[',num2str(c(1)),' ',num2str(c(2)),'],''',STAT,...
                ''',''',BnPMstr.result_dir,'Result_Part',num2str(part),'/mask.img'',[],[],''',BnPMstr.result_dir,'Result_Part',num2str(part),''',[],[],',num2str(nPiCond),')'];
        else
            commend = ['van_bnpm(''REGRESSION'',''master_flist.txt'',[',num2str(c(1)),' ',num2str(c(2)),'],''',STAT,...
                ''',''',BnPMstr.result_dir,'Result_Part',num2str(part),'/mask.img'',''',conf,''',''',BnPMstr.result_dir,'Result_Part',BnPMstr.num2str(part),''',[],[],',num2str(nPiCond),')'];
        end
    end
    fprintf(fid,'%s',commend);
    fclose(fid);
end

% save BnPM structure
fname = fullfile(BnPMstr.result_dir,'BnPMstr');   
save(fname, 'BnPMstr');

% FOR loop here to create and submit all jobs
for part = 1:BnPMstr.njob
    fp = fopen(['MTbnpm',num2str(part),'.pbs'],'wt');

    donefile = sprintf([BnPMstr.result_dir,'job%.0f.done'],rand()*100000);

    fprintf(fp,'%s\n','#!/bin/bash');
    fprintf(fp,'%s\n','#PBS -M xue.yang@vanderbilt.edu');
    fprintf(fp,'%s\n','#PBS -l nodes=1:ppn=1:x86');
    fprintf(fp,'#PBS -l walltime=%s\n',walltime);
    fprintf(fp,'%s\n','#PBS -l mem=2000mb');
    fprintf(fp,'%s\n','#PBS -o xue-cluster.txt');
    fprintf(fp,'%s\n','#PBS -j oe');

    fprintf(fp,'%s\n',['matlab -nodesktop < ',BnPMstr.result_dir,'bnpm_cluster',num2str(part),'.m >matlab.output']);
    fprintf(fp,'touch %s\n',donefile);
    fclose(fp);

    if part==1
        system(sprintf('qsub %s','-m ba MTbnpm1.pbs'));
    else
        if part == numparts
            system(sprintf('qsub %s',['-m e MTbnpm',num2str(part),'.pbs']));
        else
            system(sprintf('qsub %s',['MTbnpm',num2str(part),'.pbs']));
        end
    end
    
end

% CLOSE FOR LOOP 

% NEW FOR LOOP - wait fro all jobs to finish
while(length(dir(donefile))<1)
    disp('Looking...');
    pause(10);
end
% END FOR LOOP
return
end



%
function [fileh, swd, name] = wfu_bpm_get_any_file(title, dir)
if exist('spm_get')
    fileh = spm_get(1, '*', title);
else
    fileh = spm_select(1, 'any', title, [], dir, '.*');
end
[swd, name] = fileparts(fileh);
end


%
function [fileh, swd, name] = wfu_bpm_get_image_file(title, dir)
if exist('spm_get')
    fileh = spm_get(1, '*.img', title);
else
    fileh = spm_select(1, 'IMAGE', title, [], dir, '.*');
end
[swd, name] = fileparts(fileh);
end
