m_popul.f90

Go to the documentation of this file.

!> @file m_popul.f90
!! The Population objects for the AHA Model.
!! @author Sergey Budaev <sergey.budaev@uib.no>
!! @author Jarl Giske <jarl.giske@uib.no>
!! @date 2016-2017
 
!-------------------------------------------------------------------------------
! $Id: m_popul.f90 11210 2021-06-03 09:47:06Z sbu062 $
!-------------------------------------------------------------------------------
 
!-------------------------------------------------------------------------------
!> @brief   Define the population of agents object, its properties and
!!          functions.
!> @section the_population_module THE_POPULATION module
!! @details The population is in the simplest case an array of individual
!!          agent objects. Individual properties of the population members can
!!          be referred as, e.g. `proto_parents%individual(i)%fitness` .
module the_population
 
  use commondata      ! Global definitions of the model objects
  use the_individual  ! The individual - properties of the individuals
 
  implicit none
 
  character (len=*), parameter, private :: modname = "(THE_POPULATION)"
 
  !> @brief   Definition of individual member of a population.
  !! @details Add an additional object, a member of the population.
  !!          It is the`INDIVIDUAL_AGENT` adding an (unique) integer
  !!          ID: `person_number`. Here we also have two functions for
  !!          setting and retrieving the IDs. We do not init members
  !!          of the population using a single init function, Instead,
  !!          init normally the INDIVIDUAL_AGENT object and for the
  !!          MEMBER_POPULATION just set_id. This is because we cannot
  !!          set id for isolated subject, only as a member of population.
  type, public, extends(individual_agent) :: member_population
    integer :: person_number
    contains
      private
      !> Set integer ID number to individual member of the population object.
      !! See `the_population::set_individual_id()`.
      procedure, private :: set_id => set_individual_id
      !> Get integer ID number to individual member of the population object.
      !! See `the_population::get_individual_id()`.
      procedure, public  :: get_id => get_individual_id
      !> Places the individual agent, a member of the population, within a specific
      !! environment at random with a **uniform** distribution.
      !! See `the_population::individ_posit_in_environ_uniform()`.
      procedure, public  :: place_uniform => individ_posit_in_environ_uniform
      !> Set the individual to be **dead**.
      !! See `the_population::genome_individual_set_dead_non_pure()`.
      procedure, public  :: dies_debug => genome_individual_set_dead_non_pure
  end type member_population
 
  !> @brief   Definition of the population object.
  !! @details This is basically an array of individuals (of the type
  !!          `MEMBER_POPULATION`) plus population or generation descriptors.
  type, public :: population
    !> The size of the population.
    integer :: population_size
    !> @brief   `POPULATION` is represented by an array of objects
    !!          of the type `MEMBER_POPULATION`
    !! @details `POPULATION` is an array of objects of the type
    !!          `MEMBER_POPULATION`. It is represented as the \%individual
    !!          component of the type. The `i`-th individual of the population
    !!          `this` is accessed as `this%individual(i)`. the population also
    !!          has two descriptors: integer `pop_number` and string `pop_name`.
    type(member_population), allocatable, dimension(:) :: individual
    !> The numeric ID of the population (if we have several populations))
    integer :: pop_number
    !> The descriptive name of the population
    character (len=LABEL_LENGTH) :: pop_name
    contains
      private
      !> Initialise the population object.
      !! See `the_population::init_population_random()`.
      procedure, public :: init => init_population_random
      !> Impose selective mortality at birth of the agents.
      !! See `the_population::population_birth_mortality_init()`.
      procedure, public :: mortality_birth => population_birth_mortality_init
      !> Destroys this population and deallocates the array of individual
      !! member  objects.
      !! See `the_population::population_destroy_deallocate_objects()`.
      procedure, public :: destroy => population_destroy_deallocate_objects
      !> Get the size of this population.
      !! See `the_population::population_get_popsize()`.
      procedure, public :: get_size => population_get_popsize
      !> Get the population number ID.
      !! See `the_population::population_get_pop_number()`.
      procedure, public :: get_num_id => population_get_pop_number
      !> Get the population character label ID.
      !! See `the_population::population_get_pop_name()`.
      procedure, public :: get_name => population_get_pop_name
      !> Reset individual IDs of the population members.
      !! See `the_population::reset_population_id_random()`.
      procedure, public :: reset_id => reset_population_id_random
      !> Determine the sex for each member of the population.
      !! See `the_population::sex_initialise_from_genome()`.
      procedure, public :: sex => sex_initialise_from_genome
      !> Position each member of the population randomly within a
      !! bounding environment with the **uniform** distribution.
      !! See `the_population::position_individuals_uniform()`.
      procedure, public :: scatter_uniform => position_individuals_uniform
      !> Perform a single step of random walk by all agents, in 3D.
      !! See `the_population::population_rwalk3d_all_agents_step()`.
      procedure, public :: rwalk3d => population_rwalk3d_all_agents_step
      !> Perform a single step of random walk by all agents, in 2.5D.
      !! See `the_population::population_rwalk25d_all_agents_step()`.
      procedure, public :: rwalk25d => population_rwalk25d_all_agents_step
      !> This subroutine sorts the population `individual` object by
      !! their \%fitness components.
      !! See `the_population::sort_population_by_fitness()`.
      procedure, public :: sort_by_fitness => sort_population_by_fitness
      !> Save data for all agents within the population into a csv file.
      !! See `the_population::population_save_data_all_agents_csv()`.
      procedure, public :: save_csv => population_save_data_all_agents_csv
      !> Save the genome data of all agents in this population to a CSV file.
      !! See `the_population::population_save_data_all_genomes()`.
      procedure, public :: save_genomes_csv => population_save_data_all_genomes
      !> Load agent genome data in this population from a CSV file.
      !! See `the_population::population_load_data_all_genomes()`.
      procedure, public :: load_genomes_csv => population_load_data_all_genomes
      !> Save the perceptual and emotional memory stack data of all agents
      !! in this population to a CSV file.
      !! See `the_population::population_save_data_memory()`.
      procedure, public :: save_memory_csv => population_save_data_memory
      !> Save the latest movement history of all agents.
      !! See `the_population::population_save_data_movements()`.
      procedure, public :: save_movements_csv => population_save_data_movements
      !> Save the behaviours history the_neurobio::behaviour::history_behave
      !! for all agents.
      !! See `the_population::population_save_data_behaviours()`.
      procedure, public :: save_behaviour_csv => population_save_data_behaviours
      !> Subject the population to an attack by a specific predator.
      !! See `the_population::population_subject_predator_attack()`.
      procedure, public :: attacked => population_subject_predator_attack
      !> Subject the population to mortality caused by habitat-specific
      !! mortality risk. Each agent is affected by the risk associated with
      !! the habitat it is currently in.
      !! See `the_population::population_subject_other_risks()`.
      procedure, public :: mortality_habitat => population_subject_other_risks
      !> Subject all members of this population to their individual mortality
      !! risks.
      !!  See `the_population::population_subject_individual_risk_mortality()`.
      procedure, public :: mortality_individ =>                               &
                                  population_subject_individual_risk_mortality
      !> Calculate fitness for the pre-evolution phase of the genetic algorithm.
      !! **Pre-evolution** is based on selection for a simple criterion without
      !! explicit reproduction etc. The criterion for selection at this phase
      !! is set by the integer the_individual::individual_agent::fitness
      !! component. This procedure provides a whole-population wrapper for the
      !! the_individual::fitness_calc() function.
      !! See `the_population::population_preevol_fitness_calc()`.
      procedure, public :: fitness_calc => population_preevol_fitness_calc
      !> Determine the number of parents that have fitness higher than the
      !! minimum acceptable value.
      !! See `the_population::population_ga_reproduce_max()`.
      procedure, public :: ga_reproduce_max => population_ga_reproduce_max
      !> This function implements adaptive mutation rate that increases
      !! as the population size reduces.
      !! See `the_population::population_ga_mutation_rate_adaptive()`.
      procedure, public :: ga_mutat_adaptive =>                               &
                                          population_ga_mutation_rate_adaptive
      !> Perform a single step of the life cycle of the population.
      !! See `the_population::population_lifecycle_step_preevol()`.
      procedure, public :: lifecycle_step => population_lifecycle_step_preevol
      !> Perform a single step of the life cycle of the population. This
      !! version includes only optimal food selection and eating without
      !! the full fledged behaviour selection cascade of procedures
      !! ::do_behave().
      !! See `the_population::population_lifecycle_step_eatonly_preevol()`.
      procedure, public :: lifecycle_eatonly =>                               &
                                      population_lifecycle_step_eatonly_preevol
 
  end type population
 
  !> Global indicator variable that keeps the number of agents that have died
  !! as a consequence of predatory attacks. All other dies are therefore caused
  !! by starvation.
  !! @note Note that this variable is initialised to zero at the start of each
  !!       generation in `GENERATIONS_PREEVOL` named loop in
  !!       the_evolution::generations_loop_ga().
  integer, public :: global_ind_n_eaten_by_predators
 
