text:
/*
MPCOTool:
The Multi-Purposes Calibration and Optimization Tool. A software to perform
calibrations or optimizations of empirical parameters.

AUTHORS: Javier Burguete and Borja Latorre.

Copyright 2012-2019, AUTHORS.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

    1. Redistributions of source code must retain the above copyright notice,
        this list of conditions and the following disclaimer.

    2. Redistributions in binary form must reproduce the above copyright notice,
        this list of conditions and the following disclaimer in the
        documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY AUTHORS ``AS IS'' AND ANY EXPRESS OR IMPLIED
WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
OF SUCH DAMAGE.
*/

/**
 * \file optimize.c
 * \brief Source file to define the optimization functions.
 * \authors Javier Burguete and Borja Latorre.
 * \copyright Copyright 2012-2019, all rights reserved.
 */
#define _GNU_SOURCE
#include "config.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/param.h>
#include <gsl/gsl_rng.h>
#include <libxml/parser.h>
#include <libintl.h>
#include <glib.h>
#include <glib/gstdio.h>
#include <json-glib/json-glib.h>
#ifdef G_OS_WIN32
#include <windows.h>
#elif !defined(__BSD_VISIBLE) && !defined(NetBSD)
#include <alloca.h>
#endif
#if HAVE_MPI
#include <mpi.h>
#endif
#include "genetic/genetic.h"
#include "utils.h"
#include "experiment.h"
#include "variable.h"
#include "input.h"
#include "optimize.h"

#define DEBUG_OPTIMIZE 0        ///< Macro to debug optimize functions.

/**
 * \def RM
 * \brief Macro to define the shell remove command.
 */
#ifdef G_OS_WIN32
#define RM "del"
#else
#define RM "rm"
#endif

Optimize optimize[1];           ///< Optimization data.
unsigned int nthreads_climbing;
///< Number of threads for the hill climbing method.

static void (*optimize_algorithm) ();
///< Pointer to the function to perform a optimization algorithm step.
static double (*optimize_estimate_climbing) (unsigned int variable,
                                             unsigned int estimate);
///< Pointer to the function to estimate the climbing.
static double (*optimize_norm) (unsigned int simulation);
///< Pointer to the error norm function.

/**
 * Function to write the simulation input file.
 */
static inline void
optimize_input (unsigned int simulation,        ///< Simulation number.
                char *input,    ///< Input file name.
                GMappedFile * stencil)  ///< Template of the input file name.
{
  char buffer[32], value[32];
  GRegex *regex;
  FILE *file;
  char *buffer2, *buffer3 = NULL, *content;
  gsize length;
  unsigned int i;

#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_input: start\n");
#endif

  // Checking the file
  if (!stencil)
    goto optimize_input_end;

  // Opening stencil
  content = g_mapped_file_get_contents (stencil);
  length = g_mapped_file_get_length (stencil);
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_input: length=%lu\ncontent:\n%s", length, content);
#endif
  file = g_fopen (input, "w");

  // Parsing stencil
  for (i = 0; i < optimize->nvariables; ++i)
    {
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_input: variable=%u\n", i);
#endif
      snprintf (buffer, 32, "@variable%u@", i + 1);
      regex = g_regex_new (buffer, (GRegexCompileFlags) 0, (GRegexMatchFlags) 0,
                           NULL);
      if (i == 0)
        {
          buffer2 = g_regex_replace_literal (regex, content, length, 0,
                                             optimize->label[i],
                                             (GRegexMatchFlags) 0, NULL);
#if DEBUG_OPTIMIZE
          fprintf (stderr, "optimize_input: buffer2\n%s", buffer2);
#endif
        }
      else
        {
          length = strlen (buffer3);
          buffer2 = g_regex_replace_literal (regex, buffer3, length, 0,
                                             optimize->label[i],
                                             (GRegexMatchFlags) 0, NULL);
          g_free (buffer3);
        }
      g_regex_unref (regex);
      length = strlen (buffer2);
      snprintf (buffer, 32, "@value%u@", i + 1);
      regex = g_regex_new (buffer, (GRegexCompileFlags) 0, (GRegexMatchFlags) 0,
                           NULL);
      snprintf (value, 32, format[optimize->precision[i]],
                optimize->value[simulation * optimize->nvariables + i]);

#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_input: value=%s\n", value);
#endif
      buffer3 = g_regex_replace_literal (regex, buffer2, length, 0, value,
                                         (GRegexMatchFlags) 0, NULL);
      g_free (buffer2);
      g_regex_unref (regex);
    }

  // Saving input file
  fwrite (buffer3, strlen (buffer3), sizeof (char), file);
  g_free (buffer3);
  fclose (file);

optimize_input_end:
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_input: end\n");
#endif
  return;
}

/**
 * Function to parse input files, simulating and calculating the objective
 *   function.
 *
 * \return Objective function value.
 */
static double
optimize_parse (unsigned int simulation,        ///< Simulation number.
                unsigned int experiment)        ///< Experiment number.
{
  unsigned int i;
  double e;
  char buffer[512], input[MAX_NINPUTS][32], output[32], result[32], *buffer2,
    *buffer3, *buffer4;
  FILE *file_result;

#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_parse: start\n");
  fprintf (stderr, "optimize_parse: simulation=%u experiment=%u\n",
           simulation, experiment);
#endif

  // Opening input files
  for (i = 0; i < optimize->ninputs; ++i)
    {
      snprintf (&input[i][0], 32, "input-%u-%u-%u", i, simulation, experiment);
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_parse: i=%u input=%s\n", i, &input[i][0]);
#endif
      optimize_input (simulation, &input[i][0], optimize->file[i][experiment]);
    }
  for (; i < MAX_NINPUTS; ++i)
    strcpy (&input[i][0], "");
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_parse: parsing end\n");
#endif

  // Performing the simulation
  snprintf (output, 32, "output-%u-%u", simulation, experiment);
  buffer2 = g_path_get_dirname (optimize->simulator);
  buffer3 = g_path_get_basename (optimize->simulator);
  buffer4 = g_build_filename (buffer2, buffer3, NULL);
  snprintf (buffer, 512, "\"%s\" %s %s %s %s %s %s %s %s %s",
            buffer4, input[0], input[1], input[2], input[3], input[4],
            input[5], input[6], input[7], output);
  g_free (buffer4);
  g_free (buffer3);
  g_free (buffer2);
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_parse: %s\n", buffer);
#endif
  if (system (buffer) == -1)
    error_message = buffer;

  // Checking the objective value function
  if (optimize->evaluator)
    {
      snprintf (result, 32, "result-%u-%u", simulation, experiment);
      buffer2 = g_path_get_dirname (optimize->evaluator);
      buffer3 = g_path_get_basename (optimize->evaluator);
      buffer4 = g_build_filename (buffer2, buffer3, NULL);
      snprintf (buffer, 512, "\"%s\" %s %s %s",
                buffer4, output, optimize->experiment[experiment], result);
      g_free (buffer4);
      g_free (buffer3);
      g_free (buffer2);
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_parse: %s\n", buffer);
      fprintf (stderr, "optimize_parse: result=%s\n", result);
#endif
      if (system (buffer) == -1)
        error_message = buffer;
      file_result = g_fopen (result, "r");
      e = atof (fgets (buffer, 512, file_result));
      fclose (file_result);
    }
  else
    {
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_parse: output=%s\n", output);
#endif
      strcpy (result, "");
      file_result = g_fopen (output, "r");
      e = atof (fgets (buffer, 512, file_result));
      fclose (file_result);
    }

  // Removing files
#if !DEBUG_OPTIMIZE
  for (i = 0; i < optimize->ninputs; ++i)
    {
      if (optimize->file[i][0])
        {
          snprintf (buffer, 512, RM " %s", &input[i][0]);
          if (system (buffer) == -1)
            error_message = buffer;
        }
    }
  snprintf (buffer, 512, RM " %s %s", output, result);
  if (system (buffer) == -1)
    error_message = buffer;
#endif

  // Processing pending events
  if (show_pending)
    show_pending ();

#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_parse: end\n");
#endif

  // Returning the objective function
  return e * optimize->weight[experiment];
}

/**
 * Function to calculate the Euclidian error norm.
 *
 * \return Euclidian error norm.
 */
static double
optimize_norm_euclidian (unsigned int simulation)       ///< simulation number.
{
  double e, ei;
  unsigned int i;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_norm_euclidian: start\n");
#endif
  e = 0.;
  for (i = 0; i < optimize->nexperiments; ++i)
    {
      ei = optimize_parse (simulation, i);
      e += ei * ei;
    }
  e = sqrt (e);
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_norm_euclidian: error=%lg\n", e);
  fprintf (stderr, "optimize_norm_euclidian: end\n");
#endif
  return e;
}

/**
 * Function to calculate the maximum error norm.
 *
 * \return Maximum error norm.
 */
static double
optimize_norm_maximum (unsigned int simulation) ///< simulation number.
{
  double e, ei;
  unsigned int i;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_norm_maximum: start\n");
#endif
  e = 0.;
  for (i = 0; i < optimize->nexperiments; ++i)
    {
      ei = fabs (optimize_parse (simulation, i));
      e = fmax (e, ei);
    }
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_norm_maximum: error=%lg\n", e);
  fprintf (stderr, "optimize_norm_maximum: end\n");
#endif
  return e;
}

/**
 * Function to calculate the P error norm.
 *
 * \return P error norm.
 */
static double
optimize_norm_p (unsigned int simulation)       ///< simulation number.
{
  double e, ei;
  unsigned int i;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_norm_p: start\n");
#endif
  e = 0.;
  for (i = 0; i < optimize->nexperiments; ++i)
    {
      ei = fabs (optimize_parse (simulation, i));
      e += pow (ei, optimize->p);
    }
  e = pow (e, 1. / optimize->p);
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_norm_p: error=%lg\n", e);
  fprintf (stderr, "optimize_norm_p: end\n");
#endif
  return e;
}

/**
 * Function to calculate the taxicab error norm.
 *
 * \return Taxicab error norm.
 */
static double
optimize_norm_taxicab (unsigned int simulation) ///< simulation number.
{
  double e;
  unsigned int i;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_norm_taxicab: start\n");
#endif
  e = 0.;
  for (i = 0; i < optimize->nexperiments; ++i)
    e += fabs (optimize_parse (simulation, i));
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_norm_taxicab: error=%lg\n", e);
  fprintf (stderr, "optimize_norm_taxicab: end\n");
#endif
  return e;
}

/**
 * Function to print the results.
 */
static void
optimize_print ()
{
  unsigned int i;
  char buffer[512];
#if HAVE_MPI
  if (optimize->mpi_rank)
    return;
#endif
  printf ("%s\n", _("Best result"));
  fprintf (optimize->file_result, "%s\n", _("Best result"));
  printf ("error = %.15le\n", optimize->error_old[0]);
  fprintf (optimize->file_result, "error = %.15le\n", optimize->error_old[0]);
  for (i = 0; i < optimize->nvariables; ++i)
    {
      snprintf (buffer, 512, "%s = %s\n",
                optimize->label[i], format[optimize->precision[i]]);
      printf (buffer, optimize->value_old[i]);
      fprintf (optimize->file_result, buffer, optimize->value_old[i]);
    }
  fflush (optimize->file_result);
}

/**
 * Function to save in a file the variables and the error.
 */
