00b SimEvo V5.08

;;=====================================================================|=======|
;; SECTION A – AUTHOR IDENTIFICATION AND CODE ABSTRACT
;;=====================================================================|=======|
;;
;; File Name: SimEvo_V5.08.nlogo
;; By Orrery Software
;; Dated: 2018-11-21
;; Author contact:
;;   Garvin H Boyle
;;   orrery@rogers.com
;;   http://orrery-software.webs.com

;;---------------------------------------------------------------------|-------|
;; This simEvo ABM is based on "Simulated Evolution" by Dr Michael Palmiter
;;   as described in the Scientific American article by A.K.Dewdney of 1989.

;; In this model, energy arrives in a steady stream from the Sun and is
;;   captured in algae, which appears randomly in patches in the pond. 
;;   Blind seekers roam the pond searching for food, but they cannot sense
;;   it.  Each must develop a heuristic strategy, an informed pattern of 
;;   behaviour, to seek the algae mats currently having food, and harvest it.
;;   The strategies start as ineffective bland heuristics, and evolve to be 
;;   much more sophisticated.  Only those with the best search pattern
;;   survive to reproduce.
;;---------------------------------------------------------------------|-------|
;; This program was developed on NetLogo Version 5.0.5
;;

;;=====================================================================|=======|
;; SECTION B – INITIAL DECLARATIONS OF GLOBAL VARIABLES AND BREEDS
;;=====================================================================|=======|
;;

;;---------------------------------------------------------------------|-------|
;; Code-determined global variables
globals
[
;;---------------------------------------------------------------------|-------|
;; MODELING ENGINE
;;---------------------------------------------------------------------|-------|
;; BIOPHYSICAL SUB-SYSTEM PARAMETERS
;;---------------------------------------------------------------------|-------|

  ;; Biophyscial life function parameters, seekers, 
  ;;   derived from Dr. Michael Palmiter's model.
  g-c2-dat-parm  ;; The death age threshold
  g-c2-det-parm  ;; The death energy threshold
  g-c2-rat-parm  ;; The reproductive age threshold
  g-c2-ret-parm  ;; The reproductive energy threshold
  g-c2-epm-parm  ;; The energy per move
  g-c2-epa-parm  ;; The maximum energy an agent may hold
  
  ;; Various parameters - soft coded
  g-energy-per-deposit ;; Energy placed in a patch
  g-prob-of-mutation   ;; Probability that a seeker mutates at birth.

  ;; The global list of possible heading deltas for moves.
  gl-heading-list  ;; List of heading deltas.
  gl-base-factors  ;; List of factors, used to mutate bases
  
  ;; Energy control variables
  g-prob-of-deposit        ;; Probability that energy will be deposited. 
  g-energy-in-sunshine     ;; Daily influx of energy held in the sunshine.

  ;; List of counts of deaths, by cause.
  gl-cod                   ;; Counts of deaths, by cause.
  
  ;; Global enumeration (ge-) codes for cause of death.
  ge-cod-none
  ge-cod-fission
  ge-cod-hunger
  ge-cod-oldage

  ;; For display.
  g-energy-entropic-index  ;; For [energy] of seekers.
  gl-ave-phenotype         ;; Average values of c1-pheno
  g-pheno-entropic-index   ;; For average phenotypes
]

;;---------------------------------------------------------------------|-------|
;; DEFINING PATCHES AND BREEDS
;;---------------------------------------------------------------------|-------|
  
;;---------------------------------------------------------------------|-------|
;; Attributes of patches
patches-own 
[
  algae  ;; The amount of energy in algae in this patch.
]

;;---------------------------------------------------------------------|-------|
;; Turtles and breeds
breed [ seekers seeker ]

;;---------------------------------------------------------------------|-------|
;; Attributes of seekers
seekers-own 
[
  ;; Chromosome 1 (C1) genes are used to distinguish behaviours.
  c1-bases ;; list of 8 [B]ases for genes. 
  c1-expos ;; list of 8 [E]xponents for the genes.
  c1-stren ;; list of 8 [S]trengths             [ Si=Bi^Ei ]
  c1-pheno ;; list of 8 [P]henotypic characteristics [Pi=(Si/Sum(Si))]
  
  ;; Chromosome 2 (C2) genes are static in this model.
  c2-DAT   ;; Death Age Threshold.
  c2-DET   ;; Death Energy Threshold.
  c2-RAT   ;; Reproductive Age Threshold.
  c2-RET   ;; Reproductive Energy Threshold.
  c2-EPM   ;; Energy Per Move.
  c2-EPA   ;; Maximum Energy Per Agent.
  
  ;; Chromosome 3 (C3) genes are not present in this model.
  
  ;; Other variable characteritics.
  mas-who                  ;; Serial number of mother agent.
  age                      ;; Age of the seeker in ticks
  energy                   ;; Energy in this seeker
  cause-of-death           ;; A switch
  b-is-ready-to-move?      ;; 0 = no; 1 = ready to move
  b-is-ready-to-reproduce? ;; Mature (in age) and healthy (in energy)
  b-is-ready-to-die?       ;; Old (in age) or starved (in energy)
]

