source: observatorio/simulacion/ModuloDinamico/simulator.C

/*
  Copyright (C) 2014
  Alejandro Mujica (amujica@cenditel.gob.ve)
  José Ruiz (jruiz@cenditel.gob.ve)
  Julie Vera (jvera@cenditel.gob.ve)
 
  CENDITEL Fundación Centro Nacional de Desarrollo e Investigación en
  Tecnologías Libres
 
  Este programa es software libre; Usted puede usarlo bajo los términos de la
  licencia de software GPL versión 2.0 de la Free Software Foundation.
 
  Este programa se distribuye con la esperanza de que sea útil, pero SIN
  NINGUNA GARANTÍA; tampoco las implícitas garantías de MERCANTILIDAD o
  ADECUACIÓN A UN PROPÓSITO PARTICULAR.
  Consulte la licencia GPL para más detalles. Usted debe recibir una copia
  de la GPL junto con este programa; si no, escriba a la Free Software
  Foundation Inc. 51 Franklin Street,5 Piso, Boston, MA 02110-1301, USA.
*/

/*
  Autor:             Alejandro J. Mujica
  Fecha de creación: 14/11/2014
  Este archivo contiene la implementación de la función que ejecuta la
  simulación de una red productiva con descrita en un archivo xml. Además
  contiene funciones de modularización de la simulación. Esta función
  fue creada con la pretensión de crear el módulo para Python.
*/

# include <simulator.H>

// Cuenta la cantidad de arcos de salida de un nodo de un grafo.
size_t count_output_arcs(Graph & network, Graph::Node * node)
{
  size_t counter = 0;

  for (Graph::Node_Arc_Iterator it(node); it.has_current(); it.next())
    {
      Graph::Arc * arc = it.get_current();

      Graph::Node * src = network.get_src_node(arc);

      // Cuento solamente los arcos para los que el nodo dado es fuente.
      if (src == node)
        ++counter;
    }

  return counter;
}

/* Crea el estado inicial de simulación para un nodo de tipo producto (productos
   pertenecientes a UE registradas en SIGESIC).
*/
void init_product(Graph & network, Graph::Node * ptr_node)
{
  Exogenous_Variables & exo_var = Exogenous_Variables::get_instance();

  Good * ptr_good = ptr_node->get_info().get();

  size_t num_output_arcs = count_output_arcs(network, ptr_node);

  Product * ptr_product = static_cast<Product *>(ptr_good);

  Product_Simulation_Attributes state_0;

  state_0.state_number = 0;
  state_0.production = ptr_product->get_production();
  state_0.stock = 0.0;
  state_0.internal_sales = ptr_product->get_internal_sales();
  state_0.external_sales = num_output_arcs == 0 ?
                           ptr_product->get_production() *
                           ptr_product->get_price() :
                           ptr_product->get_external_sales();
  state_0.internal_requested_quantity =
    state_0.internal_sales / ptr_product->get_price();
  state_0.external_requested_quantity =
    state_0.external_sales / ptr_product->get_price();

  state_0.num_employees = ptr_product->get_num_employees();
  state_0.daily_wage = exo_var.get_average_salary() / DAYS_IN_A_YEAR;
  state_0.other_wage = exo_var.get_average_salary() * WAGE_PROPORTION;
  state_0.other_daily_wage = state_0.other_wage / DAYS_IN_A_YEAR;
  state_0.daily_feeding_coupon = exo_var.UT_at(0) * UT_PROPORTION;
  state_0.feeding_coupon = state_0.daily_feeding_coupon * DAYS_IN_A_YEAR;
  state_0.integral_wage = exo_var.get_average_salary() +
                          state_0.other_wage + state_0.feeding_coupon;

  state_0.administrative_staff_cost =
    state_0.integral_wage * ptr_product->get_num_administrative_staff();

  state_0.labor_cost = state_0.integral_wage * state_0.num_employees;

  state_0.input_cost = 0.0;

  for (Graph::Node_Arc_Iterator it(ptr_node); it.has_current(); it.next())
    {
      Graph::Arc * ptr_arc = it.get_current();

      Graph::Node * ptr_src_node = network.get_src_node(ptr_arc);

      if (ptr_src_node == ptr_node)
        continue;

      IP_Relationship & ip_rel = ptr_arc->get_info();

      state_0.input_cost += ip_rel.get_purchase_price() *
                            ip_rel.get_bought_quantity();
    }

  state_0.price = ptr_product->get_price();

  state_0.income = state_0.production * state_0.price;

  real partial_cost = state_0.administrative_staff_cost +
                      state_0.labor_cost + state_0.input_cost;

  state_0.other_cost = std::abs(state_0.income - partial_cost);

  state_0.total_cost = partial_cost + state_0.other_cost;

  state_0.unitarian_administrative_staff_cost =
    state_0.administrative_staff_cost / state_0.production;

  state_0.unitarian_labor_cost = state_0.labor_cost / state_0.production;

  state_0.input_cost_per_unit = state_0.input_cost / state_0.production;

  state_0.other_unitarian_cost = state_0.other_cost / state_0.production;

  state_0.economic_status = state_0.income - state_0.total_cost;

  ptr_product->add_simulation_attributes(state_0);
}