static void
optimize_save_variables (unsigned int simulation,       ///< Simulation number.
                         double error)  ///< Error value.
{
  unsigned int i;
  char buffer[64];
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_save_variables: start\n");
#endif
  for (i = 0; i < optimize->nvariables; ++i)
    {
      snprintf (buffer, 64, "%s ", format[optimize->precision[i]]);
      fprintf (optimize->file_variables, buffer,
               optimize->value[simulation * optimize->nvariables + i]);
    }
  fprintf (optimize->file_variables, "%.14le\n", error);
  fflush (optimize->file_variables);
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_save_variables: end\n");
#endif
}

/**
 * Function to save the best simulations.
 */
static void
optimize_best (unsigned int simulation, ///< Simulation number.
               double value)    ///< Objective function value.
{
  unsigned int i, j;
  double e;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_best: start\n");
  fprintf (stderr, "optimize_best: nsaveds=%u nbest=%u\n",
           optimize->nsaveds, optimize->nbest);
#endif
  if (optimize->nsaveds < optimize->nbest
      || value < optimize->error_best[optimize->nsaveds - 1])
    {
      if (optimize->nsaveds < optimize->nbest)
        ++optimize->nsaveds;
      optimize->error_best[optimize->nsaveds - 1] = value;
      optimize->simulation_best[optimize->nsaveds - 1] = simulation;
      for (i = optimize->nsaveds; --i;)
        {
          if (optimize->error_best[i] < optimize->error_best[i - 1])
            {
              j = optimize->simulation_best[i];
              e = optimize->error_best[i];
              optimize->simulation_best[i] = optimize->simulation_best[i - 1];
              optimize->error_best[i] = optimize->error_best[i - 1];
              optimize->simulation_best[i - 1] = j;
              optimize->error_best[i - 1] = e;
            }
          else
            break;
        }
    }
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_best: end\n");
#endif
}

/**
 * Function to optimize sequentially.
 */
static void
optimize_sequential ()
{
  unsigned int i;
  double e;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_sequential: start\n");
  fprintf (stderr, "optimize_sequential: nstart=%u nend=%u\n",
           optimize->nstart, optimize->nend);
#endif
  for (i = optimize->nstart; i < optimize->nend; ++i)
    {
      e = optimize_norm (i);
      optimize_best (i, e);
      optimize_save_variables (i, e);
      if (e < optimize->threshold)
        {
          optimize->stop = 1;
          break;
        }
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_sequential: i=%u e=%lg\n", i, e);
#endif
    }
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_sequential: end\n");
#endif
}

/**
 * Function to optimize on a thread.
 *
 * \return NULL.
 */
static void *
optimize_thread (ParallelData * data)   ///< Function data.
{
  unsigned int i, thread;
  double e;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_thread: start\n");
#endif
  thread = data->thread;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_thread: thread=%u start=%u end=%u\n", thread,
           optimize->thread[thread], optimize->thread[thread + 1]);
#endif
  for (i = optimize->thread[thread]; i < optimize->thread[thread + 1]; ++i)
    {
      e = optimize_norm (i);
      g_mutex_lock (mutex);
      optimize_best (i, e);
      optimize_save_variables (i, e);
      if (e < optimize->threshold)
        optimize->stop = 1;
      g_mutex_unlock (mutex);
      if (optimize->stop)
        break;
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_thread: i=%u e=%lg\n", i, e);
#endif
    }
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_thread: end\n");
#endif
  g_thread_exit (NULL);
  return NULL;
}

/**
 * Function to merge the 2 optimization results.
 */
static inline void
optimize_merge (unsigned int nsaveds,   ///< Number of saved results.
                unsigned int *simulation_best,
                ///< Array of best simulation numbers.
                double *error_best)
                ///< Array of best objective function values.
{
  unsigned int i, j, k, s[optimize->nbest];
  double e[optimize->nbest];
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_merge: start\n");
#endif
  i = j = k = 0;
  do
    {
      if (i == optimize->nsaveds)
        {
          s[k] = simulation_best[j];
          e[k] = error_best[j];
          ++j;
          ++k;
          if (j == nsaveds)
            break;
        }
      else if (j == nsaveds)
        {
          s[k] = optimize->simulation_best[i];
          e[k] = optimize->error_best[i];
          ++i;
          ++k;
          if (i == optimize->nsaveds)
            break;
        }
      else if (optimize->error_best[i] > error_best[j])
        {
          s[k] = simulation_best[j];
          e[k] = error_best[j];
          ++j;
          ++k;
        }
      else
        {
          s[k] = optimize->simulation_best[i];
          e[k] = optimize->error_best[i];
          ++i;
          ++k;
        }
    }
  while (k < optimize->nbest);
  optimize->nsaveds = k;
  memcpy (optimize->simulation_best, s, k * sizeof (unsigned int));
  memcpy (optimize->error_best, e, k * sizeof (double));
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_merge: end\n");
#endif
}

/**
 * Function to synchronise the optimization results of MPI tasks.
 */
#if HAVE_MPI
static void
optimize_synchronise ()
{
  unsigned int i, nsaveds, simulation_best[optimize->nbest], stop;
  double error_best[optimize->nbest];
  MPI_Status mpi_stat;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_synchronise: start\n");
#endif
  if (optimize->mpi_rank == 0)
    {
      for (i = 1; (int) i < ntasks; ++i)
        {
          MPI_Recv (&nsaveds, 1, MPI_INT, i, 1, MPI_COMM_WORLD, &mpi_stat);
          MPI_Recv (simulation_best, nsaveds, MPI_INT, i, 1,
                    MPI_COMM_WORLD, &mpi_stat);
          MPI_Recv (error_best, nsaveds, MPI_DOUBLE, i, 1,
                    MPI_COMM_WORLD, &mpi_stat);
          optimize_merge (nsaveds, simulation_best, error_best);
          MPI_Recv (&stop, 1, MPI_UNSIGNED, i, 1, MPI_COMM_WORLD, &mpi_stat);
          if (stop)
            optimize->stop = 1;
        }
      for (i = 1; (int) i < ntasks; ++i)
        MPI_Send (&optimize->stop, 1, MPI_UNSIGNED, i, 1, MPI_COMM_WORLD);
    }
  else
    {
      MPI_Send (&optimize->nsaveds, 1, MPI_INT, 0, 1, MPI_COMM_WORLD);
      MPI_Send (optimize->simulation_best, optimize->nsaveds, MPI_INT, 0, 1,
                MPI_COMM_WORLD);
      MPI_Send (optimize->error_best, optimize->nsaveds, MPI_DOUBLE, 0, 1,
                MPI_COMM_WORLD);
      MPI_Send (&optimize->stop, 1, MPI_UNSIGNED, 0, 1, MPI_COMM_WORLD);
      MPI_Recv (&stop, 1, MPI_UNSIGNED, 0, 1, MPI_COMM_WORLD, &mpi_stat);
      if (stop)
        optimize->stop = 1;
    }
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_synchronise: end\n");
#endif
}
#endif

/**
 * Function to optimize with the sweep algorithm.
 */
static void
optimize_sweep ()
{
  unsigned int i, j, k, l;
  double e;
  GThread *thread[nthreads];
  ParallelData data[nthreads];
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_sweep: start\n");
#endif
  for (i = 0; i < optimize->nsimulations; ++i)
    {
      k = i;
      for (j = 0; j < optimize->nvariables; ++j)
        {
          l = k % optimize->nsweeps[j];
          k /= optimize->nsweeps[j];
          e = optimize->rangemin[j];
          if (optimize->nsweeps[j] > 1)
            e += l * (optimize->rangemax[j] - optimize->rangemin[j])
              / (optimize->nsweeps[j] - 1);
          optimize->value[i * optimize->nvariables + j] = e;
        }
    }
  optimize->nsaveds = 0;
  if (nthreads <= 1)
    optimize_sequential ();
  else
    {
      for (i = 0; i < nthreads; ++i)
        {
          data[i].thread = i;
          thread[i]
            = g_thread_new (NULL, (GThreadFunc) optimize_thread, &data[i]);
        }
      for (i = 0; i < nthreads; ++i)
        g_thread_join (thread[i]);
    }
#if HAVE_MPI
  // Communicating tasks results
  optimize_synchronise ();
#endif
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_sweep: end\n");
#endif
}

/**
 * Function to optimize with the Monte-Carlo algorithm.
 */
static void
optimize_MonteCarlo ()
{
  unsigned int i, j;
  GThread *thread[nthreads];
  ParallelData data[nthreads];
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_MonteCarlo: start\n");
#endif
  for (i = 0; i < optimize->nsimulations; ++i)
    for (j = 0; j < optimize->nvariables; ++j)
      optimize->value[i * optimize->nvariables + j]
        = optimize->rangemin[j] + gsl_rng_uniform (optimize->rng)
        * (optimize->rangemax[j] - optimize->rangemin[j]);
  optimize->nsaveds = 0;
  if (nthreads <= 1)
    optimize_sequential ();
  else
    {
      for (i = 0; i < nthreads; ++i)
        {
          data[i].thread = i;
          thread[i]
            = g_thread_new (NULL, (GThreadFunc) optimize_thread, &data[i]);
        }
      for (i = 0; i < nthreads; ++i)
        g_thread_join (thread[i]);
    }
#if HAVE_MPI
  // Communicating tasks results
  optimize_synchronise ();
#endif
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_MonteCarlo: end\n");
#endif
}

/**
 * Function to optimize with the orthogonal sampling algorithm.
 */
static void
optimize_orthogonal ()
{
  unsigned int i, j, k, l;
  double e;
  GThread *thread[nthreads];
  ParallelData data[nthreads];
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_orthogonal: start\n");
#endif
  for (i = 0; i < optimize->nsimulations; ++i)
    {
      k = i;
      for (j = 0; j < optimize->nvariables; ++j)
        {
          l = k % optimize->nsweeps[j];
          k /= optimize->nsweeps[j];
          e = optimize->rangemin[j];
          if (optimize->nsweeps[j] > 1)
            e += (l + gsl_rng_uniform (optimize->rng))
              * (optimize->rangemax[j] - optimize->rangemin[j])
              / optimize->nsweeps[j];
          optimize->value[i * optimize->nvariables + j] = e;
        }
    }
  optimize->nsaveds = 0;
  if (nthreads <= 1)
    optimize_sequential ();
  else
    {
      for (i = 0; i < nthreads; ++i)
        {
          data[i].thread = i;
          thread[i]
            = g_thread_new (NULL, (GThreadFunc) optimize_thread, &data[i]);
        }
      for (i = 0; i < nthreads; ++i)
        g_thread_join (thread[i]);
    }
#if HAVE_MPI
  // Communicating tasks results
  optimize_synchronise ();
#endif
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_orthogonal: end\n");
#endif
}

/**
 * Function to save the best simulation in a hill climbing method.
 */
static void
optimize_best_climbing (unsigned int simulation,        ///< Simulation number.
                        double value)   ///< Objective function value.
{
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_best_climbing: start\n");
  fprintf (stderr,
           "optimize_best_climbing: simulation=%u value=%.14le best=%.14le\n",
           simulation, value, optimize->error_best[0]);
#endif
  if (value < optimize->error_best[0])
    {
      optimize->error_best[0] = value;
      optimize->simulation_best[0] = simulation;
#if DEBUG_OPTIMIZE
      fprintf (stderr,
               "optimize_best_climbing: BEST simulation=%u value=%.14le\n",
               simulation, value);
#endif
    }
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_best_climbing: end\n");
#endif
}

/**
 * Function to estimate the hill climbing sequentially.
 */