;;=====================================================================|=======|
;; SECTION C – INITIALIZATION OR SETUP PROCEDURE(S)
;;=====================================================================|=======|

;;---------------------------------------------------------------------|-------|
;; The 'autostart' startup routine

to startup
  ;; This routine is to be executed by the observer.
  ;; It is automatically run when the model is first loaded.
  ;; This ensures that key startup variables stored in persistent controls
  ;;   such as sliders are returned to their default values after a save.

  ;; Any globals not initialized in setup should be initialized here.
  set g-random-seed 7

  ;; Run the setup routine to initialize (other) globals.
  setup
end 

;;---------------------------------------------------------------------|-------|
;; The setup button(s)

to setup
  ;; This routine is to be executed by the observer.

  ;; NOTE: The contents of switches, sliders, and choosers seem to be 
  ;;   immune to these 'clear' commands.
  clear-ticks       ;; Clears the tick counter, must be followed by reset-ticks.
  clear-turtles     ;; Kills all turtles and resets the who counter to zero.
  clear-links       ;; Kills all links.
  clear-patches     ;; Sets all attributes of patches to default values. 
  clear-drawing     ;; Clears all lines and stamps drawn by turtles.
  clear-all-plots   ;; Clears all plots and sets attributes to default values.
  clear-output      ;; Clears the output area, if there is one.
  
  random-seed ( g-random-seed )
  
  ;; Set the energy (encoded in the variable algae) in all of the patches to zero.
  ask patches
  [
    set algae 0
    set pcolor brown  ;; Empty patches (no algae) are brown.
  ]
  
  random-seed 7      ;; Tells the PRNG (a Mersene Twister) to use this seed.
  
  ;; Various global parameters - soft-coded
  set g-energy-per-deposit 40   ;; Energy placed in a patch
  set g-prob-of-deposit 0.2
  set g-prob-of-mutation 0.2 ;; Probability that a seeker mutates at birth.
  ;; Chromosome 2 (C2) biophysical controls - borrowed from Palmiter's model.
  ;; In this model these values are static.  C2 does not mutate.
  set g-c2-dat-parm        1600  ;; Death Age Threshold
  set g-c2-det-parm           4  ;; Death Energy Threshold
  set g-c2-rat-parm         800  ;; Reproductive Age Threshold
  set g-c2-ret-parm        1000  ;; Reproductive Energy Threshold
  set g-c2-epm-parm           4  ;; Energy Per Move (Metabolic cost)
  set g-c2-epa-parm        1600  ;; Maximum Energy Per Agent

  ;; Establish the list of allowed headings, each 45 degrees from the last.
  ;;   These are the possible deltas that will be added to the current heading
  ;;   based on which of the 8 genes is expressed during a move.
  set gl-heading-list [ 0 45 90 135 180 225 270 315 ] 

  ;; The factors used to mutate the base values of the genes need to be 
  ;;   calculated.
  let prime-list [ 7 11 13 17 ]
  let factor-list ( map [ 1 + ( 1 / ? ) ] prime-list )
  let inverse-list ( map [ 1 - ( 1 / ? ) ] prime-list )
  set gl-base-factors ( sentence factor-list inverse-list )

  ;; Suppressed debug code
  ;; show ( word "     prime-list: " prime-list )
  ;; show ( word "    factor-list: " factor-list )
  ;; show ( word "   inverse-list: " inverse-list )
  ;; show ( word "gl-base-factors: " gl-base-factors )

  set gl-cod [ 0 0 0 0 ]        ;; Counts of deaths, by cause.

  ;; Global enumeration (ge-) codes for cause of death.
  set ge-cod-none     0
  set ge-cod-fission  1
  set ge-cod-hunger   2
  set ge-cod-oldage   3
  
  set-default-shape seekers "arrow" ;; pulled from shapes library

  set g-energy-entropic-index 0 ;; For [energy] of seekers.
  set gl-ave-phenotype ( n-values 8 [0] ) ;; Average values of c1-pheno
  set g-pheno-entropic-index  0 ;; For average phenotypes
  
  ;; Initialize the seekers.
  f-initialize-seekers
  
  ;; Update the statistics used in the plots.
  ;; show [energy] of seekers
  f-update-aggregates
  
  reset-ticks ;; restarts tick counter, runs setup commands within plots
  
  ;; end of setup
end 

;;---------------------------------------------------------------------|-------|
;; Initialize a population of seekers.

to f-initialize-seekers
  ;; This routine is to be executed by the observer.

  ;; Initialize the seekers.
  create-seekers 100 [ f-initialize-new-seeker ]
  
  ;; Place more energy into patches.
  ask patches
  [
    set algae 0
    set pcolor brown
  ]

  ;; End of f-initialize-seekers
end 

;;---------------------------------------------------------------------|-------|
;; Calculate the strengths and phenotypic values.

