/* irt_wrapper.c
 *
 * Copyright (C) 2005 Stephane Germain
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/**
   \file

   \brief Functions to be called from Excel, R or anything that cannot used
   gsl vector and matrix.

   The matrices are stored in row major order.

   Only pointer are used as parameter.

   The memory for all the parameters should be allocated before.

   The functions are void.

   \author Stephane Germain <germste@gmail.com>
*/

#include <libirt.h>
#include <math.h>

/**
   \brief Set the verbosity level.

   @param[in] verbose The verbosity level.
*/
void set_verbose(int *verbose)
{
  libirt_verbose = *verbose;
}

/**
   \brief Set the default GSL error handling.

   @param[in] yes Set to 1 to continue on errors and to 0 to abort on error.
*/
void continue_on_error(int *yes)
{
  if(*yes) gsl_set_error_handler_off();
  else gsl_set_error_handler(NULL);
}

/**
   \brief Estimate the parameters of a 1, 2 or 3 PLM by MMLE
   or more generally (if a prior is used) by BME, or estimate
   the response functions directly by PMMLE or kernel.

   @param[in] nbr_subject The number of subject.
   @param[in] nbr_item The number of items.
   @param[in] patterns_ptr A matrix(subjects x items) of binary responses.
   @param[in] nbr_quad The number of quadrature classes.
   @param[in] quad_from The middle point of the first quadrature class.
   @param[in] quad_to The middle point of the last quadrature class.
   @param[in] model The number of parameters to used (1, 2 or 3).
   @param[in] penalized Whether to use PMMLE.
   @param[in] kernel Whether to use kernel regression.
   @param[in] smooth_factor The penalizing factor when using PMMLE.
   The bandwidth when using kernel regression.
   @param[in] slope_prior Controls whether to use a lognormal(\em slope_mean, \em slope_dev)
   prior on the slope.
   @param[in] thresh_prior Controls whether to use a normal(\em thresh_mean, \em thresh_dev)
   prior on the threshold
   @param[in] asymp_prior Controls whether to use a beta(\em asymp_mean, \em asymp_weight)
   prior on the asymptote
   @param[in] slope_mean The mean of the lognormal prior on the slope.
   @param[in] slope_dev The standard deviation of the lognormal prior on the slope.
   @param[in] thresh_mean The mean of the normal prior on the threshold.
   @param[in] thresh_dev The standard deviation of the normal prior on the threshold.
   @param[in] asymp_mean The mean of the beta prior on the asymptote.
   The alpha and beta parameter are alpha = mean * weight + 1 and beta = (1 - mean) * weight + 1).
   @param[in] asymp_weight The weight of the beta prior on the asymptote.
   The alpha and beta parameter are alpha = mean * weight + 1 and beta = (1 - mean) * weight + 1).
   @param[in] max_em_iter The maximum number of EM iterations.
   @param[in] max_nr_iter The maximum number of Newton iterations performed
   for each item at each EM iteration.
   @param[in] prec The desired precision of each parameter estimate.
   @param[in] grouping Whether to group the identical pattern together.
   @param[out] slopes_ptr A vector(items) with the estimated slope parameters.
   @param[out] slopes_stddev_ptr A vector(items) with the standard errors of the estimated slopes.
   @param[out] thresh_ptr A vector(items) with the estimated threshold parameters.
   @param[out] thresh_stddev_ptr A vector(items) with the standard errors of
   the estimated thresholds.
   @param[out] asymp_ptr A vector(items) with the estimated asymptote parameters.
   @param[out] asymp_stddev_ptr A vector(items) with the standard errors of
   the estimated asymptotes.
   @param[out] quad_points_ptr A vector(classes) with the middle points of each quadrature class.
   @param[out] quad_weights_ptr A vector(classes) with the prior weights of each quadrature class.
   @param[out] probs_ptr A matrix(items x classes) with the estimated response functions.
   @param[out] probs_stddev_ptr A matrix(items x classes) with the standard errors of
   the estimated response functions.
   @param[out] nbr_notconverge The number of item that did not converged at the last EM iteration.
   @param[out] notconverge_ptr A vector(items) indicating which item did not converged.
   @param[out] nbr_ignore The number of item that were ignored.
   @param[out] ignore_ptr A vector(items) indicating which item were ignored.
   @param[in] initialized A flag indicating if the probs are already initialized.
   @param[in] adjust_weights Controls whether to adjust the quadrature weights after each iteration.
   @param[out] em_converge A flag set to 0 if the EM algo didn't converged.
*/
void irt_wrapper(int *nbr_subject, int *nbr_item, int *patterns_ptr,
		 int *nbr_quad, double *quad_from, double *quad_to,
		 int *model, int *penalized, int *kernel, double *smooth_factor,
		 int *slope_prior, int *thresh_prior, int *asymp_prior,
		 double *slope_mean, double *slope_dev,
		 double *thresh_mean, double *thresh_dev,
		 double *asymp_mean, double *asymp_weight,
		 int *max_em_iter, int *max_nr_iter, double *precision, int *grouping,
		 double *slopes_ptr, double *thresh_ptr, double *asymp_ptr,
		 double *slopes_stddev_ptr, double *thresh_stddev_ptr,
		 double *asymp_stddev_ptr,
		 double *quad_points_ptr, double *quad_weights_ptr,
		 double *probs_ptr, double *probs_stddev_ptr,
		 int *nbr_notconverge, int *notconverge_ptr,
		 int *nbr_ignore, int *ignore_ptr, int *initialized,
		 int *adjust_weights, int *em_converge)
{
  int i, nbr_pattern;

  /* create views on the vector and matrix to be used by the other functions in the library */
  gsl_matrix_int_view patterns_view = 
    gsl_matrix_int_view_array(patterns_ptr, *nbr_subject, *nbr_item);
  gsl_matrix_int * patterns = &patterns_view.matrix;
  gsl_vector_view slopes_view =
    gsl_vector_view_array(slopes_ptr, *nbr_item);
  gsl_vector * slopes = &slopes_view.vector;
  gsl_vector_view thresh_view =
    gsl_vector_view_array(thresh_ptr, *nbr_item);
  gsl_vector * thresholds = &thresh_view.vector;
  gsl_vector_view asymp_view =
    gsl_vector_view_array(asymp_ptr, *nbr_item);
  gsl_vector * asymptotes = &asymp_view.vector;
  gsl_vector_view slopes_stddev_view =
    gsl_vector_view_array(slopes_stddev_ptr, *nbr_item);
  gsl_vector * slopes_stddev = &slopes_stddev_view.vector;
  gsl_vector_view thresh_stddev_view =
    gsl_vector_view_array(thresh_stddev_ptr, *nbr_item);
  gsl_vector * thresh_stddev = &thresh_stddev_view.vector;
  gsl_vector_view asymp_stddev_view =
    gsl_vector_view_array(asymp_stddev_ptr, *nbr_item);
  gsl_vector * asymp_stddev = &asymp_stddev_view.vector;
  gsl_vector_view quad_points_view =
    gsl_vector_view_array(quad_points_ptr, *nbr_quad);
  gsl_vector * quad_points = &quad_points_view.vector;
  gsl_vector_view quad_weights_view =
    gsl_vector_view_array(quad_weights_ptr, *nbr_quad);
  gsl_vector * quad_weights = &quad_weights_view.vector;
  gsl_matrix_view probs_view =
    gsl_matrix_view_array(probs_ptr, *nbr_item, *nbr_quad);
  gsl_matrix * probs = &probs_view.matrix;
  gsl_matrix_view probs_stddev_view =
    gsl_matrix_view_array(probs_stddev_ptr, *nbr_item, *nbr_quad);
  gsl_matrix * probs_stddev = &probs_stddev_view.matrix;
  gsl_vector_int_view ignore_view =
    gsl_vector_int_view_array(ignore_ptr, *nbr_item);
  gsl_vector_int *ignore = &ignore_view.vector;
  gsl_vector_int_view notconverge_view =
    gsl_vector_int_view_array(notconverge_ptr, *nbr_item);
  gsl_vector_int *notconverge = &notconverge_view.vector;
  gsl_matrix_int * patterns_group;
  gsl_vector * counts;
  gsl_vector_int * index;

  /* printing the libirt version */
  if (libirt_verbose > 0)
    {
      printf ("Using version %s of libirt.\n", libirt_version);
    }

  /* group the identicals patterns together */
  if (*grouping)
    {
      patterns_counts (patterns, &index, &patterns_group, &counts);
      nbr_pattern = patterns_group->size1;
      if (libirt_verbose > 0)
	{
	  printf ("%d patterns after grouping.\n", nbr_pattern);
	}
    }
  else
    {
      if (libirt_verbose > 0 && *nbr_subject / *nbr_item > 99)
	{
	  printf ("Grouping the patterns might speed up the process.\n");
	}
      patterns_group = patterns;
      nbr_pattern = patterns_group->size1;
      counts = NULL;
    }

  /* generate the quadrature points and weights */
  quadrature (*nbr_quad, *quad_from, *quad_to, quad_points, quad_weights);

  /* initialize the values */
  if (*penalized && !*initialized)
    {
      /* estimate the response functions by kernel */
      nadaraya_watson (2.7 * pow(*nbr_subject, -0.2), patterns, 
		       quad_points, quad_weights, probs, NULL);
    }
  else if(!*kernel)
    {
      gsl_vector_set_all (slopes_stddev, 0.0);
      gsl_vector_set_all (thresh_stddev, 0.0);
      gsl_vector_set_all (asymp_stddev, 0.0);
    }

  /* check for degenerate items */
  *nbr_ignore = set_ignore(patterns, probs, thresholds, ignore);
  if (libirt_verbose > 0 && *nbr_ignore > 0)
    {
      printf ("Ignoring %d items.\n", *nbr_ignore);
      if (libirt_verbose > 1)
	{
	  printf ("Ignoring items");
	  for (i = 0; i < ignore->size; i++)
	    if(gsl_vector_int_get(ignore, i))
	      printf(" %d",i+1);
	  printf (".\n");
	}
    }

  /* reset the convergence flag */
  /* it will be the number of items that didn't converged */
  *nbr_notconverge = 0;

  /* print the description of the estimation method */
  if (libirt_verbose > 0)
    {
      if(!*kernel)
	{
	  printf ("EM iterations : %d.\n", *max_em_iter);
	  printf ("Newton iterations : %d.\n", *max_nr_iter);
	  printf ("Precision : %.5lg.\n", *precision);
	}
      if (*penalized)
	printf ("Using PMMLE with a smooth factor of %.5lg.\n",
		*smooth_factor);
      else if (*kernel)
	printf ("Using kernel regression with a bandwidth of %.5lg.\n",
		*smooth_factor);
      else
	{
	  if (*model == 1)
	    printf ("Using the Rasch model with fixed slopes.\n");
	  else if (*model == 2)
	    printf ("Using the 2PLM.\n");
	  else if (*model == 3)
	    printf ("Using the 3PLM.\n");
	  else
	    printf ("This model is not implemented yet !\n");

	  if (*model > 1 && *slope_prior)
	    printf
	      ("With a lognormal prior on the slopes (mu=%.5lg, sigma=%.5lg).\n",
	       *slope_mean, *slope_dev);

	  if (*thresh_prior)
	    printf
	      ("With a normal prior on the thresholds (mu=%.5lg, sigma=%.5lg).\n",
	       *thresh_mean, *thresh_dev);

	  if (*model > 2 && *asymp_prior)
	    printf
	      ("With a beta prior on the asymptotes (p=%.5lg, m=%.5lg).\n",
	       *asymp_mean, *asymp_weight);
	}
      fflush (stdout);
    }

  if(*penalized)
    {
      *em_converge =
	em_mple_wave (*max_em_iter, *max_nr_iter, *precision, *smooth_factor,
		      patterns_group, counts, quad_points, quad_weights,
		      probs, probs_stddev, ignore, nbr_notconverge,
		      notconverge, *adjust_weights);
    }
  else if(*kernel)
    {
      nadaraya_watson (*smooth_factor, patterns,
		       quad_points, quad_weights,
		       probs, probs_stddev);
    }
  else
    {
      /* estimate the parameters */
      *em_converge = em_bme_3plm (*model, *max_em_iter, *max_nr_iter, *precision,
				  patterns_group, counts, quad_points,
				  quad_weights, *slope_prior, *slope_mean, *slope_dev,
				  *thresh_prior, *thresh_mean, *thresh_dev,
				  *asymp_prior, *asymp_mean, *asymp_weight, slopes,
				  thresholds, asymptotes, slopes_stddev,
				  thresh_stddev, asymp_stddev, ignore, 
				  nbr_notconverge, notconverge,
				  *adjust_weights);
    }

  /* print a message if some items didn't converged */
  if (libirt_verbose > 0 && *nbr_notconverge > 0)
    {
      printf ("Warning : %d items didn't converged.\n", *nbr_notconverge);
      if (libirt_verbose > 1)
	{
	  printf ("Items that didn't converged");
	  for (i = 0; i < notconverge->size; i++)
	    if(gsl_vector_int_get(notconverge, i))
	      printf(" %d",i+1);
	  printf (".\n");
	}
    }

  if (!*penalized && !*kernel)
    /* compute the response functions */
    probs_3plm (slopes, thresholds, asymptotes, quad_points, probs);

  /* free the memory */
  if (*grouping)
    {
      gsl_matrix_int_free(patterns_group);
      gsl_vector_int_free(index);
      gsl_vector_free(counts);
    }
}