static inline void
optimize_climbing_sequential (unsigned int simulation)  ///< Simulation number.
{
  double e;
  unsigned int i, j;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_climbing_sequential: start\n");
  fprintf (stderr, "optimize_climbing_sequential: nstart_climbing=%u "
           "nend_climbing=%u\n",
           optimize->nstart_climbing, optimize->nend_climbing);
#endif
  for (i = optimize->nstart_climbing; i < optimize->nend_climbing; ++i)
    {
      j = simulation + i;
      e = optimize_norm (j);
      optimize_best_climbing (j, e);
      optimize_save_variables (j, e);
      if (e < optimize->threshold)
        {
          optimize->stop = 1;
          break;
        }
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_climbing_sequential: i=%u e=%lg\n", i, e);
#endif
    }
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_climbing_sequential: end\n");
#endif
}

/**
 * Function to estimate the hill climbing on a thread.
 *
 * \return NULL
 */
static void *
optimize_climbing_thread (ParallelData * data)  ///< Function data.
{
  unsigned int i, thread;
  double e;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_climbing_thread: start\n");
#endif
  thread = data->thread;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_climbing_thread: thread=%u start=%u end=%u\n",
           thread,
           optimize->thread_climbing[thread],
           optimize->thread_climbing[thread + 1]);
#endif
  for (i = optimize->thread_climbing[thread];
       i < optimize->thread_climbing[thread + 1]; ++i)
    {
      e = optimize_norm (i);
      g_mutex_lock (mutex);
      optimize_best_climbing (i, e);
      optimize_save_variables (i, e);
      if (e < optimize->threshold)
        optimize->stop = 1;
      g_mutex_unlock (mutex);
      if (optimize->stop)
        break;
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_climbing_thread: i=%u e=%lg\n", i, e);
#endif
    }
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_climbing_thread: end\n");
#endif
  g_thread_exit (NULL);
  return NULL;
}

/**
 * Function to estimate a component of the hill climbing vector.
 */
static double
optimize_estimate_climbing_random (unsigned int variable,
                                   ///< Variable number.
                                   unsigned int estimate
                                   __attribute__((unused)))
  ///< Estimate number.
{
  double x;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_estimate_climbing_random: start\n");
#endif
  x = optimize->climbing[variable]
    + (1. - 2. * gsl_rng_uniform (optimize->rng)) * optimize->step[variable];
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_estimate_climbing_random: climbing%u=%lg\n",
           variable, x);
  fprintf (stderr, "optimize_estimate_climbing_random: end\n");
#endif
  return x;
}

/**
 * Function to estimate a component of the hill climbing vector.
 */
static double
optimize_estimate_climbing_coordinates (unsigned int variable,
                                        ///< Variable number.
                                        unsigned int estimate)
                                        ///< Estimate number.
{
  double x;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_estimate_climbing_coordinates: start\n");
#endif
  x = optimize->climbing[variable];
  if (estimate >= (2 * variable) && estimate < (2 * variable + 2))
    {
      if (estimate & 1)
        x += optimize->step[variable];
      else
        x -= optimize->step[variable];
    }
#if DEBUG_OPTIMIZE
  fprintf (stderr,
           "optimize_estimate_climbing_coordinates: climbing%u=%lg\n",
           variable, x);
  fprintf (stderr, "optimize_estimate_climbing_coordinates: end\n");
#endif
  return x;
}

/**
 * Function to do a step of the hill climbing method.
 */
static inline void
optimize_step_climbing (unsigned int simulation)        ///< Simulation number.
{
  GThread *thread[nthreads_climbing];
  ParallelData data[nthreads_climbing];
  unsigned int i, j, k, b;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_step_climbing: start\n");
#endif
  for (i = 0; i < optimize->nestimates; ++i)
    {
      k = (simulation + i) * optimize->nvariables;
      b = optimize->simulation_best[0] * optimize->nvariables;
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_step_climbing: simulation=%u best=%u\n",
               simulation + i, optimize->simulation_best[0]);
#endif
      for (j = 0; j < optimize->nvariables; ++j, ++k, ++b)
        {
#if DEBUG_OPTIMIZE
          fprintf (stderr,
                   "optimize_step_climbing: estimate=%u best%u=%.14le\n",
                   i, j, optimize->value[b]);
#endif
          optimize->value[k]
            = optimize->value[b] + optimize_estimate_climbing (j, i);
          optimize->value[k] = fmin (fmax (optimize->value[k],
                                           optimize->rangeminabs[j]),
                                     optimize->rangemaxabs[j]);
#if DEBUG_OPTIMIZE
          fprintf (stderr,
                   "optimize_step_climbing: estimate=%u variable%u=%.14le\n",
                   i, j, optimize->value[k]);
#endif
        }
    }
  if (nthreads_climbing == 1)
    optimize_climbing_sequential (simulation);
  else
    {
      for (i = 0; i <= nthreads_climbing; ++i)
        {
          optimize->thread_climbing[i]
            = simulation + optimize->nstart_climbing
            + i * (optimize->nend_climbing - optimize->nstart_climbing)
            / nthreads_climbing;
#if DEBUG_OPTIMIZE
          fprintf (stderr,
                   "optimize_step_climbing: i=%u thread_climbing=%u\n",
                   i, optimize->thread_climbing[i]);
#endif
        }
      for (i = 0; i < nthreads_climbing; ++i)
        {
          data[i].thread = i;
          thread[i] = g_thread_new
            (NULL, (GThreadFunc) optimize_climbing_thread, &data[i]);
        }
      for (i = 0; i < nthreads_climbing; ++i)
        g_thread_join (thread[i]);
    }
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_step_climbing: end\n");
#endif
}

/**
 * Function to optimize with a hill climbing method.
 */
static inline void
optimize_climbing ()
{
  unsigned int i, j, k, b, s, adjust;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_climbing: start\n");
#endif
  for (i = 0; i < optimize->nvariables; ++i)
    optimize->climbing[i] = 0.;
  b = optimize->simulation_best[0] * optimize->nvariables;
  s = optimize->nsimulations;
  adjust = 1;
  for (i = 0; i < optimize->nsteps; ++i, s += optimize->nestimates, b = k)
    {
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_climbing: step=%u old_best=%u\n",
               i, optimize->simulation_best[0]);
#endif
      optimize_step_climbing (s);
      k = optimize->simulation_best[0] * optimize->nvariables;
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_climbing: step=%u best=%u\n",
               i, optimize->simulation_best[0]);
#endif
      if (k == b)
        {
          if (adjust)
            for (j = 0; j < optimize->nvariables; ++j)
              optimize->step[j] *= 0.5;
          for (j = 0; j < optimize->nvariables; ++j)
            optimize->climbing[j] = 0.;
          adjust = 1;
        }
      else
        {
          for (j = 0; j < optimize->nvariables; ++j)
            {
#if DEBUG_OPTIMIZE
              fprintf (stderr,
                       "optimize_climbing: best%u=%.14le old%u=%.14le\n",
                       j, optimize->value[k + j], j, optimize->value[b + j]);
#endif
              optimize->climbing[j]
                = (1. - optimize->relaxation) * optimize->climbing[j]
                + optimize->relaxation
                * (optimize->value[k + j] - optimize->value[b + j]);
#if DEBUG_OPTIMIZE
              fprintf (stderr, "optimize_climbing: climbing%u=%.14le\n",
                       j, optimize->climbing[j]);
#endif
            }
          adjust = 0;
        }
    }
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_climbing: end\n");
#endif
}

/**
 * Function to calculate the objective function of an entity.
 *
 * \return objective function value.
 */
static double
optimize_genetic_objective (Entity * entity)    ///< entity data.
{
  unsigned int j;
  double objective;
  char buffer[64];
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_genetic_objective: start\n");
#endif
  for (j = 0; j < optimize->nvariables; ++j)
    {
      optimize->value[entity->id * optimize->nvariables + j]
        = genetic_get_variable (entity, optimize->genetic_variable + j);
    }
  objective = optimize_norm (entity->id);
  g_mutex_lock (mutex);
  for (j = 0; j < optimize->nvariables; ++j)
    {
      snprintf (buffer, 64, "%s ", format[optimize->precision[j]]);
      fprintf (optimize->file_variables, buffer,
               genetic_get_variable (entity, optimize->genetic_variable + j));
    }
  fprintf (optimize->file_variables, "%.14le\n", objective);
  g_mutex_unlock (mutex);
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_genetic_objective: end\n");
#endif
  return objective;
}

/**
 * Function to optimize with the genetic algorithm.
 */
static void
optimize_genetic ()
{
  double *best_variable = NULL;
  char *best_genome = NULL;
  double best_objective = 0.;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_genetic: start\n");
  fprintf (stderr, "optimize_genetic: ntasks=%u nthreads=%u\n", ntasks,
           nthreads);
  fprintf (stderr,
           "optimize_genetic: nvariables=%u population=%u generations=%u\n",
           optimize->nvariables, optimize->nsimulations, optimize->niterations);
  fprintf (stderr,
           "optimize_genetic: mutation=%lg reproduction=%lg adaptation=%lg\n",
           optimize->mutation_ratio, optimize->reproduction_ratio,
           optimize->adaptation_ratio);
#endif
  genetic_algorithm_default (optimize->nvariables,
                             optimize->genetic_variable,
                             optimize->nsimulations,
                             optimize->niterations,
                             optimize->mutation_ratio,
                             optimize->reproduction_ratio,
                             optimize->adaptation_ratio,
                             optimize->seed,
                             optimize->threshold,
                             &optimize_genetic_objective,
                             &best_genome, &best_variable, &best_objective);
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_genetic: the best\n");
#endif
  optimize->error_old = (double *) g_malloc (sizeof (double));
  optimize->value_old
    = (double *) g_malloc (optimize->nvariables * sizeof (double));
  optimize->error_old[0] = best_objective;
  memcpy (optimize->value_old, best_variable,
          optimize->nvariables * sizeof (double));
  g_free (best_genome);
  g_free (best_variable);
  optimize_print ();
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_genetic: end\n");
#endif
}

/**
 * Function to save the best results on iterative methods.
 */
static inline void
optimize_save_old ()
{
  unsigned int i, j;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_save_old: start\n");
  fprintf (stderr, "optimize_save_old: nsaveds=%u\n", optimize->nsaveds);
#endif
  memcpy (optimize->error_old, optimize->error_best,
          optimize->nbest * sizeof (double));
  for (i = 0; i < optimize->nbest; ++i)
    {
      j = optimize->simulation_best[i];
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_save_old: i=%u j=%u\n", i, j);
#endif
      memcpy (optimize->value_old + i * optimize->nvariables,
              optimize->value + j * optimize->nvariables,
              optimize->nvariables * sizeof (double));
    }
#if DEBUG_OPTIMIZE
  for (i = 0; i < optimize->nvariables; ++i)
    fprintf (stderr, "optimize_save_old: best variable %u=%lg\n",
             i, optimize->value_old[i]);
  fprintf (stderr, "optimize_save_old: end\n");
#endif
}

/**
 * Function to merge the best results with the previous step best results on
 *   iterative methods.
 */
static inline void
optimize_merge_old ()
{
  unsigned int i, j, k;
  double v[optimize->nbest * optimize->nvariables], e[optimize->nbest],
    *enew, *eold;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_merge_old: start\n");
#endif
  enew = optimize->error_best;
  eold = optimize->error_old;
  i = j = k = 0;
  do
    {
      if (*enew < *eold)
        {
          memcpy (v + k * optimize->nvariables,
                  optimize->value
                  + optimize->simulation_best[i] * optimize->nvariables,
                  optimize->nvariables * sizeof (double));
          e[k] = *enew;
          ++k;
          ++enew;
          ++i;
        }
      else
        {
          memcpy (v + k * optimize->nvariables,
                  optimize->value_old + j * optimize->nvariables,
                  optimize->nvariables * sizeof (double));
          e[k] = *eold;
          ++k;
          ++eold;
          ++j;
        }
    }
  while (k < optimize->nbest);
  memcpy (optimize->value_old, v, k * optimize->nvariables * sizeof (double));
  memcpy (optimize->error_old, e, k * sizeof (double));
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_merge_old: end\n");
#endif
}