to f-find-strens-n-phenos
  ;; This routine is to be executed by a seeker.

  ;; It uses the second example of the map feature.
  ;; Examples of the map feature.
  ;; show (map + [1 2 3] [2 4 6])
  ;; => [3 6 9]
  ;; show (map [?1 + ?2 = ?3] 
  ;;           [1 2 3] 
  ;;           [2 4 6] 
  ;;           [3 5 9])
  ;; => [true false true]
    
  ;; Compute the strength as Si=Bi^Ei.
  set c1-stren ( map [?1 ^ ?2] 
                     c1-bases 
                     c1-expos )
  ;; Compute the phenotypic character as Pi=(Si/sum(Si)).
  set c1-pheno 
    ( map [?1 / ?2] 
          c1-stren 
          ( n-values 8 [sum c1-stren] ) )

  ;; End of f-find-strens-n-phenos
end 

;;---------------------------------------------------------------------|-------|
;; Initialize a single seeker.

to f-initialize-new-seeker
  ;; This routine is to be executed by a seeker.

  set color yellow
  
  ;; Load chromosome 1 with Palmiter genes.
  set c1-bases [ 2 2 2 2 2 2 2 2 ]  ;; 8 unbiased genes
  set c1-expos [ 0 0 0 0 0 0 0 0 ]  ;; 8 unbiased genes
  ;; Calculate the strengths and phenotypic characters.
  f-find-strens-n-phenos
      
  ;; The biophysical body function genes are static in this model.
  ;; I.e. this chromosome does not suffer mutations.
  ;; Load chromosome 2 with the parameters from sliders.
  set c2-DAT g-c2-dat-parm
  set c2-DET g-c2-det-parm
  set c2-RAT g-c2-rat-parm
  set c2-RET g-c2-ret-parm
  set c2-EPM g-c2-epm-parm
  set c2-EPA g-c2-epa-parm

  ;; Associated with seeker dynamics.
  set energy ( g-c2-rat-parm )
  set age ( random ( g-c2-ret-parm / 2 ) )
  set mas-who -1                 ;; serial number of parent seeker.

  ;; Set the logic trigger flags.
  set b-is-ready-to-move?      true     
  set b-is-ready-to-reproduce? false 
  set b-is-ready-to-die?       false 
  
  ;; end f-initialize-new-seeker
end 

;;=====================================================================|=======|
;; SECTION D – GO OR MAIN-LOOP PROCEDURE(S)
;;=====================================================================|=======|

;;---------------------------------------------------------------------|-------|
;; The go button

to go
  ;; This routine is to be executed by the observer.

  let b-should-stop-now false
  if( count turtles <= 0 ) [ set b-should-stop-now true ]
  if( b-should-stop-now = true )
  [
    stop
  ]
  ;; stop commands should only be seen in the go routine.

  ;; Major steps or functions, done once per tick, in order of execution.
  do-pre-tick   ;; Maintenance activities.
  do-energize   ;; Sun deposits energy into patches of algae.
  do-move       ;; Seekers move according to their unique stochastic pattern.
  do-feed       ;; Seekers eat any algae in the patch they occupy, if hungry.
  do-reproduce  ;; Seekers fission, if old enough and healthy enough.
  do-die        ;; Seekers die, if too old, or too sickly.
  do-post-tick  ;; Maintenance activities.

  ;; End of go.
end 

;;---------------------------------------------------------------------|-------|
;; D1 - do-pre-tick procedure( s )
;;---------------------------------------------------------------------|-------|

to do-pre-tick
  ;; This routine is to be executed by the observer.
  
  ;; Advance the ticks by one and update the plots.
  tick
  ask seekers [ set age ( age + 1 ) ]

  ;; End of do-pre-tick.
end 

;;---------------------------------------------------------------------|-------|
;; D2 – do-energize procedure(s)
;;---------------------------------------------------------------------|-------|

to do-energize
  ;; This routine is to be executed by the observer.
  
  ;; The Sun radiates a fixed amount of energy per tick into the system.
  set g-energy-in-sunshine 400
  
  ;; Make a list of the patches that are without algae (without energy)
  let empty-patch-list ( patches with [algae = 0] )

  ;; Use the probability of deposit to determine if algae is added.
  ask empty-patch-list
  [
    let random-number ( random 100000 )
    let threshold ( 100000 * g-prob-of-deposit )
    if ( ( random-number <= threshold ) and 
         ( g-energy-in-sunshine > g-energy-per-deposit ) )
    [
      set algae ( algae + g-energy-per-deposit )
      set pcolor green
      ;; Record in system energy accounts
      set g-energy-in-sunshine ( g-energy-in-sunshine - g-energy-per-deposit )
    ] ;; End if(() and ())
  ] ;; End ask
      
  ;; End of do-energize.
end 

;;---------------------------------------------------------------------|-------|
;; D3 – do-move procedure(s)
;;---------------------------------------------------------------------|-------|

to do-move
  ;; This routine is to be executed by the observer.
  
  ;; The seekers move.
  ask seekers
  [
    if( b-is-ready-to-move? = true )
    [
      f-seeker-moves
    ]  ;; End if( b-is-ready-to-move? = true )
  ]  ;; End ask seekers

  ;; End of do-move.
end 

;;---------------------------------------------------------------------|-------|
;; A seeker moves according to genes and heuristics.

