function dcm2dot(varargin)
%
% dcm2dot(dcm_files) renders the dcm models in the cell dcm_files as Graphviz .dot
% network models
%
% dcm_files is a cell of dcm model file names. These files are read
% and the model structure is written out as respective filename.dot
% If the model has been estimated the edges contain the weights as
% labels (DCM.Ep) with the estimated probability (DCM.Pp) in brakets.
%
% Example: To convert all files in the working directory to pdf, do in Matlab/octave:
% files = dir(); for n=1:size(files,1); ff{n}=files(n).name; end;
% dcm2dot(ff);
% and then in the same directory in a unix-shell with graphviz installed:
% ls *.dot | xargs -l1 -I '{}' dot {} -Tpdf -o '{}'.pdf
%
% 2012-08-29 Ingo R. Keck (ingokeck@ingokeck.de)
% Licence: GPL 3.0


if nargin==0
    dcm_files = {};
else
    dcm_files = varargin{1};
end
file_no = length(dcm_files);
for n=1:file_no
    
    dcm_model = dcm_files{n};
    %
    [path,filename,~] = fileparts(dcm_model);
    %
    try
        DCM = load(dcm_model,'DCM');
    catch
        msg = sprintf('error reading variable DCM from %s - skipping\n',dcm_model);
        warning(msg)
        continue
    end
    DCM = DCM.DCM;
    %
    % get regions
    regions = DCM.Y.name;
    % get inputs
    inputs = DCM.U.name;
    %
    dotfile = fullfile(path,[filename ,'.dot']);
    fid = fopen(dotfile,'w');
    fprintf(fid,'digraph %s { rankdir=LR;\n    charset = "UTF-8";\n',filename);
    % generate inputs
    for input = 1:length(inputs)
        fprintf(fid,'   "%s" [shape=box];\n', inputs{input});
    end
    % generate regions
    for region = 1:length(regions)
        fprintf(fid,'   "%s" [shape=oval];\n', regions{region});
    end
    % generate edges input -> region
    for input = 1:size(DCM.c,2)
        for region = 1:size(DCM.c,1)
            if DCM.c(region,input) == 1
                if isfield(DCM,'Ep')
                    fprintf(fid, '  "%s" -> "%s" [label="%.2f\\n(%.2f)"]\n',...
                        inputs{input},regions{region},DCM.Ep.C(region,input),DCM.Pp.C(region,input));
                else
                    fprintf(fid,'   "%s" -> "%s" [label=""];\n',inputs{input},regions{region});
                end
            end
        end
    end
    % connect regions
    for reg1 = 1:size(DCM.a,2)
        for reg2 = 1: size(DCM.a,1)
            if DCM.a(reg2,reg1) == 1
                if isfield(DCM,'Ep')
                    fprintf(fid,'   "%s" -> "%s" [label="%.2f\\n(%.2f)"]\n',...
                        regions{reg1},regions{reg2},DCM.Ep.A(reg2,reg1),DCM.Pp.A(reg2,reg1));
                else
                    fprintf(fid,('   "%s" -> "%s" [label=""] ;\n'),regions{reg1}, regions{reg2});
                end
            end
        end
    end
    % generate modulations
    for input = 1:size(DCM.b,3)
        for reg1 = 1:size(DCM.b,2)
            for reg2 = 1: size(DCM.b,1)
                if DCM.b(reg2,reg1,input) == 1
                    fprintf(fid,'   "%s" [shape=hexagon]; \n',[regions{reg1},' &#8594; ', regions{reg2}]);
                    if isfield(DCM,'Ep')
                        fprintf(fid,'   edge [style=dotted, label="%.2f\\n(%.2f)"]\n',...
                            DCM.Ep.B(reg2,reg1,input),DCM.Pp.B(reg2,reg1,input));
                    else
                        fprintf(fid,'   edge [style=dotted, label=""]\n');
                    end
                    fprintf(fid,'   "%s" -> "%s";\n',inputs{input},[regions{reg1},' &#8594; ', regions{reg2}]);
                end
            end
        end
    end
    
    fprintf(fid,'}\n');
    fclose(fid);
end