/**
 * Function to refine the search ranges of the variables in iterative
 *   algorithms.
 */
static inline void
optimize_refine ()
{
  unsigned int i, j;
  double d;
#if HAVE_MPI
  MPI_Status mpi_stat;
#endif
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_refine: start\n");
#endif
#if HAVE_MPI
  if (!optimize->mpi_rank)
    {
#endif
      for (j = 0; j < optimize->nvariables; ++j)
        {
          optimize->rangemin[j] = optimize->rangemax[j]
            = optimize->value_old[j];
        }
      for (i = 0; ++i < optimize->nbest;)
        {
          for (j = 0; j < optimize->nvariables; ++j)
            {
              optimize->rangemin[j]
                = fmin (optimize->rangemin[j],
                        optimize->value_old[i * optimize->nvariables + j]);
              optimize->rangemax[j]
                = fmax (optimize->rangemax[j],
                        optimize->value_old[i * optimize->nvariables + j]);
            }
        }
      for (j = 0; j < optimize->nvariables; ++j)
        {
          d = optimize->tolerance
            * (optimize->rangemax[j] - optimize->rangemin[j]);
          switch (optimize->algorithm)
            {
            case ALGORITHM_MONTE_CARLO:
              d *= 0.5;
              break;
            default:
              if (optimize->nsweeps[j] > 1)
                d /= optimize->nsweeps[j] - 1;
              else
                d = 0.;
            }
          optimize->rangemin[j] -= d;
          optimize->rangemin[j]
            = fmax (optimize->rangemin[j], optimize->rangeminabs[j]);
          optimize->rangemax[j] += d;
          optimize->rangemax[j]
            = fmin (optimize->rangemax[j], optimize->rangemaxabs[j]);
          printf ("%s min=%lg max=%lg\n", optimize->label[j],
                  optimize->rangemin[j], optimize->rangemax[j]);
          fprintf (optimize->file_result, "%s min=%lg max=%lg\n",
                   optimize->label[j], optimize->rangemin[j],
                   optimize->rangemax[j]);
        }
#if HAVE_MPI
      for (i = 1; (int) i < ntasks; ++i)
        {
          MPI_Send (optimize->rangemin, optimize->nvariables, MPI_DOUBLE, i,
                    1, MPI_COMM_WORLD);
          MPI_Send (optimize->rangemax, optimize->nvariables, MPI_DOUBLE, i,
                    1, MPI_COMM_WORLD);
        }
    }
  else
    {
      MPI_Recv (optimize->rangemin, optimize->nvariables, MPI_DOUBLE, 0, 1,
                MPI_COMM_WORLD, &mpi_stat);
      MPI_Recv (optimize->rangemax, optimize->nvariables, MPI_DOUBLE, 0, 1,
                MPI_COMM_WORLD, &mpi_stat);
    }
#endif
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_refine: end\n");
#endif
}

/**
 * Function to do a step of the iterative algorithm.
 */
static void
optimize_step ()
{
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_step: start\n");
#endif
  optimize_algorithm ();
  if (optimize->nsteps)
    optimize_climbing ();
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_step: end\n");
#endif
}

/**
 * Function to iterate the algorithm.
 */
static inline void
optimize_iterate ()
{
  unsigned int i;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_iterate: start\n");
#endif
  optimize->error_old = (double *) g_malloc (optimize->nbest * sizeof (double));
  optimize->value_old =
    (double *) g_malloc (optimize->nbest * optimize->nvariables *
                         sizeof (double));
  optimize_step ();
  optimize_save_old ();
  optimize_refine ();
  optimize_print ();
  for (i = 1; i < optimize->niterations && !optimize->stop; ++i)
    {
      optimize_step ();
      optimize_merge_old ();
      optimize_refine ();
      optimize_print ();
    }
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_iterate: end\n");
#endif
}

/**
 * Function to free the memory used by the Optimize struct.
 */
void
optimize_free ()
{
  unsigned int i, j;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_free: start\n");
#endif
  for (j = 0; j < optimize->ninputs; ++j)
    {
      for (i = 0; i < optimize->nexperiments; ++i)
        g_mapped_file_unref (optimize->file[j][i]);
      g_free (optimize->file[j]);
    }
  g_free (optimize->error_old);
  g_free (optimize->value_old);
  g_free (optimize->value);
  g_free (optimize->genetic_variable);
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_free: end\n");
#endif
}

/**
 * Function to open and perform a optimization.
 */
void
optimize_open ()
{
  GTimeZone *tz;
  GDateTime *t0, *t;
  unsigned int i, j;

#if DEBUG_OPTIMIZE
  char *buffer;
  fprintf (stderr, "optimize_open: start\n");
#endif

  // Getting initial time
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_open: getting initial time\n");
#endif
  tz = g_time_zone_new_utc ();
  t0 = g_date_time_new_now (tz);

  // Obtaining and initing the pseudo-random numbers generator seed
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_open: getting initial seed\n");
#endif
  if (optimize->seed == DEFAULT_RANDOM_SEED)
    optimize->seed = input->seed;
  gsl_rng_set (optimize->rng, optimize->seed);

  // Replacing the working directory
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_open: replacing the working directory\n");
#endif
  g_chdir (input->directory);

  // Getting results file names
  optimize->result = input->result;
  optimize->variables = input->variables;

  // Obtaining the simulator file
  optimize->simulator = input->simulator;

  // Obtaining the evaluator file
  optimize->evaluator = input->evaluator;

  // Reading the algorithm
  optimize->algorithm = input->algorithm;
  switch (optimize->algorithm)
    {
    case ALGORITHM_MONTE_CARLO:
      optimize_algorithm = optimize_MonteCarlo;
      break;
    case ALGORITHM_SWEEP:
      optimize_algorithm = optimize_sweep;
      break;
    case ALGORITHM_ORTHOGONAL:
      optimize_algorithm = optimize_orthogonal;
      break;
    default:
      optimize_algorithm = optimize_genetic;
      optimize->mutation_ratio = input->mutation_ratio;
      optimize->reproduction_ratio = input->reproduction_ratio;
      optimize->adaptation_ratio = input->adaptation_ratio;
    }
  optimize->nvariables = input->nvariables;
  optimize->nsimulations = input->nsimulations;
  optimize->niterations = input->niterations;
  optimize->nbest = input->nbest;
  optimize->tolerance = input->tolerance;
  optimize->nsteps = input->nsteps;
  optimize->nestimates = 0;
  optimize->threshold = input->threshold;
  optimize->stop = 0;
  if (input->nsteps)
    {
      optimize->relaxation = input->relaxation;
      switch (input->climbing)
        {
        case CLIMBING_METHOD_COORDINATES:
          optimize->nestimates = 2 * optimize->nvariables;
          optimize_estimate_climbing = optimize_estimate_climbing_coordinates;
          break;
        default:
          optimize->nestimates = input->nestimates;
          optimize_estimate_climbing = optimize_estimate_climbing_random;
        }
    }

#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_open: nbest=%u\n", optimize->nbest);
#endif
  optimize->simulation_best
    = (unsigned int *) alloca (optimize->nbest * sizeof (unsigned int));
  optimize->error_best = (double *) alloca (optimize->nbest * sizeof (double));

  // Reading the experimental data
#if DEBUG_OPTIMIZE
  buffer = g_get_current_dir ();
  fprintf (stderr, "optimize_open: current directory=%s\n", buffer);
  g_free (buffer);
#endif
  optimize->nexperiments = input->nexperiments;
  optimize->ninputs = input->experiment->ninputs;
  optimize->experiment
    = (char **) alloca (input->nexperiments * sizeof (char *));
  optimize->weight = (double *) alloca (input->nexperiments * sizeof (double));
  for (i = 0; i < input->experiment->ninputs; ++i)
    optimize->file[i] = (GMappedFile **)
      g_malloc (input->nexperiments * sizeof (GMappedFile *));
  for (i = 0; i < input->nexperiments; ++i)
    {
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_open: i=%u\n", i);
#endif
      optimize->experiment[i] = input->experiment[i].name;
      optimize->weight[i] = input->experiment[i].weight;
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_open: experiment=%s weight=%lg\n",
               optimize->experiment[i], optimize->weight[i]);
#endif
      for (j = 0; j < input->experiment->ninputs; ++j)
        {
#if DEBUG_OPTIMIZE
          fprintf (stderr, "optimize_open: stencil%u\n", j + 1);
#endif
          optimize->file[j][i]
            = g_mapped_file_new (input->experiment[i].stencil[j], 0, NULL);
        }
    }

  // Reading the variables data
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_open: reading variables\n");
#endif
  optimize->label = (char **) alloca (input->nvariables * sizeof (char *));
  j = input->nvariables * sizeof (double);
  optimize->rangemin = (double *) alloca (j);
  optimize->rangeminabs = (double *) alloca (j);
  optimize->rangemax = (double *) alloca (j);
  optimize->rangemaxabs = (double *) alloca (j);
  optimize->step = (double *) alloca (j);
  j = input->nvariables * sizeof (unsigned int);
  optimize->precision = (unsigned int *) alloca (j);
  optimize->nsweeps = (unsigned int *) alloca (j);
  optimize->nbits = (unsigned int *) alloca (j);
  for (i = 0; i < input->nvariables; ++i)
    {
      optimize->label[i] = input->variable[i].name;
      optimize->rangemin[i] = input->variable[i].rangemin;
      optimize->rangeminabs[i] = input->variable[i].rangeminabs;
      optimize->rangemax[i] = input->variable[i].rangemax;
      optimize->rangemaxabs[i] = input->variable[i].rangemaxabs;
      optimize->precision[i] = input->variable[i].precision;
      optimize->step[i] = input->variable[i].step;
      optimize->nsweeps[i] = input->variable[i].nsweeps;
      optimize->nbits[i] = input->variable[i].nbits;
    }
  if (input->algorithm == ALGORITHM_SWEEP
      || input->algorithm == ALGORITHM_ORTHOGONAL)
    {
      optimize->nsimulations = 1;
      for (i = 0; i < input->nvariables; ++i)
        {
          optimize->nsimulations *= optimize->nsweeps[i];
#if DEBUG_OPTIMIZE
          fprintf (stderr, "optimize_open: nsweeps=%u nsimulations=%u\n",
                   optimize->nsweeps[i], optimize->nsimulations);
#endif
        }
    }
  if (optimize->nsteps)
    optimize->climbing
      = (double *) alloca (optimize->nvariables * sizeof (double));

  // Setting error norm
  switch (input->norm)
    {
    case ERROR_NORM_EUCLIDIAN:
      optimize_norm = optimize_norm_euclidian;
      break;
    case ERROR_NORM_MAXIMUM:
      optimize_norm = optimize_norm_maximum;
      break;
    case ERROR_NORM_P:
      optimize_norm = optimize_norm_p;
      optimize->p = input->p;
      break;
    default:
      optimize_norm = optimize_norm_taxicab;
    }

  // Allocating values
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_open: allocating variables\n");
  fprintf (stderr, "optimize_open: nvariables=%u algorithm=%u\n",
           optimize->nvariables, optimize->algorithm);