to f-seeker-moves
  ;; This routine is to be executed by a seeker.

  ;; When a seeker moves it expends energy out of its pool of energy.
  ;; Determine if this seeker has sufficient energy to move.
  ifelse ( energy >= c2-EPM )
  [
    ;; Establish a heading.
    f-seeker-sets-heading
    
    ;; Step forward
    forward 1
    
    set energy ( energy - c2-EPM )
    
    set b-is-ready-to-move? true
  ] 
  ;; Else
  [
    ;; The seeker is marked for death, and energy is removed.
    ;; It will die and be removed when do-die is executed.    
    set energy 0 
    set cause-of-death ge-cod-hunger
    set b-is-ready-to-move? false
    set b-is-ready-to-die?  true
  ]
  ;; End else

  ;; End of f-seeker-moves  
end 

;;---------------------------------------------------------------------|-------|
;; A seeker sets a heading using c1-pheno.

to f-seeker-sets-heading
  ;; This routine is to be executed by a seeker.

  let old-heading heading
  let sum-of-phenos ( sum c1-pheno ) ;; Should be 1.0.
  let random-number ( random-float sum-of-phenos )

  let counter 0
  let good-index -1
  let this-pheno 0
  let this-sum 0
  let next-sum 0

  while [ counter < 8 ]  ;; Do not overshoot
  [
    set this-pheno ( item counter c1-pheno )
    set next-sum ( this-pheno + this-sum )
    if ( ( random-number >= this-sum ) and ( random-number < next-sum  ) )
    [
      set good-index counter
    ]
    set this-sum next-sum
    set counter ( counter + 1 )
  ]

  let heading-delta ( item good-index gl-heading-list )
  set heading ( heading + heading-delta )
  
  ;; End of f-seeker-sets-heading  
end 

;;---------------------------------------------------------------------|-------|
;; D4 – do-feed procedure(s)
;;---------------------------------------------------------------------|-------|

to do-feed
  ;; This routine is to be executed by the observer.
  
  ;; Agents feed on algae found in patches.
  ask seekers
  [
    let this-patch patch-here  ;; handle to the patch under the seeker.
    let energy-available ( [algae] of this-patch )
    if ( ( energy < ( c2-EPA - g-energy-per-deposit ) ) and 
         ( energy-available > 0 ) )
    [
      ;; Case of there is food to eat.
      ;; Seeker eats.
      set energy ( energy + energy-available )
      ask this-patch
      [
        set algae 0
        set pcolor brown
      ]  ;; End of ask this-patch
    ]  ;; End of if(() and ())
  ]  ;; End of ask seekers.
  
  ;; End Do-feed procedure.
end 

;;---------------------------------------------------------------------|-------|
;; D5 – do-reproduce procedure(s)
;;---------------------------------------------------------------------|-------|

to do-reproduce
  ;; This routine is to be executed by the observer.
  
  ask seekers
  [
    f-set-seeker-repro-flag        ;; A seeker routine.
    f-reproduce-seeker-by-fission  ;; A seeker routine.
  ]  ;; end of ask seekers.

  ;; End do-reproduce.  
end 

;;---------------------------------------------------------------------|-------|
;; f-set-seeker-repro-flag

to f-set-seeker-repro-flag
  ;; This routine is to be executed by a seeker.
  
  set b-is-ready-to-reproduce? true        ;; i.e. true is the default.
  if( energy < c2-RET )
    [ set b-is-ready-to-reproduce? false ] ;; due to lack of health.
  
  if( age < c2-RAT )
    [ set b-is-ready-to-reproduce? false ] ;; due to lack of maturity.
  
  ;; End f-set-seeker-repro-flag
end 

;;---------------------------------------------------------------------|-------|
;; A seeker reproduces via fission, one mother having two daughters.

to f-reproduce-seeker-by-fission
  ;; This routine is to be executed by a seeker.
  
  if( b-is-ready-to-reproduce? = true )  
  [

    let mother self
    let mothers-who who
    let mothers-patch patch-here

    let first-share-of-energy floor( energy / 2 )
    let second-share-of-energy ( energy - first-share-of-energy )
  
    let daughter-count 0
    ask mothers-patch
    [
      sprout-seekers 2
      [
        set daughter-count ( daughter-count + 1 )
        f-initialize-new-seeker
        set color ( [color] of mother )
        
        ;; Copy the C1 genetic/learned material.
        set c1-bases ( [c1-bases] of mother )
        set c1-expos ( [c1-expos] of mother )
        set c1-stren ( [c1-stren] of mother )
        set c1-pheno ( [c1-pheno] of mother )

        ;; C2 genes are static.
        set c2-DAT ( [c2-DAT] of mother )
        set c2-DET ( [c2-DET] of mother )
        set c2-RAT ( [c2-RAT] of mother )
        set c2-RET ( [c2-RET] of mother )
        set c2-EPM ( [c2-EPM] of mother )
        set c2-EPA ( [c2-EPA] of mother )
        
        ;; Note the mother of this daughter.
        set mas-who ( [who] of mother )
        set age 0

        ifelse ( daughter-count = 1 )
        [ set energy first-share-of-energy ]
        [ set energy second-share-of-energy ]
        
        set cause-of-death           ge-cod-none
        set b-is-ready-to-move?      true
        set b-is-ready-to-reproduce? false
        set b-is-ready-to-die?       false
        
        f-mutate-new-seeker
      ]  ;; End of sprout.
    ]  ;; End of ask mothers-patch
    
    ;; Set the cause of death for the mother.
    set b-is-ready-to-die? true
    set cause-of-death ge-cod-fission
    ;; The mother disappears after fission, leaving two daughters.
    ;; Death actually occurs in the Do-die step.
    
  ]  ;; End of if(ready-to-reproduce)

  ;; End f-reproduce-seeker-by-fission