/**
   \brief Compute the likelihood of each patterns.
   
   @param[in] nbr_subject The number of subject.
   @param[in] nbr_item The number of items.
   @param[in] nbr_quad The number of quadrature classes.
   @param[in] patterns A matrix (patterns x items) of binary responses.
   @param[in] probs A matrix(items x classes) with response functions.
   @param[out] like A matrix(patterns x classes) with the likelihood of each pattern
   in each quadrature class.
*/
void
likelihood_wrapper (int *nbr_subject, int *nbr_item, int *nbr_quad,
		    int * patterns_ptr, double * probs_ptr,
		    double * like_ptr)
{
  /* create views on the vector and matrix to be used by the other functions in the library */
  gsl_matrix_int_view patterns_view = 
    gsl_matrix_int_view_array(patterns_ptr, *nbr_subject, *nbr_item);
  gsl_matrix_int * patterns = &patterns_view.matrix;
  gsl_matrix_view probs_view =
    gsl_matrix_view_array(probs_ptr, *nbr_item, *nbr_quad);
  gsl_matrix * probs = &probs_view.matrix;
  gsl_matrix_view like_view =
    gsl_matrix_view_array(like_ptr, *nbr_subject, *nbr_quad);
  gsl_matrix * like = &like_view.matrix;

  /* printing the libirt version */
  if (libirt_verbose > 0)
    {
      printf ("Using version %s of libirt.\n", libirt_version);
    }

  likelihood (patterns, probs, like);
}

/**
   \brief Compute the likelihood of each multiple choice patterns.
   
   @param[in] nbr_subject The number of subject.
   @param[in] nbr_item The number of items.
   @param[in] nbr_quad The number of quadrature classes.
   @param[in] nbr_option_tot The total number of options.
   @param[in] patterns A matrix (patterns x options) of binary responses.
   @param[in] probs A matrix(items x classes) with response functions.
   @param[in] nbr_options_ptr A vector(items) with the number of option of each items.
   @param[in] items_pos_ptr A vector(items) with the position of the first option of each item
   in patterns.
   @param[out] like A matrix(patterns x classes) with the likelihood of each pattern
   in each quadrature class.
*/
void
likelihood_mc_wrapper (int *nbr_subject, int *nbr_item, int *nbr_quad, int *nbr_option_tot,
		       int * patterns_ptr, double * probs_ptr,
		       int * nbr_options_ptr, int * items_pos_ptr,
		       double * like_ptr)
{
  /* create views on the vector and matrix to be used by the other functions in the library */
  gsl_matrix_int_view patterns_view = 
    gsl_matrix_int_view_array(patterns_ptr, *nbr_subject, *nbr_option_tot);
  gsl_matrix_int * patterns = &patterns_view.matrix;
  gsl_matrix_view probs_view =
    gsl_matrix_view_array(probs_ptr, *nbr_option_tot, *nbr_quad);
  gsl_matrix * probs = &probs_view.matrix;
  gsl_matrix_view like_view =
    gsl_matrix_view_array(like_ptr, *nbr_subject, *nbr_quad);
  gsl_matrix * like = &like_view.matrix;
  gsl_vector_int_view nbr_options_view =
    gsl_vector_int_view_array(nbr_options_ptr, *nbr_item);
  gsl_vector_int * nbr_options = &nbr_options_view.vector;
  gsl_vector_int_view items_pos_view =
    gsl_vector_int_view_array(items_pos_ptr, *nbr_item);
  gsl_vector_int * items_pos = &items_pos_view.vector;

  /* printing the libirt version */
  if (libirt_verbose > 0)
    {
      printf ("Using version %s of libirt.\n", libirt_version);
    }

  likelihood_mc (patterns, probs, nbr_options, items_pos, like);
}

/**
   \brief Estimate the abilities by Bayes exact a posteriori.

   @param[in] nbr_subject The number of subject.
   @param[in] nbr_item The number of items.
   @param[in] nbr_quad The number of quadrature classes.
   @param[in] patterns A matrix(patterns x items) of binary responses.
   @param[in] probs A matrix(items x classes) with the response functions.
   @param[in] quad_points A vector(classes) with the middle points of each quadrature class.
   @param[in] quad_weights A vector(classes) with the prior weights of each quadrature class.
   @param[out] abilities A vector(patterns) with the estimated abilities.
   @param[out] abilities_stddev A vector(patterns) with the standard errors
   of the estimated abilities.

   \warning The memory for \em abilities and \em abilities_stddev should be allocated before.
*/
void eap_abilities_wrapper (int *nbr_subject, int *nbr_item, int *nbr_quad,
			    int * patterns_ptr, double * probs_ptr,
			    double * quad_points_ptr, double * quad_weights_ptr,
			    double * abilities_ptr, double * abilities_stddev_ptr)
{
  /* create views on the vector and matrix to be used by the other functions in the library */
  gsl_matrix_int_view patterns_view = 
    gsl_matrix_int_view_array(patterns_ptr, *nbr_subject, *nbr_item);
  gsl_matrix_int * patterns = &patterns_view.matrix;
  gsl_vector_view quad_weights_view =
    gsl_vector_view_array(quad_weights_ptr, *nbr_quad);
  gsl_vector * quad_weights = &quad_weights_view.vector;
  gsl_matrix_view probs_view =
    gsl_matrix_view_array(probs_ptr, *nbr_item, *nbr_quad);
  gsl_matrix * probs = &probs_view.matrix;
  gsl_vector_view quad_points_view =
    gsl_vector_view_array(quad_points_ptr, *nbr_quad);
  gsl_vector * quad_points = &quad_points_view.vector;
  gsl_vector_view abilities_view =
    gsl_vector_view_array(abilities_ptr, *nbr_subject);
  gsl_vector * abilities = &abilities_view.vector;
  gsl_vector_view abilities_stddev_view =
    gsl_vector_view_array(abilities_stddev_ptr, *nbr_subject);
  gsl_vector * abilities_stddev = &abilities_stddev_view.vector;
  gsl_matrix *post = gsl_matrix_alloc(*nbr_subject, *nbr_quad);

  /* printing the libirt version */
  if (libirt_verbose > 0)
    {
      printf ("Using version %s of libirt.\n", libirt_version);
    }

  posteriors (patterns, probs, quad_weights, post);

  eap_abilities (post, quad_points, quad_weights,
		 abilities, abilities_stddev);

  gsl_matrix_free(post);
}