#endif
  optimize->genetic_variable = NULL;
  if (optimize->algorithm == ALGORITHM_GENETIC)
    {
      optimize->genetic_variable = (GeneticVariable *)
        g_malloc (optimize->nvariables * sizeof (GeneticVariable));
      for (i = 0; i < optimize->nvariables; ++i)
        {
#if DEBUG_OPTIMIZE
          fprintf (stderr, "optimize_open: i=%u min=%lg max=%lg nbits=%u\n",
                   i, optimize->rangemin[i], optimize->rangemax[i],
                   optimize->nbits[i]);
#endif
          optimize->genetic_variable[i].minimum = optimize->rangemin[i];
          optimize->genetic_variable[i].maximum = optimize->rangemax[i];
          optimize->genetic_variable[i].nbits = optimize->nbits[i];
        }
    }
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_open: nvariables=%u nsimulations=%u\n",
           optimize->nvariables, optimize->nsimulations);
#endif
  optimize->value = (double *)
    g_malloc ((optimize->nsimulations
               + optimize->nestimates * optimize->nsteps)
              * optimize->nvariables * sizeof (double));

  // Calculating simulations to perform for each task
#if HAVE_MPI
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_open: rank=%u ntasks=%u\n",
           optimize->mpi_rank, ntasks);
#endif
  optimize->nstart = optimize->mpi_rank * optimize->nsimulations / ntasks;
  optimize->nend = (1 + optimize->mpi_rank) * optimize->nsimulations / ntasks;
  if (optimize->nsteps)
    {
      optimize->nstart_climbing
        = optimize->mpi_rank * optimize->nestimates / ntasks;
      optimize->nend_climbing
        = (1 + optimize->mpi_rank) * optimize->nestimates / ntasks;
    }
#else
  optimize->nstart = 0;
  optimize->nend = optimize->nsimulations;
  if (optimize->nsteps)
    {
      optimize->nstart_climbing = 0;
      optimize->nend_climbing = optimize->nestimates;
    }
#endif
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_open: nstart=%u nend=%u\n", optimize->nstart,
           optimize->nend);
#endif

  // Calculating simulations to perform for each thread
  optimize->thread
    = (unsigned int *) alloca ((1 + nthreads) * sizeof (unsigned int));
  for (i = 0; i <= nthreads; ++i)
    {
      optimize->thread[i] = optimize->nstart
        + i * (optimize->nend - optimize->nstart) / nthreads;
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_open: i=%u thread=%u\n", i,
               optimize->thread[i]);
#endif
    }
  if (optimize->nsteps)
    optimize->thread_climbing = (unsigned int *)
      alloca ((1 + nthreads_climbing) * sizeof (unsigned int));

  // Opening result files
  optimize->file_result = g_fopen (optimize->result, "w");
  optimize->file_variables = g_fopen (optimize->variables, "w");

  // Performing the algorithm
  switch (optimize->algorithm)
    {
      // Genetic algorithm
    case ALGORITHM_GENETIC:
      optimize_genetic ();
      break;

      // Iterative algorithm
    default:
      optimize_iterate ();
    }

  // Getting calculation time
  t = g_date_time_new_now (tz);
  optimize->calculation_time = 0.000001 * g_date_time_difference (t, t0);
  g_date_time_unref (t);
  g_date_time_unref (t0);
  g_time_zone_unref (tz);
  printf ("%s = %.6lg s\n", _("Calculation time"), optimize->calculation_time);
  fprintf (optimize->file_result, "%s = %.6lg s\n",
           _("Calculation time"), optimize->calculation_time);

  // Closing result files
  fclose (optimize->file_variables);
  fclose (optimize->file_result);

#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_open: end\n");
#endif
}

ext: c

max_stars_repo_name: jburguete/mpcotool

max_issues_repo_name: jburguete/mpcotool

max_forks_repo_name: jburguete/mpcotool


text:
/*
   LAPACKE_dgesv Example
   =====================

   The program computes the solution to the system of linear
   equations with a square matrix A and multiple
   right-hand sides B, where A is the coefficient matrix
   and b is the right-hand side matrix:

   Description
   ===========

   The routine solves for X the system of linear equations A*X = B,
   where A is an n-by-n matrix, the columns of matrix B are individual
   right-hand sides, and the columns of X are the corresponding
   solutions.

   The LU decomposition with partial pivoting and row interchanges is
   used to factor A as A = P*L*U, where P is a permutation matrix, L
   is unit lower triangular, and U is upper triangular. The factored
   form of A is then used to solve the system of equations A*X = B.

   LAPACKE Interface
   =================

   LAPACKE_dgesv (row-major, high-level) Example Program Results

  -- LAPACKE Example routine (version 3.7.0) --
  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
     December 2016

*/
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <math.h>
#include <sys/time.h>
#include <unistd.h>
#include <lapacke.h>
#include <string.h>
#ifndef _LAPACKE_EXAMPLE_AUX_
#define _LAPACKE_EXAMPLE_AUX_

void print_matrix_rowmajor( char* desc, lapack_int m, lapack_int n, double* mat, lapack_int ldm );
void print_matrix_colmajor( char* desc, lapack_int m, lapack_int n, double* mat, lapack_int ldm );
void print_vector( char* desc, lapack_int n, lapack_int* vec );

#endif /* _LAPACKE_EXAMPLE_AUX_*/

#define FLT double
//void mset(FLT **m, int n, int in) {
void mset(FLT *m, int n, int in) {
	int i,j;
    for(i=0;i<n;i++)
       for(j=0;j<n;j++) {
           if(i == j) {
               //m[i][j]=in;
               m[i+j*n]=in;
           } else {
              // m[i][j]=0.1;
               m[i+j*n]=0.1;
           }
       }

}
FLT system_clock() {
	FLT t;
	FLT six=1.0e-6;
	struct timeval tb;
	struct timezone tz;
	gettimeofday(&tb,&tz);
	t=(FLT)tb.tv_sec+((FLT)tb.tv_usec)*six;
 	return(t);
}

/* Main program */
int main(int argc, char **argv) {

        /* Locals */
        lapack_int n, nrhs, lda, ldb, info;
		int i, j;
		int DOP=0;
        /* Local arrays */
		double *A, *b;
		lapack_int *ipiv;
		int icount,jcount;
		double t1,t2,tmin,tmax,dt,tstart;
		double *tvect;
		tmin=1e30;
		tmax=0.0;

        /* Default Value */
	    n = 5; nrhs = 1;

        /* Arguments */
	    for( i = 1; i < argc; i++ ) {
	    	if( strcmp( argv[i], "-n" ) == 0 ) {
		    	n  = atoi(argv[i+1]);
			    i++;
		    }
			if( strcmp( argv[i], "-nrhs" ) == 0 ) {
				nrhs  = atoi(argv[i+1]);
				i++;
			}
		}
        scanf("%d",&n);
        /* Initialization */
        lda=n, ldb=nrhs;
		A = (double *)malloc(n*n*sizeof(double)) ;
		if (A==NULL){ printf("error of memory allocation\n"); exit(0); }
		b = (double *)malloc(n*nrhs*sizeof(double)) ;
		if (b==NULL){ printf("error of memory allocation\n"); exit(0); }
		ipiv = (lapack_int *)malloc(n*sizeof(lapack_int)) ;
		if (ipiv==NULL){ printf("error of memory allocation\n"); exit(0); }
		tstart=system_clock();
		tvect=(double*)malloc(500*sizeof(double)) ;
        for (icount=0; icount<500;icount++){
			mset(A,n,10);

			for(i=0;i<n*nrhs;i++)
				b[i] = 10.0;
				//b[i] = ((double) rand()) / ((double) RAND_MAX) - 0.5;

			/* Print Entry Matrix */
			if(DOP)print_matrix_rowmajor( "Entry Matrix A", n, n, A, lda );
			/* Print Right Rand Side */
			if(DOP)print_matrix_rowmajor( "Right Rand Side b", n, nrhs, b, ldb );
			if(DOP)printf( "\n" );
			/* Executable statements */
			if(DOP)printf( "LAPACKE_dgesv (row-major, high-level) Example Program Results\n" );
			/* Solve the equations A*X = B */
			t1=system_clock();
			info = LAPACKE_dgesv( LAPACK_ROW_MAJOR, n, nrhs, A, lda, ipiv,
							b, ldb );
			t2=system_clock();
			dt=t2-t1;
			jcount=icount+1;
			if(dt > tmax)tmax=dt;
			if(dt < tmin)tmin=dt;
			tvect[icount]=dt;
			if(t2-tstart> 120.0)icount=1000000;
			/* Check for the exact singularity */
			if( info > 0 ) {
					printf( "The diagonal element of the triangular factor of A,\n" );
					printf( "U(%i,%i) is zero, so that A is singular;\n", info, info );
					printf( "the solution could not be computed.\n" );
					exit( 1 );
			}
			if (info <0) exit( 1 );
			/* Print solution */
			if(DOP)print_matrix_rowmajor( "Solution", n, nrhs, b, ldb );
			/* Print details of LU factorization */
			if(DOP)print_matrix_rowmajor( "Details of LU factorization", n, n, A, lda );
			/* Print pivot indices */
			if(DOP)print_vector( "Pivot indices", n, ipiv );
        }
        printf("size= %d min=%g max=%g total=%g inverts=%d\n",n,tmin,tmax,t2-tstart,jcount);
        n=0;
        for(icount=0;icount<jcount;icount++){
        	printf("%f10.5 ",tvect[icount]);
        	n=n+1;
        	if(n==8){
        		printf("\n");
        		n=0;}
        }
        exit( 0 );
} /* End of LAPACKE_dgesv Example */

max_stars_repo_name: timkphd/examples

max_issues_repo_name: timkphd/examples

max_forks_repo_name: timkphd/examples


text:
/* specfunc/gsl_sf_airy.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#ifndef __GSL_SF_AIRY_H__
#define __GSL_SF_AIRY_H__

#include <gsl/gsl_mode.h>
#include <gsl/gsl_sf_result.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

/* Airy function Ai(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_airy_Ai_e(const double x, const gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_airy_Ai(const double x, gsl_mode_t mode);

/* Airy function Bi(x)
 *
 * exceptions: GSL_EOVRFLW
 */