/* Crea el estado inicial de simulación para un nodo de tipo insumo (productos
   pertenecientes a UE nacionales no registradas en SIGESIC).
*/
void init_input(Graph::Node * ptr_node)
{
  Good * ptr_good = ptr_node->get_info().get();

  Input * ptr_input = static_cast<Input *>(ptr_good);

  Input_Simulation_Attributes state_0;

  state_0.state_number = 0;

  state_0.production = (ptr_input->get_internal_sales() +
                        ptr_input->get_external_sales()) /
                        ptr_input->get_price();
  state_0.stock = 0.0;
  state_0.internal_sales = ptr_input->get_internal_sales();
  state_0.external_sales = ptr_input->get_external_sales();
  state_0.internal_requested_quantity =
    state_0.internal_sales / ptr_input->get_price();
  state_0.external_requested_quantity =
    state_0.external_sales / ptr_input->get_price();
  state_0.price = ptr_input->get_price();

  ptr_input->add_simulation_attributes(state_0);  
}

/* Crea el estado inicial de simulación para un nodo de tipo producto importado
   (productos pertenecientes a UE extranjeras).
*/
void init_imported_product(Graph::Node * ptr_node)
{
  Good * ptr_good = ptr_node->get_info().get();

  Imported_Product * ptr_imported_product =
    static_cast<Imported_Product *>(ptr_good);

  Imported_Product_Simulation_Attributes state_0;

  state_0.state_number = 0;
  state_0.imports = ptr_imported_product->get_imports();
  state_0.stock = 0.0;
  state_0.sales = ptr_imported_product->get_sale();
  state_0.requested_quantity = ptr_imported_product->get_requested_quantity();
  state_0.price = ptr_imported_product->get_price();

  ptr_imported_product->add_simulation_attributes(state_0);
}

// Asigna los estados iniciales de simulación para cada nodo y arco del grafo
void init_network(Graph & network)
{
  Exogenous_Variables & exo_var = Exogenous_Variables::get_instance();

  exo_var.set_average_salary(exo_var.get_average_salary() * NUM_MONTHS);

  for (Graph::Node_Iterator it(network); it.has_current(); it.next())
    {
      Graph::Node * ptr_node = it.get_current();

      switch (ptr_node->get_info()->get_type())
        {
        case PRODUCT:
          init_product(network, ptr_node);
          break;

        case INPUT:
          init_input(ptr_node);
          break;

        case IMPORTED_PRODUCT:
          init_imported_product(ptr_node);
          break;
        }
    }

  for (Graph::Arc_Iterator it(network); it.has_current(); it.next())
    {
      Graph::Arc * ptr_arc = it.get_current();

      IP_Relationship & ip_rel = ptr_arc->get_info();

      IP_Simulation_Attributes state_0;

      state_0.state_number = 0;

      state_0.purchase_price = ip_rel.get_purchase_price();

      state_0.bought_quantity = ip_rel.get_bought_quantity();

      ip_rel.add_simulation_attributes(state_0);
    }
}