/**
   \brief Estimate the abilities by Bayes exact a posteriori for multiple choice patterns.

   @param[in] nbr_subject The number of subject.
   @param[in] nbr_item The number of items.
   @param[in] nbr_quad The number of quadrature classes.
   @param[in] nbr_option_tot The total number of options.
   @param[in] item_pos_ptr A vector(items) with the position of each items in the outputs.
   @param[in] nbr_options_ptr A vector(items) with the number of option of each items.
   @param[in] patterns A matrix(patterns x options) of binary responses.
   @param[in] probs A matrix(options x classes) with the response functions.
   @param[in] quad_points A vector(classes) with the middle points of each quadrature class.
   @param[in] quad_weights A vector(classes) with the prior weights of each quadrature class.
   @param[out] abilities A vector(patterns) with the estimated abilities.
   @param[out] abilities_stddev A vector(patterns) with the standard errors
   of the estimated abilities.

   \warning The memory for \em abilities and \em abilities_stddev should be allocated before.
*/
void eap_abilities_mc_wrapper (int *nbr_subject, int *nbr_item, int *nbr_quad,
			       int *nbr_option_tot, int * items_pos_ptr, int * nbr_options_ptr,
			       int * patterns_ptr, double * probs_ptr,
			       double * quad_points_ptr, double * quad_weights_ptr,
			       double * abilities_ptr, double * abilities_stddev_ptr)
{
  /* create views on the vector and matrix to be used by the other functions in the library */
  gsl_vector_int_view items_pos_view =
    gsl_vector_int_view_array(items_pos_ptr, *nbr_item);
  gsl_vector_int * items_pos = &items_pos_view.vector;
  gsl_vector_int_view nbr_options_view =
    gsl_vector_int_view_array(nbr_options_ptr, *nbr_item);
  gsl_vector_int * nbr_options = &nbr_options_view.vector;
  gsl_matrix_int_view patterns_view = 
    gsl_matrix_int_view_array(patterns_ptr, *nbr_subject, *nbr_option_tot);
  gsl_matrix_int * patterns = &patterns_view.matrix;
  gsl_vector_view quad_weights_view =
    gsl_vector_view_array(quad_weights_ptr, *nbr_quad);
  gsl_vector * quad_weights = &quad_weights_view.vector;
  gsl_matrix_view probs_view =
    gsl_matrix_view_array(probs_ptr, *nbr_option_tot, *nbr_quad);
  gsl_matrix * probs = &probs_view.matrix;
  gsl_vector_view quad_points_view =
    gsl_vector_view_array(quad_points_ptr, *nbr_quad);
  gsl_vector * quad_points = &quad_points_view.vector;
  gsl_vector_view abilities_view =
    gsl_vector_view_array(abilities_ptr, *nbr_subject);
  gsl_vector * abilities = &abilities_view.vector;
  gsl_vector_view abilities_stddev_view =
    gsl_vector_view_array(abilities_stddev_ptr, *nbr_subject);
  gsl_vector * abilities_stddev = &abilities_stddev_view.vector;
  gsl_matrix *post = gsl_matrix_alloc(*nbr_subject, *nbr_quad);

  /* printing the libirt version */
  if (libirt_verbose > 0)
    {
      printf ("Using version %s of libirt.\n", libirt_version);
    }

  posteriors_mc (patterns, probs, nbr_options, items_pos, quad_weights, post);

  eap_abilities (post, quad_points, quad_weights,
		 abilities, abilities_stddev);

  gsl_matrix_free(post);
}

/**
   \brief Estimate the parameters of a multiple choice model by MMLE
   or more generally (if a prior is used) by BME, or estimate
   the response functions directly by PMMLE or kernel.

   @param[in] nbr_subject The number of subject.
   @param[in] nbr_item The number of items.
   @param[in] nbr_option_tot The total number of options.
   @param[out] item_pos_ptr A vector(items) with the position of each items in the outputs.
   @param[in] nbr_options_ptr A vector(items) with the number of option of each items.
   @param[in] patterns_ptr A matrix(subjects x items) of mc responses.
   @param[in] options_weights_ptr A vector(options) of options weights.
   @param[in] nbr_quad The number of quadrature classes.
   @param[in] quad_from The middle point of the first quadrature class.
   @param[in] quad_to The middle point of the last quadrature class.
   @param[in] penalized Whether to use PMMLE.
   @param[in] kernel Whether to use kernel regression.
   @param[in] smooth_factor The penalizing factor when using PMMLE.
   The bandwidth when using kernel regression.
   @param[in] graded Whether to use the graded model.
   @param[in] max_em_iter The maximum number of EM iterations.
   @param[in] max_nr_iter The maximum number of Newton iterations performed
   for each item at each EM iteration.
   @param[in] prec The desired precision of each parameter estimate.
   @param[in] grouping Whether to group the identical pattern together.
   @param[in,out] slopes_ptr A vector(options) with the estimated slope parameters.
   @param[out] slopes_stddev_ptr A vector(options) with the standard errors of the estimated slopes.
   @param[in,out] thresh_ptr A vector(options) with the estimated threshold parameters.
   @param[out] thresh_stddev_ptr A vector(options) with the standard errors of
   the estimated thresholds.
   @param[in] slope_init The initial slope.
   @param[in] thresh_init The initial threshold.
   @param[out] quad_points_ptr A vector(classes) with the middle points of each quadrature class.
   @param[out] quad_weights_ptr A vector(classes) with the prior weights of each quadrature class.
   @param[in,out] probs_ptr A matrix(options x classes) with the estimated response functions.
   @param[out] probs_stddev_ptr A matrix(options x classes) with the standard errors of
   the estimated response functions.
   @param[out] iccs_ptr A matrix(items+1 x classes) with the estimated icc and the tcc.
   @param[out] iccs_stddev_ptr A matrix(items+1 x classes) with the standard errors of
   the estimated icc and tcc.
   @param[out] nbr_notconverge The number of item that did not converged at the last EM iteration.
   @param[out] notconverge_ptr A vector(items) indicating which item did not converged.
   @param[out] nbr_ignore The number of item that were ignored.
   @param[out] ignore_ptr A vector(items) indicating which item were ignored.
   @param[in] initialized A flag indicating if the parameters or probs are already initialized.
   @param[out] patterns_exp_ptr A matrix(subjects x options) of binary responses.
   @param[in] adjust_weights Controls whether adjust the quadrature weights after each iteration.
   @param[out] em_converge A flag set to 0 if the EM algo didn't converged.
*/
void mirt_wrapper(int *nbr_subject, int *nbr_item, int *nbr_option_tot,
		  int *items_pos_ptr, int *nbr_options_ptr,
		  int *patterns_ptr, double *options_weights_ptr,
		  int *nbr_quad, double *quad_from, double *quad_to,
		  int *penalized, int *kernel, double *smooth_factor, int *graded,
		  int *max_em_iter, int *max_nr_iter, double *precision, int *grouping,
		  double *slopes_ptr, double *thresh_ptr,
		  double *slopes_stddev_ptr, double *thresh_stddev_ptr,
		  double *slope_init, double *thresh_init,
		  double *quad_points_ptr, double *quad_weights_ptr,
		  double *probs_ptr, double *probs_stddev_ptr,
		  double *iccs_ptr, double *iccs_stddev_ptr,
		  int *nbr_notconverge, int *notconverge_ptr,
		  int *nbr_ignore, int *ignore_ptr, int *initialized,
		  int *patterns_exp_ptr, int *adjust_weights, int *em_converge)
{
  int i, o, nbr_pattern, nbr_option, pos;
  double max_weight, min_weight, range_weight, weight_01, delta=0;

  /* create views on the vector and matrix to be used by the other functions in the library */
  gsl_vector_view options_weights_view =
    gsl_vector_view_array(options_weights_ptr, *nbr_option_tot);
  gsl_vector * options_weights = &options_weights_view.vector;
  gsl_vector_int_view items_pos_view =
    gsl_vector_int_view_array(items_pos_ptr, *nbr_item);
  gsl_vector_int * items_pos = &items_pos_view.vector;
  gsl_vector_int_view nbr_options_view =
    gsl_vector_int_view_array(nbr_options_ptr, *nbr_item);
  gsl_vector_int * nbr_options = &nbr_options_view.vector;
  gsl_matrix_int_view patterns_view = 
    gsl_matrix_int_view_array(patterns_ptr, *nbr_subject, *nbr_item);
  gsl_matrix_int * patterns = &patterns_view.matrix;
  gsl_matrix_int_view patterns_exp_view = 
    gsl_matrix_int_view_array(patterns_exp_ptr, *nbr_subject, *nbr_option_tot);
  gsl_matrix_int * patterns_exp = &patterns_exp_view.matrix;
  gsl_vector_view slopes_view =
    gsl_vector_view_array(slopes_ptr, *graded?(*nbr_item):(*nbr_option_tot));
  gsl_vector * slopes = &slopes_view.vector;
  gsl_vector_view thresh_view =
    gsl_vector_view_array(thresh_ptr, *nbr_option_tot);
  gsl_vector * thresholds = &thresh_view.vector;
  gsl_vector_view slopes_stddev_view =
    gsl_vector_view_array(slopes_stddev_ptr, *graded?(*nbr_item):(*nbr_option_tot));
  gsl_vector * slopes_stddev = &slopes_stddev_view.vector;
  gsl_vector_view thresh_stddev_view =
    gsl_vector_view_array(thresh_stddev_ptr, *nbr_option_tot);
  gsl_vector * thresh_stddev = &thresh_stddev_view.vector;
  gsl_vector_view quad_points_view =
    gsl_vector_view_array(quad_points_ptr, *nbr_quad);
  gsl_vector * quad_points = &quad_points_view.vector;
  gsl_vector_view quad_weights_view =
    gsl_vector_view_array(quad_weights_ptr, *nbr_quad);
  gsl_vector * quad_weights = &quad_weights_view.vector;
  gsl_matrix_view probs_view =
    gsl_matrix_view_array(probs_ptr, *nbr_option_tot, *nbr_quad);
  gsl_matrix * probs = &probs_view.matrix;
  gsl_matrix_view probs_stddev_view =
    gsl_matrix_view_array(probs_stddev_ptr, *nbr_option_tot, *nbr_quad);
  gsl_matrix * probs_stddev = &probs_stddev_view.matrix;
  gsl_matrix_view iccs_view =
    gsl_matrix_view_array(iccs_ptr, *nbr_item+1, *nbr_quad);
  gsl_matrix * iccs = &iccs_view.matrix;
  gsl_matrix_view iccs_stddev_view =
    gsl_matrix_view_array(iccs_stddev_ptr, *nbr_item+1, *nbr_quad);
  gsl_matrix * iccs_stddev = &iccs_stddev_view.matrix;
  gsl_vector_int_view ignore_view =
    gsl_vector_int_view_array(ignore_ptr, *nbr_item);
  gsl_vector_int *ignore = &ignore_view.vector;
  gsl_vector_int_view notconverge_view =
    gsl_vector_int_view_array(notconverge_ptr, *nbr_item);
  gsl_vector_int *notconverge = &notconverge_view.vector;
  gsl_matrix_int * patterns_group;
  gsl_vector * counts;
  gsl_vector_int * index;

  /* printing the libirt version */
  if (libirt_verbose > 0)
    {
      printf ("Using version %s of libirt.\n", libirt_version);
    }

  /* expand the options of each items */
  patterns_expand (patterns, nbr_options, patterns_exp, items_pos);

  /* group the identicals patterns together */
  if (*grouping)
    {
      patterns_counts (patterns_exp, &index, &patterns_group, &counts);
      nbr_pattern = patterns_group->size1;
      if (libirt_verbose > 0)
	{
	  printf ("%d patterns after grouping.\n", nbr_pattern);
	}
    }
  else
    {
      if (libirt_verbose > 0 && *nbr_subject / *nbr_item > 99)
	{
	  printf ("Grouping the patterns might speed up the process.\n");
	}
      patterns_group = patterns_exp;
      nbr_pattern = patterns_group->size1;
      counts = NULL;
    }

  /* generate the quadrature points and weights */
  quadrature (*nbr_quad, *quad_from, *quad_to, quad_points, quad_weights);

  /* set the initials values */
  if (*penalized && !*initialized)
    {
      /* estimate the response functions by kernel */
      nadaraya_watson_mc (2.7 * pow(*nbr_subject, -0.2), patterns_exp,
			  items_pos, nbr_options, 
			  options_weights, quad_points, quad_weights,
			  probs, NULL);
    }
  else if(!*kernel && !*initialized)
    {
      for (i = 0; i < *nbr_item; i++)
	{
	  pos = gsl_vector_int_get(items_pos, i);
	  nbr_option = gsl_vector_int_get(nbr_options, i);
	  min_weight = 1e10;
	  max_weight = 0;

	  for (o = 0; o < nbr_option; o++)
	    {
	      if(gsl_vector_get(options_weights,pos+o) < min_weight)
		min_weight = gsl_vector_get(options_weights,pos+o);
	      if(gsl_vector_get(options_weights,pos+o) > max_weight)
		max_weight = gsl_vector_get(options_weights,pos+o);
	    }
	  range_weight = max_weight - min_weight;
	  if(range_weight == 0) range_weight = 1.0;

	  if(*graded)
	    {
	      if(nbr_option > 2) delta = *thresh_init/(nbr_option-2.0);
	      gsl_vector_set(thresholds, pos, -*thresh_init);
	      gsl_vector_set(thresholds, pos+nbr_option-1, *thresh_init);
	      for (o = 1; o < nbr_option-1; o++)
		{
		  gsl_vector_set(thresholds, pos+o, -*thresh_init+(2*o-1)*delta);
		}
	      gsl_vector_set(slopes, i, *slope_init);
	    }
	  else
	    {
	      for (o = 0; o < nbr_option; o++)
		{
		  weight_01 = (gsl_vector_get(options_weights,pos+o) - min_weight)
		    / range_weight + (0.01 * o) / nbr_option;
		  gsl_vector_set(thresholds, pos+o, *thresh_init*(2*weight_01-1));
		  gsl_vector_set(slopes, pos+o, *slope_init*(2*weight_01-1));
		}
	    }
	}
    }

  /* check for degenerate items */
  *nbr_ignore = set_ignore_mc (patterns_exp, nbr_options, items_pos,
			      probs, thresholds, ignore);
  if (libirt_verbose > 0 && *nbr_ignore > 0)
    {
      printf ("Ignoring %d items.\n", *nbr_ignore);
      if (libirt_verbose > 1)
	{
	  printf ("Ignoring items");
	  for (i = 0; i < ignore->size; i++)
	    if(gsl_vector_int_get(ignore, i))
	      printf(" %d",i+1);
	  printf (".\n");
	}
    }

  /* reset the convergence flag */
  /* it will be the number of items that didn't converged */
  *nbr_notconverge = 0;

  /* print the description of the estimation method */
  if (libirt_verbose > 0)
    {
      if(!*kernel)
	{
	  printf ("EM iterations : %d.\n", *max_em_iter);
	  printf ("Newton iterations : %d.\n", *max_nr_iter);
	  printf ("Precision : %.5lg.\n", *precision);
	}
      if (*penalized)
	printf ("Using PMMLE with a smooth factor of %.5lg.\n",
		*smooth_factor);
      else if (*kernel)
	printf ("Using kernel regression with a bandwidth of %.5lg.\n",
		*smooth_factor);
      else
	printf ("Using the multivariate logistic model.\n");
      fflush (stdout);
    }

  /* estimate the parameters */
  if(*penalized)
    {
      *em_converge =
	em_mple_wave_mc (*max_em_iter, *max_nr_iter, *precision, *smooth_factor,
			 patterns_group, counts, quad_points, quad_weights,
			 items_pos, nbr_options, 
			 probs, probs_stddev, ignore, 
			 nbr_notconverge, notconverge, *adjust_weights);
    }
  else if(*kernel)
    {
      nadaraya_watson_mc (*smooth_factor, patterns_exp, 
			  items_pos, nbr_options, 
			  options_weights, quad_points, quad_weights, 
			  probs, probs_stddev);
    }
  else if(*graded)
    {
      *em_converge = mmle_2plm_grad (*max_em_iter, *max_nr_iter, *precision,
				     patterns_group, counts, quad_points, quad_weights, 
				     items_pos, nbr_options,
				     thresholds, thresh_stddev, slopes, slopes_stddev,
				     ignore, nbr_notconverge, notconverge, *adjust_weights);
    }
  else
    {
      *em_converge = mmle_2plm_mc (*max_em_iter, *max_nr_iter, *precision,
				   patterns_group, counts, quad_points, quad_weights, 
				   items_pos, nbr_options,
				   thresholds, thresh_stddev, slopes, slopes_stddev,
				   ignore, nbr_notconverge, notconverge, *adjust_weights);
    }

  /* print a message if some items didn't converged */
  if (libirt_verbose > 0 && *nbr_notconverge > 0)
    {
      printf ("Warning : %d items didn't converged.\n", *nbr_notconverge);
      if (libirt_verbose > 1)
	{
	  printf ("Items that didn't converged");
	  for (i = 0; i < notconverge->size; i++)
	    if(gsl_vector_int_get(notconverge, i))
	      printf(" %d",i+1);
	  printf (".\n");
	}
    }

  /* compute the response functions */
  if (!*penalized && !*kernel)
    {
      if(*graded) probs_2plm_grad (slopes, thresholds,
				   nbr_options, items_pos,
				   quad_points, probs, NULL);
      else probs_2plm_mc (slopes, thresholds, nbr_options,
			  items_pos, quad_points, probs);
    }

  /* compute the icc */
  icc_from_probs (probs, probs_stddev, options_weights,
		  nbr_options, items_pos,
		  iccs, iccs_stddev);

  /* free the memory */
  if (*grouping)
    {
      gsl_matrix_int_free(patterns_group);
      gsl_vector_int_free(index);
      gsl_vector_free(counts);
    }
}