end 

;;---------------------------------------------------------------------|-------|
;; A new seeker mutates, changing the genetic basis of strategies.

to f-mutate-new-seeker
  ;; This routine is to be executed by a seeker.
  
  ;; Decide if a genetic mutation is to happen.
  let random-number ( random-float 1 )
  let threshold ( g-prob-of-mutation )
  
  ;; show ( word "f-m-n-c random-number: " random-number )
  ;; show ( word "f-m-n-c     threshold: " threshold )
  
  if ( random-number <= threshold )
  [
    ;; Case of mutation to be done.
    ;; Select the gene to be mutated.  These genes control 
    ;;   the search strategy.
    let gene-to-be-mutated ( random 8 )
    
    ;; show ( word "Base Mutated." )
    ;; Mutate the gene base.
    let oldbase ( item gene-to-be-mutated c1-bases )
    ;; Choose a factor for the base.
    let base-factor ( item (random 8) gl-base-factors )
    ;; Mutate it
    let newbase ( oldbase * base-factor )
    set c1-bases ( replace-item gene-to-be-mutated c1-bases newbase )
    
    ;; show ( word "Gene Mutated." )
    ;; Mutate the gene.  The gene value is an integer of either sign.
    let oldgene ( item gene-to-be-mutated c1-expos )
    ;; Decide whether it will increase or decrease in value.
    let delta ( -1 + ( 2 * ( random 2 ) ) )  ;; Either a -1 or a 1.
    let newgene ( oldgene + delta )
    set c1-expos ( replace-item gene-to-be-mutated c1-expos newgene )

    ;; Calculate the strengths and phenotypic characters.
    f-find-strens-n-phenos
  ]
  
  ;; End of f-mutate-new-seeker
end 
  
;;---------------------------------------------------------------------|-------|
;; D6 – do-die procedure(s)
;;---------------------------------------------------------------------|-------|

to do-die
  ;; This routine is to be executed by the observer.
  
  if( ( count seekers ) > 0 )
  [
    ask seekers
    [
      f-set-seeker-death-flag
      f-seeker-dies
    ]
  ]
end 

;;---------------------------------------------------------------------|-------|
;; f-set-seeker-death-flag

to f-set-seeker-death-flag
  ;; This routine is to be executed by a seeker.
  
  ;; If a cause of death has already been noted, it dies.
  ifelse( cause-of-death > ge-cod-none )
  [
    ;; A cause of death has been previously flagged.
    ;; This is either due to hunger (in do-move) or fission (in do-repro).
    ;; In both cases energy has been stripped out already.
    ;; In the cases of c2-DET and c2-DAT, the flag is not yet set, 
    ;;   and the energy remains.
    set b-is-ready-to-die? true 
  ]  ;; End if.
  ;; Else
  [
    ;; No cause of death has been set yet.  Check basic vital signs.
    if( energy <= c2-DET )  ;; Death Energy Threshold.
    [ 
      set b-is-ready-to-die? true 
      set cause-of-death ge-cod-hunger 
      set energy 0
    ]
    
    if( age > c2-DAT )   ;; Death Age Threshold.
    [ 
      set b-is-ready-to-die? true 
      set cause-of-death ge-cod-oldage 
      set energy 0
    ]
  ]  ;; End else.

  ;; End f-set-seeker-death-flag
end 

;;---------------------------------------------------------------------|-------|

to f-seeker-dies
  ;; This routine is to be executed by a seeker.
  
  if( b-is-ready-to-die? = true )
  [
    ;; show ( word "who died: " who )
    set gl-cod 
      ( replace-item cause-of-death gl-cod ( 1 + item cause-of-death gl-cod ) )  
    ;; show ( word "cod: " cause-of-death )
    ;; show gl-cod
    die  ;; The seeker disappears from the system.  
  ]

  ;; End f-seeker-dies
end 

;;---------------------------------------------------------------------|-------|
;; D7 - do-post-tick procedure(s)
;;---------------------------------------------------------------------|-------|

to do-post-tick
  ;; This routine is to be executed by the observer.
   
  f-update-aggregates  ;; Calculates entropic indices.
  display

  ;; End do-post-tick.
end 

;;=====================================================================|=======|
;; SECTION E – DRAWING AND MAINTENANCE PROCEDURE(S)
;;=====================================================================|=======|