/* Una red productiva es realmente un digrafo, la red fue modelada con grafo
   para poder efectuar navegación hacia los nodos de atrás. Para esta
   construcción se asume una convención de que si el nodo s va detrás de t,
   en un digrafo el arco sería s -> t y en el grafo se respeta ese orden
   teniendo los arcos s -- t.

   Esta función toma el grafo y lo convierte en un digrafo para aplicar sobre
   éste el algoritmo de construcción de rangos topológicos.

   En los nodos del digrafo se almacenan los punteros a los nodos del grafo y
   análogamente se hace con los arcos, pues los rangos topológicos son
   requeridos para hacer los recorridos por el grafo por cada nivel topológico
   desde el más alto hasta el más bajo, pero finalmente el proceso de simulación
   se efectúa sobre los nodos del grafo (y no los del digrafo) porque se
   requiere navegación hacia atrás.
*/
void graph_to_digraph(Graph & graph, Digraph & digraph)
{
  Map<Graph::Node *, Digraph::Node *> map;

  for (Graph::Node_Iterator it(graph); it.has_current(); it.next())
    {
      Graph::Node * p = it.get_current();

      Digraph::Node * q = digraph.insert_node(p);

      map.insert(p, q);
    }

  for (Graph::Arc_Iterator it(graph); it.has_current(); it.next())
    {
      Graph::Arc * a = it.get_current();

      Graph::Node * gs = graph.get_src_node(a);
      Graph::Node * gt = graph.get_tgt_node(a);

      Digraph::Node * ds = map.find(gs);
      Digraph::Node * dt = map.find(gt);

      digraph.insert_arc(ds, dt, a);
    }
}

// Libera las listas que almacenan cada nivel topológico.
void free_ranks(Rank_Type & ranks)
{
  for (Rank_Type::Iterator it(ranks); it.has_current(); it.next())
    delete it.get_current();
}


/* Trunca la producción en el nodo según lo que permitan los insumos, la
   cantidad requeridad de cada insumo se añade a la demanda interna del insumo,
   una vez  que se tiene la cantidad posible de producción, se actualiza la
   cantidad comprada en la relación a la cantidad que realmente necesita, lo que
   sobre de cada arco se añade al stock del proveedor.
*/
void truncate_production_by_inputs(Graph & network, Graph::Node * node,
                                   const size_t & t)
{
  Good * ptr_good = node->get_info().get();

  // Inicializo la máxima producción posible en tendiendo a infinito
  real max_possible_production = std::numeric_limits<real>::max();

  // Producción deseada por el nodo.
  real & production = ptr_good->production_at(t);

  // Recorro los arcos de entrada
  for (Graph::Node_Arc_Iterator it(node); it.has_current(); it.next())
    {
      Graph::Arc * arc = it.get_current();

      Graph::Node * src_node = network.get_src_node(arc);

      // Ignoro los arcos de salida
      if (src_node == node)
        continue;

      IP_Relationship & ip_rel = arc->get_info();

      IP_Simulation_Attributes & ip_sim_attr =
        ip_rel.get_simulation_attributes_at(t);

      // Cantidad requerida del proveedor.
      real requested_quantity = production * ip_rel.get_requested_quantity();

      // Incrementa la demanda interna al proveedor
      src_node->get_info()->internal_requested_quantity_at(t) +=
        requested_quantity;

      I(ip_sim_attr.bought_quantity > 0.0);

      real possible_production = ip_rel.get_requested_quantity() *
                                 ip_sim_attr.bought_quantity;

      if (possible_production == 0.0)
        {
          std::cout << "RQ: " << ip_rel.get_requested_quantity() << " - "
                    << "BQ: " << ip_sim_attr.bought_quantity << '\n';
        }

      // Se actualiza la máxima producción posible.
      if (possible_production < max_possible_production)
        max_possible_production = possible_production;
    }

  /* Aquí trunco la producción con el mínimo entre máxima producción posible y
     la producción actual.
  */
  production = std::min<real>(production, max_possible_production);

  /* Nuevamente recorro los arcos de entrada para asignar la cantidad que
     realmente comprará del insumo.
  */
  for (Graph::Node_Arc_Iterator it(node); it.has_current(); it.next())
    {
      Graph::Arc * arc = it.get_current();

      Graph::Node * src_node = network.get_src_node(arc);

      // Ignoro los arcos de salida
      if (src_node == node)
        continue;

      IP_Relationship & ip_rel = arc->get_info();

      IP_Simulation_Attributes & ip_sim_attr =
        ip_rel.get_simulation_attributes_at(t);

      // Cantidad que comprará.
      real requested_quantity = production * ip_rel.get_requested_quantity();

      /* La diferencia entre la cantidad en el arco y la cantidad que comprará
         se almacena en el stock del proveedor. cantidad en arco >= cantidad
         que comprará.
      */
      real diff = ip_sim_attr.bought_quantity - requested_quantity;

      ip_sim_attr.bought_quantity = requested_quantity;

      src_node->get_info()->stock_at(t) += diff;

      src_node->get_info()->internal_sales_at(t) += ip_sim_attr.bought_quantity;
  }
}