/**
   \brief Estimate the parameters of a multiple choice model by MMLE
   or more generally (if a prior is used) by BME, or estimate
   the response functions directly by PMMLE or kernel.

   @param[in] nbr_subject The number of subject.
   @param[in] nbr_item The number of items.
   @param[in] nbr_option_tot The total number of options.
   @param[in] items_models_ptr A vector(items) with the model of each items.
   @param[in] items_methods_ptr A vector(items) with the estimation methods of each items.
   @param[out] item_pos_ptr A vector(items) with the position of each items in the outputs.
   @param[in] nbr_options_ptr A vector(items) with the number of option of each items.
   @param[in] patterns_ptr A matrix(subjects x items) of mc responses.
   @param[in] options_weights_ptr A vector(options) of options weights.
   @param[in] nbr_quad The number of quadrature classes.
   @param[in] quad_from The middle point of the first quadrature class.
   @param[in] quad_to The middle point of the last quadrature class.
   @param[in] smooth_factor The penalizing factor when using PMMLE.
   The bandwidth when using kernel regression.
   @param[in] max_em_iter The maximum number of EM iterations.
   @param[in] max_nr_iter The maximum number of Newton iterations performed
   for each item at each EM iteration.
   @param[in] prec The desired precision of each parameter estimate.
   @param[in] grouping Whether to group the identical pattern together.
   @param[in,out] slopes_ptr A vector(options) with the estimated slope parameters.
   @param[out] slopes_stddev_ptr A vector(options) with the standard errors of the estimated slopes.
   @param[in,out] thresh_ptr A vector(options) with the estimated threshold parameters.
   @param[out] thresh_stddev_ptr A vector(options) with the standard errors of
   the estimated thresholds.
   @param[in] slope_init The initial slope.
   @param[in] thresh_init The initial threshold.
   @param[out] quad_points_ptr A vector(classes) with the middle points of each quadrature class.
   @param[out] quad_weights_ptr A vector(classes) with the prior weights of each quadrature class.
   @param[in,out] probs_ptr A matrix(options x classes) with the estimated response functions.
   @param[out] probs_stddev_ptr A matrix(options x classes) with the standard errors of
   the estimated response functions.
   @param[out] iccs_ptr A matrix(items+1 x classes) with the estimated icc and the tcc.
   @param[out] iccs_stddev_ptr A matrix(items+1 x classes) with the standard errors of
   the estimated icc and tcc.
   @param[out] nbr_notconverge The number of item that did not converged at the last EM iteration.
   @param[out] notconverge_ptr A vector(items) indicating which item did not converged.
   @param[out] nbr_ignore The number of item that were ignored.
   @param[out] ignore_ptr A vector(items) indicating which item were ignored.
   @param[in] initialized A flag indicating if the parameters or probs are already initialized.
   @param[out] patterns_exp_ptr A matrix(subjects x options) of binary responses.
   @param[out] em_converge A flag set to 0 if the EM algo didn't converged.
*/
void mixed_wrapper(int *nbr_subject, int *nbr_item, int *nbr_option_tot,
		   int *items_models_ptr, int *items_methods_ptr,
		   int *items_pos_ptr, int *nbr_options_ptr,
		   int *patterns_ptr, double *options_weights_ptr,
		   int *nbr_quad, double *quad_from, double *quad_to,
		   double *smooth_factor,
		   int *max_em_iter, int *max_nr_iter, double *precision, int *grouping,
		   double *slopes_ptr, double *thresh_ptr,
		   double *slopes_stddev_ptr, double *thresh_stddev_ptr,
		   double *slope_init, double *thresh_init,
		   double *quad_points_ptr, double *quad_weights_ptr,
		   double *probs_ptr, double *probs_stddev_ptr,
		   double *iccs_ptr, double *iccs_stddev_ptr,
		   int *nbr_notconverge, int *notconverge_ptr,
		   int *nbr_ignore, int *ignore_ptr, int *initialized,
		   int *patterns_exp_ptr, int *em_converge)
{
  int i, o, nbr_pattern, nbr_option, pos, model;
  double max_weight, min_weight, range_weight, weight_01, delta=0;

  /* create views on the vector and matrix to be used by the other functions in the library */
  gsl_vector_view options_weights_view =
    gsl_vector_view_array(options_weights_ptr, *nbr_option_tot);
  gsl_vector * options_weights = &options_weights_view.vector;
  gsl_vector_int_view items_models_view =
    gsl_vector_int_view_array(items_models_ptr, *nbr_item);
  gsl_vector_int * items_models = &items_models_view.vector;
  gsl_vector_int_view items_methods_view =
    gsl_vector_int_view_array(items_methods_ptr, *nbr_item);
  gsl_vector_int * items_methods = &items_methods_view.vector;
  gsl_vector_int_view items_pos_view =
    gsl_vector_int_view_array(items_pos_ptr, *nbr_item);
  gsl_vector_int * items_pos = &items_pos_view.vector;
  gsl_vector_int_view nbr_options_view =
    gsl_vector_int_view_array(nbr_options_ptr, *nbr_item);
  gsl_vector_int * nbr_options = &nbr_options_view.vector;
  gsl_matrix_int_view patterns_view = 
    gsl_matrix_int_view_array(patterns_ptr, *nbr_subject, *nbr_item);
  gsl_matrix_int * patterns = &patterns_view.matrix;
  gsl_matrix_int_view patterns_exp_view = 
    gsl_matrix_int_view_array(patterns_exp_ptr, *nbr_subject, *nbr_option_tot);
  gsl_matrix_int * patterns_exp = &patterns_exp_view.matrix;
  gsl_vector_view slopes_view =
    gsl_vector_view_array(slopes_ptr, *nbr_option_tot);
  gsl_vector * slopes = &slopes_view.vector;
  gsl_vector_view thresh_view =
    gsl_vector_view_array(thresh_ptr, *nbr_option_tot);
  gsl_vector * thresholds = &thresh_view.vector;
  gsl_vector_view slopes_stddev_view =
    gsl_vector_view_array(slopes_stddev_ptr, *nbr_option_tot);
  gsl_vector * slopes_stddev = &slopes_stddev_view.vector;
  gsl_vector_view thresh_stddev_view =
    gsl_vector_view_array(thresh_stddev_ptr, *nbr_option_tot);
  gsl_vector * thresh_stddev = &thresh_stddev_view.vector;
  gsl_vector_view quad_points_view =
    gsl_vector_view_array(quad_points_ptr, *nbr_quad);
  gsl_vector * quad_points = &quad_points_view.vector;
  gsl_vector_view quad_weights_view =
    gsl_vector_view_array(quad_weights_ptr, *nbr_quad);
  gsl_vector * quad_weights = &quad_weights_view.vector;
  gsl_matrix_view probs_view =
    gsl_matrix_view_array(probs_ptr, *nbr_option_tot, *nbr_quad);
  gsl_matrix * probs = &probs_view.matrix;
  gsl_matrix_view probs_stddev_view =
    gsl_matrix_view_array(probs_stddev_ptr, *nbr_option_tot, *nbr_quad);
  gsl_matrix * probs_stddev = &probs_stddev_view.matrix;
  gsl_matrix_view iccs_view =
    gsl_matrix_view_array(iccs_ptr, *nbr_item+1, *nbr_quad);
  gsl_matrix * iccs = &iccs_view.matrix;
  gsl_matrix_view iccs_stddev_view =
    gsl_matrix_view_array(iccs_stddev_ptr, *nbr_item+1, *nbr_quad);
  gsl_matrix * iccs_stddev = &iccs_stddev_view.matrix;
  gsl_vector_int_view ignore_view =
    gsl_vector_int_view_array(ignore_ptr, *nbr_item);
  gsl_vector_int *ignore = &ignore_view.vector;
  gsl_vector_int_view notconverge_view =
    gsl_vector_int_view_array(notconverge_ptr, *nbr_item);
  gsl_vector_int *notconverge = &notconverge_view.vector;
  gsl_matrix_int * patterns_group;
  gsl_vector * counts;
  gsl_vector_int * index;

  /* printing the libirt version */
  if (libirt_verbose > 0)
    {
      printf ("Using version %s of libirt.\n", libirt_version);
    }

  /* expand the options of each items */
  patterns_expand (patterns, nbr_options, patterns_exp, items_pos);

  /* group the identicals patterns together */
  if (*grouping)
    {
      patterns_counts (patterns_exp, &index, &patterns_group, &counts);
      nbr_pattern = patterns_group->size1;
      if (libirt_verbose > 0)
	{
	  printf ("%d patterns after grouping.\n", nbr_pattern);
	}
    }
  else
    {
      if (libirt_verbose > 0 && *nbr_subject / *nbr_item > 99)
	{
	  printf ("Grouping the patterns might speed up the process.\n");
	}
      patterns_group = patterns_exp;
      nbr_pattern = patterns_group->size1;
      counts = NULL;
    }

  /* generate the quadrature points and weights */
  quadrature (*nbr_quad, *quad_from, *quad_to, quad_points, quad_weights);

  /* set the initials values */
  if(!*initialized)
    {
      for (i = 0; i < *nbr_item; i++)
	{
	  model = gsl_vector_int_get(items_models, i);
	  pos = gsl_vector_int_get(items_pos, i);
	  nbr_option = gsl_vector_int_get(nbr_options, i);
	  min_weight = 1e10;
	  max_weight = 0;

	  for (o = 0; o < nbr_option; o++)
	    {
	      if(gsl_vector_get(options_weights,pos+o) < min_weight)
		min_weight = gsl_vector_get(options_weights,pos+o);
	      if(gsl_vector_get(options_weights,pos+o) > max_weight)
		max_weight = gsl_vector_get(options_weights,pos+o);
	    }
	  range_weight = max_weight - min_weight;
	  if(range_weight == 0) range_weight = 1.0;

	  switch(model) {
	  case CODE_NOMINAL:
	    for (o = 0; o < nbr_option; o++)
	      {
		weight_01 = (gsl_vector_get(options_weights,pos+o) - min_weight)
		  / range_weight + (0.01 * o) / nbr_option;
		gsl_vector_set(thresholds, pos+o, *thresh_init*(2*weight_01-1));
		gsl_vector_set(slopes, pos+o, *slope_init*(2*weight_01-1));
	      }
	    break;
	  case CODE_GRADED:
	    if(nbr_option > 2) delta = *thresh_init/(nbr_option-2.0);
	    gsl_vector_set(thresholds, pos, -*thresh_init);
	    gsl_vector_set(thresholds, pos+nbr_option-1, *thresh_init);
	    for (o = 1; o < nbr_option-1; o++)
	      {
		gsl_vector_set(thresholds, pos+o, -*thresh_init+(2*o-1)*delta);
		gsl_vector_set(slopes, pos+o, *slope_init);
	      }
	    break;
	  }
	}
    }

  /* check for degenerate items */
  *nbr_ignore = set_ignore_mc (patterns_exp, nbr_options, items_pos,
			      probs, thresholds, ignore);
  if (libirt_verbose > 0 && *nbr_ignore > 0)
    {
      printf ("Ignoring %d items.\n", *nbr_ignore);
      if (libirt_verbose > 1)
	{
	  printf ("Ignoring items");
	  for (i = 0; i < ignore->size; i++)
	    if(gsl_vector_int_get(ignore, i))
	      printf(" %d",i+1);
	  printf (".\n");
	}
    }

  /* reset the convergence flag */
  /* it will be the number of items that didn't converged */
  *nbr_notconverge = 0;

  /* print the description of the estimation method */
  if (libirt_verbose > 0)
    {
      printf ("EM iterations : %d.\n", *max_em_iter);
      printf ("Newton iterations : %d.\n", *max_nr_iter);
      printf ("Precision : %.5lg.\n", *precision);
      printf ("Using the multivariate logistic model.\n");
      fflush (stdout);
    }

  /* estimate the parameters */
  *em_converge = mmle_mixed (*max_em_iter, *max_nr_iter, *precision,
			     patterns_group, counts, quad_points, quad_weights, 
			     items_models, items_pos, nbr_options,
			     thresholds, thresh_stddev, slopes, slopes_stddev,
			     ignore, nbr_notconverge, notconverge);

  /* print a message if some items didn't converged */
  if (libirt_verbose > 0 && *nbr_notconverge > 0)
    {
      printf ("Warning : %d items didn't converged.\n", *nbr_notconverge);
      if (libirt_verbose > 1)
	{
	  printf ("Items that didn't converged");
	  for (i = 0; i < notconverge->size; i++)
	    if(gsl_vector_int_get(notconverge, i))
	      printf(" %d",i+1);
	  printf (".\n");
	}
    }

  /* compute the response functions */
  probs_2plm_mixed (slopes, thresholds, items_models, nbr_options, items_pos, quad_points, probs, NULL);

  /* compute the icc */
  icc_from_probs (probs, probs_stddev, options_weights,
		  nbr_options, items_pos,
		  iccs, iccs_stddev);

  /* free the memory */
  if (*grouping)
    {
      gsl_matrix_int_free(patterns_group);
      gsl_vector_int_free(index);
      gsl_vector_free(counts);
    }
}

