function generate_mylw_C(bnet, fid, nnodes, onodes)
% GENERATE_MYLW_C Create the model-specific file needed by C_lw.c etc.
%
% GENERATE_MYLW_C(BNET, FID, NNODES) will allow an arbitrary pattern of evidence to be passed in
% as needed by C_lw.c
% GENERATE_MYLW_C(BNET, FID, NNODES, ONODES) hard-codes which nodes are observed, as needed by
% C_lw1.c and C_lw2.c

if ~exist('onodes')
  fixed_obs = 0;
else
  fixed_obs = 1;
  obs = zeros(1, nnodes);
  obs(onodes) = 1;
end

ns = bnet.node_sizes;
if fixed_obs
  fprintf(fid, 'void mylw(const int ev[], int nsamples, \n');
  fprintf(fid, '        double *samples, double weights[], const double *rnd) {\n\n');
else
  fprintf(fid, 'void mylw(const int ev[], const int obs[], int nsamples, \n');
  fprintf(fid, '        double *samples, double weights[], const double *rnd) {\n\n');
end

fprintf(fid, '#define NNODES %d\n', nnodes);
for i=1:nnodes
  j = i-1;
  fprintf(fid, '#define ns%d %d\n', j, ns(i));
end
fprintf(fid, '\n\n');

fprintf(fid, '#define SAMPLES(s,i) samples[(s)+nsamples*(i)]\n');
fprintf(fid, '#define RND(i,s) rnd[(i)+NNODES*(s)]\n');
fprintf(fid, '\n\n');


CPT_to_C(bnet, fid);
fprintf(fid, '\n\n');

fprintf(fid, 'int x[NNODES];\n');

fprintf(fid, 'double ');
for i=1:nnodes-1
  j = i - 1;
  fprintf(fid, ' cum%d[ns%d], ', j, j);
end
j = nnodes -1;
fprintf(fid, ' cum%d[ns%d];\n', j, j);
  
fprintf(fid, 'int s, i, nrnd;\n');
fprintf(fid, 'double w;\n');
fprintf(fid, '\n\n');

%fprintf(fid, 'printf("max_node= %%d\\n", max_node);\n');

if fixed_obs
  for i=1:nnodes
    j = i-1;
    if obs(i)
      fprintf(fid, 'x[%d] = ev[%d];\n', j, j);
    end
  end
  fprintf(fid, '\n\n');
else
  fprintf(fid, 'for (i = 0; i < NNODES; i++) {\n');
  fprintf(fid, '  if (obs[i])\n');
  fprintf(fid, '    x[i] = ev[i];\n');
  fprintf(fid, '  else\n');
  fprintf(fid, '    x[i] = -1;\n');
  fprintf(fid, '}\n');
end


fprintf(fid, 'nrnd = 0;\n');
fprintf(fid, 'for (s = 0; s < nsamples; s++) {\n');
%fprintf(fid, 'printf("sample %%d\\n", s);\n');
fprintf(fid, '    w = 1;\n');

fprintf(fid, '\n\n');
for i=1:nnodes
  j = i-1;
  ps = parents(bnet.adj_mat, i) - 1;
  pvals = '';
  for k=1:length(ps)
    pvals = sprintf('%s[x[%d]]', pvals, ps(k));
  end
  %fprintf(fid, 'if (%d > max_node) { goto end_of_net; }\n', j);
  if fixed_obs
    if obs(i)
      fprintf(fid, 'w *= CPT%d%s[x[%d]];\n\n', j, pvals, j);
    else
      fprintf(fid, 'x[%d] = sample_discrete(ns%d, CPT%d%s, cum%d, rnd[nrnd++]);\n\n', j, j, j, pvals, j);
    end
  else
    fprintf(fid, 'if (obs[%d])\n', j);
    fprintf(fid, '  w *= CPT%d%s[x[%d]];\n', j, pvals, j);
    fprintf(fid, 'else\n');
    fprintf(fid, ' x[%d] = sample_discrete(ns%d, CPT%d%s, cum%d, rnd[nrnd++]);\n\n', j, j, j, pvals, j);
  end
end
fprintf(fid, '\n\n');

%fprintf(fid, 'end_of_net:\n');
fprintf(fid, 'for (i = 0; i < NNODES; i++) SAMPLES(s,i) = x[i];\n');
fprintf(fid, 'weights[s] = w;\n');
fprintf(fid, '} }\n');
fprintf(fid, '\n\n');


%%%%%%%

function CPT_to_C(bnet, fid)
% CPT_TO_C Create C arrays containing all the CPTs
% CPT_to_C(bnet, fid)
%
% Example - for the file examples/discrete1.m this function creates
% double CPT0[2] = {0.500, 0.500};
% double CPT1[2][2] = {0.500, 0.500, 0.900, 0.100};
% double CPT2[2][2] = {0.800, 0.200, 0.200, 0.800};
% double CPT3[2][2][2] = {1.000, 0.000, 0.100, 0.900, 0.100, 0.900, 0.010, 0.990};


n = length(bnet.adj_mat);
ns = bnet.node_sizes;
for i=1:n
  e = bnet.equiv_class(i);
  assert(strcmp(bnet.bnodes{e}.type, 'tabular'));
  ps = bnet.bnodes{e}.parents;
  CPT = bnet.bnodes{e}.CPT;
  if length(ps)>=1
    CPT = permute(CPT, length(ps)+1:-1:1); % so last dim becomes first (C order)
  end
  CPT = CPT(:);
  j = i-1; % C numbers nodes from 0
  switch length(ps)
    case 0, fprintf(fid, 'double CPT%d[%d] = {', j, ns(i));
    case 1, fprintf(fid, 'double CPT%d[%d][%d] = {', j, ns(ps(1)), ns(i));
    case 2, fprintf(fid, 'double CPT%d[%d][%d][%d] = {', j, ns(ps(1)), ns(ps(2)), ns(i));
    case 3, fprintf(fid, 'double CPT%d[%d][%d][%d][%d] = {', j, ns(ps(1)), ns(ps(2)), ns(ps(3)), ns(i));
    case 4, fprintf(fid, 'double CPT%d[%d][%d][%d][%d][%d] = {', ...
	  j, ns(ps(1)), ns(ps(2)), ns(ps(3)), ns(ps(4)), ns(i));
    case 5, fprintf(fid, 'double CPT%d[%d][%d][%d][%d][%d][%d] = {', ...
	  j, ns(ps(1)), ns(ps(2)), ns(ps(3)), ns(ps(4)), ns(ps(5)), ns(i));
    case 6, fprintf(fid, 'double CPT%d[%d][%d][%d][%d][%d][%d][%d] = {', ...
	  j, ns(ps(1)), ns(ps(2)), ns(ps(3)), ns(ps(4)), ns(ps(5)), ns(ps(6)), ns(i));
    otherwise, error('too many parents');
  end
  for j=1:length(CPT)-1
    fprintf(fid, '%5.3f, ', CPT(j));
  end
  j = length(CPT);
  fprintf(fid, '%5.3f};\n', CPT(j));
end