/* Distribuye la cantidad de producto que el nodo tenga para vender en la red
   hacia sus arcos de salida. El valor se asigna en la variable que almacena la
   cantidad comprada, sin embargo, cuando el nodo comprador 'decida' cuánto
   producto generará, actualizará estos valores a la cantidad que realmente
   necesite comprar.

   La distribución de la cantidad se hace proporcional a los montos pagados por
   cada comprador en el estado de simulación anterior.
*/
void distribute_internal_production(Graph & network, Graph::Node * node,
                                    const real & internal_quantity,
                                    const size_t & t)
{
  real total_amount = 0.0;

  for (Graph::Node_Arc_Iterator it(node); it.has_current(); it.next())
    {
      Graph::Arc * arc = it.get_current();

      Graph::Node * tgt_node = network.get_tgt_node(arc);

      if (tgt_node == node)
        continue;

      // Añado atributos de simulación para el tiempo t.
      IP_Simulation_Attributes & ip_sim_attr_t =
        arc->get_info().add_simulation_attributes(IP_Simulation_Attributes());

      ip_sim_attr_t.state_number = t;

      IP_Simulation_Attributes & ip_sim_attr_t_1 =
        arc->get_info().get_simulation_attributes_at(t - 1);

      total_amount += (ip_sim_attr_t_1.bought_quantity *
                       ip_sim_attr_t_1.purchase_price);
    }

  for (Graph::Node_Arc_Iterator it(node); it.has_current(); it.next())
    {
      Graph::Arc * arc = it.get_current();

      Graph::Node * tgt_node = network.get_tgt_node(arc);

      if (tgt_node == node)
        continue;

      IP_Simulation_Attributes & ip_sim_attr_t =
        arc->get_info().get_simulation_attributes_at(t);

      IP_Simulation_Attributes & ip_sim_attr_t_1 =
        arc->get_info().get_simulation_attributes_at(t - 1);

      real amount = ip_sim_attr_t_1.bought_quantity *
                    ip_sim_attr_t_1.purchase_price;

      I(internal_quantity > 0.0);

      I(total_amount > 0.0);

      real proportion = total_amount == 0.0 ? 0.0 : amount / total_amount;

      ip_sim_attr_t.bought_quantity = internal_quantity * proportion;

      I(ip_sim_attr_t.bought_quantity > 0.0);
    }
}

/* Función de redondeo especial.

   Un número real se puede dividir en parte entera + parte decimal, si la parte
   decimal es menor que 0.5 entonces se redondea a 0.5, de lo contrario se
   redondea a 1.

   Para el número de empleados requeridos por una UE, tener 0.5 empleados quiere
   decir que tiene un empleado contratado a medio tiempo.
*/
real special_round(const real & num)
{
  long int_part = (long) num;

  real dec_part = num - (real) int_part;

  I(dec_part < 1);

  if (dec_part == 0.0)
    return (real) int_part;
  else if (dec_part < 0.5)
    return (real) int_part + 0.5;
  else
    return (real) int_part + 1.0;
}

/* Functor para calcular la demanda externa para nodos de cualquier nivel
   topológico en la red excepto los del último nivel.
*/
struct Any_Level_Node
{
  void operator () (Product_Simulation_Attributes & sat,
                    Product_Simulation_Attributes & sat_1, gsl_rng * rng,
                    const size_t &)
  {
    sat.external_requested_quantity =
      sat_1.external_requested_quantity * RANDOM_UNIF(0.0, 2.0);
  }
};

/* Functor para calcular la demanda externa para nodos de cualquier del último
   nivel topológico en la red.
*/
struct Last_Level_Node
{
  void operator () (Product_Simulation_Attributes & sat,
                    Product_Simulation_Attributes & sat_1, gsl_rng *,
                    const size_t & t)
  {
    Exogenous_Variables & exo_var = Exogenous_Variables::get_instance();

    sat.external_requested_quantity =
      sat_1.external_requested_quantity *
      (1.0 + exo_var.rate_of_change_in_final_demand_at(t) / 100.0);
  }
};