/**
   \brief Compute the items' information curves from the response functions.

   @param[in] nbr_item The number of items.
   @param[in] nbr_quad The number of quadrature classes.
   @param[in] probs_ptr A matrix (items x classes) of response functions.
   @param[in] quad_points_ptr A vector(classes) with the middle points of each class.
   @param[out] infos_ptr A matrix (items x classes) of items' information curves.
   @param[out] test_info_ptr A vector (classes) with test information curve.

   \warning The memory for the outputs must be allocated before.
*/
void info_from_probs_wrapper (int *nbr_item, int *nbr_quad,
			      double *probs_ptr, double *quad_points_ptr,
			      double *infos_ptr, double *test_info_ptr)
{
  gsl_matrix_view probs_view =
    gsl_matrix_view_array(probs_ptr, *nbr_item, *nbr_quad);
  gsl_matrix * probs = &probs_view.matrix;
  gsl_vector_view quad_points_view =
    gsl_vector_view_array(quad_points_ptr, *nbr_quad);
  gsl_vector * quad_points = &quad_points_view.vector;
  gsl_vector_view test_info_view =
    gsl_vector_view_array(test_info_ptr, *nbr_quad);
  gsl_vector * test_info = &test_info_view.vector;
  gsl_matrix_view infos_view =
    gsl_matrix_view_array(infos_ptr, *nbr_item, *nbr_quad);
  gsl_matrix * infos = &infos_view.matrix;

  info_from_probs (probs, quad_points, infos, test_info);
}