to f-update-aggregates
  ;; this routine is to be executed by the observer.
  
  ;; show "----------"
  ;; show "f-update-aggregates"
  
  ;; show ( word "g-energy-entropic-index before: " g-energy-entropic-index )
  ;; Calculate the entropic index of the energy distribution.
  set g-energy-entropic-index 
    ( fr-get-energy-entropic-index ( [energy] of seekers ) )
  ;; show ( word "g-energy-entropic-index  after: " g-energy-entropic-index )
  print ""
  
  ;; show ""
  ;; show ""
  
  ;; show ( word "g-pheno-entropic-index before: " g-pheno-entropic-index )
  ;; Calculate the entropic index of the phenotype probability distribution.
  let l-pheno-partition ( fr-make-pheno-partition )
  set g-pheno-entropic-index ( fr-calc-pheno-entropic-index l-pheno-partition )
  ;; show ( word "g-pheno-entropic-index after:  " g-pheno-entropic-index )
  print ""

  ;; end of update-aggregates.
end 

;;---------------------------------------------------------------------|-------|
;;  Partitions in ABMs.  Entropic index in ABMs.  GammaLn() function in ABMs.
;;---------------------------------------------------------------------|-------|
;;
;; REFERENCES: The diary notes including (or more recent versions):
;; A.  150527 PPR - Definition of Ei R17.PDF
;; B.  190927 NTF Shannon Vs Boltzmann R4.PDF
;; C.  180214 NTF Entropy and Units of Measure R5.PDF
;; D.  180822 NTF FactLn() and GammaLn() R3.PDF
;; E.  190927 NTF Entropy in a Histogram R4.PDF
;;

;; NOTE:  These routines create a partition (i.e. a histogram) independent of 
;;   the capabilities of the plot controls in the interface tab.

;;---------------------------------------------------------------------|-------|

to-report fr-get-energy-entropic-index [ l-data-input ]
  ;; This routine is to be executed by the observer.
  
  ;; show "----------"
  ;; show "fr-get-energy-entropic-index"
  
  let l-energy-partition ( fr-make-energy-partition l-data-input )
  ;; show ( word "fr-nrg-ei: l-energy-partition - " l-energy-partition )

  let nrg-ei ( fr-calc-energy-entropic-index l-energy-partition )
  ;; show ( word "nrg-ei: " nrg-ei )

  report nrg-ei
  
  ;; End of fr-get-energy-entropic-index
end 

;;---------------------------------------------------------------------|-------|

to-report fr-make-energy-partition [ l-data-input ]
  ;; This routine is to be executed by the observer.
  
  ;; show "----------"
  ;; show "fr-make-energy-partition"
  ;; show ( word "l-data-input:" l-data-input )
  
  ;; It creates binned data (a partition) that can be used to plot a histogram.
  ;; A typical call might look like this:
  ;;  set l-energy-partition 
  ;;              ( fr-make-energy-partition ( [energy] in seekers ) )
  ;;  reports [ l-partition ]
  
  ;; Allocate memory.
  let a-no-of-agents length l-data-input
  let k-no-of-bins ( floor ( sqrt a-no-of-agents ) )
  let ave-data ( mean l-data-input )
  let min-data ( ( min l-data-input ) - 0.001 )
  let max-data ( ( max l-data-input ) + 0.001 )
  let bin-delta 0
    ;; show ( word "fmep: ave-data(1): " ave-data )
    ;; show ( word "fmep: min-data(1): " min-data )
    ;; show ( word "fmep: max-data(1): " max-data )
    ;; show ( word "fmep: bin-delta(1): " bin-delta )
    ;; show ""
  
  ifelse ( ( max-data - min-data ) < k-no-of-bins )
  [
    ;; Case of all data-points being the same.
    set min-data ( ave-data - ( k-no-of-bins / 2 ) )
    set max-data ( ave-data + ( k-no-of-bins / 2 ) )
    set bin-delta ( ( max-data - min-data ) / k-no-of-bins )
    ;; show ( word "fmep: min-data(2): " min-data )
    ;; show ( word "fmep: max-data(2): " max-data )
    ;; show ( word "fmep: bin-delta(2): " bin-delta )
  ]
  ;; Else
  [
    set bin-delta ( ( max-data - min-data ) / k-no-of-bins )
  ]
    
  ;; Initialize the output histogram with zeros.
  let l-partition ( n-values k-no-of-bins [0] )
  ;; show ( word "fmep: l-partition: " l-partition )
  
  let counter 0
  let data-point 0
  let bin-no 0
  let agents-in-bin 0
  let limit ( ( length l-data-input ) - 1 )
  ;; show ( word "fmep: limit: " limit )

  while [counter <= limit ] 
  [
    ;; show ( word "fmep: counter(3): " counter )
    set data-point ( item counter l-data-input )
    ;; We need to convert this data point into a bin-number for the partition.
    set bin-no ( ( ( data-point - min-data ) / bin-delta ) )
    ;; show ( word "fmep: bin-no(3): " bin-no )
    set bin-no ( floor ( ( data-point - min-data ) / bin-delta ) )
    ;; show ( word "fmep: bin-no(3): " bin-no )
    ;; bin-no is a whole number which determines which bin the agent is counted in.
    ;; show ( word "fmep: data-point(3): " data-point )
    ;; show ( word "fmep: min-data(3): " min-data )
    ;; show ( word "fmep: bin-delta(3): " bin-delta )
    
    set agents-in-bin ( item bin-no l-partition )
    set agents-in-bin ( agents-in-bin + 1 )
    set l-partition ( replace-item bin-no l-partition agents-in-bin )
    
    set counter ( counter + 1 )
  ]
  
  report l-partition
  
  ;; End fr-make-energy-partition.