/* Calcula la nueva producción para un nodo de tipo producto.

   La función es de tipo plantilla para pasar como parámetro el functor de
   cálculo de la demanda externa, éste cálculo cambia cuando se trata de un nodo
   del último nivel topológico respecto a lo de los otros niveles.
*/
template <class Update_External_Request>
void compute_product_production(Graph & network, Graph::Node * node,
                                const size_t & t, gsl_rng * rng)
{
  Exogenous_Variables & exo_var = Exogenous_Variables::get_instance();

  Good * ptr_good = node->get_info().get();

  Product * ptr_product = static_cast<Product *>(ptr_good);

  Product_Simulation_Attributes & sim_attr_t =
    ptr_product->add_simulation_attributes(Product_Simulation_Attributes());

  sim_attr_t.state_number = t;

  Product_Simulation_Attributes & sim_attr_t_1 =
    ptr_product->get_simulation_attributes_at(t - 1);

  real total_requested_quantity = sim_attr_t_1.internal_requested_quantity +
                                  sim_attr_t_1.external_requested_quantity;

  // Media para la nueva producción: cantidad total demandada - stock.
  real mean = total_requested_quantity - sim_attr_t_1.stock;

  // Varianza para la nueva producción.
  const real sigma = mean * exo_var.get_sigma() / 100.0;

  // La nueva producción viene dada por el valor aleatorio.
  real new_production = std::max<real>(0.0, RANDOM_NORMAL(mean, sigma));

  sim_attr_t.production =
    std::min<real>(new_production, ptr_product->get_production_capacity());

  truncate_production_by_inputs(network, node, t);

  // Cantidad de producto disponible para venta.
  real available_quantity = sim_attr_t.production + sim_attr_t_1.stock;

  // Proporción de demanda interna.
  real internal_request_proportion =
    sim_attr_t_1.internal_requested_quantity / total_requested_quantity;

  // Cantidad para repartir a los arcos.
  real internal_quantity = available_quantity * internal_request_proportion;

  distribute_internal_production(network, node, internal_quantity, t);

  // Cantidad para satisfacer demanda externa.
  real external_quantity = available_quantity - internal_quantity;

  /* Asigno a las ventas externas el valor mínimo entre demanda externa en
     (t- 1)y lo que tenga de producción para satisfacer la demanda externa.
  */
  sim_attr_t.external_sales =
    std::min<real>(external_quantity, sim_attr_t_1.external_requested_quantity);

  // Almaceno en stock lo que quede de las ventas externas.
  real diff =  external_quantity - sim_attr_t.external_sales;
  sim_attr_t.stock += diff;

  // Calculo la nueva demanda externa.
  Update_External_Request()(sim_attr_t, sim_attr_t_1, rng, t);

  /* Calcular número de personal

     cantidad de horas(t) = producción(t) * relación horas-hombre.

     cantidad de empleados = cantidad de horas(t) / (jornada * 231)

     La cantidad de empleados es un número real el cual puede ser dividido en
     su parte entera y decimal, si la parte decimal es mayor que 0 y menor que
     0.5 entonces ésta se redondea a 0.5, de lo contrario se redondea a 1.

     La parte decimal en 0.5 representa un empleado a medio tiempo.
   */
  real requested_num_hours = sim_attr_t.production *
                             ptr_product->get_relationship_manhours();

  real num_employees = requested_num_hours /
                       (ptr_product->get_workday() * WORKED_DAYS_IN_A_YEAR);

  sim_attr_t.num_employees = special_round(num_employees);
}

// Calcula la nueva producción para un nodo de tipo insumo.
void compute_input_production(Graph & network, Graph::Node * node,
                              const size_t & t, gsl_rng * rng)
{
  Exogenous_Variables & exo_var = Exogenous_Variables::get_instance();

  Good * ptr_good = node->get_info().get();

  Input * ptr_input = static_cast<Input *>(ptr_good);

  Input_Simulation_Attributes & sim_attr_t =
    ptr_input->add_simulation_attributes(Input_Simulation_Attributes());

  sim_attr_t.state_number = t;

  Input_Simulation_Attributes & sim_attr_t_1 =
    ptr_input->get_simulation_attributes_at(t - 1);

  real total_requested_quantity = sim_attr_t_1.internal_requested_quantity +
                                  sim_attr_t_1.external_requested_quantity;

  // Media para la nueva producción: cantidad total demandada - stock.
  real mean = total_requested_quantity - sim_attr_t_1.stock;

  // Varianza para la nueva producción.
  const real sigma = mean * exo_var.get_sigma() / 100.0;

  // La nueva producción viene dada por el valor aleatorio.
  sim_attr_t.production = std::max<real>(0.0, RANDOM_NORMAL(mean, sigma));

  // Cantidad de producto disponible para venta.
  real available_quantity = sim_attr_t.production + sim_attr_t_1.stock;

  // Proporción de demanda interna.
  real internal_request_proportion =
    sim_attr_t_1.internal_requested_quantity / total_requested_quantity;

  // Cantidad para repartir a los arcos.
  real internal_quantity = available_quantity * internal_request_proportion;

  distribute_internal_production(network, node, internal_quantity, t);

  // Cantidad para satisfacer demanda externa.
  real external_quantity = available_quantity - internal_quantity;

  /* Asigno a las ventas externas el valor mínimo entre demanda externa en
     (t- 1)y lo que tenga de producción para satisfacer la demanda externa.
  */
  sim_attr_t.external_sales =
    std::min<real>(external_quantity, sim_attr_t_1.external_requested_quantity);

  // Calculo la nueva demanda externa.
  sim_attr_t.external_requested_quantity =
    sim_attr_t_1.external_requested_quantity * RANDOM_UNIF(0.0, 2.0);

  // Almaceno en stock lo que quede de las ventas externas.
  real diff =  external_quantity - sim_attr_t.external_sales;

  sim_attr_t.stock += diff;
}