/**
   \brief Compute the options' information curves from the response functions
   for a multivariate model.

   @param[in] nbr_item The number of items.
   @param[in] nbr_option_tot The total number of options.
   @param[in] nbr_quad The number of quadrature classes.
   @param[in] probs_ptr A matrix (options x classes) of response functions.
   @param[in] quad_points_ptr A vector(classes) with the middle points of each class.
   @param[in] nbr_options_ptr A vector(items) with the number of option of each items.
   @param[in] items_pos_ptr A vector(items) with the position of the first option of each item.
   @param[out] options_infos_ptr A matrix (options x classes) of options' information curves.
   @param[out] infos_ptr A matrix (items x classes) of items' information curves.
   @param[out] test_info_ptr A vector (classes) with test information curve.

   \warning The memory for the outputs must be allocated before.
*/
void info_from_probs_mc_wrapper (int *nbr_item, int *nbr_option_tot, int *nbr_quad,
				 double *probs_ptr, double *quad_points_ptr,
				 int *nbr_options_ptr, int *items_pos_ptr,
				 double *options_infos_ptr, double *infos_ptr, 
				 double *test_info_ptr)
{
  gsl_matrix_view probs_view =
    gsl_matrix_view_array(probs_ptr, *nbr_option_tot, *nbr_quad);
  gsl_matrix * probs = &probs_view.matrix;
  gsl_vector_view quad_points_view =
    gsl_vector_view_array(quad_points_ptr, *nbr_quad);
  gsl_vector * quad_points = &quad_points_view.vector;
  gsl_vector_view test_info_view =
    gsl_vector_view_array(test_info_ptr, *nbr_quad);
  gsl_vector * test_info = &test_info_view.vector;
  gsl_vector_int_view items_pos_view =
    gsl_vector_int_view_array(items_pos_ptr, *nbr_item);
  gsl_vector_int * items_pos = &items_pos_view.vector;
  gsl_vector_int_view nbr_options_view =
    gsl_vector_int_view_array(nbr_options_ptr, *nbr_item);
  gsl_vector_int * nbr_options = &nbr_options_view.vector;
  gsl_matrix_view options_infos_view =
    gsl_matrix_view_array(options_infos_ptr, *nbr_option_tot, *nbr_quad);
  gsl_matrix * options_infos = &options_infos_view.matrix;
  gsl_matrix_view infos_view =
    gsl_matrix_view_array(infos_ptr, *nbr_item, *nbr_quad);
  gsl_matrix * infos = &infos_view.matrix;

  info_from_probs_mc (probs, quad_points, nbr_options, items_pos,
		      options_infos, infos, test_info);
}

/**
   \brief Compute the log likelihood ratio statistics of each items.

   @param[in] nbr_item The number of item.
   @param[in] nbr_par The number of quadrature class.
   @param[in] patterns A matrix(patterns x items) of binary responses.
   @param[in] quad_weights_ptr A vector(classes) with the normal weights
   of each class.
   @param[in] probs_ptr A matrix(items x classes) with response functions.
   @param[in] nbr_inter The number of intervals to use.
   @param[out] chi2_ptr A vector(items+1) with the statistics of each items
   and for the overall fit.
   @param[out] df A vector(items+1) with the degrees of freedom.
   @param[out] p_value_ptr A vector(items+1) with the p-values.
   
   \warning The memory for the outputs should be allocated before.
*/
void
llk_ratio_fit_test_wrapper (int *nbr_subject, int *nbr_item, int *nbr_quad,
			    int *patterns_ptr,
			    double *quad_weights_ptr, double *probs_ptr,
			    int *nbr_inter, double *chi2_ptr, 
			    int *df_ptr, double *p_value_ptr)
{
  gsl_vector * quad_sizes = gsl_vector_alloc(*nbr_quad);
  gsl_matrix * quad_freqs = gsl_matrix_alloc(*nbr_item, *nbr_quad);
  gsl_matrix * post = gsl_matrix_alloc(*nbr_subject, *nbr_quad);
  gsl_matrix_int_view patterns_view = 
    gsl_matrix_int_view_array(patterns_ptr, *nbr_subject, *nbr_item);
  gsl_matrix_int * patterns = &patterns_view.matrix;
  gsl_vector_view quad_weights_view =
    gsl_vector_view_array(quad_weights_ptr, *nbr_quad);
  gsl_vector * quad_weights = &quad_weights_view.vector;
  gsl_matrix_view probs_view =
    gsl_matrix_view_array(probs_ptr, *nbr_item, *nbr_quad);
  gsl_matrix * probs = &probs_view.matrix;
  gsl_vector_view chi2_view =
    gsl_vector_view_array(chi2_ptr, *nbr_item+1);
  gsl_vector * chi2 = &chi2_view.vector;
  gsl_vector_int_view df_view =
    gsl_vector_int_view_array(df_ptr, *nbr_item+1);
  gsl_vector_int * df = &df_view.vector;
  gsl_vector_view p_value_view =
    gsl_vector_view_array(p_value_ptr, *nbr_item+1);
  gsl_vector * p_value = &p_value_view.vector;

  posteriors (patterns, probs, quad_weights, post);

  frequencies (patterns, NULL, post, probs, quad_sizes, quad_freqs);

  llk_ratio_fit_test (quad_sizes, quad_freqs, quad_weights, probs,
		      *nbr_inter, chi2, df, p_value);
}

/**
   \brief Compute the log likelihood ratio statistics of mc each items.

   @param[in] nbr_subject The number of subject.
   @param[in] nbr_item The number of item.
   @param[in] nbr_par The number of quadrature class.
   @param[in] nbr_option_tot The total number of options.
   @param[in] patterns A matrix(patterns x options) of binary responses.
   @param[in] quad_weights_ptr A vector(classes) with the normal weights
   of each class.
   @param[in] probs_ptr A matrix(items x classes) with response functions.
   @param[in] nbr_options_ptr A vector(items) with the number of option of each items.
   @param[in] items_pos_ptr A vector(items) with the position of the first option of each item.
   @param[in] nbr_inter The number of intervals to use.
   @param[out] chi2_ptr A vector(items+1) with the statistics of each items
   and for the overall fit.
   @param[out] df A vector(items+1) with the degrees of freedom.
   @param[out] p_value_ptr A vector(items+1) with the p-values.
   
   \warning The memory for the outputs should be allocated before.
*/
void
llk_ratio_fit_test_mc_wrapper (int *nbr_subject, int *nbr_item, int *nbr_quad,
			       int *nbr_option_tot, int *patterns_ptr,
			       double *quad_weights_ptr, double *probs_ptr,
			       int * nbr_options_ptr, int * items_pos_ptr,
			       int *nbr_inter, double *chi2_ptr, 
			       int *df_ptr, double *p_value_ptr)
{
  gsl_vector * quad_sizes = gsl_vector_alloc(*nbr_quad);
  gsl_matrix * quad_freqs = gsl_matrix_alloc(*nbr_option_tot, *nbr_quad);
  gsl_matrix * post = gsl_matrix_alloc(*nbr_subject, *nbr_quad);
  gsl_matrix_int_view patterns_view = 
    gsl_matrix_int_view_array(patterns_ptr, *nbr_subject, *nbr_option_tot);
  gsl_matrix_int * patterns = &patterns_view.matrix;
  gsl_vector_view quad_weights_view =
    gsl_vector_view_array(quad_weights_ptr, *nbr_quad);
  gsl_vector * quad_weights = &quad_weights_view.vector;
  gsl_matrix_view probs_view =
    gsl_matrix_view_array(probs_ptr, *nbr_option_tot, *nbr_quad);
  gsl_matrix * probs = &probs_view.matrix;
  gsl_vector_view chi2_view =
    gsl_vector_view_array(chi2_ptr, *nbr_item+1);
  gsl_vector * chi2 = &chi2_view.vector;
  gsl_vector_int_view df_view =
    gsl_vector_int_view_array(df_ptr, *nbr_item+1);
  gsl_vector_int * df = &df_view.vector;
  gsl_vector_view p_value_view =
    gsl_vector_view_array(p_value_ptr, *nbr_item+1);
  gsl_vector * p_value = &p_value_view.vector;
  gsl_vector_int_view nbr_options_view =
    gsl_vector_int_view_array(nbr_options_ptr, *nbr_item);
  gsl_vector_int * nbr_options = &nbr_options_view.vector;
  gsl_vector_int_view items_pos_view =
    gsl_vector_int_view_array(items_pos_ptr, *nbr_item);
  gsl_vector_int * items_pos = &items_pos_view.vector;

  posteriors_mc (patterns, probs, nbr_options, items_pos, 
		 quad_weights, post);

  frequencies (patterns, NULL, post, probs, quad_sizes, quad_freqs);

  llk_ratio_fit_test_mc (quad_sizes, quad_freqs, quad_weights, probs,
			 nbr_options, items_pos, *nbr_inter, chi2, df, 
			 p_value);
}