end 

;;---------------------------------------------------------------------|-------|

to-report fr-calc-energy-entropic-index [ l-partition ]
  ;; This routine can be called by any agent.
  
  ;; show "----------"
  ;; show "fr-calc-energy-entropic-index"
  
  ;; It uses the formula that was developed in a series of diary notes
  ;;   that investigated the most effective formula for calculating
  ;;   entropy in an agent-based model such as "Model I" of EiLab.

  ;; I hypothesize that it is equally valid in most agent-based models, but 
  ;;   that hypothesis needs to be expored and examined and tested.

  ;; THIS IS AN EXAMPLE OF AN AS-YET UNEXAMINED USAGE.

  ;;
  ;; REFERENCES: The diary notes including (or more recent versions):
  ;; A.  150527 PPR - Definition of Ei R17.PDF
  ;; B.  190927 NTF Shannon Vs Boltzmann R4.PDF
  ;; C.  180214 NTF Entropy and Units of Measure R5.PDF
  ;; D.  180822 NTF FactLn() and GammaLn() R3.PDF
  ;; E.  190927 NTF Entropy in a Histogram R4.PDF
  ;;

  ;; The formula used is equation 30 in the Ref E document.

  ;; It can be written like this:
  ;;   I = (GammaLn(A+1) - Sum(GammaLn(ai+1))) / (A*Ln(K))
  ;;   where:
  ;;          K is the number of bins in a histogram of some conserved quantity.   
  ;;          A is the number of agents in the model.
  ;;          ai is the count of agents in bin i of the histogram.
  
  ;; The expected input is a partition constructed from a list of energy values.
  ;; The sum of all elements of all bins equals the number of agents.
  let a-no-of-agents ( sum l-partition )
  let k-no-of-bins ( length l-partition )
  
  ;; Note that GammaLn( ) is not native to NetLogo.  A version is implemented
  ;;   below.

  ;; show ( word "l-partition: " l-partition )
  
  let entropic-index fr-gammaln ( 1 + a-no-of-agents )
  foreach l-partition
     [ set entropic-index ( entropic-index - fr-gammaln ( 1 + ? ) ) ]

  let alpha ( a-no-of-agents / k-no-of-bins )
  let boltzmann-max ( fr-gammaln ( 1 + a-no-of-agents ) )
  set boltzmann-max 
    ( boltzmann-max - ( k-no-of-bins * fr-gammaln ( alpha + 1 ) ) )

  set entropic-index ( entropic-index / boltzmann-max )
  
  report entropic-index 
  
  ;; End of fr-calc-entropic-index.
end 

;;---------------------------------------------------------------------|-------|

to-report fr-make-pheno-partition 
  ;; This routine is to be executed by the observer.
  
  ;; show "----------"
  ;; show "fr-make-pheno-partition"
  
  ;; It creates binned data (a partition) that can be used to plot a histogram.
  ;; A typical call might look like this:
  ;;  set l-pheno-partition ( fr-make-pheno-partition )
  ;;  reports [ l-partition ]
  
  ;; Allocate memory.
  let a-no-of-agents count turtles
  let k-no-of-bins 8
  let new-value 0
  let l-pheno-list 0
    
  ;; Initialize the output histogram with zeros.
  let l-partition ( n-values k-no-of-bins [0] )
  ;; show ( word "fmpp: l-partition: " l-partition )
  
  ;; One-by-one, calculate the average phenotype value, and insert.
  let counter 0
  let limit ( length l-partition )
  while [ counter < limit ]
  [
    set l-pheno-list ( [ item counter c1-pheno ] of seekers )
    ;; show ( word "fmpp: l-pheno-list: " l-pheno-list )
    set new-value ( mean l-pheno-list )
    ;; show ( word "fmpp: new-value: " new-value )
    set l-partition ( replace-item counter l-partition new-value )
    ;; show ( word "fmpp: l-partition:" l-partition )
    set counter ( counter + 1 )
  ]

  ;; show ( word "fmpp: l-partition: " l-partition )

  report l-partition
  
  ;; End fr-make-pheno-partition.
end 

;;---------------------------------------------------------------------|-------|