// Calcula la nueva producción para un nodo de tipo producto importado.
void compute_imported_product_production(Graph & network, Graph::Node * node,
                                         const size_t & t, gsl_rng * rng)
{
  Exogenous_Variables & exo_var = Exogenous_Variables::get_instance();

  Good * ptr_good = node->get_info().get();

  Imported_Product * ptr_imported_product =
    static_cast<Imported_Product *>(ptr_good);

  Imported_Product_Simulation_Attributes & sim_attr_t =
    ptr_imported_product->add_simulation_attributes(
      Imported_Product_Simulation_Attributes()
    );

  sim_attr_t.state_number = t;

  Imported_Product_Simulation_Attributes & sim_attr_t_1 =
    ptr_imported_product->get_simulation_attributes_at(t - 1);

  // Media para la nueva producción: cantidad total demandada - stock.
  real mean = sim_attr_t_1.requested_quantity - sim_attr_t_1.stock;

  // Varianza para la nueva producción.
  const real sigma = mean * exo_var.get_sigma() / 100.0;

  // La nueva producción viene dada por el valor aleatorio.
  sim_attr_t.imports = std::max<real>(0.0, RANDOM_NORMAL(mean, sigma));

  // Cantidad de producto disponible para venta.
  real available_quantity = sim_attr_t.imports + sim_attr_t_1.stock;

  distribute_internal_production(network, node, available_quantity, t);
}

/* Recorre los rangos topológicos de la red productiva para calcular las nuevas
   producciones en cada nodo.

   Para los nodos de tipo producto se cuentan los arcos de salida, si no tiene,
   entonces se llama a la función de cálculo para productos con el functor de
   cálculo de demanda externa para nodos de último nivel. En caso contrario se
   llama a la función con el functor que calcula la demanda externa para nodos
   que no son del último nivel.
*/
void compute_productions(Graph & network, Rank_Type & ranks, const size_t & t,
                         gsl_rng * rng)
{
  for (Rank_Type::Iterator it(ranks); it.has_current(); it.next())
    {
      Level * level = it.get_current();

      for (Level::Iterator lit(*level); lit.has_current(); lit.next())
        {
          Digraph::Node * dnode = lit.get_current();

          Graph::Node * node = dnode->get_info();

          switch (node->get_info()->get_type())
            {
            case PRODUCT:
              if (count_output_arcs(network, node) == 0)
                compute_product_production<Last_Level_Node>(network, node, t,
                                                            rng);              
              else
                compute_product_production<Any_Level_Node>(network, node, t,
                                                           rng);
              break;
            case INPUT:
              compute_input_production(network, node, t, rng);
              break;
            case IMPORTED_PRODUCT:
              compute_imported_product_production(network, node, t, rng);
              break;
            }
        }
    }
}

/* Actualiza los precios de venta de los productos en los arcos de salida. El
   nuevo precio de venta es el precio unitario del producto (costo total de
   producción) agregándole el porcentaje de ganancia.
*/
void update_price_in_output_arcs(Graph & network, Graph::Node * node,
                                 const size_t & t)
{
  const real & new_price = node->get_info()->price_at(t);

  for (Graph::Node_Arc_Iterator it(node); it.has_current(); it.next())
    {
      Graph::Arc * arc = it.get_current();

      Graph::Node * tgt = network.get_tgt_node(arc);

      // Ignoro arcos de entrada
      if (node == tgt)
        continue;

      IP_Relationship & ip_rel = arc->get_info();

      IP_Simulation_Attributes & sim_attr =
        ip_rel.get_simulation_attributes_at(t);

      // Al precio de venta se le incrementa la tasa de ganancia.
      sim_attr.purchase_price = new_price;
    }
}