int gsl_sf_airy_Bi_e(const double x, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_airy_Bi(const double x, gsl_mode_t mode);

/* scaled Ai(x):
 *                     Ai(x)   x < 0
 *   exp(+2/3 x^{3/2}) Ai(x)   x > 0
 *
 * exceptions: none
 */
int gsl_sf_airy_Ai_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_airy_Ai_scaled(const double x, gsl_mode_t mode);

/* scaled Bi(x):
 *                     Bi(x)   x < 0
 *   exp(-2/3 x^{3/2}) Bi(x)   x > 0
 *
 * exceptions: none
 */
int gsl_sf_airy_Bi_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_airy_Bi_scaled(const double x, gsl_mode_t mode);

/* derivative Ai'(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_airy_Ai_deriv_e(const double x, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_airy_Ai_deriv(const double x, gsl_mode_t mode);

/* derivative Bi'(x)
 *
 * exceptions: GSL_EOVRFLW
 */
int gsl_sf_airy_Bi_deriv_e(const double x, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_airy_Bi_deriv(const double x, gsl_mode_t mode);

/* scaled derivative Ai'(x):
 *                     Ai'(x)   x < 0
 *   exp(+2/3 x^{3/2}) Ai'(x)   x > 0
 *
 * exceptions: none
 */
int gsl_sf_airy_Ai_deriv_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_airy_Ai_deriv_scaled(const double x, gsl_mode_t mode);

/* scaled derivative:
 *                     Bi'(x)   x < 0
 *   exp(-2/3 x^{3/2}) Bi'(x)   x > 0
 *
 * exceptions: none
 */
int gsl_sf_airy_Bi_deriv_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_airy_Bi_deriv_scaled(const double x, gsl_mode_t mode);

/* Zeros of Ai(x)
 */
int gsl_sf_airy_zero_Ai_e(unsigned int s, gsl_sf_result * result);
double gsl_sf_airy_zero_Ai(unsigned int s);

/* Zeros of Bi(x)
 */
int gsl_sf_airy_zero_Bi_e(unsigned int s, gsl_sf_result * result);
double gsl_sf_airy_zero_Bi(unsigned int s);

/* Zeros of Ai'(x)
 */
int gsl_sf_airy_zero_Ai_deriv_e(unsigned int s, gsl_sf_result * result);
double gsl_sf_airy_zero_Ai_deriv(unsigned int s);

/* Zeros of Bi'(x)
 */
int gsl_sf_airy_zero_Bi_deriv_e(unsigned int s, gsl_sf_result * result);
double gsl_sf_airy_zero_Bi_deriv(unsigned int s);

__END_DECLS

#endif /* __GSL_SF_AIRY_H__ */

ext: h


text:
#pragma once

#include <bgfx/bgfx.h>
#include <napi/napi.h>
#include <gsl/gsl>

#include <optional>

namespace Babylon
{
    class IndexBuffer final
    {
    public:
        IndexBuffer(gsl::span<uint8_t> bytes, uint16_t flags, bool dynamic);
        ~IndexBuffer();

        void Dispose();

        void Update(Napi::Env env, gsl::span<uint8_t> bytes, uint32_t startIndex);
        bool CreateHandle();
        void Set(bgfx::Encoder* encoder, uint32_t firstIndex, uint32_t numIndices);

    private:
        std::optional<std::vector<uint8_t>> m_bytes{};
        uint16_t m_flags{};
        bool m_dynamic{};

        union
        {
            bgfx::IndexBufferHandle m_handle{bgfx::kInvalidHandle};
            bgfx::DynamicIndexBufferHandle m_dynamicHandle;
        };

        bool m_disposed{};
    };
}

max_stars_repo_name: HarveyLijh/ReactNativeBabylon

max_issues_repo_name: HarveyLijh/ReactNativeBabylon

max_forks_repo_name: HarveyLijh/ReactNativeBabylon


text:
/*
 *            Copyright 2009-2017 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      Licensed under the Apache License, Version 2.0 (the "License")
 *
 * You may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *              http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

#ifndef _VOTCA_XTP_DIIS__H
#define _VOTCA_XTP_DIIS__H

#include <gsl/gsl_multimin.h>
#include <gsl/gsl_vector.h>

#include <votca/tools/linalg.h>
#include <votca/xtp/aomatrix.h>
#include <votca/xtp/orbitals.h>
#include <votca/ctp/logger.h>

using namespace votca::tools;

namespace votca { namespace xtp {
 namespace ub = boost::numeric::ublas;

 class Diis{
public:

    Diis() {_maxerrorindex=0;
                _maxerror=0.0; };
   ~Diis() {
     for (std::vector< ub::matrix<double>* >::iterator it = _mathist.begin() ; it !=_mathist.end(); ++it){
         delete *it;
     }
     _mathist.clear();
      for (std::vector< ub::matrix<double>* >::iterator it = _dmathist.begin() ; it !=_dmathist.end(); ++it){
         delete *it;
     }
     _dmathist.clear();
     for (std::vector< ub::matrix<double>* >::iterator it = _errormatrixhist.begin() ; it !=_errormatrixhist.end(); ++it){
         delete *it;
     }
    _errormatrixhist.clear();

    for (std::vector< std::vector<double>* >::iterator it = _Diis_Bs.begin() ; it !=_Diis_Bs.end(); ++it){
         delete *it;
     }
    _Diis_Bs.clear();

   }

   void Configure(bool usediis,bool noisy, unsigned histlength, bool maxout, string diismethod, double adiis_start,double diis_start,double levelshift,double levelshiftend,unsigned nocclevels){

       _usediis=usediis;
       _noisy=noisy;
       _histlength=histlength;
       _maxout=maxout;
       _diismethod=diismethod;
       _adiis_start=adiis_start;
       _diis_start=diis_start;
       _levelshift=levelshift;
       _levelshiftend=levelshiftend;

       _nocclevels=nocclevels;

   }

   void setOverlap(ub::matrix<double>* _S){
       S=_S;
   }
   void setSqrtOverlap(ub::matrix<double>* _Sminusahalf){
       Sminusahalf=_Sminusahalf;
   }
    void setLogger(ctp::Logger *pLog){_pLog=pLog;}
    double Evolve(const ub::matrix<double>& dmat,const ub::matrix<double>& H,ub::vector<double> &MOenergies,ub::matrix<double> &MOs, int this_iter,double totE);
    void SolveFockmatrix(ub::vector<double>& MOenergies,ub::matrix<double>& MOs,ub::matrix<double>&H);
    void Levelshift(ub::matrix<double>& H,const ub::matrix<double>&MOs);
    unsigned gethistlength(){return _mathist.size();}

    double get_E_adiis(const gsl_vector * x) const;

    void get_dEdx_adiis(const gsl_vector * x, gsl_vector * dEdx) const;
    void get_E_dEdx_adiis(const gsl_vector * x, double * Eval, gsl_vector * dEdx) const;

 private:

    ctp::Logger *_pLog;
    ub::matrix<double>* S;
    ub::matrix<double>* Sminusahalf;
    bool                              _usediis;
    bool                              _noisy;
    unsigned                          _histlength;
    bool                              _maxout;
    string                            _diismethod;
    ub::matrix<double>                _Sminusonehalf;
    double                              _maxerror;
    double                              _adiis_start;
    double                              _diis_start;
    double                              _levelshiftend;
    unsigned                            _maxerrorindex;
    std::vector< ub::matrix<double>* >   _mathist;
    std::vector< ub::matrix<double>* >   _dmathist;
    std::vector< ub::matrix<double>* >   _errormatrixhist;
    std::vector< std::vector<double>* >  _Diis_Bs;

    std::vector<double>                 _totE;
    ub::vector<double>                  _DiF;
    ub::matrix<double>                  _DiFj;

    ub::vector<double> ADIIsCoeff();

 ub::vector<double> compute_c(const gsl_vector * x);
 /// Compute jacobian
 ub::matrix<double> compute_jac(const gsl_vector * x);
 /// Compute energy
 double min_f(const gsl_vector * x, void * params);
 /// Compute derivative
 void min_df(const gsl_vector * x, void * params, gsl_vector * g);
 /// Compute energy and derivative
void min_fdf(const gsl_vector * x, void * params, double * f, gsl_vector * g);

ub::vector<double> DIIsCoeff();

    unsigned _nocclevels;
    double _levelshift;

 };

 namespace adiis {
  /// Compute weights
  ub::vector<double> compute_c(const gsl_vector * x);
  /// Compute jacobian
  ub::matrix<double> compute_jac(const gsl_vector * x);

  /// Compute energy
  double min_f(const gsl_vector * x, void * params);
  /// Compute derivative
  void min_df(const gsl_vector * x, void * params, gsl_vector * g);
  /// Compute energy and derivative
  void min_fdf(const gsl_vector * x, void * params, double * f, gsl_vector * g);
};

}}

#endif

max_stars_repo_name: choudarykvsp/xtp

max_issues_repo_name: choudarykvsp/xtp

max_forks_repo_name: choudarykvsp/xtp


text:
/* rb.c
 *
 * Copyright (C) 1998-2002, 2004 Free Software Foundation, Inc.
 * Copyright (C) 2018 Patrick Alken
 *
 * 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 3 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.
 */

/* This code is originally from GNU libavl, with some modifications */

#include <config.h>
#include <stdlib.h>
#include <stddef.h>
#include <string.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_bst.h>
#include <gsl/gsl_errno.h>

typedef struct gsl_bst_rb_node rb_node;
typedef gsl_bst_rb_table rb_table;
typedef gsl_bst_rb_traverser rb_traverser;

enum rb_color
{
  RB_BLACK, /* black */
  RB_RED    /* red */
};

#ifndef RB_MAX_HEIGHT
#define RB_MAX_HEIGHT GSL_BST_RB_MAX_HEIGHT
#endif

/* tree functions */
static int rb_init(const gsl_bst_allocator * allocator, gsl_bst_cmp_function * compare,
                   void * params, void * vtable);
static size_t rb_nodes (const void * vtable);
static int rb_empty (void * vtable);
static void ** rb_probe (void * item, rb_table * table);
static void * rb_insert (void * item, void * vtable);
static void * rb_find (const void * item, const void * vtable);
static void * rb_remove (const void * item, void * vtable);

/* traverser functions */
static int rb_t_init (void * vtrav, const void * vtable);
static void * rb_t_first (void * vtrav, const void * vtable);
static void * rb_t_last (void * vtrav, const void * vtable);
static void * rb_t_find (const void * item, void * vtrav, const void * vtable);
static void * rb_t_insert (void * item, void * vtrav, void * vtable);
static void * rb_t_copy (void * vtrav, const void * vsrc);
static void * rb_t_next (void * vtrav);
static void * rb_t_prev (void * vtrav);
static void * rb_t_cur (const void * vtrav);
static void * rb_t_replace (void * vtrav, void * new_item);
static void rb_trav_refresh (rb_traverser *trav);

static int
rb_init(const gsl_bst_allocator * allocator, gsl_bst_cmp_function * compare,
        void * params, void * vtable)
{
  rb_table * table = (rb_table *) vtable;

  table->rb_alloc = allocator;
  table->rb_compare = compare;
  table->rb_param = params;

  table->rb_root = NULL;
  table->rb_count = 0;
  table->rb_generation = 0;

  return GSL_SUCCESS;
}

static size_t
rb_nodes (const void * vtable)
{
  const rb_table * table = (const rb_table *) vtable;
  return table->rb_count;
}

/* empty tree (delete all nodes) but do not free the tree itself */
static int
rb_empty (void * vtable)
{
  rb_table * table = (rb_table *) vtable;
  rb_node *p, *q;

  for (p = table->rb_root; p != NULL; p = q)
    {
      if (p->rb_link[0] == NULL)
        {
          q = p->rb_link[1];
          table->rb_alloc->free (p, table->rb_param);
        }
      else
        {
          q = p->rb_link[0];
          p->rb_link[0] = q->rb_link[1];
          q->rb_link[1] = p;
        }
    }

  table->rb_root = NULL;
  table->rb_count = 0;
  table->rb_generation = 0;

  return GSL_SUCCESS;
}

/* Inserts |item| into |table| and returns a pointer to |item|'s address.
   If a duplicate item is found in the tree,
   returns a pointer to the duplicate without inserting |item|.
   Returns |NULL| in case of memory allocation failure. */

static void **
rb_probe (void * item, rb_table * table)
{
  rb_node *pa[RB_MAX_HEIGHT];      /* nodes on stack */
  unsigned char da[RB_MAX_HEIGHT]; /* directions moved from stack nodes */
  int k;                           /* stack height */

  rb_node *p; /* traverses tree looking for insertion point */
  rb_node *n; /* newly inserted node */

  pa[0] = (rb_node *) &table->rb_root;
  da[0] = 0;
  k = 1;
  for (p = table->rb_root; p != NULL; p = p->rb_link[da[k - 1]])
    {
      int cmp = table->rb_compare (item, p->rb_data, table->rb_param);
      if (cmp == 0)
        return &p->rb_data;

      pa[k] = p;
      da[k++] = cmp > 0;
    }

  n = pa[k - 1]->rb_link[da[k - 1]] =
    table->rb_alloc->alloc (sizeof *n, table->rb_param);
  if (n == NULL)
    return NULL;

  n->rb_data = item;
  n->rb_link[0] = n->rb_link[1] = NULL;
  n->rb_color = RB_RED;
  table->rb_count++;
  table->rb_generation++;

  while (k >= 3 && pa[k - 1]->rb_color == RB_RED)
    {
      if (da[k - 2] == 0)
        {
          rb_node *y = pa[k - 2]->rb_link[1];
          if (y != NULL && y->rb_color == RB_RED)
            {
              pa[k - 1]->rb_color = y->rb_color = RB_BLACK;
              pa[k - 2]->rb_color = RB_RED;
              k -= 2;
            }
          else
            {
              rb_node *x;

              if (da[k - 1] == 0)
                y = pa[k - 1];
              else
                {
                  x = pa[k - 1];
                  y = x->rb_link[1];
                  x->rb_link[1] = y->rb_link[0];
                  y->rb_link[0] = x;
                  pa[k - 2]->rb_link[0] = y;
                }

              x = pa[k - 2];
              x->rb_color = RB_RED;
              y->rb_color = RB_BLACK;

              x->rb_link[0] = y->rb_link[1];
              y->rb_link[1] = x;
              pa[k - 3]->rb_link[da[k - 3]] = y;
              break;
            }
        }
      else
        {
          rb_node *y = pa[k - 2]->rb_link[0];
          if (y != NULL && y->rb_color == RB_RED)
            {
              pa[k - 1]->rb_color = y->rb_color = RB_BLACK;
              pa[k - 2]->rb_color = RB_RED;
              k -= 2;
            }
          else
            {
              rb_node *x;

              if (da[k - 1] == 1)
                y = pa[k - 1];
              else
                {
                  x = pa[k - 1];
                  y = x->rb_link[0];
                  x->rb_link[0] = y->rb_link[1];
                  y->rb_link[1] = x;
                  pa[k - 2]->rb_link[1] = y;
                }

              x = pa[k - 2];
              x->rb_color = RB_RED;
              y->rb_color = RB_BLACK;

              x->rb_link[1] = y->rb_link[0];
              y->rb_link[0] = x;
              pa[k - 3]->rb_link[da[k - 3]] = y;
              break;
            }
        }
    }

  table->rb_root->rb_color = RB_BLACK;

  return &n->rb_data;
}

/* Inserts |item| into |table|.
   Returns |NULL| if |item| was successfully inserted
   or if a memory allocation error occurred.
   Otherwise, returns the duplicate item. */

static void *
rb_insert (void * item, void * vtable)
{
  void **p = rb_probe (item, vtable);
  return p == NULL || *p == item ? NULL : *p;
}

/* Search |table| for an item matching |item|, and return it if found.
   Otherwise return |NULL|. */
static void *
rb_find (const void * item, const void * vtable)
{
  const rb_table * table = (const rb_table *) vtable;
  const rb_node *p;

  for (p = table->rb_root; p != NULL; )
    {
      int cmp = table->rb_compare (item, p->rb_data, table->rb_param);

      if (cmp < 0)
        p = p->rb_link[0];
      else if (cmp > 0)
        p = p->rb_link[1];
      else /* |cmp == 0| */
        return p->rb_data;
    }

  return NULL;
}

#if 0 /*XXX*/

/* Inserts |item| into |table|, replacing any duplicate item.
   Returns |NULL| if |item| was inserted without replacing a duplicate,
   or if a memory allocation error occurred.
   Otherwise, returns the item that was replaced. */
void *
rb_replace (struct rb_table *table, void *item)
{
  void **p = rb_probe (table, item);
  if (p == NULL || *p == item)
    return NULL;
  else
    {
      void *r = *p;
      *p = item;
      return r;
    }
}

#endif

/* Deletes from |table| and returns an item matching |item|.
   Returns a null pointer if no matching item found. */

static void *
rb_remove (const void * item, void * vtable)
{
  rb_table * table = (rb_table *) vtable;
  rb_node *pa[RB_MAX_HEIGHT];      /* nodes on stack */
  unsigned char da[RB_MAX_HEIGHT]; /* directions moved from stack nodes */
  int k;                           /* stack height */

  rb_node *p;    /* the node to delete, or a node part way to it */
  int cmp;       /* result of comparison between |item| and |p| */

  k = 0;
  p = (rb_node *) &table->rb_root;
  for (cmp = -1; cmp != 0;
       cmp = table->rb_compare (item, p->rb_data, table->rb_param))
    {
      int dir = cmp > 0;

      pa[k] = p;
      da[k++] = dir;

      p = p->rb_link[dir];
      if (p == NULL)
        return NULL;
    }
  item = p->rb_data;

  if (p->rb_link[1] == NULL)
    pa[k - 1]->rb_link[da[k - 1]] = p->rb_link[0];
  else
    {
      enum rb_color t;
      rb_node *r = p->rb_link[1];

      if (r->rb_link[0] == NULL)
        {
          r->rb_link[0] = p->rb_link[0];
          t = r->rb_color;
          r->rb_color = p->rb_color;
          p->rb_color = t;
          pa[k - 1]->rb_link[da[k - 1]] = r;
          da[k] = 1;
          pa[k++] = r;
        }
      else
        {
          rb_node *s;
          int j = k++;

          for (;;)
            {
              da[k] = 0;
              pa[k++] = r;
              s = r->rb_link[0];
              if (s->rb_link[0] == NULL)
                break;

              r = s;
            }

          da[j] = 1;
          pa[j] = s;
          pa[j - 1]->rb_link[da[j - 1]] = s;

          s->rb_link[0] = p->rb_link[0];
          r->rb_link[0] = s->rb_link[1];
          s->rb_link[1] = p->rb_link[1];

          t = s->rb_color;
          s->rb_color = p->rb_color;
          p->rb_color = t;
        }
    }

  if (p->rb_color == RB_BLACK)
    {
      for (;;)
        {
          rb_node *x = pa[k - 1]->rb_link[da[k - 1]];
          if (x != NULL && x->rb_color == RB_RED)
            {
              x->rb_color = RB_BLACK;
              break;
            }
          if (k < 2)
            break;

          if (da[k - 1] == 0)
            {
              rb_node *w = pa[k - 1]->rb_link[1];

              if (w->rb_color == RB_RED)
                {
                  w->rb_color = RB_BLACK;
                  pa[k - 1]->rb_color = RB_RED;

                  pa[k - 1]->rb_link[1] = w->rb_link[0];
                  w->rb_link[0] = pa[k - 1];
                  pa[k - 2]->rb_link[da[k - 2]] = w;

                  pa[k] = pa[k - 1];
                  da[k] = 0;
                  pa[k - 1] = w;
                  k++;

                  w = pa[k - 1]->rb_link[1];
                }

              if ((w->rb_link[0] == NULL
                   || w->rb_link[0]->rb_color == RB_BLACK)
                  && (w->rb_link[1] == NULL
                      || w->rb_link[1]->rb_color == RB_BLACK))
                w->rb_color = RB_RED;
              else
                {
                  if (w->rb_link[1] == NULL
                      || w->rb_link[1]->rb_color == RB_BLACK)
                    {
                      rb_node *y = w->rb_link[0];
                      y->rb_color = RB_BLACK;
                      w->rb_color = RB_RED;
                      w->rb_link[0] = y->rb_link[1];
                      y->rb_link[1] = w;
                      w = pa[k - 1]->rb_link[1] = y;
                    }

                  w->rb_color = pa[k - 1]->rb_color;
                  pa[k - 1]->rb_color = RB_BLACK;
                  w->rb_link[1]->rb_color = RB_BLACK;

                  pa[k - 1]->rb_link[1] = w->rb_link[0];
                  w->rb_link[0] = pa[k - 1];
                  pa[k - 2]->rb_link[da[k - 2]] = w;
                  break;
                }
            }
          else
            {
              rb_node *w = pa[k - 1]->rb_link[0];

              if (w->rb_color == RB_RED)
                {
                  w->rb_color = RB_BLACK;
                  pa[k - 1]->rb_color = RB_RED;

                  pa[k - 1]->rb_link[0] = w->rb_link[1];
                  w->rb_link[1] = pa[k - 1];
                  pa[k - 2]->rb_link[da[k - 2]] = w;

                  pa[k] = pa[k - 1];
                  da[k] = 1;
                  pa[k - 1] = w;
                  k++;

                  w = pa[k - 1]->rb_link[0];
                }

              if ((w->rb_link[0] == NULL
                   || w->rb_link[0]->rb_color == RB_BLACK)
                  && (w->rb_link[1] == NULL
                      || w->rb_link[1]->rb_color == RB_BLACK))
                w->rb_color = RB_RED;
              else
                {
                  if (w->rb_link[0] == NULL
                      || w->rb_link[0]->rb_color == RB_BLACK)
                    {
                      rb_node *y = w->rb_link[1];
                      y->rb_color = RB_BLACK;
                      w->rb_color = RB_RED;
                      w->rb_link[1] = y->rb_link[0];
                      y->rb_link[0] = w;
                      w = pa[k - 1]->rb_link[0] = y;
                    }

                  w->rb_color = pa[k - 1]->rb_color;
                  pa[k - 1]->rb_color = RB_BLACK;
                  w->rb_link[0]->rb_color = RB_BLACK;

                  pa[k - 1]->rb_link[0] = w->rb_link[1];
                  w->rb_link[1] = pa[k - 1];
                  pa[k - 2]->rb_link[da[k - 2]] = w;
                  break;
                }
            }

          k--;
        }

    }

  table->rb_alloc->free (p, table->rb_param);
  table->rb_count--;
  table->rb_generation++;

  return (void *) item;
}

/* Initializes |trav| for use with |tree| and selects the null node. */
static int
rb_t_init (void * vtrav, const void * vtable)
{
  rb_traverser * trav = (rb_traverser *) vtrav;
  const rb_table * table = (const rb_table *) vtable;

  trav->rb_table = table;
  trav->rb_node = NULL;
  trav->rb_height = 0;
  trav->rb_generation = table->rb_generation;

  return GSL_SUCCESS;
}

/* Initializes |trav| for |table|
   and selects and returns a pointer to its least-valued item.
   Returns |NULL| if |table| contains no nodes. */
static void *
rb_t_first (void * vtrav, const void * vtable)
{
  const rb_table * table = (const rb_table *) vtable;
  rb_traverser * trav = (rb_traverser *) vtrav;
  rb_node *x;

  trav->rb_table = table;
  trav->rb_height = 0;
  trav->rb_generation = table->rb_generation;

  x = table->rb_root;
  if (x != NULL)
    {
      while (x->rb_link[0] != NULL)
        {
          if (trav->rb_height >= RB_MAX_HEIGHT)
            {
              GSL_ERROR_NULL ("traverser height exceeds maximum", GSL_ETABLE);
            }

          trav->rb_stack[trav->rb_height++] = x;
          x = x->rb_link[0];
        }
    }

  trav->rb_node = x;

  return x != NULL ? x->rb_data : NULL;
}

/* Initializes |trav| for |table|
   and selects and returns a pointer to its greatest-valued item.
   Returns |NULL| if |table| contains no nodes. */
static void *
rb_t_last (void * vtrav, const void * vtable)
{
  const rb_table * table = (const rb_table *) vtable;
  rb_traverser * trav = (rb_traverser *) vtrav;
  rb_node *x;

  trav->rb_table = table;
  trav->rb_height = 0;
  trav->rb_generation = table->rb_generation;

  x = table->rb_root;
  if (x != NULL)
    {
      while (x->rb_link[1] != NULL)
        {
          if (trav->rb_height >= RB_MAX_HEIGHT)
            {
              GSL_ERROR_NULL ("traverser height exceeds maximum", GSL_ETABLE);
            }

          trav->rb_stack[trav->rb_height++] = x;
          x = x->rb_link[1];
        }
    }

  trav->rb_node = x;

  return x != NULL ? x->rb_data : NULL;
}

/* Searches for |item| in |table|.
   If found, initializes |trav| to the item found and returns the item
   as well.
   If there is no matching item, initializes |trav| to the null item
   and returns |NULL|. */
static void *
rb_t_find (const void * item, void * vtrav, const void * vtable)
{
  const rb_table * table = (const rb_table *) vtable;
  rb_traverser * trav = (rb_traverser *) vtrav;
  rb_node *p, *q;

  trav->rb_table = table;
  trav->rb_height = 0;
  trav->rb_generation = table->rb_generation;
  for (p = table->rb_root; p != NULL; p = q)
    {
      int cmp = table->rb_compare (item, p->rb_data, table->rb_param);

      if (cmp < 0)
        q = p->rb_link[0];
      else if (cmp > 0)
        q = p->rb_link[1];
      else /* |cmp == 0| */
        {
          trav->rb_node = p;
          return p->rb_data;
        }

      if (trav->rb_height >= RB_MAX_HEIGHT)
        {
          GSL_ERROR_NULL ("traverser height exceeds maximum", GSL_ETABLE);
        }

      trav->rb_stack[trav->rb_height++] = p;
    }

  trav->rb_height = 0;
  trav->rb_node = NULL;

  return NULL;
}

/* Attempts to insert |item| into |table|.
   If |item| is inserted successfully, it is returned and |trav| is
   initialized to its location.
   If a duplicate is found, it is returned and |trav| is initialized to
   its location.  No replacement of the item occurs.
   If a memory allocation failure occurs, |NULL| is returned and |trav|
   is initialized to the null item. */
static void *
rb_t_insert (void * item, void * vtrav, void * vtable)
{
  rb_table * table = (rb_table *) vtable;
  rb_traverser * trav = (rb_traverser *) vtrav;
  void **p;

  p = rb_probe (item, table);
  if (p != NULL)
    {
      trav->rb_table = table;
      trav->rb_node = ((rb_node *) ((char *) p - offsetof (rb_node, rb_data)));
      trav->rb_generation = table->rb_generation - 1;
      return *p;
    }
  else
    {
      rb_t_init (vtrav, vtable);
      return NULL;
    }
}

/* Initializes |trav| to have the same current node as |src|. */
static void *
rb_t_copy (void * vtrav, const void * vsrc)
{
  const rb_traverser * src = (const rb_traverser *) vsrc;
  rb_traverser * trav = (rb_traverser *) vtrav;

  if (trav != src)
    {
      trav->rb_table = src->rb_table;
      trav->rb_node = src->rb_node;
      trav->rb_generation = src->rb_generation;
      if (trav->rb_generation == trav->rb_table->rb_generation)
        {
          trav->rb_height = src->rb_height;
          memcpy (trav->rb_stack, (const void *) src->rb_stack,
                  sizeof *trav->rb_stack * trav->rb_height);
        }
    }

  return trav->rb_node != NULL ? trav->rb_node->rb_data : NULL;
}

/* Returns the next data item in inorder
   within the tree being traversed with |trav|,
   or if there are no more data items returns |NULL|. */
static void *
rb_t_next (void * vtrav)
{
  rb_traverser * trav = (rb_traverser *) vtrav;
  rb_node *x;

  if (trav->rb_generation != trav->rb_table->rb_generation)
    rb_trav_refresh (trav);

  x = trav->rb_node;
  if (x == NULL)
    {
      return rb_t_first (vtrav, trav->rb_table);
    }
  else if (x->rb_link[1] != NULL)
    {
      if (trav->rb_height >= RB_MAX_HEIGHT)
        {
          GSL_ERROR_NULL ("traverser height exceeds maximum", GSL_ETABLE);
        }

      trav->rb_stack[trav->rb_height++] = x;
      x = x->rb_link[1];

      while (x->rb_link[0] != NULL)
        {
          if (trav->rb_height >= RB_MAX_HEIGHT)
            {
              GSL_ERROR_NULL ("traverser height exceeds maximum", GSL_ETABLE);
            }

          trav->rb_stack[trav->rb_height++] = x;
          x = x->rb_link[0];
        }
    }
  else
    {
      rb_node *y;

      do
        {
          if (trav->rb_height == 0)
            {
              trav->rb_node = NULL;
              return NULL;
            }

          y = x;
          x = trav->rb_stack[--trav->rb_height];
        }
      while (y == x->rb_link[1]);
    }

  trav->rb_node = x;

  return x->rb_data;
}

/* Returns the previous data item in inorder
   within the tree being traversed with |trav|,
   or if there are no more data items returns |NULL|. */
static void *
rb_t_prev (void * vtrav)
{
  rb_traverser * trav = (rb_traverser *) vtrav;
  rb_node *x;

  if (trav->rb_generation != trav->rb_table->rb_generation)
    rb_trav_refresh (trav);

  x = trav->rb_node;
  if (x == NULL)
    {
      return rb_t_last (vtrav, trav->rb_table);
    }
  else if (x->rb_link[0] != NULL)
    {
      if (trav->rb_height >= RB_MAX_HEIGHT)
        {
          GSL_ERROR_NULL ("traverser height exceeds maximum", GSL_ETABLE);
        }

      trav->rb_stack[trav->rb_height++] = x;
      x = x->rb_link[0];

      while (x->rb_link[1] != NULL)
        {
          if (trav->rb_height >= RB_MAX_HEIGHT)
            {
              GSL_ERROR_NULL ("traverser height exceeds maximum", GSL_ETABLE);
            }

          trav->rb_stack[trav->rb_height++] = x;
          x = x->rb_link[1];
        }
    }
  else
    {
      rb_node *y;

      do
        {
          if (trav->rb_height == 0)
            {
              trav->rb_node = NULL;
              return NULL;
            }

          y = x;
          x = trav->rb_stack[--trav->rb_height];
        }
      while (y == x->rb_link[0]);
    }

  trav->rb_node = x;

  return x->rb_data;
}

/* Returns |trav|'s current item. */
static void *
rb_t_cur (const void * vtrav)
{
  const rb_traverser * trav = (const rb_traverser *) vtrav;
  return trav->rb_node != NULL ? trav->rb_node->rb_data : NULL;
}

/* Replaces the current item in |trav| by |new| and returns the item replaced.
   |trav| must not have the null item selected.
   The new item must not upset the ordering of the tree. */
static void *
rb_t_replace (void * vtrav, void * new_item)
{
  rb_traverser * trav = (rb_traverser *) vtrav;
  void *old;

  old = trav->rb_node->rb_data;
  trav->rb_node->rb_data = new_item;

  return old;
}

#if 0 /*XXX*/
/* Destroys |new| with |rb_destroy (new, destroy)|,
   first setting right links of nodes in |stack| within |new|
   to null pointers to avoid touching uninitialized data. */
static void
copy_error_recovery (rb_node **stack, int height,
                     struct rb_table *new, rb_item_func *destroy)
{
  assert (stack != NULL && height >= 0 && new != NULL);

  for (; height > 2; height -= 2)
    stack[height - 1]->rb_link[1] = NULL;
  rb_destroy (new, destroy);
}

/* Copies |org| to a newly created tree, which is returned.
   If |copy != NULL|, each data item in |org| is first passed to |copy|,
   and the return values are inserted into the tree,
   with |NULL| return values taken as indications of failure.
   On failure, destroys the partially created new tree,
   applying |destroy|, if non-null, to each item in the new tree so far,
   and returns |NULL|.
   If |allocator != NULL|, it is used for allocation in the new tree.
   Otherwise, the same allocator used for |org| is used. */
struct rb_table *
rb_copy (const struct rb_table *org, rb_copy_func *copy,
          rb_item_func *destroy, struct libavl_allocator *allocator)
{
  rb_node *stack[2 * (RB_MAX_HEIGHT + 1)];
  int height = 0;

  struct rb_table *new;
  const rb_node *x;
  rb_node *y;

  assert (org != NULL);
  new = rb_create (org->rb_compare, org->rb_param,
                    allocator != NULL ? allocator : org->rb_alloc);
  if (new == NULL)
    return NULL;
  new->rb_count = org->rb_count;
  if (new->rb_count == 0)
    return new;

  x = (const rb_node *) &org->rb_root;
  y = (rb_node *) &new->rb_root;
  for (;;)
    {
      while (x->rb_link[0] != NULL)
        {
          assert (height < 2 * (RB_MAX_HEIGHT + 1));

          y->rb_link[0] =
            new->rb_alloc->libavl_malloc (new->rb_alloc,
                                           sizeof *y->rb_link[0]);
          if (y->rb_link[0] == NULL)
            {
              if (y != (rb_node *) &new->rb_root)
                {
                  y->rb_data = NULL;
                  y->rb_link[1] = NULL;
                }

              copy_error_recovery (stack, height, new, destroy);
              return NULL;
            }

          stack[height++] = (rb_node *) x;
          stack[height++] = y;
          x = x->rb_link[0];
          y = y->rb_link[0];
        }
      y->rb_link[0] = NULL;

      for (;;)
        {
          y->rb_color = x->rb_color;
          if (copy == NULL)
            y->rb_data = x->rb_data;
          else
            {
              y->rb_data = copy (x->rb_data, org->rb_param);
              if (y->rb_data == NULL)
                {
                  y->rb_link[1] = NULL;
                  copy_error_recovery (stack, height, new, destroy);
                  return NULL;
                }
            }

          if (x->rb_link[1] != NULL)
            {
              y->rb_link[1] =
                new->rb_alloc->libavl_malloc (new->rb_alloc,
                                               sizeof *y->rb_link[1]);
              if (y->rb_link[1] == NULL)
                {
                  copy_error_recovery (stack, height, new, destroy);
                  return NULL;
                }

              x = x->rb_link[1];
              y = y->rb_link[1];
              break;
            }
          else
            y->rb_link[1] = NULL;

          if (height <= 2)
            return new;

          y = stack[--height];
          x = stack[--height];
        }
    }
}

#endif

/* Refreshes the stack of parent pointers in |trav|
   and updates its generation number. */
static void
rb_trav_refresh (rb_traverser *trav)
{
  trav->rb_generation = trav->rb_table->rb_generation;

  if (trav->rb_node != NULL)
    {
      gsl_bst_cmp_function *cmp = trav->rb_table->rb_compare;
      void *param = trav->rb_table->rb_param;
      rb_node *node = trav->rb_node;
      rb_node *i;

      trav->rb_height = 0;
      for (i = trav->rb_table->rb_root; i != node; )
        {
          if (trav->rb_height >= RB_MAX_HEIGHT)
            {
              GSL_ERROR_VOID ("traverser height exceeds maximum", GSL_ETABLE);
            }

          trav->rb_stack[trav->rb_height++] = i;
          i = i->rb_link[cmp (node->rb_data, i->rb_data, param) > 0];
        }
    }
}

static const gsl_bst_type rb_tree_type =
{
  "red-black",
  sizeof(rb_node),
  rb_init,
  rb_nodes,
  rb_insert,
  rb_find,
  rb_remove,
  rb_empty,

  rb_t_init,
  rb_t_first,
  rb_t_last,
  rb_t_find,
  rb_t_insert,
  rb_t_copy,
  rb_t_next,
  rb_t_prev,
  rb_t_cur,
  rb_t_replace
};

const gsl_bst_type * gsl_bst_rb = &rb_tree_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense