function dcmsetup = sweetDCMvariate( varargin )
%SWEETDCMVARIATE   Create DCM models within a given model space based on a
%   given template file. Template files should be created with the
%   official SPM functions or compatible implementation. This tool will
%   then create all variations of connectivity priors as defined in the
%   template model. 
%
%   Connectivity priors can either be 0 (no connection), 1 (connection in
%   all models) or 999 (connections without a priori knowledge). 999 is
%   used to indicate that models with and without this connection should be
%   created.
%
%   Example: DCM_tmpl.mat:
%            DCM.a = [   1  0   
%                      999  1   ]
%            sweetDCMvariate will create two models:
%            DCM_m0001.mat:
%            DCM.a = [   1  0   
%                        0  1   ]
%            DCM_m0002.mat:
%            DCM.a = [   1  0   
%                        1  1   ]
%
%   sweetDCMvariate( dcmtemplatedir, dcmtemplatefile )
%   dcmtemplatedir    - Directory of template file
%                       e.g. pwd
%   dcmtemplatefile   - Valid DCM.mat-file with template information
%                       e.g. 'DCM_tmpl.mat'
%                       
%   sweetDCMvariate( ..., dcmtemplateprefix )
%   dcmtemplateprefix - Prefix used for template creation
%                       Default: ''
%
%   sweetDCMvariate( dcmsetup )
%   dcmsetup          - Structure containing
%                       dcmsetup.dcmtemplatedir
%                       dcmsetup.dcmtemplatefile
%                       dcmsetup.dcmtemplateprefix
%
%
%Copyright Ronald Sladky - sweetneuron.at
%Last edited: 2016/12/27



    if nargin == 1 && isstruct( varargin{1} )
        dcmtemplatedir    = dcmsetup.dcmtemplatedir;
        dcmtemplatefile   = dcmsetup.dcmtemplatefile;
        dcmtemplateprefix = dcmsetup.dcmtemplateprefix;
    elseif nargin == 2
        dcmtemplatedir    = char( varargin{1} );
        dcmtemplatefile   = char( varargin{2} );
        dcmtemplateprefix = '';      
    elseif nargin == 3
        dcmtemplatedir    = char( varargin{1} );
        dcmtemplatefile   = char( varargin{2} );
        dcmtemplateprefix = char( varargin{3} );        
    else
        help sweetDCMvariate
        return;
    end
    
    
    dcmtemplatedir    = [dcmtemplatedir '/' ];
    dcmdir            = [dcmtemplatedir '/' strrep( dcmtemplatefile, '.mat', '/' ) ];
    dcmsetup.dcmdir   = dcmdir;
    
    if ~isdir( dcmdir ) 
        mkdir( dcmdir )
    else
        display( 'WARNING DCM Template directory not empty.'  )
    end
    
    load( [ dcmtemplatedir dcmtemplatefile] );
    DCM_tmpl = DCM;
    
    sa = size( DCM.a ); 
    sb = size( DCM.b ); 
    
    sc = size( DCM.c );
    sd = size( DCM.d );
    
    tmpl = [ reshape(DCM.a,1,[]) reshape(DCM.b,1,[]) reshape(DCM.c,1,[]) reshape(DCM.d,1,[]) ];
    dcms = do_dcmvariate( tmpl );

    j = 1;
    for i = 1:size(dcms,1)       

        
        DCM = DCM_tmpl;
        
        DCM.a = reshape( dcms(i, 1:sa(1)*sa(2)), sa(1), sa(2));
        part = sa(1)*sa(2);
        if length(sb) == 2
            DCM.b = reshape( dcms(i, part+1:part+sb(1)*sb(2) ), sb(1), sb(2));
            part = part + sb(1)*sb(2);
        elseif length(sb) == 3
            DCM.b = reshape( dcms(i, part+1:part+sb(1)*sb(2)*sb(3) ), sb(1), sb(2), sb(3) );
            part = part + sb(1)*sb(2)*sb(3);
        end
            
        
        DCM.c = reshape( dcms(i, part+1:part+sc(1)*sc(2)), sc(1), sc(2));
        part = part + sc(1)*sc(2);
        
        DCM.d = reshape( dcms(i, part+1:part+sd(1)*sd(2)*sd(3)), sd(1), sd(2), sd(3) );

        
        if isfield( dcmsetup, 'requires_postprocessing' ) && dcmsetup.requires_postprocessing
            for row = dcmsetup.A'
                row = row';
                voi1 = char( row{2} );
                voi2 = char( row{3} );
                
                if row{5}>1 && row{5}<999
                    for lookup_row = dcmsetup.A'
                        if lookup_row{1} == row{5}
                            lookup_voi1 = char( lookup_row{2} );
                            lookup_voi2 = char( lookup_row{3} );
                            v = DCM.a( find( strcmp( lookup_voi2, dcmsetup.xY ) ),  find( strcmp( lookup_voi1, dcmsetup.xY ) ) );
                            DCM.a( find( strcmp( voi2, dcmsetup.xY ) ),  find( strcmp( voi1, dcmsetup.xY ) ) ) = v;
                        end
                    end
                end
            end
            
            for u = 1:length( dcmsetup.u_names )
                uname = char( dcmsetup.u_names{u} );
                is_input = 0;
                
                if isfield( dcmsetup.B, uname )
                    if dcmsetup.B.(uname) == 'A'
                        DCM.b(:,:,u) = DCM.a - eye( size(DCM.a,1), size(DCM.a,2) );
                    end
                end
            end
            
        end
        
        if sum( sum( sum(DCM.c) ) ~= 0 )
            newfilename = sprintf( '%s/DCM%s_m%04d.mat', dcmdir, dcmtemplateprefix, j );
            save( newfilename, 'DCM' );
            j = j+1;
        end
    end
    
    

    
function dcms = do_dcmvariate( tmpl )
    
    dcms = tmpl;
    
    for i=1:length( tmpl )
        for j=1:size( dcms, 1 )
            if dcms(j,i) == 999
                dcms(j,i) = 1;
                dcms(end+1,:) = dcms(j,:);
                dcms(end,i) = 0;
                
            end
        end
    end
    
    