// Cálcula el costo de insumos para un producto.
real compute_input_cost(Graph & network, Graph::Node * node, const size_t & t)
{
  real input_cost = 0.0;

  for (Graph::Node_Arc_Iterator it(node); it.has_current(); it.next())
    {
      Graph::Arc * arc = it.get_current();

      Graph::Node * src = network.get_src_node(arc);

      // Ignoro arcos de salida
      if (node == src)
        continue;

      IP_Relationship & ip_rel = arc->get_info();

      IP_Simulation_Attributes & sim_attr =
        ip_rel.get_simulation_attributes_at(t);

      input_cost += (sim_attr.purchase_price * sim_attr.bought_quantity);
    }

  return input_cost;
}

/* Calcula los nuevos costos para un nodo de tipo producto, los nuevos salarios
   y el costo unitario de producción.
*/
void compute_product_price(Graph & network, Graph::Node * node,
                           const size_t & t)
{
  Exogenous_Variables & exo_var = Exogenous_Variables::get_instance();

  Good * ptr_good = node->get_info().get();

  Product * ptr_product = static_cast<Product *>(ptr_good);

  Product_Simulation_Attributes & sim_attr_t =
    ptr_product->get_simulation_attributes_at(t);

  Product_Simulation_Attributes & sim_attr_t_1 =
    ptr_product->get_simulation_attributes_at(t - 1);

  // Actualización de salarios

  sim_attr_t.daily_wage = sim_attr_t_1.daily_wage *
                         (1.0 + exo_var.rate_of_change_of_salary_at(t) / 100.0);

  // Nuevo salario anual
  real wage = sim_attr_t.daily_wage * DAYS_IN_A_YEAR;

  sim_attr_t.other_wage = wage * WAGE_PROPORTION;
  sim_attr_t.other_daily_wage = sim_attr_t.other_wage / DAYS_IN_A_YEAR;

  sim_attr_t.daily_feeding_coupon = exo_var.UT_at(t) * UT_PROPORTION;

  sim_attr_t.feeding_coupon = sim_attr_t.daily_feeding_coupon * DAYS_IN_A_YEAR;

  sim_attr_t.integral_wage = wage + sim_attr_t.other_wage +
                             sim_attr_t.feeding_coupon;

  // Actualización de costos

  // Cálculo de costos administrativos
  sim_attr_t.administrative_staff_cost =
    ptr_product->get_num_administrative_staff() * sim_attr_t.integral_wage;
  sim_attr_t.unitarian_administrative_staff_cost =
    sim_attr_t.administrative_staff_cost / sim_attr_t.production;

  // Cálculo de costos en mano de obra
  sim_attr_t.labor_cost = sim_attr_t.num_employees * sim_attr_t.integral_wage;
  sim_attr_t.unitarian_labor_cost =
    sim_attr_t.labor_cost / sim_attr_t.production;

    // Actualización de costo de insumos.
  real input_cost = compute_input_cost(network, node, t);

  if (input_cost == 0.0) // No tenía insumos, cambia por tasa de precios
    sim_attr_t.input_cost = sim_attr_t_1.input_cost *
                            (1.0 + exo_var.rate_of_price_change_at(t) / 100.0);
  else
    sim_attr_t.input_cost = input_cost;

  sim_attr_t.input_cost_per_unit =
    sim_attr_t.input_cost / sim_attr_t.production;

  // Actualización de otros costos
  sim_attr_t.other_cost = sim_attr_t_1.other_cost *
                          (1.0 + exo_var.rate_of_price_change_at(t) / 100.0);
  sim_attr_t.other_unitarian_cost =
    sim_attr_t.other_cost / sim_attr_t.production;

  sim_attr_t.total_cost = sim_attr_t.administrative_staff_cost +
                          sim_attr_t.labor_cost + sim_attr_t.input_cost +
                          sim_attr_t.other_cost;

  sim_attr_t.price = (sim_attr_t.unitarian_administrative_staff_cost +
                     sim_attr_t.unitarian_labor_cost +
                     sim_attr_t.input_cost_per_unit +
                     sim_attr_t.other_unitarian_cost)  *
                     (1.0 + exo_var.rate_of_gain_at(t) / 100.0);

  sim_attr_t.internal_sales *= sim_attr_t.price;

  sim_attr_t.external_sales *= sim_attr_t.price;

  update_price_in_output_arcs(network, node, t);

  sim_attr_t.economic_status = sim_attr_t.income - sim_attr_t.total_cost;
}