contains ! ........ implementation of procedures for this level ................
 
  !-----------------------------------------------------------------------------
  !> Set integer ID number to individual member of the population object.
  !! @note Note that this subroutine is private as setting individual IDs
  !!       makes sense only during the initialisation phase of the population.
  subroutine set_individual_id (this, idnumber)
    !> @param class, member of population.
    class(member_population), intent(inout) :: this
    !> @param id, integer id number assigned.
    integer, intent(in) :: idnumber
 
    this%person_number = idnumber    !> `person_number` is assigned the idnumber
 
  end subroutine set_individual_id
 
  !-----------------------------------------------------------------------------
  !> Get integer ID number to individual member of the population object.
  !! @note Note that this is public, so we can retrieve individual IDs but not
  !!       set them (id's are always set at initialisation)
  function get_individual_id (this) result (idnumber)
    !> @param class, member of population.
    class(member_population), intent(in) :: this
    !> @param id, integer id number retreived.
    integer :: idnumber
 
    idnumber = this%person_number    !> `person_number` is assigned the idnumber
 
  end function get_individual_id
 
  !-----------------------------------------------------------------------------
  !> Places the individual agent, a member of the population, within a specific
  !! environment at random with a **uniform** distribution.
  !! The agents can be positioned with respect to their initial depth
  !! - fixed depth, see commondata::init_agents_depth_is_fixed
  !! - Gaussian depth, see commondata::init_agents_depth_is_gauss
  !! - fully uniform distribution.
  !! .
  subroutine individ_posit_in_environ_uniform(this, environ)
    class(member_population), intent(inout) :: this
    !> @param environ sets the environment object where the agent is
    !!        placed randomly
    class(environment), intent(in) :: environ
 
    if (init_agents_depth_is_fixed) then
      call this%position( environ%uniform( within( init_agents_depth,         &
                                                   environ%depth_min(),       &
                                                   environ%depth_max()  ) ) )
    else if (init_agents_depth_is_gauss) then
      call this%position(                                                     &
          environ%uniform( within( rnorm( init_agents_depth,                  &
                                         cv2variance(init_agents_depth_cv,    &
                                             init_agents_depth)),             &
                                   environ%depth_min(),                       &
                                   environ%depth_max()                  ) ) )
    else
      call this%position( environ%uniform() )
    end if
 
  end subroutine individ_posit_in_environ_uniform
 
  !-----------------------------------------------------------------------------
  !> Set the individual to be **dead**. Note that this function does not
  !! deallocate the individual agent object, this may be a separate destructor
  !! function.
  !! @note    This is a non-pure function logging some extended diagnostics
  !!          about the dying agent. See `dies` procedure from @ref the_genome
  !!          and all its overrides:
  !!            - the_genome::individual_genome::dies();
  !!            - the_neurobio::appraisal::dies();
  !!            - the_neurobio::gos_global::dies();
  !!            - the_individual::individual_agent::dies().
  !!            .
  !! @warning This function should be normally used only while **debugging**.
  !!          In normal runs use the pure subroutine `dies` from
  !!          `ÌNDIVIDUAL_GENOME`.
  subroutine genome_individual_set_dead_non_pure(this, non_debug_log)
    class(member_population), intent(inout) :: this
    logical, optional, intent(in) :: non_debug_log
 
    !> Local flag to do show log.
    logical :: do_show_log
 
    !> Local parameter showing how many latest memory values to log out.
    integer, parameter :: HIST_N=10
 
    do_show_log = .false. !> Don't show log by default.
 
    call this%dies() !> This is a pure function from `THE_GENOME` level.
 
    !> Turn on output logging only if in the DEBUG mode or if `non_debug_log`
    !! is explicitly set to TRUE.
    if (is_debug) then
      do_show_log = .true.
    else
      if (present(non_debug_log)) then
        if (non_debug_log) do_show_log = .true.
      end if
    end if
 
    if (.not. do_show_log) return   !> Just exit back if no logging.
 
    call log_msg( "Agent # " //  tostr(this%get_id()) // " with name " //     &
                  trim(this%individ_label()) // " dies." )
    call log_msg( "Agent properties: " )
    call log_msg( "  body mass: " // tostr(this%get_mass()) //                &
                  ", body mass at birth: " // tostr(this%body_mass_birth) //  &
                  ", max body mass: " // tostr(this%body_mass_maximum) //     &
                  ", body length: " // tostr(this%get_length()) //            &
                  ", energy reserves: " // tostr(this%get_energy()) //        &
                  ", SMR: " // tostr(this%get_smr()) //                       &
                  ", stomach content: " // tostr(this%get_stom_content()) )
    call log_msg( "  Latest body mass (" // tostr(hist_n) // ") history: " // &
                  tostr(this%body_mass_history(                               &
                    history_size_agent_prop-hist_n+1:history_size_agent_prop)) )
    call log_msg( "  Latest body length (" // tostr(hist_n) // ") history: "  &
                  // tostr(this%body_length_history(                          &
                    history_size_agent_prop-hist_n+1:history_size_agent_prop)) )
    call log_msg( "Agent's GOS is " // this%gos_label() //                    &
                  ", arousal: " // tostr(this%arousal()) // "." )
    !> @note Note that `TOSTR` accepts arrays including (concatenated)
    !!       character arrays.
    call log_msg( "  Latest GOS states: " //                                  &
                  tostr(this%memory_motivations%gos_main(                     &
                    history_size_motivation-hist_n+1:history_size_motivation)) )
 
  end subroutine genome_individual_set_dead_non_pure
 
  !-----------------------------------------------------------------------------
  !> @brief   Initialise the population object.
  !! @details Initialise the population object, init it with random individuals
  !!          (function init on individuals), and assign sequential
  !!          person_number`s.
  subroutine init_population_random(this, pop_size, pop_number_here,          &
                                                                  pop_name_here)
    ! Parameters for this subroutine:
    !> @param class, This -- member of population class (this)).
    class(population), intent(inout) :: this
    !> @param pop_size, The size of the population.
    integer, intent(in) :: pop_size
    !> @param id, Optional numeric population ID.
    integer, optional, intent(in) :: pop_number_here
    !> @param descriptor, Optional population string descriptor.
    character (len=*), optional, intent(in) :: pop_name_here
 
    ! Local variables:
    integer :: i                !> local variable `i` is counter
 
    ! PROCNAME is the procedure name for logging and debugging
    character(len=*), parameter :: PROCNAME = "(init_population_random)"
 
    !> Set population size from input parameter.
    this%population_size = pop_size
 
    !>  Allocate the population
    if (.not. allocated(this%individual))                                     &
                                allocate(this%individual(this%population_size))
 
    !> Initialise all individuals of the population
    do i=1, this%population_size
      !> first, call `init to object individual(i)
      call this%individual(i)%init()
      !! second, ! call `set_id(i)` to identify agents within the population.
      call this%individual(i)%set_id(i)
    end do
 
    !> Set optional descriptors for the whole population. Numeric ID and
    !! a short text description.
    if (present(pop_number_here)) then
      this%pop_number = pop_number_here
    else
      !> If optional ID is absent, ID is set to a random integer value from 1
      !! to the maximum integer allowed for the pop_number type (minus 1).
      this%pop_number = rand_i(1,huge(this%pop_number-1))
    end if
    if (present(pop_name_here)) then
      this%pop_name = pop_name_here
    else
      !> If optional text description string of the population is absent,
      !! it is set to the string representation of its numeric ID.
      this%pop_name = tostr(this%pop_number)
    end if
 
    call log_msg( ltag_info // "Initialised population " // this%pop_name //  &
                  " # " //  tostr(this%pop_number) // " with " //             &
                  tostr(this%population_size) // " agents." )
 
    !> Log the initial location of the agents.
    !! @warning The logic of the logger constructs here must coincide with
    !!          that in the population::individ_posit_in_environ_uniform().
    if (init_agents_depth_is_fixed) then
      call log_msg( ltag_info // "Initial location is FIXED at " //           &
                    tostr(init_agents_depth) )
    else if (init_agents_depth_is_gauss) then
      call log_msg( ltag_info // "Initial location is GAUSSIAN at " //        &
                    tostr(init_agents_depth) //  ", with CV " //              &
                    tostr(init_agents_depth_cv) )
    else
      call log_msg( ltag_info // "Initial location is fully uniform." )
    end if
 
  end subroutine init_population_random
 
  !-----------------------------------------------------------------------------
  !> Destroys this population and deallocates the array of individual member
  !!  objects.
  subroutine population_destroy_deallocate_objects(this)
    class(population), intent(inout) :: this
 
    this%population_size = 0
    this%pop_number = unknown
    this%pop_name = ""
    if (allocated(this%individual)) deallocate(this%individual)
 
  end subroutine population_destroy_deallocate_objects
 
  !-----------------------------------------------------------------------------
  !> Impose selective mortality at birth on the agents. Selective mortality
  !! sets a fixed limit on uncontrolled evolution of the *energy reserves*
  !! in newborn agents. If some newborn has too high energy at birth
  !! (genetically fixed), such a deviating agent is killed at once.
  !!
  !! The values of the risk are chosen such that it is zero at the
  !! fixed population mean (agent does not deviate) and reaches 1.0 if
  !! the agent's energy reserves are 3.0 standard deviations higher than
  !! the fixed mean value (agent strongly deviates).
  subroutine population_birth_mortality_init(this, energy_mean, energy_sd)
    class(population), intent(inout) :: this
    !> @param[in] energy_mean optional mean energy at birth, if absent, is
    !!            calculated from the population data.
    real(SRP), optional, intent(in) :: energy_mean
    !> @param[in] energy_sd optional mean energy at birth, if absent, is
    !!            calculated from the population data.
    real(SRP), optional, intent(in) :: energy_sd
 
    ! Local copies of optionals
    real(SRP) :: energy_mean_loc, energy_sd_loc
 
    real(SRP) :: mortality
    integer :: ind
 
    !> - `MORTALITY_BIRTH_INIT_ENERG_ABSCISSA` is the baseline abscissa for the
    !!    nonparametric interpolation function grid, in units of the standard
    !!    deviation (sigma) of the energy reserves in the newborn population,
    !!    i.e. the_body::condition::energy_birth.
    ! htintrpl.exe [1.4 4 5] [0 0.5 1]
    ! htintrpl.exe [0.7 1 2 3] [0 0.005 0.1 1]
    real(SRP), parameter, dimension(*) :: MORTALITY_BIRTH_INIT_ENERG_ABSCISSA &
              = [ 0.7_srp,   1.0_srp,  1.5_srp, 2.0_srp, 3.0_srp ]
 
    !> - `MORTALITY_BIRTH_INIT_ENERG_ORDINATE` is the ordinate for the
    !!    nonparametric interpolation function grid: sets the probability of
    !!    the agent's death given its energy reserves deviate from the fixed
    !!    mean by the number of standard deviations set by
    !!    `MORTALITY_BIRTH_INIT_ENERG_ABSCISSA`.
    !! .
    real(SRP), parameter, dimension(*) :: MORTALITY_BIRTH_INIT_ENERG_ORDINATE &
              = [ 0.0_srp, 0.002_srp, 0.01_srp, 0.1_srp, 1.0_srp ]
 
    if (present(energy_mean)) then
      energy_mean_loc = energy_mean
    else
      energy_mean_loc = average( this%individual%get_energ_birth() )
    end if
 
    if (present(energy_sd)) then
      energy_sd_loc = energy_sd
    else
      energy_sd_loc = std_dev( this%individual%get_energ_birth() )
    end if
 
    inds: do ind=1, this%population_size
 
      !> Selective mortality sets a fixed limit on uncontrolled evolution of
      !! the energy reserves in newborn agents. All newborn agents that have
      !! the energy reserves exceeding some point may die with the probability
      !! determined by the grid arrays:
      !! - `MORTALITY_BIRTH_INIT_ENERG_ABSCISSA`;
      !! - `MORTALITY_BIRTH_INIT_ENERG_ORDINATE`.
      !! .
      !!
      !! The probability of death (mortality risk) is determined by the
      !! nonparametric nonlinear function defined by `DDPINTERPOL`, where
      !! the actual grid abscissa is the energy reserve of the agent in
      !! and grid ordinate is the mortality risk.
      !!
      !! The values of the risk are chosen such that it is zero at the
      !! fixed population mean (agent does not deviate) and reaches 1.0 if
      !! the agent's energy reserves are 3.0 standard deviations higher than
      !! the mean value.
      mortality = within(ddpinterpol(energy_mean_loc +                        &
                                       mortality_birth_init_energ_abscissa *  &
                                       energy_sd_loc,                         &
                                     mortality_birth_init_energ_ordinate,     &
                                     this%individual(ind)%get_energ_birth()), &
                         0.0_srp, 1.0_srp )
 
      !> Interpolation plots can be saved in the @ref intro_debug_mode
      !! "debug mode" using this plotting command:
      !! `commondata::debug_interpolate_plot_save()`.
      !! @warning Involves **huge** number of plots, should normally be
      !!          disabled.
      call debug_interpolate_plot_save(                                       &
            grid_xx=energy_mean_loc + mortality_birth_init_energ_abscissa *   &
                          energy_sd_loc,                                      &
            grid_yy=mortality_birth_init_energ_ordinate,                      &
            ipol_value=this%individual(ind)%get_energ_birth(),                &
            algstr="DDPINTERPOL",                                             &
            output_file="plot_debug_mortality_birth_" //                      &
                        "gen_" // tostr(global_generation_number_current) //  &
                        "_step_"// tostr(global_time_step_model_current) //   &
                        mmdd // "_a_"//                                       &
                        trim(this%individual(ind)%individ_label()) //         &
                        "_" // rand_string(label_length, label_cst,label_cen) &
                        // ps )
 
      if ( rand_r4() < mortality ) then
        call this%individual(ind)%dies()
      end if
    end do inds
 
  end subroutine population_birth_mortality_init
 
  !-----------------------------------------------------------------------------
  !> Accessor get-function for the size of this population.
  function population_get_popsize(this) result (pop_size_output)
    class(population), intent(in) :: this
    integer :: pop_size_output
 
    pop_size_output = this%population_size
 
  end function population_get_popsize
 
  !-----------------------------------------------------------------------------
  !> Accessor get-function for the population number ID.
  function population_get_pop_number(this) result(pop_number_output)
    class(population), intent(in) :: this
    integer :: pop_number_output
 
    pop_number_output = this%pop_number
 
  end function population_get_pop_number
 
  !-----------------------------------------------------------------------------
  !> Accessor get-function for the population character label ID.
  function population_get_pop_name(this) result (pop_name_string_out)
    class(population), intent(in) :: this
    character(len=LABEL_LENGTH) :: pop_name_string_out
 
    pop_name_string_out = this%pop_name
 
  end function population_get_pop_name
 
  !-----------------------------------------------------------------------------
  !> @brief   Reset individual IDs of the population members.
  !! @details Makes new random individual IDs for the population members.
  subroutine reset_population_id_random(this, pop_number_here, pop_name_here)
  ! Parameters for this subroutine:
    !> @param class, This -- member of population class (this)).
    class(population), intent(inout) :: this
    !> @param id, Optional numeric population ID.
    integer, optional, intent(in) :: pop_number_here
    !> @param descriptor, Optional population string descriptor.
    character (len=*), optional, intent(in) :: pop_name_here
 
    ! Local variables:
    integer :: i                !> local variable `i` is counter
 
    !> Reset all individual IDs of the population members
    do i=1, this%population_size
      call this%individual(i)%set_id(i) ! call set_id(i) to the same object
    end do
 
    !> Set optional descriptors for the whole population. Numeric ID and
    !! a short text description.
    if (present(pop_number_here)) then
      this%pop_number = pop_number_here
    else
      !> If optional ID is absent, ID is set to a random integer value from 1
      !! to the maximum integer allowed for the pop_number type (minus 1).
      this%pop_number = rand_i(1,huge(this%pop_number-1))
    end if
    if (present(pop_name_here)) then
      this%pop_name = pop_name_here
    else
      !> If optional text description string of the population is absent,
      !! it is set to the string representation of its numeric ID.
      this%pop_name = tostr(this%pop_number)
    end if
 
  end subroutine reset_population_id_random
 
  !-----------------------------------------------------------------------------
  !> @brief Determine the sex for each member of the population.
  subroutine sex_initialise_from_genome(this)
    class(population), intent(inout) :: this
 
    integer :: i
 
    do i=1, this%population_size
      call this%individual(i)%sex_init()
    end do
 
  end subroutine sex_initialise_from_genome
 
  !-----------------------------------------------------------------------------
  !> @brief Position each member of the population randomly within a
  !!        bounding environment.
  !! @note  Moved the positioning procedure into a separate procedure as
  !!        initialising population may involve different spatial positioning,
  !!        e.g. uniform within the bounding environment or gaussian localised.
  subroutine position_individuals_uniform(this, environ)
    class(population), intent(inout) :: this
 
    !> @param   environ the environment where we place the population
    !! @warning Even though this parameter is optional, the bounding
    !!          environment should in most cases (almost?) be fixed, and
    !!          provided as a parameter. In most cases, unlimited environment
    !!          is useful for debugging only.
    !! TODO: convert to class
    class(environment), optional, intent(in) :: environ
 
    !> Local object representing the bounding environment for this population
    type(environment) :: environ_here
 
    !> Local counter
    integer :: i
 
    !> Check if the bounding environment is provided, if not, place agents
    !! without limits
    !! @warning The bounding environment should in most cases be fixed, and
    !!          provided as a parameter.
    if (present(environ)) then
      call environ_here%build( environ%lim_min(), environ%lim_max() )
    else
      call environ_here%build_unlimited()
    end if
 
    !> Position agents randomly (uniform distribution) within the
    !! bounding environment.
    do i=1, this%population_size
      call this%individual(i)%place_uniform(environ_here)
    end do
 
  end subroutine position_individuals_uniform
 
  !-----------------------------------------------------------------------------
  !> @brief   This subroutine sorts the population `individual` object by
  !!          their \%fitness components.
  !! @details The two subroutines `qsort` and `qs_partition_fitness` are a
  !! variant of the recursive quick sort algorithm adapted for
  !! `MEMBER_POPULATION` integer fitness component
  elemental subroutine sort_population_by_fitness(this)
    !> @param class, This -- member of population class (this)).
    class(population), intent(inout) :: this
 
    call qsort(this%individual) !> This is the array component we sort.
 
    contains
 
    !...........................................................................
    !> `qsort` is a recursive frontend for `MEMBER_POPULATION` objects
    ! Based on:
    ! - Juli Rew, SCD Consulting (juliana@ucar.edu), 9/03
    ! - Based on algorithm from Cormen et al., Introduction to Algorithms, 1997
    ! - Made F conformant by Walt Brainerd
    ! - Source: http://www.fortran.com/qsort_c.f95
    recursive pure subroutine qsort(A, is_reverse)
 
      !> @param `A` has the same type as the individual component objects
      !!         of the array-object that we are going to sort.
      type(member_population), intent(in out), dimension(:) :: a
      !> Logical flag for reverse (descending) sorting.
      logical, optional, intent(in) :: is_reverse
      integer :: iq
 
      if(size(a) > 1) then
        call qs_partition_fitness(a, iq) ! partition
        call qsort(a(:iq-1))
        call qsort(a(iq:))
      endif
 
      ! Check if reverse sorted array is requested
      if (present(is_reverse)) then
        if (is_reverse) a = a( size(a):1:-1 )
      end if
 
    end subroutine qsort
 
    !...........................................................................
    !> partition is a pivot backend for `fitness`
    pure subroutine qs_partition_fitness(A, marker)
 
      type(member_population), intent(in out), dimension(:) :: a
      integer, intent(out) :: marker
      integer :: i, j
      type(member_population) :: temp
      !> @note Pivot point `x`, has the same type **as
      !!       the sorted object component**.
      integer :: x
 
      !> Fitness is hardwired in this partition subroutine, but it can be used
      !! as a model for similar othert sorting functions. Note that here integer
      !! array is sorted.
      x = a(1)%fitness
      i= 0
      j= size(a) + 1
 
      do
        j = j-1
        do
            if (a(j)%fitness <= x) exit
            j = j-1
        end do
        i = i+1
        do
            if (a(i)%fitness >= x) exit
            i = i+1
        end do
        if (i < j) then
            ! exchange A(i) and A(j)
            temp = a(i)
            a(i) = a(j)
            a(j) = temp
        elseif (i == j) then
            marker = i+1
            return
        else
            marker = i
            return
        endif
      end do
 
    end subroutine qs_partition_fitness
 
  end subroutine sort_population_by_fitness
 
  !-----------------------------------------------------------------------------
  !> Perform one or several steps of random walk by all agents.
  !! @note This procedure was used for debugging.
  subroutine population_rwalk3d_all_agents_step( this, dist_array, cv_array,  &
                                                 dist_all, cv_all,            &
                                                 environment_limits, n_walks )
    class(population), intent(inout) :: this
 
    !> @param[in] step_size_array an array of step sizes for each individual.
    real(SRP), optional, dimension(:), intent(in) :: dist_array
    !> @param[in]cv_array Coefficients of variation for the walk.
    real(SRP), optional, dimension(:), intent(in) :: cv_array
    !> @param[in] dist_all the value of the walk step size that is identical in
    !!            all agents within the population.
    real(SRP), optional, intent(in) :: dist_all
    !> @param[in] cv_all the value of the walk coefficient of variation that is
    !!            identical in all agents within the population.
    real(SRP), optional, intent(in) :: cv_all
    !> @param environment_limits Limits of the environment area available for
    !!        the random walk. The moving object cannot get beyond this limit.
    !!        If this parameter is not provided, the environmental limits are
    !!        obtained automatically from the global array
    !!        the_environment::global_habitats_available.
    class(environment), intent(in), optional :: environment_limits
 
    !> @param[in] n_walk optional number of walk steps that should be
    !!            performed, default just one.
    integer, optional, intent(in) :: n_walks
 
    ! Local variables, copies of optionals.
    real(SRP), dimension(this%population_size) :: dist_array_here, cv_array_here
    integer :: n_walks_here
 
    ! Local params.
    real(SRP), dimension(this%population_size) :: step_size_walk
    integer :: j, i, ind, pop_n
    integer, dimension(this%population_size) :: pop_permutation
 
    ! Default walk step CV.
    real(SRP), parameter ::  CV_DEFAULT = 0.5_srp
 
    !> ### Implementation details ###
    !> - Calculate the distance array size.
    pop_n = this%population_size
 
    if (present(dist_array)) then
      dist_array_here = dist_array
    else
      dist_array_here = this%individual%get_length()
    end if
 
    if (present(cv_array)) then
      cv_array_here = cv_array
    else
      cv_array_here = cv_default
    end if
 
    if (present(dist_all)) then
      dist_array_here = dist_all
    else
      dist_array_here = this%individual%get_length()
    end if
 
    if (present(cv_all)) then
      cv_array_here = cv_all
    else
      cv_array_here = cv_default
    end if
 
    if (present(n_walks)) then
      n_walks_here = n_walks
    else
      n_walks_here = 1
    end if
 
    !> - calculate the step size along the axes from the distance array.
    step_size_walk = dist2step(dist_array_here)
 
    !> - Calculate the random permutation of individual indices.
    !!   @warning Random order here is a prototype for testing for use in
    !!          behaviour selection by population members.
    pop_permutation = permute_random(pop_n)
 
    !> - Perform Gaussian random walks for each of the individuals in a random
    !!   order that is set by the `pop_permutation` array.
    !! .
    environ_restrict: if (present(environment_limits)) then
      do j=1, n_walks_here
        do i=1, pop_n
          ind = pop_permutation(i)
          if (this%individual(ind)%is_alive())                                &
                    call this%individual(ind)%rwalk( step_size_walk(ind),     &
                                                     cv_array_here(ind),      &
                                                     environment_limits  )
        end do
      end do
    else environ_restrict
      do j=1, n_walks_here
        do i=1, pop_n
          ind = pop_permutation(i)
          if (this%individual(ind)%is_alive())                                &
                    call this%individual(ind)%rwalk(                          &
                                step_size_walk(ind),                          &
                                cv_array_here(ind),                           &
                                global_habitats_available(                    &
                                    this%individual(ind)%find_environment(    &
                                              global_habitats_available) ) )
        end do
      end do
    end if environ_restrict
 
  end subroutine population_rwalk3d_all_agents_step
 
  !-----------------------------------------------------------------------------
  !> Perform one or several steps of random walk by all agents.
  !! @note This procedure was used for debugging.
  subroutine population_rwalk25d_all_agents_step ( this,                      &
                                            dist_array_xy, cv_array_xy,       &
                                            dist_array_depth, cv_array_depth, &
                                            dist_all_xy, cv_all_xy,           &
                                            dist_all_depth, cv_all_depth,     &
                                            environment_limits, n_walks )
    class(population), intent(inout) :: this
 
    !> @param[in] dist_array_xy an array of step sizes for each individual.
    real(SRP), optional, dimension(:), intent(in) :: dist_array_xy
    !> @param[in]cv_array_xy Coefficients of variation for the walk.
    real(SRP), optional, dimension(:), intent(in) :: cv_array_xy
 
    !> @param[in] dist_array_depth an array of step sizes for each individual.
    real(SRP), optional, dimension(:), intent(in) :: dist_array_depth
    !> @param[in]cv_array_depth Coefficients of variation for the walk.
    real(SRP), optional, dimension(:), intent(in) :: cv_array_depth
 
    !> @param[in] dist_all_xy the value of the walk step size for horizontal
    !!            plane that is identical in all agents within the population.
    real(SRP), optional, intent(in) :: dist_all_xy
    !> @param[in] cv_all_xy the value of the walk coefficient of variation in
    !!            the horizontal plane that is identical in all agents within
    !!            the population.
    real(SRP), optional, intent(in) :: cv_all_xy
 
    !> @param[in] dist_all_depth the value of the walk step size for the depth
    !!            plane that is identical in all agents within the population.
    real(SRP), optional, intent(in) :: dist_all_depth
    !> @param[in] cv_all_depth the value of the walk coefficient of variation
    !!            in the depth plane that is identical in all agents within
    !!            the population.
    real(SRP), optional, intent(in) :: cv_all_depth
    !> @param environment_limits Limits of the environment area available for
    !!        the random walk. The moving object cannot get beyond this limit.
    !!        If this parameter is not provided, the environmental limits are
    !!        obtained automatically from the global array
    !!        the_environment::global_habitats_available.
    class(environment), intent(in), optional :: environment_limits
    !> @param[in] n_walk optional number of walk steps that should be
    !!            performed, default just one.
    integer, optional, intent(in) :: n_walks
 
    ! Local variables, copies of optionals.
    real(SRP), dimension(this%population_size) ::                             &
                                      dist_array_xy_here, cv_array_xy_here
    real(SRP), dimension(this%population_size) ::                             &
                                      dist_array_depth_here, cv_array_depth_here
    integer :: n_walks_here
 
    ! Local params.
    real(SRP), dimension(this%population_size) :: step_size_walk_xy,          &
                                                  step_size_walk_depth
    integer :: j, i, ind, pop_n
    integer, dimension(this%population_size) :: pop_permutation
 
    ! Default walk step CV.
    real(SRP), parameter ::  CV_DEFAULT = 0.5_srp
 
    !> ### Implementation details ###
    !> - Calculate the distance array size.
    pop_n = this%population_size
 
    if (present(dist_array_xy)) then
      dist_array_xy_here = dist_array_xy
    else
      dist_array_xy_here = this%individual%get_length()
    end if
 
    if (present(cv_array_xy)) then
      cv_array_xy_here = cv_array_xy
    else
      cv_array_xy_here = cv_default
    end if
 
    !> - If the depth walk step distance is not provided as a parameter,
    !!   1/2 of the agent body size is used as the default value. Thus,
    !!   it is assumed that the extent of random movements of the agents
    !!   in the horizontal plane is greater than vertical movements.
    if (present(dist_array_depth)) then
      dist_array_depth_here = dist_array_depth
    else
      dist_array_depth_here = this%individual%get_length() / 2.0_srp
    end if
 
    if (present(cv_array_depth)) then
      cv_array_depth_here = cv_array_depth
    else
      cv_array_depth_here = cv_default
    end if
 
    if (present(dist_all_xy)) then
      dist_array_xy_here = dist_all_xy
    else
      dist_array_xy_here = this%individual%get_length()
    end if
 
    if (present(cv_all_xy)) then
      cv_array_xy_here = cv_all_xy
    else
      cv_array_xy_here = cv_default
    end if
 
    if (present(dist_all_depth)) then
      dist_array_depth_here = dist_all_depth
    else
      dist_array_depth_here = this%individual%get_length() / 2.0_srp
    end if
 
    if (present(cv_all_depth)) then
      cv_array_depth_here = cv_all_depth
    else
      cv_array_depth_here = cv_default
    end if
 
    if (present(n_walks)) then
      n_walks_here = n_walks
    else
      n_walks_here = 1
    end if
 
    !> - Calculate the step size along the axes from the distance array.
    step_size_walk_xy = dist2step(dist_array_xy_here)
    step_size_walk_depth = dist2step(dist_array_depth_here)
 
    !> - Calculate the random permutation of individual indices.
    !!   @warning Random order here is a prototype for testing for use in
    !!          behaviour selection by population members.
    pop_permutation = permute_random(pop_n)
 
    !> - Perform Gaussian random walks for each of the individuals in a random
    !!   order that is set by the `pop_permutation` array.
    !! .
    environ_restrict: if (present(environment_limits)) then
      do j=1, n_walks_here
        do i=1, pop_n
          ind = pop_permutation(i)
          if (this%individual(ind)%is_alive())                                &
              call this%individual(ind)%rwalk25d                              &
                    ( meanshift_xy = step_size_walk_xy(ind),                  &
                      cv_shift_xy = cv_array_xy_here(ind),                    &
                      meanshift_depth = step_size_walk_depth(ind),            &
                      cv_shift_depth = cv_array_depth_here(ind),              &
                      environment_limits = environment_limits  )
        end do
      end do
    else environ_restrict
      do j=1, n_walks_here
        do i=1, pop_n
          ind = pop_permutation(i)
          if (this%individual(ind)%is_alive())                                &
              call this%individual(ind)%rwalk25d                              &
                    ( meanshift_xy = step_size_walk_xy(ind),                  &
                      cv_shift_xy = cv_array_xy_here(ind),                    &
                      meanshift_depth = step_size_walk_depth(ind),            &
                      cv_shift_depth = cv_array_depth_here(ind),              &
                      environment_limits=global_habitats_available(           &
                                      this%individual(ind)%find_environment(  &
                                              global_habitats_available) ) )
        end do
      end do
    end if environ_restrict
 
  end subroutine population_rwalk25d_all_agents_step
 
  !-----------------------------------------------------------------------------
  !> Subject the population to an attack by a specific predator. The predator
  !! acts on agents in its proximity and takes account of the predation
  !! confusion and dilution effects (see
  !! the_environment::predator::risk_fish_group()).
  subroutine population_subject_predator_attack(this, this_predator,          &
                                                                time_step_model)
    class(population), intent(inout) :: this
    class(predator), intent(in) :: this_predator
    integer, optional, intent(in) :: time_step_model
 
    ! Local copie sof optionals
    integer :: time_step_model_here
 
    !> ### Notable variables ###
    !> - p_risk is an array with the size equal to the agent population size,
    !!   that keeps all the attack probabilities calculated by the
    !!   the_environment::predator::risk_fish_group() function.
    real(SRP), dimension(size(this%individual)) :: p_risk
 
    !> - prey_index is the partial index of the prey agents that are in proximity of
    !!   this predator.
    !! .
    integer, dimension(size(this%individual)) :: prey_index
 
    ! This is a temporary array of the the_environment::spatial type to keep
    ! the location of the prey agent data.
    ! COMPILER BUG REPORT:
    ! @warning It is necessary as a work around the bug in GNU gfortran
    !          where the array is passed into the (neighbours) subroutine
    !          and suddenly loses the standard Fortran array bounds starting
    !          from 1 and gets C-specific bounds starting from 0, keeping
    !          the overall array size intact.
    type(spatial), dimension(size(this%individual)) :: tmp_location
 
    ! Local counter
    integer :: i
 
    ! PROCNAME is the procedure name for logging and debugging
    character(len=*), parameter :: PROCNAME =                                 &
                                        "(population_subject_predator_attack)"
 
    !> ### Implementation notes ###
    !> **Preparations:** Check optional time step parameter. If not provided,
    !! use global commondata::global_time_step_model_current parameter value.
    if (present(time_step_model)) then
      time_step_model_here = time_step_model
    else
      time_step_model_here = global_time_step_model_current
    end if
 
    !> **First,** calculate the risk of predation from this specific predator
    !! to each of the agents in the population using the
    !! the_environment::predator::risk_fish_group() method. The "raw" indexed
    !! output scheme is used here to avoid multiple cycling over the whole
    !! large population of agents, only the agents that are closest to the
    !! predator are processed and attacked (maximum number is equal to the
    !! index limit commondata::predator_risk_group_select_index_partial).
    tmp_location = this%individual%location()
    call this_predator%risk_fish_group(                                       &
                          prey_spatial = tmp_location,                        &
                          prey_length = this%individual%get_length(),         &
                          is_freezing = this%individual%freeze%is_executed(), &
                          time_step_model = time_step_model_here,             &
                          risk_indexed = p_risk,                              &
                          index_dist = prey_index )
 
    !> **Second,** cycle through all the agents in close proximity of the
    !! predator, up to the maximum size of partial indexing parameter
    !! commondata::predator_risk_group_select_index_partial. The predator then
    !! stochastically attacks each of these agents with the probability equal
    !! to the risk of predation. If the attack is successful, the agent
    !! the_genome::individual_genome::dies() and loop is exited because the
    !! predator is assumed to catch only one agent at a time.
    do i=1, predator_risk_group_select_index_partial
      if ( this%individual(prey_index(i))%is_alive() ) then
        if ( rand_r4() < p_risk(i) ) then
          ! Increment global counter of agents that are eaten by predators
          global_ind_n_eaten_by_predators = global_ind_n_eaten_by_predators + 1
          call this%individual(prey_index(i))%dies()
          call log_dbg( ltag_info // "The agent # " // tostr(prey_index(i)) //&
                        " (" // tostr(i) // ")" //                            &
                        " is caught by the predator and dies.",               &
                        procname, modname)
          call log_dbg( ltag_info // "Agent status is: " //                   &
                        tostr(this%individual(prey_index(i))%is_alive()) )
          exit ! the predator catches only one agent at a time
        end if
      end if
    end do
 
  end subroutine population_subject_predator_attack
 
  !-----------------------------------------------------------------------------
  !> Subject the population to mortality caused by habitat-specific
  !! mortality risk. Each agent is affected by the risk associated with
  !! the habitat it is currently in.
  !> @note Note that there is no such a function for a single agent as it does
  !!       not seem to be necessary.
  subroutine population_subject_other_risks(this)
    class(population), intent(inout) :: this
 
    ! Local variable that sets the habitat (number) the agent is currently in
    ! within the global array the_environment::global_habitats_available.
    integer :: agent_in
 
    ! Local counter.
    integer :: i
 
    ! PROCNAME is the procedure name for logging and debugging
    character(len=*), parameter :: PROCNAME = "(population_subject_other_risks)"
 
    !> ### Implementation notes ###
    !> All agents in the population are randomly subjected to the mortality
    !! risk  the_environment::habitat::risk_mortality that is linked to the
    !! habitat object the agent is currently in in a loop.
    do i=1, this%population_size
      agent_in = this%individual(i)%find_environment()
      if (rand_r4() < global_habitats_available(agent_in)%get_mortality()) then
        !> If the agent is unhappy and is subjected to mortality event, it
        !! immediately the_genome::individual_genome::dies().
        call this%individual(i)%dies()
        call log_dbg( ltag_info // "Agent " //                                &
                  tostr(this%individual(i)%get_id()) // " dies due to " //    &
                  "habitat linked mortality risk " //                         &
                  tostr(global_habitats_available(agent_in)%get_mortality()), &
                  procname, modname  )
      end if
    end do
 
  end subroutine population_subject_other_risks
 
  !-----------------------------------------------------------------------------
  !> Subject all members of this population to their individual mortality
  !! risks.
  subroutine population_subject_individual_risk_mortality(this)
    class(population), intent(inout) :: this
 
    ! Local counter.
    integer :: i
 
    !> ### Implementation notes ###
    !> The procedure is simple, loop over all individual agents in the
    !! population and stochastically call the_genome::individual_genome::dies()
    !! method with probability equal to individual mortality of the agent.
    do i=1, this%population_size
      if ( rand_r4() < this%individual(i)%get_mortality() )                   &
                                                call this%individual(i)%dies()
    end do
 
  end subroutine population_subject_individual_risk_mortality
 
  !-----------------------------------------------------------------------------
  !> This procedure performs a **single step** of the life cycle of the whole
  !! population, the agents for the step are selected in a random order.
  subroutine population_lifecycle_step_preevol(this)
    class(population), intent(inout) :: this
 
    !> ### Notable variables ###
    !> - `inds_order` is an array that sets the order in which the agents are
    !!   being drawn out of the population to perform the step.
    integer, dimension(this%population_size) :: inds_order
 
    !> - `ind_seq` is the sequential number of the agent as drawn from the
    !!   population; it is the loop control counter variable.
    integer :: ind_seq
 
    !> - `ind_real` is the real sequential number of the agent in the
    !!   population; this real id is obtained from the order array
    !!   `inds_order`.
    !! .
    integer :: ind_real
 
    ! Local variable that sets the habitat (number) the agent is currently in
    ! within the global array the_environment::global_habitats_available.
    integer :: agent_in
 
    ! PROCNAME is the procedure name for logging and debugging
    character(len=*), parameter :: PROCNAME =                                 &
                                          "(population_lifecycle_step_preevol)"
 
    !> ### Implementation notes ###
    !> First, an ordering array `inds_order` is calculated that sets the order
    !! in which the agents are drawn from the population. In the simplest case,
    !! the order of the agents is random, so this array is actually an array of
    !! random integers. The procedure
    !! [PERMUTE_RANDOM](http://ahamodel.uib.no/doc/ar01s09.html#_random_permutation_permute_random_function)
    !! from HEDTOOLS is used here then.
    inds_order = permute_random(this%population_size)
 
    !> The agents can also be processed in any **non-random** order. This would
    !! require invoking an array indexing procedure
    !! [ARRAY_INDEX](http://ahamodel.uib.no/doc/ar01s07.html#_subroutines_array_index_and_array_rank)
    !! instead of `PERMUTE_RANDOM`.
    !!
    !! For example, to process the agents in the order of the body mass,
    !! the ordering array can be obtained obtained as below:
    !! @code
    !! call ARRAY_INDEX(this%individual%get_mass(), inds_order)
    !! @endcode
    !> But this code would rank the agents in an *increasing* order of their
    !! body mass. If this is not what is expected, e.g. if the
    !! non-random selection is used to mimic a competitive advantage for bigger
    !! and heavier agents, a reverse of the ordering is obtained like this:
    !! @code
    !! inds_order = inds_order(this%population_size:1:-1)
    !! @endcode
 
    !> Second, the agents are drawn from the population, one by one, in the
    !! named loop construct `INDIVIDUALS`.
    individuals: do ind_seq = 1, popsize
 
      !> - One individual is then drawn from the `inds_order` ordering array.
      ind_real = inds_order(ind_seq)
 
      associate( agent => this%individual(ind_real) )
 
        !> #### Initialisations ####
        !> - First, a check is done if the agent is dead or starved to death;
        !!   if yes, no further processing is done. Also, in the former case
        !!   the individual_genome::dies() method is called.
        if (agent%is_dead()) then
          cycle individuals
        end if
        if (agent%starved_death()) then
          call agent%dies()
          cycle individuals
        end if
 
        !> - `agent_in` index calculates the habitat number where the selected
        !!   agent is currently in, calling the
        !!   the_environment::spatial::find_environment() method.
        agent_in = agent%find_environment()
        if (agent_in<1 .or. agent_in>size(global_habitats_available) ) then
          call log_msg( ltag_error // "agent_in zero in " // procname //      &
                        tostr([agent%xpos(),agent%ypos(),agent%dpos()]) )
          call log_msg( ltag_error // "Agent: " // tostr(ind_real) //         &
                        " in " // tostr(ind_seq) )
        end if
 
        !> #### Get perceptions ####
        !! - Inner perceptions are obtained from the agent's "organism":
        !!   stomach contents, bodymass, energy, age, reproductive factor:
        !!   the_neurobio::perception::perceptions_inner().
        call agent%perceptions_inner()
 
        !> - Simple environmental perceptions are obtained: light, depth:
        !!   the_neurobio::perception::perceptions_environ().
        call agent%perceptions_environ()
 
        !> - Spatial perceptions are obtained for food items, conspecifics and
        !!   predators in proximity of this agent:
        !!   - the_neurobio::perception::see_food()
        !!   - the_neurobio::perception::see_consp()
        !!   - the_neurobio::perception::see_pred()
        !!   .
        call agent%see_food( global_habitats_available( agent_in )%food )
        call agent%see_consp( this%individual )
        call agent%see_pred( global_habitats_available( agent_in )%predators )
 
        !> - The above perceptions are added into the memory stack calling
        !!   the_neurobio::perception::perception_to_memory().
        !! .
        call agent%perception_to_memory()
 
        !> #### Appraisal: Produce motivations ####
        !> - Perception components are calculated for each of the motivational
        !!   states of the agent via the neuronal response function(s):
        !!   the_neurobio::appraisal::motivations_percept_components().
        call agent%motivations_percept_components()
        !> - Then, primary motivation values are calculated from the perception
        !!   components: the_neurobio::appraisal::motivations_primary_calc().
        call agent%motivations_primary_calc()
        !> - The values of the primary motivations are subjected to modulation
        !!   the_neurobio::appraisal::modulation().
        call agent%modulation()
 
        !> - Final motivations of the agent are saved into the emotional
        !!   memory stack: the_neurobio::appraisal::motivations_to_memory().
        !! .
        call agent%motivations_to_memory()
 
        !> #### Determine the Global Organismic State ####
        !> - The Global Organismic State of the agent is calculated using the
        !!   the_neurobio::gos_global::gos_find() method.
        call agent%gos_find()
        call log_dbg( ltag_info //  "Agent " // tostr(agent%get_id()) //      &
                      ", GOS is: " // agent%gos_label() //                    &
                      ", GOS arousal :" // tostr(agent%arousal()),            &
                      procname, modname )
 
        !> - The population-wise maximum motivation value that is used for
        !!   rescaling is calculated now.
        global_rescale_maximum_motivation =                                   &
                maxval( this%individual%motivations%max_perception() )
 
        !> #### Determine and execute the optimal behaviour ####
        !> - Once the GOS is determined for this agent, it selects the
        !!   individually optimal behaviour that minimises its expected
        !!   arousal; this behaviour is then executed. Both steps are
        !!   implemented in the the_neurobio::behaviour::do_behave()
        !!   method.
        !! .
        call agent%do_behave(                                                 &
                      rescale_max_motivation=global_rescale_maximum_motivation )
        call log_dbg( ltag_info // "Agent " // tostr(agent%get_id()) //       &
                      " executed behaviour: " // agent%behaviour_is() //      &
                      " (at global time_step " //                             &
                      tostr(global_time_step_model_current) //                &
                      ")", procname, modname )
 
        !> #### Update characteristics of the agent ####
        !> - Finally, characteristics of this agent are updated depending on
        !!   the consequences of the behaviour. For example, hormone levels
        !!   are updated the_body::condition::sex_steroids_update() and the
        !!   living cost of the agent is subtracted
        !!   (the_body::condition::subtract_living_cost()).
        !!   Digestion also occurs by emptying the stomach by a fixed
        !!   fraction (the_body::condition::stomach_empify(). Also, the
        !!   energy reserves of the agent are updated based on the current
        !!   mass and length (the_body::condition::energy_update().
        !!   Note that the agent characteristics that directly follow from
        !!   the behaviour unit that has been executed (e.g. food gain if the
        !!   agent  decided to eat a food item or travel cost if the agent
        !!   migrated) are processed and updated in the respective behaviour
        !!   execution method.
        call agent%sex_steroids_update()
        call agent%subtract_living_cost()
        call agent%stomach_empify()
        call agent%energy_update()
 
        !> - Finally, the age of the agent is incremented to one time step.
        call agent%age_increment()
 
        !> - Finally, another check is done if the agent is starved to death,
        !!   if yes, the agent individual_genome::dies().
        !! .
        if (agent%starved_death()) then
          call agent%dies()
        end if
 
      end associate
 
    end do individuals
 
  end subroutine population_lifecycle_step_preevol
 
  !-----------------------------------------------------------------------------
  !> This procedure performs a **single step** of the life cycle of the whole
  !! population, the agents for the step are selected in a random order.
  !! @warning This version of the life cycle step includes only optimal food
  !!          selection and eating and does not include the full fledged
  !!          behaviour selection cascade of procedures ::do_behave().
  subroutine population_lifecycle_step_eatonly_preevol(this)
    class(population), intent(inout) :: this
 
    !> ### Notable variables ###
    !> - `inds_order` is an array that sets the order in which the agents are
    !!   being drawn out of the population to perform the step.
    integer, dimension(this%population_size) :: inds_order
 
    !> - `ind_seq` is the sequential number of the agent as drawn from the
    !!   population; it is the loop control counter variable.
    integer :: ind_seq
 
    !> - `ind_real` is the real sequential number of the agent in the
    !!   population; this real id is obtained from the order array
    !!   `inds_order`.
    !! .
    integer :: ind_real
 
    ! Local variable that sets the habitat (number) the agent is currently in
    ! within the global array the_environment::global_habitats_available.
    integer :: agent_in
 
    ! Local variable that keeps the optimal food item
    integer :: food_item_selected
 
    ! PROCNAME is the procedure name for logging and debugging
    character(len=*), parameter :: PROCNAME =                                 &
                                  "(population_lifecycle_step_eatonly_preevol)"
 
    !> ### Implementation notes ###
    !> First, an ordering array `inds_order` is calculated that sets the order
    !! in which the agents are drawn from the population. In the simplest case,
    !! the order of the agents is random, so this array is actually an array of
    !! random integers. The procedure
    !! [PERMUTE_RANDOM](http://ahamodel.uib.no/doc/ar01s09.html#_random_permutation_permute_random_function)
    !! from HEDTOOLS is used here then.
    inds_order = permute_random(this%population_size)
 
    !> The agents can also be processed in any **non-random** order. This would
    !! require invoking an array indexing procedure
    !! [ARRAY_INDEX](http://ahamodel.uib.no/doc/ar01s07.html#_subroutines_array_index_and_array_rank)
    !! instead of `PERMUTE_RANDOM`.
    !!
    !! For example, to process the agents in the order of the body mass,
    !! the ordering array can be obtained obtained as below:
    !! @code
    !! call ARRAY_INDEX(this%individual%get_mass(), inds_order)
    !! @endcode
    !> But this code would rank the agents in an *increasing *order of their
    !! of their body mass. If this is not what is expected, e.g. if the
    !! non-random selection is used to mimic a competitive advantage for bigger
    !! and heavier agents, a reverse of the ordering is obtained like this:
    !! @code
    !! inds_order = inds_order(this%population_size:1:-1)
    !! @endcode
 
    !> Second, the agents are drawn from the population, one by one, in the
    !! named loop construct `INDIVIDUALS`.
    individuals: do ind_seq = 1, popsize
 
      !> - One individual is then drawn from the `inds_order` ordering array.
      ind_real = inds_order(ind_seq)
 
      associate( agent => this%individual(ind_real) )
 
        !> #### Initialisations ####
        !> - First, a check is done if the agent is dead or starved to death;
        !!   if yes, no further processing is done. Also, in the former case
        !!   the individual_genome::dies() method is called.
        if (agent%is_dead()) then
          cycle individuals
        end if
        if (agent%starved_death()) then
          call agent%dies()
          cycle individuals
        end if
 
        !> - `agent_in` index calculates the habitat number where the selected
        !!   agent is currently in, calling the
        !!   the_environment::spatial::find_environment() method.
        agent_in = agent%find_environment()
        if (agent_in<1 .or. agent_in>size(global_habitats_available) ) then
          call log_msg( ltag_error // "agent_in zero in " // procname //      &
                        tostr([agent%xpos(),agent%ypos(),agent%dpos()]) )
          call log_msg( ltag_error // "Agent: " // tostr(ind_real) //         &
                        " in " // tostr(ind_seq) )
        end if
 
        !> #### Get perceptions ####
        !! - Inner perceptions are obtained from the agent's "organism":
        !!   stomach contents, bodymass, energy, age, reproductive factor:
        !!   the_neurobio::perception::perceptions_inner().
        call agent%perceptions_inner()
 
        !> - Simple environmental perceptions are obtained: light, depth:
        !!   the_neurobio::perception::perceptions_environ().
        call agent%perceptions_environ()
 
        !> - Spatial perceptions are obtained for food items, conspecifics and
        !!   predators in proximity of this agent:
        !!   - the_neurobio::perception::see_food()
        !!   - the_neurobio::perception::see_consp()
        !!   - the_neurobio::perception::see_pred()
        !!   .
        call agent%see_food( global_habitats_available( agent_in )%food )
        call agent%see_consp( this%individual )
        call agent%see_pred( global_habitats_available( agent_in )%predators )
 
        !> - The above perceptions are added into the memory stack calling
        !!   the_neurobio::perception::perception_to_memory().
        !! .
        call agent%perception_to_memory()
 
        !> #### Appraisal: Produce motivations ####
        !> - Perception components are calculated for each of the motivational
        !!   states of the agent via the neuronal response function(s):
        !!   the_neurobio::appraisal::motivations_percept_components().
        call agent%motivations_percept_components()
        !> - Then, primary motivation values are calculated from the perception
        !!   components: the_neurobio::appraisal::motivations_primary_calc().
        call agent%motivations_primary_calc()
        !> - The values of the primary motivations are subjected to modulation
        !!   the_neurobio::appraisal::modulation().
        call agent%modulation()
 
        !> - Final motivations of the agent are saved into the emotional
        !!   memory stack: the_neurobio::appraisal::motivations_to_memory().
        !! .
        call agent%motivations_to_memory()
 
        !> #### Determine the Global Organismic State ####
        !> - The Global Organismic State of the agent is calculated using the
        !!   the_neurobio::gos_global::gos_find() method.
        call agent%gos_find()
        call log_dbg( ltag_info //  "Agent " // tostr(agent%get_id()) //      &
                      ", GOS is: " // agent%gos_label() //                    &
                      ", GOS arousal :" // tostr(agent%arousal()),            &
                      procname, modname )
 
        !> - The population-wise maximum motivation value that is used for
        !!   rescaling is calculated now.
        global_rescale_maximum_motivation =                                   &
                maxval( this%individual%motivations%max_perception() )
 
        !> #### Determine and execute the optimal behaviour ####
        !> - Once the GOS is determined for this agent, it selects the optimal
        !!   food item that minimises its expected arousal and then tries to
        !!   eat this food item. However, if the agent has no food in its
        !!   perception a default random walk is executed.
        !! .
        do_behave: if ( agent%has_food() ) then
          food_item_selected = agent%food_item_select(                        &
                      rescale_max_motivation=global_rescale_maximum_motivation)
          ! ++++ EAT FOOD ITEM
          call agent%do_eat_food_item(                                        &
                food_item_selected, global_habitats_available( agent_in )%food )
        else do_behave
          call agent%do_walk()
        end if do_behave
 
        call log_dbg( ltag_info // "Agent " // tostr(agent%get_id()) //       &
                      " executed behaviour: " // agent%behaviour_is() //      &
                      " (at global time_step " //                             &
                      tostr(global_time_step_model_current) //                &
                      ")", procname, modname )
 
        !> #### Update characteristics of the agent ####
        !> - Finally, characteristics of this agent are updated depending on
        !!   the consequences of the behaviour. For example, hormone levels
        !!   are updated and the living cost of the agent is subtracted.
        !!   Note that the agent characteristics that directly follow from
        !!   the behaviour unit that has been executed (e.g. food gain if the
        !!   agent  decided to eat a food item or travel cost if the agent
        !!   migrated) are processed and updated in the respective behaviour
        !!   execution method.
        call agent%sex_steroids_update()
        call agent%subtract_living_cost()
        call agent%energy_update()
 
        !> - Finally, the age of the agent is incremented to one time step.
        call agent%age_increment()
 
        !> - Finally, another check is done if the agent is starved to death,
        !!   if yes, the agent individual_genome::dies().
        !! .
        if (agent%starved_death()) then
          call agent%dies()
        end if
 
      end associate
 
    end do individuals
 
  end subroutine population_lifecycle_step_eatonly_preevol
 
  !-----------------------------------------------------------------------------
  !> Save data for all agents within the population into a CSV file.
  !> @note Note that this procedure is not using the @ref file_io wrappers.
  subroutine population_save_data_all_agents_csv(this, csv_file_name,         &
                                        save_header, is_logging, is_success)
    use csv_io
    class(population), intent(in) :: this
    !> @param[in] csv_file_name is the name of the CSV output file.
    character(len=*), intent(in) :: csv_file_name
    !> @param[in] save_header turn ON/OFF of the descriptive file header.
    !!            Header is saved into the first row of the CSV output file
    !!            If not present, default is FALSE.
    logical, optional, intent(in) :: save_header
    !> @param[in] is_logging turn ON/OFF writing the file name and data into
    !!            the logger. If not present, default is TRUE if it is the
    !!            debug mode.
    logical, optional, intent(in) :: is_logging
    !> @param[out] is_success Flag showing that data save was successful
    !!             (if TRUE).
    logical, optional, intent(out) :: is_success
 
    ! Local copies of optionals.
    logical :: logging_enabled
 
    ! Counter
    integer :: ind
 
    !> ### Implementation notes ###
    !> #### Local variables for CSV backend ####
    !> - `handle_csv` is the CSV file handle object defining the file name,
    !!   Fortran unit and error descriptor, see HEDTOOLS manual for details.
    type(csv_file) :: handle_csv
    !> - `csv_record_tmp` is the temporary character string that keeps the
    !!   whole record of the file, i.e. the whole row of the spreadsheet table.
    character(len=:), allocatable :: csv_record_tmp
    !> - `COLUMNS` is a parameter array that keeps all column headers; its
    !!   size is equal to the total number of variables (columns) in the data
    !!   spreadsheet file.
    !! .
    character(len=LABEL_LENGTH), dimension(*),                       &
      parameter :: COLUMNS =  [ character(len=label_length) ::       &  ! COLS
                                                      "ID_NUM   ",   &  !  1
                                                      "PERS_NAME",   &  !  2
                                                      "ALIVE    ",   &  !  3
                                                      "SEX_MALE ",   &  !  4
                                                      "BODY_LEN ",   &  !  5
                                                      "BIRTH_LEN",   &  !  6
                                                      "CTRL_RND ",   &  !  7
                                                      "BODY_MASS",   &  !  8
                                                      "BIRTHMASS",   &  !  9
                                                      "ENERGY   ",   &  ! 10
                                                      "BIRT_ENER",   &  ! 11
                                                      "STOMACH  ",   &  ! 12
                                                      "MAXSTOMCP",   &  ! 13
                                                      "SMR      ",   &  ! 14
                                                      "S_COST_SW",   &  ! 15
                                                      "LIV_COST ",   &  ! 16
                                                      "MORTALITY",   &  ! 17
                                                      "HORM_GROW",   &  ! 18
                                                      "HORM_THYR",   &  ! 19
                                                      "HORM_ADRE",   &  ! 20
                                                      "HORM_CORT",   &  ! 21
                                                      "HORM_TEST",   &  ! 22
                                                      "HORM_ESTR",   &  ! 23
                                                      "HTEST_BAS",   &  ! 24
                                                      "HESTR_BAS",   &  ! 25
                                                      "REPR_FAC ",   &  ! 26
                                                      "P_REPROD ",   &  ! 27
                                                      "N_REPROD ",   &  ! 28
                                                      "N_OFFSPNG",   &  ! 29
                                                      "AGE      ",   &  ! 30
                                                      "GOS_MAIN ",   &  ! 31
                                                      "GOS_ARUSL",   &  ! 32
                                                      "GOS_REPET",   &  ! 33
                                                      "EAT_ATMPT",   &  ! 34
                                                      "N_EATEN  ",   &  ! 35
                                                      "MASSEATEN",   &  ! 36
                                                      "PERC_FOOD",   &  ! 37
                                                      "PERC_CONS",   &  ! 38
                                                      "PERC_PRED",   &  ! 39
                                                      "POS_X    ",   &  ! 40
                                                      "POS_Y    ",   &  ! 41
                                                      "POS_DEPTH",   &  ! 42
                                                      "HABITAT  ",   &  ! 43
                                                      "FITNNESS "  ]    ! 44
 
    ! PROCNAME is the procedure name for logging and debugging
    character(len=*), parameter ::                                            &
                            PROCNAME = "(population_save_data_all_agents_csv)"
 
    ! Local, the name of the habitat the agent is currently in, out of the
    ! commondata::global_habitats_available array.
    character(len=LABEL_LENGTH) :: habitat_in
 
    ! Local timer object for file write.
    type(timer_cpu) :: file_write_timing
 
    if (present(is_logging)) then
      logging_enabled = is_logging
    else
      logging_enabled = is_debug
    end if
 
    if (logging_enabled) then
      call log_msg (ltag_info // "Saving all individuals in population # " // &
              tostr(this%pop_number) // " (name '" // trim(this%pop_name) //  &
              "'), " //                                                       &
              "generation # " // tostr(global_generation_number_current) //   &
              ", time step " // tostr(global_time_step_model_current) //      &
              " to file: " // csv_file_name )
      call file_write_timing%start(                                           &
                  "Writing population data for population '" //               &
                  trim(this%pop_name) // "' " // tostr(size(this%individual)) &
                  // " individuals" )
    end if
 
    !> #### Save data in CSV file ####
    !> @note Note that this subroutine does not use the object oriented
    !!       wrappers from the @ref file_io module.
    !!
    !> - Define the file name \%name component of the CSV file handle. The
    !!   file handle object `handle_csv` is now used as the file identifier.
    handle_csv%name = csv_file_name
    !> - Open the output file defined by the `handle_csv` handle object for
    !!   writing.
    call csv_open_write( handle_csv )
 
    !> - Possible error status of the latest file operation is obtained by the
    !!   \%status component of the file handle. Check if there were any errors
    !!   opening the file and report in the logger with the error tag.
    if ( .not. handle_csv%status ) then
        call log_msg( ltag_error // "Opening output CSV file FAILED: " //     &
                      csv_file_name // ", in " // procname )
        call log_msg( ltag_error // "Data file " // csv_file_name //          &
                      " is not written in " // procname )
        if (present(is_success)) is_success = handle_csv%status
        return
    end if
 
    !> - If the `save_header` flag is set to TRUE, save the CSV file header.
    if (present(save_header)) then
      if (save_header) then
        call csv_header_write( "Population: " // this%pop_name, handle_csv )
        if ( .not. handle_csv%status ) then    ! treat possible write error.
          if (present(is_success)) is_success = .false.
          call csv_close( handle_csv )
          return
        end if
      end if
    end if
 
    !> - Prepare the character string variable `csv_record_tmp` that keeps the
    !!   whole record (row) of data in the output CSV data file. The length of
    !!   this string should be enough to fit all the record data, otherwise
    !!   the record is truncated.
    csv_record_tmp = repeat(" ", size(columns) * len(columns(1)) )
 
    !> - Produce the first record containing the column headers (variable
    !!   names). Note that
    !!   [CSV_RECORD_APPEND()](http://ahamodel.uib.no/doc/ar01s08.html#_subroutine_csv_record_append)
    !!   accepts both arrays and scalar values for appending. Also, write the
    !!   first record physically to the file.
    call csv_record_append( csv_record_tmp, columns )
    call csv_record_write ( csv_record_tmp, handle_csv )
    if ( .not. handle_csv%status ) then    ! treat possible write error.
      if (present(is_success)) is_success = .false.
      call csv_close( handle_csv )
      return
    end if
 
 
    !> - The actual data are written to the CSV file in a loop over all the
    !!   individual members of the population. One record (row) of the data
    !!   file then represents a single individual.
    do ind = 1, size(this%individual)
      !>   - the `csv_record_tmp` character string variable is produced such
      !!     that it can fit the whole record;
      csv_record_tmp = repeat(" ",                                            &
                      max( csv_guess_record_length(size(columns) + 1,0.0_srp),&
                            len(this%individual(ind)%genome_label) ) )
 
      !>   - the actual data for the individual is appended to the current
      !!     record one by one. Note that logical values are converted to
      !!     integers using commondata::conv_l2r() function.
      associate( agent => this%individual(ind) )                          !COLS
        call csv_record_append(csv_record_tmp,agent%person_number        ) ! 1
        call csv_record_append(csv_record_tmp,agent%genome_label         ) ! 2
        call csv_record_append(csv_record_tmp,conv_l2r(agent%alive)      ) ! 3*
        call csv_record_append(csv_record_tmp,conv_l2r(agent%sex_is_male)) ! 4*
        call csv_record_append(csv_record_tmp,agent%body_length          ) ! 5
        call csv_record_append(csv_record_tmp,agent%body_length_birth    ) ! 6
        call csv_record_append(csv_record_tmp,agent%control_unselected   ) ! 7
        call csv_record_append(csv_record_tmp,agent%body_mass            ) ! 8
        call csv_record_append(csv_record_tmp,agent%body_mass_birth      ) ! 9
        call csv_record_append(csv_record_tmp,agent%energy_current       ) ! 10
        call csv_record_append(csv_record_tmp,agent%energy_birth         ) ! 11
        call csv_record_append(csv_record_tmp,agent%stomach_content_mass ) ! 12
        call csv_record_append(csv_record_tmp,agent%maxstomcap           ) ! 13
        call csv_record_append(csv_record_tmp,agent%smr                  ) ! 14
        call csv_record_append(csv_record_tmp,agent%cost_swim_std()      ) ! 15
        call csv_record_append(csv_record_tmp,agent%living_cost()        ) ! 16
        call csv_record_append(csv_record_tmp,agent%ind_mortality        ) ! 17
        call csv_record_append(csv_record_tmp,agent%growhorm_level       ) ! 18
        call csv_record_append(csv_record_tmp,agent%thyroid_level        ) ! 19
        call csv_record_append(csv_record_tmp,agent%adrenaline_level     ) ! 20
        call csv_record_append(csv_record_tmp,agent%cortisol_level       ) ! 21
        call csv_record_append(csv_record_tmp,agent%testosterone_level   ) ! 22
        call csv_record_append(csv_record_tmp,agent%estrogen_level       ) ! 23
        call csv_record_append(csv_record_tmp,agent%testosterone_baseline) ! 24
        call csv_record_append(csv_record_tmp,agent%estrogen_baseline    ) ! 25
        if ( agent%is_male() ) then  ! repr. factor in males and females
          call csv_record_append(csv_record_tmp,agent%testosterone_level )
        else
          call csv_record_append(csv_record_tmp,agent%estrogen_baseline  ) ! 26
        endif
        call csv_record_append(csv_record_tmp,                           &
                                        agent%probability_reproduction() ) ! 27
        call csv_record_append(csv_record_tmp,agent%n_reproductions      ) ! 28
        call csv_record_append(csv_record_tmp,agent%n_offspring          ) ! 29
        call csv_record_append(csv_record_tmp,agent%age                  ) ! 30
        call csv_record_append(csv_record_tmp,agent%gos_main             ) ! 31
        call csv_record_append(csv_record_tmp,agent%gos_arousal          ) ! 32
        call csv_record_append(csv_record_tmp,agent%gos_repeated         ) ! 33
        call csv_record_append(csv_record_tmp,agent%n_eats_all_indicator ) ! 34
        call csv_record_append(csv_record_tmp,agent%n_eaten_indicator    ) ! 35
        call csv_record_append(csv_record_tmp,agent%mass_eaten_indicator ) ! 36
        call csv_record_append(csv_record_tmp,                           &
                                  agent%memory_stack%get_food_mean_n()   ) ! 37
        call csv_record_append(csv_record_tmp,                           &
                                  agent%memory_stack%get_consp_mean_n()  ) ! 38
        call csv_record_append(csv_record_tmp,                           &
                                  agent%memory_stack%get_pred_mean()     ) ! 39
        call csv_record_append(csv_record_tmp,agent%x                    ) ! 40
        call csv_record_append(csv_record_tmp,agent%y                    ) ! 41
        call csv_record_append(csv_record_tmp,agent%depth                ) ! 42
        if ( agent%is_alive() ) then
          habitat_in = global_habitats_available(                        &
                                    agent%find_environment())%get_label()
        else
          ! If the agent is dead and nullified, its habitat cannot be
          ! determined.
          habitat_in = "agent_dead"
        end if
        call csv_record_append(csv_record_tmp,habitat_in                 ) ! 43
        call csv_record_append(csv_record_tmp,agent%fitness              ) ! 44
      end associate
 
      !>   - after all data are appended to the record, this record is
      !!     physically written to the disk using
      !!     [CSV_RECORD_WRITE()](http://ahamodel.uib.no/doc/ar01s08.html#_subroutine_csv_record_write).
      !!   .
      call csv_record_write( csv_record_tmp, handle_csv )
      if ( .not. handle_csv%status ) then    ! treat possible write error.
        if (present(is_success)) is_success = .false.
        call csv_close( handle_csv )
        return
      end if
 
    end do
 
    !> - When all the records are saved, the CSV file is closed with
    !!   [CSV_CLOSE()](http://ahamodel.uib.no/doc/ar01s08.html#_subroutine_csv_close).
    call csv_close( handle_csv )
 
    !> - This is finally sent to the logger (if `logging_enabled` is TRUE).
    !! .
    if (logging_enabled) then
      call log_msg (ltag_info // "Individual data saved, population size " // &
              tostr(this%population_size) //                                  &
              ", number of columns " // tostr(size(columns)) )
      if ( .not. handle_csv%status ) call log_msg( ltag_error //              &
                                              "File write operation FAILED." )
      call log_msg( file_write_timing%log() )
    end if
 
    if (present(is_success)) is_success = handle_csv%status
 
    !> The CSV output data file can be optionally compressed with the
    !! commondata::cmd_zip_output command if commondata::is_zip_outputs is set
    !! to TRUE.
    if ( is_zip_outputs ) then
      call call_external(command=cmd_zip_output // " " // csv_file_name,      &
                         suppress_output=.true.,                              &
                         is_background_task=zip_outputs_background )
    end if
 
  end subroutine population_save_data_all_agents_csv
 
  !-----------------------------------------------------------------------------
  !> Save the genome data of all agents in this population to a CSV file.
  subroutine population_save_data_all_genomes(this, csv_file_name, is_success)
    use file_io
    use csv_io
    class(population), intent(in) :: this
    !> @param[in] csv_file_name is the name of the CSV output file.
    character(len=*), intent(in) :: csv_file_name
    !> @param[out] is_success Flag showing that data save was successful
    !!             (if TRUE).
    logical, optional, intent(out) :: is_success
 
    ! Local counters.
    integer :: i, j, k, l, m
 
    ! File handle object.
    type(file_handle) :: genome_file
 
    ! **Column names**
    ! The column names are built from the components below with the
    ! consecutive digits indicating the order of each component.
    character(len=*), parameter :: TAG_CRO = "CRO_"   ! Chromosome
    character(len=*), parameter :: TAG_GAP = "_"      ! a single gap
    character(len=*), parameter :: TAG_ALE = "_ALE_"  ! Allele
    character(len=*), parameter :: TAG_ACO = "_AC_"   ! Allele component
 
    ! The length of the digit part of the column name, assumed two for the
    ! normal range 1:99
    integer, parameter :: DIG_LEN = 2
 
    ! Column length hard limit
    integer, parameter :: COL_LEN = len(tag_cro) + dig_len +                  &
                                    len(tag_gap) + dig_len +                  &
                                    len(tag_ale) + dig_len +                  &
                                    len(tag_aco) + dig_len
 
    ! Column (variable) headers.
    character(len=COL_LEN), allocatable, dimension(:) :: colname
 
    ! The number of columns in the output data file.
    integer :: n_columns
 
    ! *Record* that keeps the whole row of data. See @ref file_io.
    character(len=:), allocatable :: record_csv
 
    ! Separate variables were necessary to work around ifort 17 compiler bug:
    ! it crashed when everything was placed inline into the CSV_RECORD_APPEND():
    ! HEDG2_04.f90(45689): internal error: Please visit
    ! 'http://www.intel.com/software/products/support' for assistance.
    !       record_csv = repeat(" ", len(TOSTR(ALLELERANGE_MAX)) * n_columns +
    ! ^
    ! [ Aborting due to internal error. ]
    integer :: record_csv_max_length
    integer :: allele_val_append
 
    !> ### Implementation notes ###
    !> #### Build column names ####
    !> First, determine the number of columns in the data file. The number of
    !! columns is calculated by unwinding the whole data structure
    !! - `chromosome(j,k)%allele(l)%allele_value(m)`
    !! .
    !! See @ref the_genome for more details on the data structure.
    n_columns=0
    do j=1, n_chromosomes
      do k=1, chromosome_ploidy
        do l=1, len_chromosomes(j)
          do m=1, additive_comps
            n_columns = n_columns + 1
          end do
        end do
      end do
    end do
 
    !> The array of the column names is then allocated to the above number.
    allocate(colname(n_columns))
 
    !> The column names are built again by unwinding the whole genome data
    !! structure into a linear sequence. The column names are like this:
    !!
    !! `CRO_1_1_ALE_01_AC_1, CRO_1_1_ALE_01_AC_2, CRO_1_1_ALE_01_AC_3, ... `
    i=1
    do j=1, n_chromosomes
      do k=1, chromosome_ploidy
        do l=1, len_chromosomes(j)
          do m=1, additive_comps
            ! CRO_1_1_ALE_01_AC_1, CRO_1_1_ALE_01_AC_2, CRO_1_1_ALE_01_AC_3, ...
            colname(i) = tag_cro // tostr(j) //                               &
                          tag_gap // tostr(k) //                              &
                            tag_ale // tostr(l,10) //                         &
                              tag_aco // tostr(m)
            i=i+1
          end do
        end do
      end do
    end do
 
    !> #### Write data to the file ####
    !> First, the file is opened for writing.
    call genome_file%open_write( csv_file_name, format_csv )
    if ( .not. genome_file%is_success() ) then
      if (present(is_success)) is_success = .false.
      call genome_file%close()
      return
    end if
 
    !> The first record of the data that contains the column names is then
    !! "appended" into the complete record and written to the file.
    !! The length of this record is calculated based on the length of the
    !! columns and their number.
    record_csv = repeat( " ", label_length * 2 + 2 * 3 +                      &
                              col_len * n_columns + n_columns * 3  )
    !! - The first two columns contain the identifiers for each agent:
    !!   - numeric ID of the agent
    !!   - test string label ("name") of the agent
    !!   .
    call csv_record_append( record_csv,                                       &
                    [ character(len=LABEL_LENGTH) :: "ID_NUM", "AGENT_NAME"] )
 
    !> - The remaining columns contain the chromosome and gene labels
    !!   (see above);
    !! .
    call csv_record_append( record_csv, colname )
 
    !> This first line consisting of column names is then written to the
    !! output file.
    call genome_file%record_write( record_csv )
    if ( .not. genome_file%is_success() ) then
      if (present(is_success)) is_success = .false.
      call genome_file%close()
      return
    end if
 
    !> The maximum length of the data record is calculated as the maximum
    !! string length of a single data value multiplied by the number of
    !! columns. Because the record also adds separators, the number of columns
    !! multiplied by three is added to this value.
    record_csv_max_length = label_length * 2 + 2 * 3 +                        &
                            len(tostr(allelerange_max)) * n_columns +         &
                                                                n_columns * 3
 
    !> Finally, cycle over all individuals in this population and
    !! save the genome data. The first two columns are
    !! the_genome::individual_genome::person_number and
    !! the_genome::individual_genome::genome_label. The other columns
    !! "unwind" the genome data structure over the inner loops for
    !! chromosomes, homologues, alleles and allele components.
    !! - - `chromosome(j,k)%allele(l)%allele_value(m)`
    !! .
    !! See @ref the_genome for more details on the data structure.
    inds: do i = 1, this%population_size
 
      record_csv = repeat(" ", record_csv_max_length )
 
      call csv_record_append( record_csv, this%individual(i)%person_number )
      call csv_record_append( record_csv, this%individual(i)%genome_label )
 
      genome: do j=1, n_chromosomes
                do k=1, chromosome_ploidy
                  do l=1, len_chromosomes(j)
                    do m=1, additive_comps
                      allele_val_append = this%individual(i)%                 &
                                  chromosome(j,k)%allele(l)%allele_value(m)
                      call csv_record_append( record_csv, allele_val_append )
                    end do
                  end do
                end do
              end do genome
      call genome_file%record_write( record_csv )
      if ( .not. genome_file%is_success() ) then
        if (present(is_success)) is_success = .false.
        call genome_file%close()
        return
      end if
 
    end do inds
 
    !> Once all individuals are saved, the file is closed.
    call genome_file%close()
    if ( .not. genome_file%is_success() ) then
      if (present(is_success)) is_success = .false.
    else
      if (present(is_success)) is_success = .true.
    end if
 
    !> The CSV output data file can be optionally compressed with the
    !! commondata::cmd_zip_output command if commondata::is_zip_outputs is set
    !! to TRUE.
    if ( is_zip_outputs ) then
      call call_external(command=cmd_zip_output // " " // csv_file_name,      &
                         suppress_output=.true.,                              &
                         is_background_task=zip_outputs_background )
    end if
 
  end subroutine population_save_data_all_genomes
 
  !-----------------------------------------------------------------------------
  !> Load the genome data of all agents in this population from a CSV file.
  !! Note that the procedure implements several error correcting measures,
  !! e.g. checks for minimum number of rows in the file and minimum row length.
  !! The input CSV file therefore can include short text notes that are then
  !! ignored when reading data.
  subroutine population_load_data_all_genomes( this, pop_size,                &
    pop_number_here, pop_name_here, csv_file_name, missing_random, is_success )
 
    use csv_io
    use base_strings, only : value, parse, split, compact, is_numeric, delall
 
    class(population), intent(inout) :: this
    !> @param pop_size, The size of the population.
    integer, intent(in) :: pop_size
    !> @param id, Optional numeric population ID.
    integer, optional, intent(in) :: pop_number_here
    !> @param descriptor, Optional population string descriptor.
    character (len=*), optional, intent(in) :: pop_name_here
    !> @param[in] csv_file_name is the name of the CSV output file.
    character(len=*), intent(in) :: csv_file_name
    !> @param[in] missing_random Flag dictating to generate random genomes
    !!            if the commondata::popsize parameter exceeds the population
    !!            size in the CSV genome data size. TRUE = generate random
    !!            genome, FALSE = leave empty data commondata::unknown
    logical, optional, intent(in) :: missing_random
    !> @param[out] is_success Flag showing that data save was successful
    !!             (if TRUE).
    logical, optional, intent(out) :: is_success
 
    ! Local copies of optionals
    integer :: pop_number_here_loc
    character (len=:), allocatable :: pop_name_here_loc
    logical :: missing_random_loc
 
    ! PROCNAME is the procedure name for logging and debugging
    character(len=*), parameter ::                                            &
                            PROCNAME = "(population_load_data_all_genomes)"
 
    !> - `handle_csv` is the CSV file handle object defining the file name,
    !!   Fortran unit and error descriptor, see HEDTOOLS manual for details.
    type(csv_file) :: handle_csv
 
    !> - `line_data_buff` is the input data buffer that is being read
    !!    via the `CSV_IO` procedure `READLINE()`
    character(len=:), allocatable :: line_data_buff
 
    !> - `n_rows_in` is the number of rows in the input CSV file excluding
    !!    the first row containing the variable names.
    integer :: n_rows_in
 
    !> - `N_ROWS_MIN` is the minimum number of readable rows in the CSV
    !!    genome data file, if there are less in the file, data are considered
    !!    invalid.
    integer, parameter :: N_ROWS_MIN = 10
 
    !> - `MIN_FILE_INP_LEN` is the minimum length of the data row to be
    !!    included in the read genome.
    integer, parameter :: MIN_FILE_INP_LEN = label_length +                   &
                    n_chromosomes * chromosome_ploidy * 2 * additive_comps * 2
 
    ! Local logical flags
    logical :: success_flag, is_new_population, is_pop_id_reset
 
    ! Delimiters for data fields:
    character, parameter :: TAB =achar(9)
    character, parameter :: DQT =achar(34) ! double quote "
    character, parameter :: SQT =achar(39) ! single quote '
    character(len=*), parameter :: SDELIM = " ," // tab
 
    ! Length of a single field within record. Must be big enough to fit
    ! ID label and all numeric data fields.
    integer, parameter :: LEN_CSV_FIELD = label_length * 2
    ! Minimum length of a data field.
    integer, parameter :: MIN_FIELD = 1
 
    !> - `line_data_substrings` - an array of row string data values after
    !!    parding the whole input raw data.
    character(len=LEN_CSV_FIELD), dimension(:), allocatable  ::               &
                                                          line_data_substrings
 
    ! Local counters
    integer :: icase, nvars, jfield, ivalue, error_iflag,                     &
               line_data_buff_length, line_data_nflds, n_ignored_lns
 
    ! Local counters for loops
    integer :: i, j, k, l, m
 
    !> - `matrix_row` -- array representing the genome data for a single row
    !!    after parsing the input data file row.
    integer, allocatable, dimension(:) :: matrix_row
 
    !> ### Implementation notes ###
    !! Genome data are obtained from the CSV data file provided by the
    !! `intent(in)` parameter `csv_file_name`.
    call log_delimiter(log_level_chapter)
    call log_msg( ltag_stage // "Reading genome data from file " //           &
                  csv_file_name // ", file rows less than " //                &
                  tostr(min_file_inp_len) // " ignored." )
 
    missing_random_loc = .false.
    if (present(missing_random)) then
      if (missing_random .eqv. .true. ) then
        missing_random_loc = .true.
      end if
    end if
 
    !> - First, calculate the number of records in the input CSV file and check
    !!   that this number is greater than the minimum value set by the
    !!   local parameter `N_ROWS_MIN`.  Also check if the file can be read
    !!   at all. If either of these conditions hold, exit from the procedure
    !!   with
    n_rows_in = csv_file_lines_count( csv_file_name, numeric_only=.false.,    &
                                      csv_file_status=success_flag) - 1
    if ( n_rows_in < n_rows_min .or. (success_flag .eqv. .false.) ) then
      call log_write_error("Input CSV file has too few rows: "//              &
                           tostr(n_rows_in) // " or file access error, " //   &
                           " success flag is " // tostr(success_flag) )
      call log_msg(ltag_warn // "Check N_ROWS_MIN (=" // tostr(n_rows_min) // &
                   ") in " // procname )
      return
    end if
    call log_msg( ltag_info // "Genome data " // csv_file_name // " contain " &
                  // tostr(n_rows_in) // " rows" )
    if ( n_rows_in < pop_size ) call log_msg( ltag_warn //                    &
               "Number of rows in genome data " // csv_file_name // "," //    &
               tostr(n_rows_in) // " is less than input pop_size " //         &
               tostr(pop_size) )
 
    !> - Second, open the input data file for reading
    !>    - Define the file name \%name component of the CSV file handle. The
    !!      file handle object `handle_csv` is now used as the file identifier.
    handle_csv%name = csv_file_name
    !>    - Open the output file defined by the `handle_csv` handle object for
    !!      reading.
    call csv_open_read( handle_csv )
 
    !>    - Possible error status of the latest file operation is obtained by
    !!      the \%status component of the file handle. Check if there were
    !!      any errors opening the file and report in the logger with the
    !!      error tag.
    if ( .not. handle_csv%status ) then
      call log_write_error("Opening input CSV file FAILED")
      return
    end if
 
    !> - Third, read the variable names line of the CSV file. Note that
    !!   they are not used here.
    if ( readline(handle_csv%unit, line_data_buff, .true.) .eqv. .false. ) then
      call log_write_error("Reading input CSV file FAILED")
      call csv_close( handle_csv )
      return
    end if
 
    !> - **HAVE_POP** check if the population is allocated and the population
    !!   size is set, if not, allocate and set to commondata::popsize.
    have_pop: if (.not. allocated(this%individual)) then
      is_new_population = .true.
      !>    - Check and reset population IDs if present
      if (present(pop_number_here)) then
        pop_number_here_loc = pop_number_here
      else
        pop_number_here_loc = rand_i(1,huge(this%pop_number-1))
      end if
      if (present(pop_name_here)) then
        pop_name_here_loc = pop_name_here
      else
        pop_name_here_loc = tostr(pop_number_here_loc)
      end if
      !>    - Initialise population with the size from input parameter.
      call this%init(pop_size,pop_number_here_loc,pop_name_here_loc)
      !>    - Note that if the population is not allocated, all individuals
      !!      that are not obtained from the data file are initialised random.
      missing_random_loc = .true.
      call log_msg( ltag_info // "Initialised population to POPSIZE " //      &
                    tostr(this%population_size) // ", new population: " //    &
                    tostr(is_new_population) )
    else
      is_new_population = .false.
      call log_msg( ltag_info // "Population is already allocated POPSIZE="   &
                    // tostr(this%population_size) )
      !>    - Check and reset population IDs if present
      is_pop_id_reset = .false.
      if (present(pop_number_here)) then
        this%pop_number = pop_number_here
        is_pop_id_reset = .true.
      end if
      if (present(pop_name_here)) then
        this%pop_name = pop_name_here
        is_pop_id_reset = .true.
      end if
      if (is_pop_id_reset)                                                    &
        call log_msg( ltag_info // "Reset population identifiers " //         &
                      "from parameters: " // tostr(this%pop_number) //        &
                      ", " // trim(this%pop_name) // ", new population: " //  &
                      tostr(is_new_population) )
    end if have_pop
 
    call log_msg(ltag_info // "Start reading genome data for "//              &
                 tostr(this%pop_number) // ", " // trim(this%pop_name) //     &
                 ", population array size " // tostr(size(this%individual)) )
 
    !> - **MAIN_READ:** Cycle through the data lines using the non-advancing
    !!   `READLINE` function and get the numerical genome data.
    icase = 0
    n_ignored_lns = 0
    main_read: do while ( readline(handle_csv%unit, line_data_buff, .true.)   &
                          .and. icase < this%population_size                  &
                          .and. icase < pop_size )
 
      !>     - All input file lines that contain data shorter than
      !!       `MIN_FILE_INP_LEN` characters are ignored.
      if ( len_trim(line_data_buff) < min_file_inp_len ) then
        n_ignored_lns = n_ignored_lns + 1
        call log_dbg( ltag_info // "Ignored line " // tostr(icase) // ": " // &
                      trim(line_data_buff), procname, modname )
        cycle main_read
      end if
 
      icase = icase + 1
 
      !>    - Stripe away single and double quotes if any.
      call delall(line_data_buff, sqt)
      call delall(line_data_buff, dqt)
 
      line_data_buff_length = len_trim(line_data_buff)
 
      !>    - An approximate upper bound for the number of data fields in the
      !!      current record is `len_trim(line_data_buff)/MIN_FIELD`), so we
      !!      can now allocate the temporary array of substrings. This is
      !!      also the total number of variables including the first two
      !!      non-data records.
      if ( .not. allocated(line_data_substrings) )                            &
        allocate( line_data_substrings(line_data_buff_length/min_field) )
 
      !>    - Split the input string buffer `line_data_buff` into an array of
      !!      individual strings representing each file data field. The number
      !!      of such non-empty strings (fields) is given by `line_data_nflds`.
      call parse(line_data_buff, sdelim, line_data_substrings, line_data_nflds)
 
      !>    - **N_GCOLS:** Check if `line_data_nflds` obtained from CSV data
      !!      coincides with the model (`this` data structure) as defined by
      !!      the number of chromosomes, ploidy, chromosome length and N of
      !!      additive components:
      !!        - commondata::n_chromosomes,
      !!        - commondata::chromosome_ploidy
      !!        - commondata::len_chromosomes
      !!        - commondata::additive_comps
      !!        .
      nvars = 0
      n_gcols: do j=1, n_chromosomes
                 do k=1, chromosome_ploidy
                   do l=1, len_chromosomes(j)
                     do m=1, additive_comps
                       nvars = nvars + 1
                     end do
                   end do
                 end do
               end do n_gcols
      if ( line_data_nflds - 2 /= nvars ) then
        call log_write_error( ltag_error // "Number of fields in data file "  &
                          // tostr(line_data_nflds) //                        &
                        " is unequal to N of model data: " // tostr(nvars) // &
                        " for case (row) " // tostr(icase) )
        if (present(is_success)) is_success = .false.
        return
      end if
 
      !>     - We have to allocate the array of row values and initialise
      !!       it for correct processing by the `VALUE` subroutine.
      if (.not. allocated(matrix_row)) allocate(matrix_row(line_data_nflds-2))
      matrix_row = unknown
 
      !>     - first value from the file row is `ID_NUM`
      if (is_numeric(line_data_substrings(1))) then
        call value( trim(line_data_substrings(1)), ivalue, error_iflag)
      else
        call log_write_error("Error reading ID_NUM field, row " // tostr(icase))
        ivalue = unknown
      end if
      this%individual(icase)%person_number = ivalue
 
      !>     - second value from the file row is `AGENT_NAME`.
      this%individual(icase)%genome_label = trim(line_data_substrings(2))
 
      !>     - now read a single row of the genome data for this (`icase`)
      !!       individual agent, note that if any data in the file are
      !!       abridged, random values are generated within the appropriate
      !!       genome limits (commondata::allelerange_min,
      !!       commondata::allelerange_max).
      do jfield=3, line_data_nflds
        if ( trim(line_data_substrings(jfield))=="" )  then
          !matrix_row(jfield-2:line_data_nflds-2) = UNKNOWN
          call rand_array( matrix_row(jfield-2:line_data_nflds-2),            &
                           allelerange_min, allelerange_max )
          call log_msg( ltag_warn // "Incomplete data in field " //           &
                        tostr(jfield) // " , row " // tostr(icase) //         &
                        "; set random array: " //                             &
                        tostr(matrix_row(jfield-2:line_data_nflds-2))  )
          cycle
        end if
        if (is_numeric(line_data_substrings(jfield))) then
          call value(trim(line_data_substrings(jfield)), matrix_row(jfield-2),&
                     error_iflag)
        else
          !matrix_row(jfield-2) = UNKNOWN
          matrix_row(jfield-2) = rand_i(allelerange_min, allelerange_max)
          call log_msg( ltag_warn // "Wrong non-numeric data in row " //      &
                        tostr(jfield) // " , row " // tostr(icase) //         &
                        ": set random " // tostr(matrix_row(jfield-2)) )
        end if
      end do
 
      !>     - Finally, split (parse) the single line array for the `icase`-s
      !!       individual into the genome data structure.
      jfield = 0
      genome: do j=1, n_chromosomes
                do k=1, chromosome_ploidy
                  do l=1, len_chromosomes(j)
                    do m=1, additive_comps
                      jfield = jfield + 1
                      this%individual(icase)%chromosome(j,k)%allele(l)%allele_value(m) =  matrix_row(jfield)
                    end do
                  end do
                end do
              end do genome
      call log_dbg( ltag_info // "Read genome data case " // tostr(icase) )
 
      deallocate( line_data_substrings )
 
    end do main_read
 
    call log_msg( ltag_info // "All " // tostr(icase) //                      &
                  " rows read from data file " // handle_csv%name )
 
    !> - **FILL_EXTRA:** If the number of rows in the data file read is smaller
    !!   than the population size commondata::popsize, the remaining data may
    !!   need to be filled with some default values. The flag `missing_random`
    !!   determines to initialise all the remaining individuals in the
    !!   population as random.
    fill_extra: if ( icase < this%population_size .and.                       &
                     (is_new_population .eqv. .false.) ) then
 
      call log_msg( ltag_info // "Processing extra data" )
      if ( missing_random_loc .eqv. .true. ) then
          !>    Initialise all individuals of the population
          do i=icase+1, this%population_size
            !>       - first, call `init to object individual(i)
            call this%individual(i)%init()
            !!       - second, ! call `set_id(i)` to identify agents
            !!         within the population.
            !!       .
            call this%individual(i)%set_id(i)
          end do
          call log_msg( ltag_info // "Initialised random cases from " //      &
                    tostr(icase+1) // " to " // tostr(this%population_size) )
      else
        call log_msg(ltag_info // "The extra data are left intact" )
      end if
 
    end if fill_extra
 
    call csv_close( handle_csv )
    if (present(is_success)) is_success = .true.
    call log_msg( ltag_info // "Read genome data from " //                    &
                  trim(handle_csv%name) // " completed, file closed. " //     &
                  tostr(n_ignored_lns) // " rows ignored." )
    call log_delimiter(log_level_chapter)
 
    contains  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      !> This subroutine writes error message to the main log file.
      subroutine log_write_error(message)
        character(len=*), intent(in) :: message
 
        call log_msg( ltag_error // message // ": " //                        &
                      csv_file_name // ", in " // procname )
        call log_msg( ltag_error // "Data file " // csv_file_name //          &
                      " cannot be read in " // procname )
        if (present(is_success)) is_success = .false.
      end subroutine log_write_error
 
  end subroutine population_load_data_all_genomes
 
  !-----------------------------------------------------------------------------
  !> Save the perceptual and emotional memory stack data of all agents in this
  !! population to a CSV file.
  subroutine population_save_data_memory(this, csv_file_name, is_success)
    use file_io
    use csv_io
    class(population), intent(in) :: this
    !> @param[in] csv_file_name is the name of the CSV output file.
    character(len=*), intent(in) :: csv_file_name
    !> @param[out] is_success Flag showing that data save was successful
    !!             (if TRUE).
    logical, optional, intent(out) :: is_success
 
    ! Local counters.
    integer :: i, j, agent
 
    ! File handle object.
    type(file_handle) :: memory_file
 
    ! *Record* that keeps the whole row of data. See @ref file_io.
    character(len=:), allocatable :: record_csv
 
    ! The length of a single record of data.
    integer :: record_csv_max_length
 
    !> ### Implementation details ###
    !> #### Notable variables ####
    !> - **COLUMNS_PERC** defines the column names for all components of the
    !!   perceptual memory stack. They must agree with the components of
    !!   perceptual memory: the_neurobio::memory_perceptual
    character(len=LABEL_LENGTH), dimension(*), parameter ::                   &
                  COLUMNS_PERC =  [ character(len=label_length) ::            &
                                      ! Perception memory:
                                          "PERC_LIGHT",               & !  1
                                          "PERC_DEPTH",               & !  2
                                          "PERC_FOOD",                & !  3
                                          "PERC_FOODSIZ",             & !  4
                                          "PERC_FOODIST",             & !  5
                                          "PERC_CONSP",               & !  6
                                          "PERC_PRED",                & !  7
                                          "PERC_STOM",                & !  8
                                          "PERC_BDMASS",              & !  9
                                          "PERC_ENERG",               & ! 10
                                          "PERC_REPRFAC"  ]             ! 11
 
    !> - **COLUMNS_EMOT** defines the column name for all components of the
    !!   emotional memory stack. They must agree with the components of
    !!   emotional memory: the_neurobio::memory_emotional.
    ! .
    character(len=LABEL_LENGTH), dimension(*), parameter ::                   &
                  COLUMNS_EMOT = [ character(len=label_length) ::             &
                                      ! Emotional memory:
                                          "MOTIV_HUNGER",             & ! 1
                                          "MOTIV_AVOIDACT",           & ! 2
                                          "MOTIV_REPROD",             & ! 3
                                      ! GOS memory:
                                          "GOS_MAIN",                 & ! 4
                                          "GOS_AROUSAL",              & ! 5
                                          "GOS_REPEATED"  ]             ! 6
 
    !> #### Preliminary ####
    !> First, the file `csv_file_name` is opened for writing.
    call memory_file%open_write( csv_file_name, format_csv )
    if ( .not. memory_file%is_success() ) then
      if (present(is_success)) is_success = .false.
      call memory_file%close()
      return
    end if
 
    !> #### Build column names ####
    !> The maximum length of the data record is calculated from three components
    !! - individual IDs: numeric ID and label;
    !! - perceptual memory components from `COLUMNS_PERC` array for
    !!   commondata::history_size_perception steps;
    !! - emotional memory components from `COLUMNS_EMOT` array for
    !!   commondata::history_size_motivation steps.
    !! .
    !! Note that, because the record also adds separators, such as comma and
    !! possibly double quotes, the number of columns multiplied by three is
    !! added to this value.
    record_csv_max_length = 50 +                                              &
          label_length * 2 + 2 * 6 +                                          &
          (label_length * size(columns_perc) + size(columns_perc) * 6) *      &
                  history_size_perception +                                   &
          (label_length * size(columns_emot) + size(columns_emot) * 6) *      &
                  history_size_motivation
 
    !> The first record of the data that contains the column names is now being
    !! "appended" into the complete record and written to the file.
    !! The maximum length of this record is calculated above.
    record_csv = repeat( " ", record_csv_max_length )
 
    !! The first two columns contain the identifiers for each agent:
    !! - numeric ID of the agent
    !! - test string label ("name") of the agent
    !! .
    call csv_record_append( record_csv,                                       &
                    [ character(len=LABEL_LENGTH) :: "ID_NUM", "AGENT_NAME"] )
 
    !> The next portion is composed of the perceptual memory columns from
    !! `COLUMNS_PERC` array for each of the commondata::history_size_perception
    !! steps in the memory.
    !! @note Note that the order of data is: (1) perception component,
    !!       (2) history steps:
    !!       `PERC_LIGHT01, PERC_LIGHT02, PERC_LIGHT03, ...
    !!        PERC_DEPTH01, PERC_DEPTH02, PERC_DEPTH03, ...`
    do j = 1, size(columns_perc)
      do i = 1, history_size_perception
        call csv_record_append( record_csv,                                   &
                                trim(columns_perc(j)) //                      &
                                    tostr(i, history_size_perception) )
      end do
    end do
 
    !> The third portion is composed of the emotional memory columns from
    !! `COLUMNS_EMOT` array for each of the commondata::history_size_motivation
    !! steps in the memory.
    !! @note Note that the order of data is: (1) motivation component,
    !!       (2) history steps:
    !!       `MOTIV_HUNGER01, MOTIV_HUNGER02, MOTIV_HUNGER03, ...
    !!        MOTIV_AVOIDPAS01, MOTIV_AVOIDPAS02, MOTIV_AVOIDPAS03, ...`.
    do j = 1, size(columns_emot)
      do i = 1, history_size_motivation
        call csv_record_append( record_csv,                                   &
                                trim(columns_emot(j)) //                      &
                                    tostr(i, history_size_motivation) )
      end do
    end do
 
    !> After this step, the first record of column names is ready to be
    !! written to the file.
    call memory_file%record_write( record_csv )
    if ( .not. memory_file%is_success() ) then
      if (present(is_success)) is_success = .false.
      call memory_file%close()
      return
    end if
 
    !> #### Write the numerical data ####
    !> The actual data are written in the same order as above, looping over
    !! all individual agents in this population.
    !!
    !! The maximum record length `record_csv_max_length` is here the same as
    !! for writing the column headers, it is assumed that any numeric value
    !! in the data matrix occupies less than commondata::label_length
    !! characters.
    !!
    !! So, for each agent, the following data are written with full history:
    inds: do agent=1, this%population_size
      record_csv = repeat(" ", record_csv_max_length )
      associate( agent => this%individual(agent) )
        !> - ID data: numeric ID and the string label;
        call csv_record_append(record_csv,agent%person_number)
        call csv_record_append(record_csv,agent%genome_label )
        !> - Perceptual memory components;
        call csv_record_append(record_csv,agent%memory_stack%memory_light  ) !1
        call csv_record_append(record_csv,agent%memory_stack%memory_depth  ) !2
        call csv_record_append(record_csv,agent%memory_stack%memory_food   ) !3
        call csv_record_append(record_csv,agent%memory_stack%memory_foodsiz) !4
        call csv_record_append(record_csv,agent%memory_stack%memory_foodist) !5
        call csv_record_append(record_csv,agent%memory_stack%memory_consp  ) !6
        call csv_record_append(record_csv,agent%memory_stack%memory_pred   ) !7
        call csv_record_append(record_csv,agent%memory_stack%memory_stom   ) !8
        call csv_record_append(record_csv,agent%memory_stack%memory_bdmass ) !9
        call csv_record_append(record_csv,agent%memory_stack%memory_energ  ) !10
        call csv_record_append(record_csv,agent%memory_stack%memory_reprfac) !11
        !> - Emotional memory components;
        call csv_record_append(record_csv,                                 & !1
                                  agent%memory_motivations%hunger)
        call csv_record_append(record_csv,                                 & !2
                                  agent%memory_motivations%defence_fear)
        call csv_record_append(record_csv,                                 & !3
                                  agent%memory_motivations%reproduction)
        !> - GOS memory components, note here that `gos_main` is a text value
        !!   that is undefined (empty string) at the initialisation stage;
        !! .
        call csv_record_append(record_csv,                                 & !4
                                  agent%memory_motivations%gos_main)
        call csv_record_append(record_csv,                                 & !5
                                  agent%memory_motivations%gos_arousal)
        call csv_record_append(record_csv,                                 & !6
                                  agent%memory_motivations%gos_repeated)
      end associate
 
      !> Each complete record is written to the file as it is built.
      call memory_file%record_write( record_csv )
      if ( .not. memory_file%is_success() ) then
        if (present(is_success)) is_success = .false.
        call memory_file%close()
        return
      end if
 
    end do inds
 
    !> Once all individuals are saved, the file is closed.
    call memory_file%close()
    if ( .not. memory_file%is_success() ) then
      if (present(is_success)) is_success = .false.
    else
      if (present(is_success)) is_success = .true.
    end if
 
    !> The CSV output data file can be optionally compressed with the
    !! commondata::cmd_zip_output command if commondata::is_zip_outputs is set
    !! to TRUE.
    if ( is_zip_outputs ) then
      call call_external(command=cmd_zip_output // " " // csv_file_name,      &
                         suppress_output=.true.,                              &
                         is_background_task=zip_outputs_background )
    end if
 
  end subroutine population_save_data_memory
 
  !-----------------------------------------------------------------------------
  !> Save the latest movement history of all agents.
  !! This method makes use of the the_environment::spatial_moving::history
  !! structure that saves latest movements of each agent.
  subroutine population_save_data_movements(this, csv_file_name, is_success)
    use file_io
    use csv_io
    class(population), intent(in) :: this
    !> @param[in] csv_file_name is the name of the CSV output file.
    character(len=*), intent(in) :: csv_file_name
    !> @param[out] is_success Flag showing that data save was successful
    !!             (if TRUE).
    logical, optional, intent(out) :: is_success
 
    ! Local counters.
    integer :: i, agent
 
    ! File handle object.
    type(file_handle) :: history_file
 
    ! *Record* that keeps the whole row of data. See @ref file_io.
    character(len=:), allocatable :: record_csv
 
    ! The length of a single record of data.
    integer :: record_csv_max_length
 
    !> ### Implementation notes ###
    !> The maximum length of the data record is calculated from three
    !! components: (1) commondata::history_size_spatial * 3 columns of X, Y,
    !! and depth coordinates plus (2) the same number of separators for these
    !! data columns (assuming 3 characters) plus (3) two additional columns
    !! that contain the numeric ID of the agent and its string label.
    record_csv_max_length =  history_size_spatial * label_length * 3 +        &
                             history_size_spatial * 3 * 3  +                  &
                             label_length * 2 + 2 * 3
 
    !> First, the file `csv_file_name` is opened for writing.
    call history_file%open_write( csv_file_name, format_csv )
    if ( .not. history_file%is_success() ) then
      if (present(is_success)) is_success = .false.
      call history_file%close()
      return
    end if
 
    !> #### Build column names ####
    !> The first record of the data that contains the column names is now being
    !! built. The maximum length of this record is calculated above.
    !! First thing to do is to cleanup the record string.
    record_csv = repeat( " ", record_csv_max_length )
 
    !> After this, the first record is built by appending components.
    !!
    !! The first two columns contain the identifiers for each agent:
    !! - numeric ID of the agent
    !! - test string label ("name") of the agent
    !! .
    call csv_record_append( record_csv,                                       &
                    [ character(len=LABEL_LENGTH) :: "ID_NUM", "AGENT_NAME"] )
 
    !> All remaining columns are built in a loop, by triplets: "X", "Y",
    !! "Depth" for each step of the history, up to
    !! commondata::history_size_spatial triplets.
    do i=1, history_size_spatial
      call csv_record_append( record_csv,                                     &
                              "X_" // tostr(i, history_size_spatial),         &
                              "Y_" // tostr(i, history_size_spatial),         &
                              "D_" // tostr(i, history_size_spatial)  )
    end do
 
    !> Now the first record of column names is ready to be written to the file.
    call history_file%record_write( record_csv )
    if ( .not. history_file%is_success() ) then
      if (present(is_success)) is_success = .false.
      call history_file%close()
      return
    end if
 
    !> #### Write the numerical data ####
    !> The actual data are written in the same order as above, looping over
    !! all individual agents in this population.
    inds: do agent=1, this%population_size
 
      !> - The record string is cleaned;
      record_csv = repeat(" ", record_csv_max_length )
 
      associate( agent => this%individual(agent) )
        !> - ID data appended: numeric ID and the string label;
        call csv_record_append( record_csv, agent%person_number )
        call csv_record_append( record_csv, agent%genome_label  )
        !> - Movement history triplets (x, y, depth) for all
        !!   commondata::history_size_spatial are appended to the
        !!   record string in a loop.
        do i = 1, history_size_spatial
          call csv_record_append( record_csv, agent%history(i)%x ,            &
                                              agent%history(i)%y ,            &
                                              agent%history(i)%depth  )
        end do
      end associate
 
      !> Each complete record is written to the file as it is built.
      call history_file%record_write( record_csv )
      if ( .not. history_file%is_success() ) then
        if (present(is_success)) is_success = .false.
        call history_file%close()
        return
      end if
 
    end do inds
 
    !> Once all individuals are saved, the file is closed.
    call history_file%close()
    if ( .not. history_file%is_success() ) then
      if (present(is_success)) is_success = .false.
    else
      if (present(is_success)) is_success = .true.
    end if
 
    !> The CSV output data file can be optionally compressed with the
    !! commondata::cmd_zip_output command if commondata::is_zip_outputs is set
    !! to TRUE.
    if ( is_zip_outputs ) then
      call call_external(command=cmd_zip_output // " " // csv_file_name,      &
                         suppress_output=.true.,                              &
                         is_background_task=zip_outputs_background )
    end if
 
  end subroutine population_save_data_movements
 
  !-----------------------------------------------------------------------------
  !> Save the behaviours history stack the_neurobio::behaviour::history_behave
  !! for all agents.
  subroutine population_save_data_behaviours(this, csv_file_name, is_success)
    use file_io
    use csv_io
    class(population), intent(in) :: this
    !> @param[in] csv_file_name is the name of the CSV output file.
    character(len=*), intent(in) :: csv_file_name
    !> @param[out] is_success Flag showing that data save was successful
    !!             (if TRUE).
    logical, optional, intent(out) :: is_success
 
    ! Local counters.
    integer :: i, agent
 
    ! File handle object.
    type(file_handle) :: history_file
 
    ! *Record* that keeps the whole row of data. See @ref file_io.
    character(len=:), allocatable :: record_csv
 
    ! The length of a single record of data.
    integer :: record_csv_max_length
 
    !> ### Implementation notes ###
    !> The maximum length of the data record is calculated from three
    !! components: (1) commondata::history_size_behaviours labels of behaviours
    !! plus (2) the same number of separators for these data columns (assuming
    !! 3 characters) plus (3) two additional columns that contain the numeric
    !! ID of the agent and its string label.
    record_csv_max_length =  history_size_behaviours * label_length +        &
                             history_size_behaviours * 3  +                  &
                             label_length * 2 + 2 * 3
 
    !> First, the file `csv_file_name` is opened for writing.
    call history_file%open_write( csv_file_name, format_csv )
    if ( .not. history_file%is_success() ) then
      if (present(is_success)) is_success = .false.
      call history_file%close()
      return
    end if
 
    !> #### Build column names ####
    !> The first record of the data that contains the column names is now being
    !! built. The maximum length of this record is calculated above.
    !! First thing to do is to cleanup the record string.
    record_csv = repeat( " ", record_csv_max_length )
 
    !> After this, the first record is built by appending components.
    !!
    !! The first two columns contain the identifiers for each agent:
    !! - numeric ID of the agent
    !! - test string label ("name") of the agent
    !! .
    call csv_record_append( record_csv,                                       &
                    [ character(len=LABEL_LENGTH) :: "ID_NUM", "AGENT_NAME"] )
 
    !> All remaining columns are built in a loop for each step of the history,
    !! up to commondata::history_size_behaviours steps back in history.
    do i=1, history_size_behaviours
      call csv_record_append( record_csv,                                     &
                              "BEHAV_" // tostr(i, history_size_behaviours) )
    end do
 
    !> Now the first record of column names is ready to be written to the file.
    call history_file%record_write( record_csv )
    if ( .not. history_file%is_success() ) then
      if (present(is_success)) is_success = .false.
      call history_file%close()
      return
    end if
 
    !> #### Write the numerical data ####
    !> The actual data are written in the same order as above, looping over
    !! all individual agents in this population.
    inds: do agent=1, this%population_size
 
      !> - The record string is cleaned;
      record_csv = repeat(" ", record_csv_max_length )
 
      associate( agent => this%individual(agent) )
        !> - ID data appended: numeric ID and the string label;
        call csv_record_append( record_csv, agent%person_number )
        call csv_record_append( record_csv, agent%genome_label  )
        !> - Behaviour history for all commondata::history_size_behaviours
        !!   steps is appended to the record string in a loop.
        do i = 1, history_size_behaviours
          call csv_record_append( record_csv, agent%history_behave(i) )
        end do
      end associate
 
      !> Each complete record is written to the file as it is built.
      call history_file%record_write( record_csv )
      if ( .not. history_file%is_success() ) then
        if (present(is_success)) is_success = .false.
        call history_file%close()
        return
      end if
 
    end do inds
 
    !> Once all individuals are saved, the file is closed.
    call history_file%close()
    if ( .not. history_file%is_success() ) then
      if (present(is_success)) is_success = .false.
    else
      if (present(is_success)) is_success = .true.
    end if
 
    !> The CSV output data file can be optionally compressed with the
    !! commondata::cmd_zip_output command if commondata::is_zip_outputs is set
    !! to TRUE.
    if ( is_zip_outputs ) then
      call call_external(command=cmd_zip_output // " " // csv_file_name,      &
                         suppress_output=.true.,                              &
                         is_background_task=zip_outputs_background )
    end if
 
  end subroutine population_save_data_behaviours
 
  !-----------------------------------------------------------------------------
  !> Calculate fitness for the pre-evolution phase of the genetic algorithm.
  !! **Pre-evolution** is based on selection for a simple criterion without
  !! explicit reproduction etc. The criterion for selection at this phase
  !! is set by the integer the_individual::individual_agent::fitness component.
  !! This procedure provides a whole-population wrapper for the
  !! the_individual::individual_agent::fitness_calc() function.
  !! @warning Note that fitness here is actually an inverse of the fitness: the
  !!          higher its value, the worse fitting is the agent.
  pure subroutine population_preevol_fitness_calc(this)
    class(population), intent(inout) :: this
 
    call this%individual%fitness_calc()
 
  end subroutine population_preevol_fitness_calc
 
  !-----------------------------------------------------------------------------
  !> Determine the number of parents that have fitness higher than the
  !! minimum acceptable value. Also, only alive agents are included into
  !! the reproducing number.
  !! @note This procedure is used in the fixed explicit fitness genetic
  !!       algorithm, see the_evolution::generations_loop_ga().
  pure function population_ga_reproduce_max (this) result (val_out)
    class(population), intent(in) :: this
    integer :: val_out
 
    !> - `MIN_FITNESS` is the normal limit of the fitness value for inclusion
    !!    into the reproducing elite group. See also
    !!    the_individual::individual_agent::individual_preevol_fitness_calc().
    integer, parameter :: min_fitness = ga_fitness_select
 
    !> - `MIN_GA_REPRODUCE` is the minimum GA_REPRODUCE, the final value
    !!    cannot be smaller than. It is set as the minimum proportion
    !!    commondata::ga_reproduce_min_prop of the commondata::popsize.
    !!    However, it cannot be smaller than the absolute minimum
    !!    commondata::ga_reproduce_n_min.
    !! .
    integer, parameter :: min_ga_reproduce =                                  &
                max( ga_reproduce_n_min, nint(ga_reproduce_min_prop*popsize) )
 
    ! - n_alive is the number of agents alive.
    integer :: n_alive
 
    val_out = within( count( this%individual%fitness < min_fitness            &
                               .and. this%individual%fitness >= 0  ),         &
                      min_ga_reproduce, ga_reproduce_n  )
 
    ! Only alive agents are reproducing, so the output val_out is limited
    ! by the number of alive agents.
    n_alive = count(this%individual%alive)
    if ( val_out > n_alive )  val_out = max( min_ga_reproduce, n_alive )
 
  end function population_ga_reproduce_max
 
  !-----------------------------------------------------------------------------
  !> This function implements adaptive mutation rate that *increases*
  !! as the population size *reduces*.
  function population_ga_mutation_rate_adaptive(this, baseline, maxvalue)     &
                                                        result(mutat_rate_out)
    class(population), intent(in) :: this
    !> @param[in] baseline baseline mutation rate
    real(srp), intent(in) :: baseline
    !> @param[in] maxvalue maximum mutation rate
    real(srp), optional, intent(in) :: maxvalue
    !> @return The adjusted adaptive mutation rate.
    real(srp) :: mutat_rate_out
 
    !> ### Implementation notes ###
    !> #### Notable variables and parameters ####
    !> - `mutationrate_max` -- the maximum limit to the mutation rate. Can be
    !!    obtained from the optional parameter `maxvalue`. The function
    !!    returns its value at the lowest population size.
    !!     @remark It is actually local copy of optional `maxvalue` parameter.
    real(srp) :: mutationrate_max
 
    !> - `MUTATIONRATE_MAX_DEF` -- the default maximum limit to the mutation
    !!    rate `mutationrate_max` if `maxvalue` is not provided.
    real(srp), parameter :: mutationrate_max_def = 0.4_srp
 
    real(srp) :: step
 
    !> - `n_base_point` -- this is the base value for calculation of the
    !!    adaptive mutation rate. It can be the number of agents alive or
    !!    the number of agents that have grown.
    integer :: n_base_point
 
    !> - `mutation_grid_abscissa` and `mutation_grid_ordinate` are the arrays
    !!    that define the interpolation grid for the adaptive mutation rate.
    real(srp), dimension(3) :: mutation_grid_abscissa, mutation_grid_ordinate
 
    !> - `MIN_GROWING` a parameter setting the minimum number of growing
    !!    agents in the population. If their actual number is smaller,
    !!    the mutation rate further incremented by a factor set by the
    !!    parameter
    integer, parameter :: min_growing = 4
 
    !> - `NON_GROW_INCREMENT` is an increment factor for the mutation rate
    !!    in case the number of growing agents is below the lower limit
    !!    `MIN_GROWING`.
    !! .
    real(srp), parameter :: non_grow_increment = 1.3_srp
 
    if (present(maxvalue)) then
      mutationrate_max = maxvalue
    else
      mutationrate_max = mutationrate_max_def
    end if
 
    !> #### Procedure ####
    !> First, calculate the base point `n_base_point`. This value is the
    !! base for the calculation of the adaptive mutation rate.
    !! - If this number reduces, mutation rate increases up to
    !!  `mutationrate_max`
    !! .
    !! The `n_base_point` can be:
    !! - Number of agents that are **alive**:
    !!       n_base_point = count( this%individual%is_alive() )
    !! - Number of agents that have **grown**:
    !!       n_base_point = count( this%individual%get_mass() >               &
    !!                             this%individual%get_mass_birth() )
    !! .
    n_base_point = count( this%individual%is_alive() )
    !n_base_point = count( this%individual%get_mass() >                       &
    !                      this%individual%get_mass_birth() )
 
    step = ( mutationrate_max - baseline ) / 4.0_srp
 
    !> - If the mutation rate is based on the number of alive agents, the grid
    !!   abscissa `mutation_grid_abscissa` is an array with three
    !!   elements:
    !!   - minimum full population size
    !!   - 1/2 of the full population size commondata::popsize
    !!   - full commondata::popsize
    !!   .
    mutation_grid_abscissa = [ 0.0_srp, popsize/2.0_srp, real(popsize, srp) ]
 
    !> - If, on the other hand, adaptive mutation rate is based on the number of
    !!   agents that have grown, abscissa is set to specific arbitrary numbers
    !!   that seem more or less optimal for the model performance.
    !! .
    !mutation_grid_abscissa = [ 0.0_SRP, 50.0_SRP, 100.0_SRP ]
 
    !> The grid ordinate is defined as
    !! - the maximum mutation rate limit `mutationrate_max`
    !! - a small middle value that is calculated as 1/4 of the range between
    !!   the maximum (`MUTATIOlNRATE_MAX`) and the minimum (`baseline`)
    !!   mutation values; the latter value increments the minimum `baseline`.
    !! - the lowest baseline value that is set by the `baseline` dummy
    !!   parameter.
    !! .
    mutation_grid_ordinate = [ mutationrate_max, baseline + step, baseline  ]
 
    !> Adaptive mutation rate is calculated based on the DDPINTERPOL()
    !! procedure with the grid array set by `mutation_grid_abscissa` and
    !! `mutation_grid_ordinate` and the interpolation value set by the
    !! number of agents that are the_genome::individual_genome::is_alive()
    !! in the population.
    !! An example pattern of the adaptive mutation rate function is plotted
    !! below. Here @f$ P_{max} @f$ is the maximum mutation rate defined by
    !! `mutationrate_max` and @f$ P_{b} @f$ is the baseline (normal, low)
    !! mutation rate defined by the `baseline` parameter.
    !! @image html img_doxygen_adapt_mutation.svg  "Adaptive mutation rate"
    !! @image latex img_doxygen_adapt_mutation.eps "Adaptive mutation rate" width=14cm
    ! plotting commands:
    ! htintrpl.exe [0 500 1000] [0.5 0.275 0.2]
    ! htintrpl.exe [0 500 1000] [0.5   0.2 0.1]
    mutat_rate_out = within(  ddpinterpol(                                    &
                                        mutation_grid_abscissa,               &
                                        mutation_grid_ordinate,               &
                                        real(n_base_point, srp) ),            &
                              baseline,                                       &
                              mutationrate_max )
 
    !> If the number of agents that had grown from their birth mass is less
    !! than the minimum number (`MIN_GROWING`), mutation rate is incremented
    !! by a fixed factor `NON_GROW_INCREMENT`. However, it is still forced to
    !! be within the range `[ baseline, mutationrate_max ]`.
    if ( count( this%individual%get_mass() >                                  &
                this%individual%get_mass_birth() ) < min_growing ) then
      mutat_rate_out = within( mutat_rate_out * non_grow_increment,           &
                               baseline, mutationrate_max  )
    end if
 
    !> Interpolation plots can be saved in the @ref intro_debug_mode
    !! "debug mode" using this plotting command:
    !! commondata::debug_interpolate_plot_save().
    !! @warning Involves **huge** number of plots, should normally be
    !!          disabled.
    call debug_interpolate_plot_save(                                         &
            grid_xx=mutation_grid_abscissa, grid_yy=mutation_grid_ordinate,   &
            ipol_value=real(n_base_point, srp), algstr="DDPINTERPOL",         &
            output_file="plot_debug_adaptive_mutation_rate_" // mmdd // "_g_" &
                        // tostr(global_generation_number_current) // ps )
 
  end function population_ga_mutation_rate_adaptive
 
end module the_population