%  test generation de donn�s
clear all
close all

% Metaparams
hyp = 'mcar';
ntests = 1;
nessais = 1;
m = 5000;

% INIT
missing = 2;
misv = -9999;
test = zeros(nessais,6);
chi2=[0.10,0.05,0.02,0.01,0.005,0.0025,0.001,0.0005];
nchi = 8;
ntests = ntests+1;

%%%%%%%%%% TESTS %%%%%%%%%%%
for essai=1:nessais

  N = ceil(rand*11) + 2
  base_proba = .15 + rand/4
  
  fan_in = ceil(rand*4)
  dag=mk_rnd_dag(N,fan_in)
  discrete=1:N;
  for i=1:N, node_sizes(i) = ceil(rand*4)+1; end
  bnet = mk_bnet(dag, node_sizes, discrete);
  node_sizes
  for i=1:N, bnet.CPD{i} = tabular_CPD(bnet, i); end

  if hyp == 'mcar',
    bnet_miss = gener_MCAR_net(bnet, base_proba);
    disp(' - MCAR net generated');
  elseif hyp == 'mar',
    bnet_miss = gener_MAR_net(bnet, base_proba);
    disp(' - MAR net generated');
  end
  
  [data, comp_data, bnet_miss, taux, bnet_orig, notok] = gener_data_from_bnet_miss(bnet_miss, m, base_proba ,1);
  
  test(essai,1,nchi)=notok;
  node0=[]; param0=[];
  for i=2:ntests
    % choose a node
    node=[]; tmp=0;
    while isempty(myintersect(node,node0)) & tmp<N-1,
      node = ceil(rand*N); tmp=tmp+1;
    end
    node0 = myunion(node0,node);

    % choose a param
    par = parents(bnet_miss.dag, node+N);
    sorted=[]; ord=[];
    for par1=1:length(par)
      ord(par1)=find(par(par1)==bnet_miss.order);
    end
    [tmp, sorted] = sort(ord);
    family = par(sorted); family(end+1)=node+N;

    state = []; stateM = [];
    for par1=1:length(family)-1
      state(par1) = ceil(rand*bnet_miss.node_sizes(family(par1)));
      if state(par1) == bnet_miss.node_sizes(family(par1)), stateM(par1) = misv; else, stateM(par1) = state(par1); end 
    end
    state(end+1) = missing; stateM(end+1)=misv;

    CPT = CPT_from_bnet(bnet_miss);
    value = CPT{node+N}(subv2ind(bnet_miss.node_sizes(family), state));

    if hyp == 'mcar',
      family = family-N;
    elseif hyp == 'mar',
      family = family-N;    1
    end

    % test this param
    datamat = bnt_to_mat(data(family,:), misv);
    indj=1:m;
    for par1 = 1:(length(family)-1)
      [tmp, indj] = find(datamat(par1,indj)==stateM(par1));
    end
    mtest = length(indj);
    [tmp, indj] = find(datamat(length(family),indj)==stateM(length(family)));

    toto=mtest; nbr_miss=length(indj);
    d = ((nbr_miss-base_proba*toto)^2)/(base_proba*toto) + (((toto-nbr_miss)-(1-base_proba)*toto)^2)/((1-base_proba)*toto);

    %df = prod(bnet_miss.node_sizes); 
    df = prod(bnet.node_sizes-1); 
    %for chiv=1:nchi,
    %  if d>chi2_table(chi2(chiv),df), test(essai,i,chiv) = 1; end
    %end
    test(essai,i) = chi2cdf(d,df);

    %fprintf('%2.1f%% OF MISSING DATA (%2.1f%%, Khi2 : %2.1f > %2.1f)\n', round(base_proba*10000)/100, round(taux*10000)/100, d, chi2(nchi)); 

  end
  clear bnet_miss bnet datamat data node0 family node_sizes CPT indj tmp 
end

test
ERREURS = length(find(test))/length(test(:))

