#include <apop.h>

/** Convert a crosstab to a PMF model, which you can use to make random draws, as a prior
    for Bayesian updating (via \ref apop_update), and so on.

\param  d  An input crosstab, expressed either in the matrix of an apop_data set.
\return    A PMF model which has a single line for each nonzero line of the crosstab. 
           The model points to a data set generated here, which you will need to 
           apop_data_free whenever you apop_model_free this model.
 */
apop_model *apop_crosstab_to_pmf (apop_data *d){
    Apop_stopif(!d, return NULL, 0, "NULL input.");
    Apop_stopif(!d->matrix, return NULL, 0, "You gave me an input set with no matrix data.");
    int tsize = d->matrix->size1 * d->matrix->size2;
    apop_data *outdata =apop_data_alloc(tsize, 2);
    outdata->weights = gsl_vector_alloc(tsize);

    size_t ctr = 0;
    double x;
    for (size_t i=0; i < d->matrix->size1; i++)
        for (size_t j=0; j < d->matrix->size2; j++)
            if ((x = gsl_matrix_get(d->matrix, i, j))) {
                apop_data_set(outdata, ctr, 0, i);
                apop_data_set(outdata, ctr, 1, j);
                gsl_vector_set(outdata->weights, ctr++, x);
            }
    return apop_estimate(outdata, apop_pmf);
}

int main(){ //test script with magic square.
   apop_model *pmf = apop_crosstab_to_pmf(
                      apop_data_fill(apop_data_alloc(3,3),
                                    8, 1, 6,
                                    3, 5, 7,
                                    4, 9, 2));
   int draw_ct=1e5;
   apop_data *draws= apop_data_alloc(draw_ct, 2);
   gsl_rng *r = apop_rng_alloc(213);
   for (int i=0; i < draw_ct; i++){
        Apop_row(draws, i, onerow);
        apop_draw(onerow->data, r, pmf);
   }
   apop_data *sum = apop_data_summarize(draws);
   assert(fabs(apop_data_get(sum, 0, .colname="mean")-1) < 1e-3);
   assert(fabs(apop_data_get(sum, 1, .colname="mean")-1) < 1e-3);
}