to-report fr-calc-pheno-entropic-index [ l-partition ]
  ;; This routine can be called by any agent.
  
  ;; show "----------"
  ;; show "fr-calc-pheno-entropic-index"
  
  ;; It uses the formula that was developed in a series of diary notes
  ;;   that investigated the most effective formula for calculating
  ;;   entropy in an agent-based model such as "Model I" of EiLab.

  ;; I hypothesize that it is equally valid in most agent-based models, but 
  ;;   that hypothesis needs to be expored and examined and tested.

  ;; THIS IS AN EXAMPLE OF AN AS-YET UNEXAMINED USAGE.

  ;;
  ;; REFERENCES: The diary notes including (or more recent versions):
  ;; A.  150527 PPR - Definition of Ei R17.PDF
  ;; B.  190927 NTF Shannon Vs Boltzmann R4.PDF
  ;; C.  180214 NTF Entropy and Units of Measure R5.PDF
  ;; D.  180822 NTF FactLn() and GammaLn() R3.PDF
  ;; E.  190927 NTF Entropy in a Histogram R4.PDF
  ;;

  ;; The formula used is equation 30 in the Ref E document.

  ;; It can be written like this:
  ;;   I = (GammaLn(A+1) - Sum(GammaLn(ai+1))) / (A*Ln(K))
  ;;   where:
  ;;          K is the number of bins in a histogram of some conserved quantity.   
  ;;          A is the number of agents in the model.
  ;;          ai is the count of agents in bin i of the histogram.
  
  ;; The expected input is a partition constructed from a list of energy values.
  ;; The sum of all elements of all bins equals the number of agents.
  let a-total ( sum l-partition )
  let k-no-of-bins ( length l-partition )
  
  ;; Note that GammaLn( ) is not native to NetLogo.  A version is implemented
  ;;   below.

  ;; show ( word "frcpei: l-partition: " l-partition )
  
  let entropic-index fr-gammaln ( 1 + a-total )
  foreach l-partition
     [ set entropic-index ( entropic-index - fr-gammaln ( 1 + ? ) ) ]

  let alpha ( a-total / k-no-of-bins )
  let boltzmann-max ( fr-gammaln ( 1 + a-total ) )
  set boltzmann-max 
    ( boltzmann-max - ( k-no-of-bins * fr-gammaln ( alpha + 1 ) ) )

  set entropic-index ( entropic-index / boltzmann-max )
  
  report entropic-index 
  
  ;; End of fr-calc-pheno-entropic-index.
end 

;;---------------------------------------------------------------------|-------|

to-report fr-gammaln [ input-value ]
  ;; This routine can be called by any agent.
  
  ;; It reproduces a routine found on page 207 of this book:
  ;; Press, Teukolsky, Vetterling and Flanery (1986) "Numerical Recipes in 
  ;;   Fortran", Cambridge, pp 206-209.
  ;; It uses Lanczos' approximation of the Ln(Gamma(x)) function.
  ;; Note that, when input-value is an integer, GammaLn(x+1) approximates Ln(x!).
  ;; When x = 160, the error seems to be about 1 in a billion.  That's pretty good.

  ;; This follows the Fortran code rather closely here.
  let x input-value
  let y x
  let tmp ( x + 5.5 )
  ;; Note that here I use ln( ) instead of log( ) as shown in the text.
  ;; The same problem comes up down below also.  With ln() it seems
  ;;   to work correctly.
  set tmp ( ( ( x + 0.5 ) * ( ln tmp ) ) - tmp )
  let ser 1.000000000190015

  ;; Next, I have unscrolled the six iterations of the loop.
  ;; cof(6) is now six different parameters.

  ;; Iteration 1
  set y ( y + 1 )
  set ser ( ser + ( 76.18009172947146 / y ) )

  ;; Iteration 2
  set y ( y + 1 )
  set ser ( ser + ( -86.50532032941677 / y ) )

  ;; Iteration 3
  set y ( y + 1 )
  set ser ( ser + ( 24.01409824083091 / y ) )

  ;; Iteration 4
  set y ( y + 1 )
  set ser ( ser + ( -1.231739572450155 / y ) )

  ;; Iteration 5
  set y ( y + 1 )
  set ser ( ser + ( 0.001208650973866179 / y ) )

  ;; Iteration 6
  set y ( y + 1 )
  set ser ( ser + ( -0.000005395239384953 / y ) )

  ;; Now, I follow the Fortran code closely again.
  let stp 2.5066282746310005

  ;; Note that here I use ln( ) instead of log( ) as shown in the text.
  let gammaln ( tmp + ( ln ( stp * ser / x ) ) )

  report precision gammaln 15

  ;; End of fr-gammaln.
end 

;;---------------------------------------------------------------------|-------|

to-report fr-factorial [ input-value ]
  ;; This routine can be called by any agent.
  
  ;; It uses recursion to calculate n!.
  let n input-value
  ifelse ( n > 0 ) [ report n * fr-factorial ( n - 1 ) ]
                   [ report 1 ]

  ;; End of fr-factorial.
end 

;;---------------------------------------------------------------------|-------|

to-report fr-factorialln [ input-value ]
  ;; This routine can be called by any agent.
  
  ;; It uses the recursive function to calculate n!, then takes the logarithm.
  let n input-value
  let factorialln fr-factorial n
  
  report precision ( ln factorialln ) 15
  
  ;; End of fr-factorialln.
end