/**
   \brief Compute the classical test theory statistics.

   @param[in] nbr_subject The number of subject.
   @param[in] nbr_item The number of item.
   @param[in] patterns_ptr A matrix(subjects x items) with the binary patterns.
   @param[out] items_mean_ptr A vector(items) with the means of each items.
   @param[out] items_sd_ptr A vector(items) with the sd of each items.
   @param[out] items_corr_ptr A vector(items) with the correlations 
   between each items and the score of each subject without the item.
   @param[out] items_bis_corr_ptr A vector(items) with the biserial correlations 
   between each items and the score of each subject without the item.
   @param[out] items_nbr_ptr A vector(items) with the number of non missing response
   for each items.
   @param[out] subjects_score_ptr A vector(subjects) with the scores of each subjects.
   @param[out] subjects_nbr_ptr A vector(subjects) with the number of non missing
   response of each subjects.
   @param[out] nbr The total number of non missing response.
   @param[out] mean The score mean.
   @param[out] sd The score sd.
   @param[out] alpha The fiability coefficient alpha.
   @param[out] pairs_corr A matrix(items x items) of correlation for each pair of items

   \warning The memory for the outputs should be allocated before.
*/
void 
classical_statistics_wrapper(int * nbr_subject, int * nbr_item, int * patterns_ptr, 
			     double * items_mean_ptr, double * items_sd_ptr, 
			     double * items_corr_ptr, double * items_bis_corr_ptr, 
			     int * items_nbr_ptr,
			     double * subjects_score_ptr, int * subjects_nbr_ptr,
			     int * nbr, double * mean, double * sd, double * alpha,
			     double * pairs_corr_ptr)
{
  gsl_matrix_view pairs_corr_view = 
    gsl_matrix_view_array(pairs_corr_ptr, *nbr_item, *nbr_item);
  gsl_matrix * pairs_corr = &pairs_corr_view.matrix;

  gsl_matrix_int_view patterns_view = 
    gsl_matrix_int_view_array(patterns_ptr, *nbr_subject, *nbr_item);
  gsl_matrix_int * patterns = &patterns_view.matrix;

  gsl_vector_view items_mean_view =
    gsl_vector_view_array(items_mean_ptr, *nbr_item);
  gsl_vector * items_mean = &items_mean_view.vector;

  gsl_vector_view items_sd_view =
    gsl_vector_view_array(items_sd_ptr, *nbr_item);
  gsl_vector * items_sd = &items_sd_view.vector;

  gsl_vector_view items_corr_view =
    gsl_vector_view_array(items_corr_ptr, *nbr_item);
  gsl_vector * items_corr = &items_corr_view.vector;

  gsl_vector_view items_bis_corr_view =
    gsl_vector_view_array(items_bis_corr_ptr, *nbr_item);
  gsl_vector * items_bis_corr = &items_bis_corr_view.vector;

  gsl_vector_int_view items_nbr_view =
    gsl_vector_int_view_array(items_nbr_ptr, *nbr_item);
  gsl_vector_int * items_nbr = &items_nbr_view.vector;

  gsl_vector_view subjects_score_view =
    gsl_vector_view_array(subjects_score_ptr, *nbr_subject);
  gsl_vector * subjects_score = &subjects_score_view.vector;

  gsl_vector_int_view subjects_nbr_view =
    gsl_vector_int_view_array(subjects_nbr_ptr, *nbr_subject);
  gsl_vector_int * subjects_nbr = &subjects_nbr_view.vector;

  classical_statistics(patterns, items_mean, items_sd, items_corr, items_bis_corr, items_nbr,
		       subjects_score, subjects_nbr, nbr, mean, sd, alpha, pairs_corr);
}

/**
   \brief Compute the classical test theory statistics.

   @param[in] nbr_option_tot The total number of options.
   @param[in] nbr_subject The number of subject.
   @param[in] nbr_item The number of item.
   @param[in] patterns_ptr A matrix(subjects x items) with the binary patterns.
   @param[in] options_weights A vector(options) of weights to compute the score.
   @param[in] items_pos A vector(items) with the position of the first option of each item
   in patterns (and probs).
   @param[in] nbr_options A vector(items) with the number of option of each item
   in patterns (and probs).
   @param[out] items_mean_ptr A vector(items) with the means of each items.
   @param[out] items_sd_ptr A vector(items) with the sd of each items.
   @param[out] items_corr_ptr A vector(items) with the correlations 
   between each items and the score of each subject without the item.
   @param[out] items_poly_corr_ptr A vector(items) with the polyserial correlations 
   between each items and the score of each subject without the item.
   @param[out] items_nbr_ptr A vector(items) with the number of non missing response
   for each items.
   @param[out] subjects_score_ptr A vector(subjects) with the scores of each subjects.
   @param[out] subjects_nbr_ptr A vector(subjects) with the number of non missing
   response of each subjects.
   @param[out] nbr The total number of non missing response.
   @param[out] mean The score mean.
   @param[out] sd The score sd.
   @param[out] alpha The fiability coefficient alpha.
   @param[out] pairs_corr A matrix(items x items) of correlation for each pair of items

   \warning The memory for the outputs should be allocated before.
*/
void 
classical_statistics_mc_wrapper(int * nbr_option_tot, 
				int * nbr_subject, int * nbr_item, int * patterns_ptr, 
				double * options_weights_ptr,
				int * items_pos_ptr, int * nbr_options_ptr,
				double * items_mean_ptr, double * items_sd_ptr, 
				double * items_corr_ptr, double * items_poly_corr_ptr, 
				int * items_nbr_ptr,
				double * subjects_score_ptr, int * subjects_nbr_ptr,
				int * nbr, double * mean, double * sd, double * alpha,
				double * pairs_corr_ptr)
{
  gsl_matrix_view pairs_corr_view = 
    gsl_matrix_view_array(pairs_corr_ptr, *nbr_item, *nbr_item);
  gsl_matrix * pairs_corr = &pairs_corr_view.matrix;
  gsl_matrix_int_view patterns_view = 
    gsl_matrix_int_view_array(patterns_ptr, *nbr_subject, *nbr_item);
  gsl_matrix_int * patterns = &patterns_view.matrix;
  gsl_vector_view items_mean_view =
    gsl_vector_view_array(items_mean_ptr, *nbr_item);
  gsl_vector * items_mean = &items_mean_view.vector;
  gsl_vector_view items_sd_view =
    gsl_vector_view_array(items_sd_ptr, *nbr_item);
  gsl_vector * items_sd = &items_sd_view.vector;
  gsl_vector_view items_corr_view =
    gsl_vector_view_array(items_corr_ptr, *nbr_item);
  gsl_vector * items_corr = &items_corr_view.vector;
  gsl_vector_view items_poly_corr_view =
    gsl_vector_view_array(items_poly_corr_ptr, *nbr_item);
  gsl_vector * items_poly_corr = &items_poly_corr_view.vector;
  gsl_vector_int_view items_nbr_view =
    gsl_vector_int_view_array(items_nbr_ptr, *nbr_item);
  gsl_vector_int * items_nbr = &items_nbr_view.vector;
  gsl_vector_view subjects_score_view =
    gsl_vector_view_array(subjects_score_ptr, *nbr_subject);
  gsl_vector * subjects_score = &subjects_score_view.vector;
  gsl_vector_int_view subjects_nbr_view =
    gsl_vector_int_view_array(subjects_nbr_ptr, *nbr_subject);
  gsl_vector_int * subjects_nbr = &subjects_nbr_view.vector;
  gsl_vector_view options_weights_view =
    gsl_vector_view_array(options_weights_ptr, *nbr_option_tot);
  gsl_vector * options_weights = &options_weights_view.vector;
  gsl_vector_int_view items_pos_view =
    gsl_vector_int_view_array(items_pos_ptr, *nbr_item);
  gsl_vector_int * items_pos = &items_pos_view.vector;
  gsl_vector_int_view nbr_options_view =
    gsl_vector_int_view_array(nbr_options_ptr, *nbr_item);
  gsl_vector_int * nbr_options = &nbr_options_view.vector;

  classical_statistics_mc(patterns,  options_weights, items_pos, nbr_options,
			  items_mean, items_sd, items_corr, items_poly_corr, items_nbr,
			  subjects_score, subjects_nbr, nbr, mean, sd, alpha, pairs_corr);
}

/**
   \brief Compute the log likelihood ratio statistics of each items pairs
   for local dependance.

   @param[in] nbr_subject The number of subject.
   @param[in] nbr_item The number of item.
   @param[in] nbr_quad The number of quadrature classes.
   @param[in] patterns_ptr A matrix(subjects x items) with binary response.
   @param[in] probs_ptr A matrix(items x classes) with response functions.
   @param[in] quad_weights_ptr A vector(classes) with the normal weights of each class.
   @param[out] chi2_ptr A matrix(items x items) with the statistics of each items pairs.
   @param[out] df_ptr A matrix(items x items) with the degrees of freedom of each items pairs.
   @param[out] p_value_ptr A matrix(items x items) with the p-values of each items pairs.

   \warning The memory for the outputs should be allocated before.
*/
void
llk_ratio_ld_test_wrapper (int * nbr_subject, int * nbr_item, int * nbr_quad,
			   int * patterns_ptr, double * probs_ptr,
			   double * quad_weights_ptr,
			   double * chi2_ptr, int * df_ptr, double * p_value_ptr)
{
  gsl_matrix_int_view patterns_view = 
    gsl_matrix_int_view_array(patterns_ptr, *nbr_subject, *nbr_item);
  gsl_matrix_int * patterns = &patterns_view.matrix;
  gsl_vector_view quad_weights_view =
    gsl_vector_view_array(quad_weights_ptr, *nbr_quad);
  gsl_vector * quad_weights = &quad_weights_view.vector;
  gsl_matrix_view probs_view =
    gsl_matrix_view_array(probs_ptr, *nbr_item, *nbr_quad);
  gsl_matrix * probs = &probs_view.matrix;
  gsl_matrix_view chi2_view =
    gsl_matrix_view_array(chi2_ptr, *nbr_item, *nbr_item);
  gsl_matrix * chi2 = &chi2_view.matrix;
  gsl_matrix_int_view df_view =
    gsl_matrix_int_view_array(df_ptr, *nbr_item, *nbr_item);
  gsl_matrix_int * df = &df_view.matrix;
  gsl_matrix_view p_value_view =
    gsl_matrix_view_array(p_value_ptr, *nbr_item, *nbr_item);
  gsl_matrix * p_value = &p_value_view.matrix;

  llk_ratio_ld_test (patterns, probs,quad_weights, chi2, df, p_value);
}

