function generate_mylw_C(bnet, fid)

ns = bnet.node_sizes;
nnodes = length(bnet.adj_mat);

fprintf(fid, 'void mylw(const int ev[], const int obs[], int nsamples, \n');
fprintf(fid, '        double *samples, double weights[], const double *rnd) {\n\n');


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, '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');

fprintf(fid, 'nrnd = 0;\n');
fprintf(fid, 'for (s = 0; s < nsamples; 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 (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
fprintf(fid, '\n\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');