// Calcula el costo unitario de producción para un nodo de tipo insumo.
void compute_input_price(Graph & network, Graph::Node * node,
                         const size_t & t)
{
  Exogenous_Variables & exo_var = Exogenous_Variables::get_instance();

  Good * ptr_good = node->get_info().get();

  Input * ptr_input = static_cast<Input *>(ptr_good);

  Input_Simulation_Attributes & sim_attr_t =
    ptr_input->get_simulation_attributes_at(t);

  Input_Simulation_Attributes & sim_attr_t_1 =
    ptr_input->get_simulation_attributes_at(t - 1);

  sim_attr_t.price = sim_attr_t_1.price *
                     (1.0 + exo_var.rate_of_price_change_at(t) / 100.0);

  sim_attr_t.internal_sales *= sim_attr_t.price;

  sim_attr_t.external_sales *= sim_attr_t.price;

  update_price_in_output_arcs(network, node, t);
}

/* Calcula el costo unitario de producción para un nodo de tipo producto
   importado.
*/
void compute_imported_product_price(Graph & network, Graph::Node * node,
                                    const size_t & t)
{
  Exogenous_Variables & exo_var = Exogenous_Variables::get_instance();

  Good * ptr_good = node->get_info().get();

  Imported_Product * ptr_imported_product =
    static_cast<Imported_Product *>(ptr_good);

  Imported_Product_Simulation_Attributes & sim_attr_t =
    ptr_imported_product->get_simulation_attributes_at(t);

  Imported_Product_Simulation_Attributes & sim_attr_t_1 =
    ptr_imported_product->get_simulation_attributes_at(t - 1);

  sim_attr_t.price = sim_attr_t_1.price *
                     (1.0 + exo_var.nominal_exchange_rate_at(t) / 100.0);

  update_price_in_output_arcs(network, node, t);
}

/* Recorre los rangos topológicos de la red productiva para calcular los nuevos
   precios en cada nodo.
*/
void compute_new_prices(Graph & network, Rank_Type & ranks,
                        const size_t & t)
{
  for (Rank_Type::Iterator it(ranks); it.has_current(); it.next())
    {
      Level * level = it.get_current();

      for (Level::Iterator lit(*level); lit.has_current(); lit.next())
        {
          Digraph::Node * dnode = lit.get_current();

          Graph::Node * node = dnode->get_info();

          switch (node->get_info()->get_type())
            {
            case PRODUCT:
              compute_product_price(network, node, t);
              break;
            case INPUT:
              compute_input_price(network, node, t);
              break;
            case IMPORTED_PRODUCT:
              compute_imported_product_price(network, node, t);
              break;
            }
        }
    }
}

/* Función que efectúa la simulación sobre la red productiva descrita en un
   archivo xml.
*/
int simulate(const char * const input_file_name,
             const char * const output_file_name)
{
  // Inicialización del generador de números aleatorios.
  gsl_rng * rng = gsl_rng_alloc(gsl_rng_mt19937);
  gsl_rng_set(rng, time(0));

  Graph network;

  Digraph aux_network;

  Rank_Type ranks;

  IO_Manager io_manager;

  // Lectura del grafo descrito en el archivo xml.
  io_manager.read_graph_from_xml(input_file_name, network);

  // Lectura de las variables exógenas del sistema descritas en el archivo xml.
  io_manager.load_exogenous_variables_from_xml(input_file_name);

  // Asignación de valores al estado inicial de la simulación.
  init_network(network);

  // Conversión de grafo a digrafo.
  graph_to_digraph(network, aux_network);

  // Cálculo de los rangos topológicos.
  Q_Topological_Sort<Digraph>()(aux_network, ranks);

  const size_t & max_t = Exogenous_Variables::get_instance().get_num_it();

  // Bucle de la simulación.
  for (size_t t = 1; t <= max_t; ++t)
    {
      compute_productions(network, ranks, t, rng);

      compute_new_prices(network, ranks, t);
    }

  // Libera la memoria rangos topológicos.
  free_ranks(ranks);

  // Libera la memoria del generador de números aleatorios.
  gsl_rng_free(rng);  

  // Escribe la red
  io_manager.add_sim_graph_to_xml(network, input_file_name, output_file_name);

  return 0;
}