/**
   \brief Compute the log likelihood ratio statistics of each items pairs
   for local dependance.

   @param[in] nbr_subject The number of subject.
   @param[in] nbr_item The number of item.
   @param[in] nbr_quad The number of quadrature classes.
   @param[in] nbr_option_tot The total number of options.
   @param[in] patterns_ptr A matrix(subjects x items) with binary response.
   @param[in] probs_ptr A matrix(probs x classes) with response functions.
   @param[in] quad_weights_ptr A vector(classes) with the normal weights of each class.
   @param[in] nbr_options_ptr A vector(items) with the number of option of each items.
   @param[in] items_pos_ptr A vector(items) with the position of the first option of each item.
   @param[out] chi2_ptr A matrix(items x items) with the statistics of each items pairs.
   @param[out] df_ptr A matrix(items x items) with the degrees of freedom of each items pairs.
   @param[out] p_value_ptr A matrix(items x items) with the p-values of each items pairs.

   \warning The memory for the outputs should be allocated before.
*/
void
llk_ratio_ld_test_mc_wrapper (int * nbr_subject, int * nbr_item, int * nbr_quad,
			      int * nbr_option_tot, int * patterns_ptr, 
			      double * probs_ptr, double * quad_weights_ptr, 
			      int *nbr_options_ptr, int *items_pos_ptr,
			      double * chi2_ptr, int *df_ptr, double * p_value_ptr)
{
  gsl_matrix_int_view patterns_view = 
    gsl_matrix_int_view_array(patterns_ptr, *nbr_subject, *nbr_item);
  gsl_matrix_int * patterns = &patterns_view.matrix;
  gsl_vector_view quad_weights_view =
    gsl_vector_view_array(quad_weights_ptr, *nbr_quad);
  gsl_vector * quad_weights = &quad_weights_view.vector;
  gsl_matrix_view probs_view =
    gsl_matrix_view_array(probs_ptr, *nbr_option_tot, *nbr_quad);
  gsl_matrix * probs = &probs_view.matrix;
  gsl_matrix_view chi2_view =
    gsl_matrix_view_array(chi2_ptr, *nbr_item, *nbr_item);
  gsl_matrix * chi2 = &chi2_view.matrix;
  gsl_matrix_int_view df_view =
    gsl_matrix_int_view_array(df_ptr, *nbr_item, *nbr_item);
  gsl_matrix_int * df = &df_view.matrix;
  gsl_matrix_view p_value_view =
    gsl_matrix_view_array(p_value_ptr, *nbr_item, *nbr_item);
  gsl_matrix * p_value = &p_value_view.matrix;
  gsl_vector_int_view items_pos_view =
    gsl_vector_int_view_array(items_pos_ptr, *nbr_item);
  gsl_vector_int * items_pos = &items_pos_view.vector;
  gsl_vector_int_view nbr_options_view =
    gsl_vector_int_view_array(nbr_options_ptr, *nbr_item);
  gsl_vector_int * nbr_options = &nbr_options_view.vector;

  llk_ratio_ld_test_mc (patterns, probs, quad_weights, 
			nbr_options, items_pos,
			chi2, df, p_value);
}

/**
   \brief Estimate the abilities by WMLE.

   @param[in] max_iter The maximum number of Newton iterations performed.
   @param[in] prec The desired precision.
   @param[in] nbr_subject The number of subject.
   @param[in] nbr_item The number of items.
   @param[in] nbr_quad The number of quadrature classes.
   @param[in] patterns A matrix(patterns x items) of binary responses.
   @param[in] probs A matrix(items x classes) with the response functions.
   @param[in] quad_points A vector(classes) with the middle points of each quadrature class.
   @param[out] abilities A vector(patterns) with the estimated abilities.
   @param[out] abilities_stddev A vector(patterns) with the standard errors
   of the estimated abilities.

   \warning The memory for \em abilities and \em abilities_stddev should be allocated before.
*/
void wmle_abilities_wrapper (int *max_iter, double *prec, 
			     int *nbr_subject, int *nbr_item, int *nbr_quad,
			     int * patterns_ptr, double * probs_ptr,
			     double * quad_points_ptr,
			     double * abilities_ptr, double * abilities_stddev_ptr)
{
  /* create views on the vector and matrix to be used by the other functions in the library */
  gsl_matrix_int_view patterns_view = 
    gsl_matrix_int_view_array(patterns_ptr, *nbr_subject, *nbr_item);
  gsl_matrix_int * patterns = &patterns_view.matrix;
  gsl_matrix_view probs_view =
    gsl_matrix_view_array(probs_ptr, *nbr_item, *nbr_quad);
  gsl_matrix * probs = &probs_view.matrix;
  gsl_vector_view quad_points_view =
    gsl_vector_view_array(quad_points_ptr, *nbr_quad);
  gsl_vector * quad_points = &quad_points_view.vector;
  gsl_vector_view abilities_view =
    gsl_vector_view_array(abilities_ptr, *nbr_subject);
  gsl_vector * abilities = &abilities_view.vector;
  gsl_vector_view abilities_stddev_view =
    gsl_vector_view_array(abilities_stddev_ptr, *nbr_subject);
  gsl_vector * abilities_stddev = &abilities_stddev_view.vector;

  gsl_matrix *like = gsl_matrix_alloc(*nbr_subject, *nbr_quad);
  gsl_matrix *infos = gsl_matrix_alloc(*nbr_item, *nbr_quad);
  gsl_vector *test_info = gsl_vector_alloc(*nbr_quad);

  /* printing the libirt version */
  if (libirt_verbose > 0)
    {
      printf ("Using version %s of libirt.\n", libirt_version);
    }

  likelihood (patterns, probs, like);

  info_from_probs (probs, quad_points, infos, test_info);

  wmle_abilities (*max_iter, *prec, like, test_info, quad_points,
		  abilities, abilities_stddev);

  gsl_matrix_free(like);
  gsl_matrix_free(infos);
  gsl_vector_free(test_info);
}

/**
   \brief Estimate the abilities by WMLE for multiple choice patterns.

   @param[in] max_iter The maximum number of Newton iterations performed.
   @param[in] prec The desired precision.
   @param[in] nbr_subject The number of subject.
   @param[in] nbr_item The number of items.
   @param[in] nbr_quad The number of quadrature classes.
   @param[in] nbr_option_tot The total number of options.
   @param[in] nbr_options_ptr A vector(items) with the number of option of each items.
   @param[in] items_pos_ptr A vector(items) with the position of the first option of each item
   in patterns.
   @param[in] patterns A matrix(patterns x options) of binary responses.
   @param[in] probs A matrix(items x classes) with the response functions.
   @param[in] quad_points A vector(classes) with the middle points of each quadrature class.
   @param[out] abilities A vector(patterns) with the estimated abilities.
   @param[out] abilities_stddev A vector(patterns) with the standard errors
   of the estimated abilities.

   \warning The memory for \em abilities and \em abilities_stddev should be allocated before.
*/
void wmle_abilities_mc_wrapper (int *max_iter, double *prec, 
				int *nbr_subject, int *nbr_item, int *nbr_quad,
				int *nbr_option_tot, int * nbr_options_ptr, int * items_pos_ptr,
				int * patterns_ptr, double * probs_ptr, double * quad_points_ptr,
				double * abilities_ptr, double * abilities_stddev_ptr)
{
  /* create views on the vector and matrix to be used by the other functions in the library */
  gsl_matrix_int_view patterns_view = 
    gsl_matrix_int_view_array(patterns_ptr, *nbr_subject, *nbr_option_tot);
  gsl_matrix_int * patterns = &patterns_view.matrix;
  gsl_matrix_view probs_view =
    gsl_matrix_view_array(probs_ptr, *nbr_option_tot, *nbr_quad);
  gsl_matrix * probs = &probs_view.matrix;
  gsl_vector_view quad_points_view =
    gsl_vector_view_array(quad_points_ptr, *nbr_quad);
  gsl_vector * quad_points = &quad_points_view.vector;
  gsl_vector_view abilities_view =
    gsl_vector_view_array(abilities_ptr, *nbr_subject);
  gsl_vector * abilities = &abilities_view.vector;
  gsl_vector_view abilities_stddev_view =
    gsl_vector_view_array(abilities_stddev_ptr, *nbr_subject);
  gsl_vector * abilities_stddev = &abilities_stddev_view.vector;
  gsl_vector_int_view nbr_options_view =
    gsl_vector_int_view_array(nbr_options_ptr, *nbr_item);
  gsl_vector_int * nbr_options = &nbr_options_view.vector;
  gsl_vector_int_view items_pos_view =
    gsl_vector_int_view_array(items_pos_ptr, *nbr_item);
  gsl_vector_int * items_pos = &items_pos_view.vector;

  gsl_matrix *like = gsl_matrix_alloc(*nbr_subject, *nbr_quad);
  gsl_matrix *options_infos = gsl_matrix_alloc(*nbr_option_tot, *nbr_quad);
  gsl_matrix *infos = gsl_matrix_alloc(*nbr_item, *nbr_quad);
  gsl_vector *test_info = gsl_vector_alloc(*nbr_quad);

  /* printing the libirt version */
  if (libirt_verbose > 0)
    {
      printf ("Using version %s of libirt.\n", libirt_version);
    }

  likelihood_mc (patterns, probs, nbr_options, items_pos, like);

  info_from_probs_mc (probs, quad_points, nbr_options, items_pos,
		      options_infos, infos, test_info);

  wmle_abilities (*max_iter, *prec, like, test_info, quad_points,
		  abilities, abilities_stddev);

  gsl_matrix_free(like);
  gsl_matrix_free(options_infos);
  gsl_matrix_free(infos);
  gsl_vector_free(test_info);
}
