natas_hw
/
jidt
mirror of https://github.com/jlizier/jidt
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Adding a code suite to run active information storage and transfer entropy analysis of positional / heading data from swarms/flocks/schools, including a demonstration of how to use the suite with a NetLogo example. Documentation on a wiki page to follow later ...

This commit is contained in:
Joseph Lizier 2019-08-29 23:08:16 +10:00
parent c22e37d1ab
commit 4966fdabc2
11 changed files with 1969 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

740
demos/octave/FlockingAnalysis/NetLogoExample/Flocking.nlogo Executable file
View File

@ -0,0 +1,740 @@
;; Adapted by J.T. Lizier from the original CC BY NC SA 3.0 licensed Flocking model from NetLogo
;; and released under the same license.
;; See license info on the Info tab (and below if viewing text file)
globals [
headings-file
positionsx-file
positionsy-file
sorted-agents
]
turtles-own [
flockmates ;; agentset of nearby turtles
nearest-neighbor ;; closest one of our flockmates
]
to setup
clear-all
create-turtles population
[ set color yellow - 2 + random 7 ;; random shades look nice
set size 1.5 ;; easier to see
setxy random-xcor random-ycor
set flockmates no-turtles ]
reset-ticks
;; Set up for file output of data:
set headings-file "headings.txt"
set positionsx-file "positionsx.txt"
set positionsy-file "positionsy.txt"
if file-exists? headings-file [
file-delete headings-file ;; clear file
]
if file-exists? positionsx-file [
file-delete positionsx-file ;; clear file
]
if file-exists? positionsy-file [
file-delete positionsy-file ;; clear file
]
set sorted-agents sort turtles
end
to go
ask turtles [ flock ]
;; the following line is used to make the turtles
;; animate more smoothly.
repeat 5 [ ask turtles [ fd 0.2 ] display ]
;; for greater efficiency, at the expense of smooth
;; animation, substitute the following line instead:
;; ask turtles [ fd 1 ]
capture-data ;; write the current position+heading data of turtles to file
tick
end
to flock ;; turtle procedure
find-flockmates
if any? flockmates
[ find-nearest-neighbor
ifelse distance nearest-neighbor < minimum-separation
[ separate ]
[ align
cohere ] ]
end
to find-flockmates ;; turtle procedure
set flockmates other turtles in-radius vision
end
to find-nearest-neighbor ;; turtle procedure
set nearest-neighbor min-one-of flockmates [distance myself]
end
;;; SEPARATE
to separate ;; turtle procedure
turn-away ([heading] of nearest-neighbor) max-separate-turn
end
;;; ALIGN
to align ;; turtle procedure
turn-towards average-flockmate-heading max-align-turn
end
to-report average-flockmate-heading ;; turtle procedure
;; We can't just average the heading variables here.
;; For example, the average of 1 and 359 should be 0,
;; not 180. So we have to use trigonometry.
let x-component sum [dx] of flockmates
let y-component sum [dy] of flockmates
ifelse x-component = 0 and y-component = 0
[ report heading ]
[ report atan x-component y-component ]
end
;;; COHERE
to cohere ;; turtle procedure
turn-towards average-heading-towards-flockmates max-cohere-turn
end
to-report average-heading-towards-flockmates ;; turtle procedure
;; "towards myself" gives us the heading from the other turtle
;; to me, but we want the heading from me to the other turtle,
;; so we add 180
let x-component mean [sin (towards myself + 180)] of flockmates
let y-component mean [cos (towards myself + 180)] of flockmates
ifelse x-component = 0 and y-component = 0
[ report heading ]
[ report atan x-component y-component ]
end
;;; HELPER PROCEDURES
to turn-towards [new-heading max-turn] ;; turtle procedure
turn-at-most (subtract-headings new-heading heading) max-turn
end
to turn-away [new-heading max-turn] ;; turtle procedure
turn-at-most (subtract-headings heading new-heading) max-turn
end
;; turn right by "turn" degrees (or left if "turn" is negative),
;; but never turn more than "max-turn" degrees
to turn-at-most [turn max-turn] ;; turtle procedure
ifelse abs turn > max-turn
[ ifelse turn > 0
[ rt max-turn ]
[ lt max-turn ] ]
[ rt turn ]
end
to capture-data
;; To do a one-off:
;; file-open "headings.txt"
;; ;; Does not select turtles in order:
;; ask turtles [ file-print heading ]
;; file-close
if ticks > 4000
[ stop ]
;; Write the current headings
file-open headings-file ;; Opening file for writing
foreach sorted-agents [
[the-turtle] ->
ask the-turtle [
file-write heading
]
]
file-print " " ;; Terminate the line
file-close
;; Write the x positions
file-open positionsx-file ;; Opening file for writing
foreach sorted-agents [
[the-turtle] ->
ask the-turtle [
file-write xcor
]
]
file-print " " ;; Terminate the line
file-close
;; Write the y positions
file-open positionsy-file ;; Opening file for writing
foreach sorted-agents [
[the-turtle] ->
ask the-turtle [
file-write ycor
]
]
file-print " " ;; Terminate the line
file-close
end
; Copyright 1998 Uri Wilensky.
; See Info tab for full copyright and license.
@#$#@#$#@
GRAPHICS-WINDOW
250
10
755
516
-1
-1
7.0
1
10
1
1
1
0
1
1
1
-35
35
-35
35
1
1
1
ticks
30.0
BUTTON
39
93
116
126
NIL
setup
NIL
1
T
OBSERVER
NIL
NIL
NIL
NIL
1
BUTTON
122
93
203
126
NIL
go
T
1
T
OBSERVER
NIL
NIL
NIL
NIL
0
SLIDER
9
51
232
84
population
population
1.0
1000.0
300.0
1.0
1
NIL
HORIZONTAL
SLIDER
4
217
237
250
max-align-turn
max-align-turn
0.0
20.0
5.0
0.25
1
degrees
HORIZONTAL
SLIDER
4
251
237
284
max-cohere-turn
max-cohere-turn
0.0
20.0
3.0
0.25
1
degrees
HORIZONTAL
SLIDER
4
285
237
318
max-separate-turn
max-separate-turn
0.0
20.0
1.5
0.25
1
degrees
HORIZONTAL
SLIDER
9
135
232
168
vision
vision
0.0
10.0
3.0
0.5
1
patches
HORIZONTAL
SLIDER
9
169
232
202
minimum-separation
minimum-separation
0.0
5.0
1.0
0.25
1
patches
HORIZONTAL
@#$#@#$#@
## WHAT IS IT?
This model is an attempt to mimic the flocking of birds. (The resulting motion also resembles schools of fish.) The flocks that appear in this model are not created or led in any way by special leader birds. Rather, each bird is following exactly the same set of rules, from which flocks emerge.
## HOW IT WORKS
The birds follow three rules: "alignment", "separation", and "cohesion".
"Alignment" means that a bird tends to turn so that it is moving in the same direction that nearby birds are moving.
"Separation" means that a bird will turn to avoid another bird which gets too close.
"Cohesion" means that a bird will move towards other nearby birds (unless another bird is too close).
When two birds are too close, the "separation" rule overrides the other two, which are deactivated until the minimum separation is achieved.
The three rules affect only the bird's heading. Each bird always moves forward at the same constant speed.
## HOW TO USE IT
First, determine the number of birds you want in the simulation and set the POPULATION slider to that value. Press SETUP to create the birds, and press GO to have them start flying around.
The default settings for the sliders will produce reasonably good flocking behavior. However, you can play with them to get variations:
Three TURN-ANGLE sliders control the maximum angle a bird can turn as a result of each rule.
VISION is the distance that each bird can see 360 degrees around it.
## THINGS TO NOTICE
Central to the model is the observation that flocks form without a leader.
There are no random numbers used in this model, except to position the birds initially. The fluid, lifelike behavior of the birds is produced entirely by deterministic rules.
Also, notice that each flock is dynamic. A flock, once together, is not guaranteed to keep all of its members. Why do you think this is?
After running the model for a while, all of the birds have approximately the same heading. Why?
Sometimes a bird breaks away from its flock. How does this happen? You may need to slow down the model or run it step by step in order to observe this phenomenon.
## THINGS TO TRY
Play with the sliders to see if you can get tighter flocks, looser flocks, fewer flocks, more flocks, more or less splitting and joining of flocks, more or less rearranging of birds within flocks, etc.
You can turn off a rule entirely by setting that rule's angle slider to zero. Is one rule by itself enough to produce at least some flocking? What about two rules? What's missing from the resulting behavior when you leave out each rule?
Will running the model for a long time produce a static flock? Or will the birds never settle down to an unchanging formation? Remember, there are no random numbers used in this model.
## EXTENDING THE MODEL
Currently the birds can "see" all around them. What happens if birds can only see in front of them? The `in-cone` primitive can be used for this.
Is there some way to get V-shaped flocks, like migrating geese?
What happens if you put walls around the edges of the world that the birds can't fly into?
Can you get the birds to fly around obstacles in the middle of the world?
What would happen if you gave the birds different velocities? For example, you could make birds that are not near other birds fly faster to catch up to the flock. Or, you could simulate the diminished air resistance that birds experience when flying together by making them fly faster when in a group.
Are there other interesting ways you can make the birds different from each other? There could be random variation in the population, or you could have distinct "species" of bird.
## NETLOGO FEATURES
Notice the need for the `subtract-headings` primitive and special procedure for averaging groups of headings. Just subtracting the numbers, or averaging the numbers, doesn't give you the results you'd expect, because of the discontinuity where headings wrap back to 0 once they reach 360.
## RELATED MODELS
* Moths
* Flocking Vee Formation
* Flocking - Alternative Visualizations
## CREDITS AND REFERENCES
(Note: This is an adaptation by J.T. Lizier of the original Flocking model distributed in the NetLogo Models Library, under CC BY NC SA license (see below))
This model is inspired by the Boids simulation invented by Craig Reynolds. The algorithm we use here is roughly similar to the original Boids algorithm, but it is not the same. The exact details of the algorithm tend not to matter very much -- as long as you have alignment, separation, and cohesion, you will usually get flocking behavior resembling that produced by Reynolds' original model. Information on Boids is available at http://www.red3d.com/cwr/boids/.
## HOW TO CITE
If you mention this model or the NetLogo software in a publication, we ask that you include the citations below.
For the model itself:
* Wilensky, U. (1998). NetLogo Flocking model. http://ccl.northwestern.edu/netlogo/models/Flocking. Center for Connected Learning and Computer-Based Modeling, Northwestern University, Evanston, IL.
Please cite the NetLogo software as:
* Wilensky, U. (1999). NetLogo. http://ccl.northwestern.edu/netlogo/. Center for Connected Learning and Computer-Based Modeling, Northwestern University, Evanston, IL.
## COPYRIGHT AND LICENSE
Copyright 1998 Uri Wilensky.
![CC BY-NC-SA 3.0](http://ccl.northwestern.edu/images/creativecommons/byncsa.png)
This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 License. To view a copy of this license, visit https://creativecommons.org/licenses/by-nc-sa/3.0/ or send a letter to Creative Commons, 559 Nathan Abbott Way, Stanford, California 94305, USA.
Commercial licenses are also available. To inquire about commercial licenses, please contact Uri Wilensky at uri@northwestern.edu.
This model was created as part of the project: CONNECTED MATHEMATICS: MAKING SENSE OF COMPLEX PHENOMENA THROUGH BUILDING OBJECT-BASED PARALLEL MODELS (OBPML). The project gratefully acknowledges the support of the National Science Foundation (Applications of Advanced Technologies Program) -- grant numbers RED #9552950 and REC #9632612.
This model was converted to NetLogo as part of the projects: PARTICIPATORY SIMULATIONS: NETWORK-BASED DESIGN FOR SYSTEMS LEARNING IN CLASSROOMS and/or INTEGRATED SIMULATION AND MODELING ENVIRONMENT. The project gratefully acknowledges the support of the National Science Foundation (REPP & ROLE programs) -- grant numbers REC #9814682 and REC-0126227. Converted from StarLogoT to NetLogo, 2002.
<!-- 1998 2002 -->
@#$#@#$#@
default
true
0
Polygon -7500403 true true 150 5 40 250 150 205 260 250
airplane
true
0
Polygon -7500403 true true 150 0 135 15 120 60 120 105 15 165 15 195 120 180 135 240 105 270 120 285 150 270 180 285 210 270 165 240 180 180 285 195 285 165 180 105 180 60 165 15
arrow
true
0
Polygon -7500403 true true 150 0 0 150 105 150 105 293 195 293 195 150 300 150
box
false
0
Polygon -7500403 true true 150 285 285 225 285 75 150 135
Polygon -7500403 true true 150 135 15 75 150 15 285 75
Polygon -7500403 true true 15 75 15 225 150 285 150 135
Line -16777216 false 150 285 150 135
Line -16777216 false 150 135 15 75
Line -16777216 false 150 135 285 75
bug
true
0
Circle -7500403 true true 96 182 108
Circle -7500403 true true 110 127 80
Circle -7500403 true true 110 75 80
Line -7500403 true 150 100 80 30
Line -7500403 true 150 100 220 30
butterfly
true
0
Polygon -7500403 true true 150 165 209 199 225 225 225 255 195 270 165 255 150 240
Polygon -7500403 true true 150 165 89 198 75 225 75 255 105 270 135 255 150 240
Polygon -7500403 true true 139 148 100 105 55 90 25 90 10 105 10 135 25 180 40 195 85 194 139 163
Polygon -7500403 true true 162 150 200 105 245 90 275 90 290 105 290 135 275 180 260 195 215 195 162 165
Polygon -16777216 true false 150 255 135 225 120 150 135 120 150 105 165 120 180 150 165 225
Circle -16777216 true false 135 90 30
Line -16777216 false 150 105 195 60
Line -16777216 false 150 105 105 60
car
false
0
Polygon -7500403 true true 300 180 279 164 261 144 240 135 226 132 213 106 203 84 185 63 159 50 135 50 75 60 0 150 0 165 0 225 300 225 300 180
Circle -16777216 true false 180 180 90
Circle -16777216 true false 30 180 90
Polygon -16777216 true false 162 80 132 78 134 135 209 135 194 105 189 96 180 89
Circle -7500403 true true 47 195 58
Circle -7500403 true true 195 195 58
circle
false
0
Circle -7500403 true true 0 0 300
circle 2
false
0
Circle -7500403 true true 0 0 300
Circle -16777216 true false 30 30 240
cow
false
0
Polygon -7500403 true true 200 193 197 249 179 249 177 196 166 187 140 189 93 191 78 179 72 211 49 209 48 181 37 149 25 120 25 89 45 72 103 84 179 75 198 76 252 64 272 81 293 103 285 121 255 121 242 118 224 167
Polygon -7500403 true true 73 210 86 251 62 249 48 208
Polygon -7500403 true true 25 114 16 195 9 204 23 213 25 200 39 123
cylinder
false
0
Circle -7500403 true true 0 0 300
dot
false
0
Circle -7500403 true true 90 90 120
face happy
false
0
Circle -7500403 true true 8 8 285
Circle -16777216 true false 60 75 60
Circle -16777216 true false 180 75 60
Polygon -16777216 true false 150 255 90 239 62 213 47 191 67 179 90 203 109 218 150 225 192 218 210 203 227 181 251 194 236 217 212 240
face neutral
false
0
Circle -7500403 true true 8 7 285
Circle -16777216 true false 60 75 60
Circle -16777216 true false 180 75 60
Rectangle -16777216 true false 60 195 240 225
face sad
false
0
Circle -7500403 true true 8 8 285
Circle -16777216 true false 60 75 60
Circle -16777216 true false 180 75 60
Polygon -16777216 true false 150 168 90 184 62 210 47 232 67 244 90 220 109 205 150 198 192 205 210 220 227 242 251 229 236 206 212 183
fish
false
0
Polygon -1 true false 44 131 21 87 15 86 0 120 15 150 0 180 13 214 20 212 45 166
Polygon -1 true false 135 195 119 235 95 218 76 210 46 204 60 165
Polygon -1 true false 75 45 83 77 71 103 86 114 166 78 135 60
Polygon -7500403 true true 30 136 151 77 226 81 280 119 292 146 292 160 287 170 270 195 195 210 151 212 30 166
Circle -16777216 true false 215 106 30
flag
false
0
Rectangle -7500403 true true 60 15 75 300
Polygon -7500403 true true 90 150 270 90 90 30
Line -7500403 true 75 135 90 135
Line -7500403 true 75 45 90 45
flower
false
0
Polygon -10899396 true false 135 120 165 165 180 210 180 240 150 300 165 300 195 240 195 195 165 135
Circle -7500403 true true 85 132 38
Circle -7500403 true true 130 147 38
Circle -7500403 true true 192 85 38
Circle -7500403 true true 85 40 38
Circle -7500403 true true 177 40 38
Circle -7500403 true true 177 132 38
Circle -7500403 true true 70 85 38
Circle -7500403 true true 130 25 38
Circle -7500403 true true 96 51 108
Circle -16777216 true false 113 68 74
Polygon -10899396 true false 189 233 219 188 249 173 279 188 234 218
Polygon -10899396 true false 180 255 150 210 105 210 75 240 135 240
house
false
0
Rectangle -7500403 true true 45 120 255 285
Rectangle -16777216 true false 120 210 180 285
Polygon -7500403 true true 15 120 150 15 285 120
Line -16777216 false 30 120 270 120
leaf
false
0
Polygon -7500403 true true 150 210 135 195 120 210 60 210 30 195 60 180 60 165 15 135 30 120 15 105 40 104 45 90 60 90 90 105 105 120 120 120 105 60 120 60 135 30 150 15 165 30 180 60 195 60 180 120 195 120 210 105 240 90 255 90 263 104 285 105 270 120 285 135 240 165 240 180 270 195 240 210 180 210 165 195
Polygon -7500403 true true 135 195 135 240 120 255 105 255 105 285 135 285 165 240 165 195
line
true
0
Line -7500403 true 150 0 150 300
line half
true
0
Line -7500403 true 150 0 150 150
pentagon
false
0
Polygon -7500403 true true 150 15 15 120 60 285 240 285 285 120
person
false
0
Circle -7500403 true true 110 5 80
Polygon -7500403 true true 105 90 120 195 90 285 105 300 135 300 150 225 165 300 195 300 210 285 180 195 195 90
Rectangle -7500403 true true 127 79 172 94
Polygon -7500403 true true 195 90 240 150 225 180 165 105
Polygon -7500403 true true 105 90 60 150 75 180 135 105
plant
false
0
Rectangle -7500403 true true 135 90 165 300
Polygon -7500403 true true 135 255 90 210 45 195 75 255 135 285
Polygon -7500403 true true 165 255 210 210 255 195 225 255 165 285
Polygon -7500403 true true 135 180 90 135 45 120 75 180 135 210
Polygon -7500403 true true 165 180 165 210 225 180 255 120 210 135
Polygon -7500403 true true 135 105 90 60 45 45 75 105 135 135
Polygon -7500403 true true 165 105 165 135 225 105 255 45 210 60
Polygon -7500403 true true 135 90 120 45 150 15 180 45 165 90
square
false
0
Rectangle -7500403 true true 30 30 270 270
square 2
false
0
Rectangle -7500403 true true 30 30 270 270
Rectangle -16777216 true false 60 60 240 240
star
false
0
Polygon -7500403 true true 151 1 185 108 298 108 207 175 242 282 151 216 59 282 94 175 3 108 116 108
target
false
0
Circle -7500403 true true 0 0 300
Circle -16777216 true false 30 30 240
Circle -7500403 true true 60 60 180
Circle -16777216 true false 90 90 120
Circle -7500403 true true 120 120 60
tree
false
0
Circle -7500403 true true 118 3 94
Rectangle -6459832 true false 120 195 180 300
Circle -7500403 true true 65 21 108
Circle -7500403 true true 116 41 127
Circle -7500403 true true 45 90 120
Circle -7500403 true true 104 74 152
triangle
false
0
Polygon -7500403 true true 150 30 15 255 285 255
triangle 2
false
0
Polygon -7500403 true true 150 30 15 255 285 255
Polygon -16777216 true false 151 99 225 223 75 224
truck
false
0
Rectangle -7500403 true true 4 45 195 187
Polygon -7500403 true true 296 193 296 150 259 134 244 104 208 104 207 194
Rectangle -1 true false 195 60 195 105
Polygon -16777216 true false 238 112 252 141 219 141 218 112
Circle -16777216 true false 234 174 42
Rectangle -7500403 true true 181 185 214 194
Circle -16777216 true false 144 174 42
Circle -16777216 true false 24 174 42
Circle -7500403 false true 24 174 42
Circle -7500403 false true 144 174 42
Circle -7500403 false true 234 174 42
turtle
true
0
Polygon -10899396 true false 215 204 240 233 246 254 228 266 215 252 193 210
Polygon -10899396 true false 195 90 225 75 245 75 260 89 269 108 261 124 240 105 225 105 210 105
Polygon -10899396 true false 105 90 75 75 55 75 40 89 31 108 39 124 60 105 75 105 90 105
Polygon -10899396 true false 132 85 134 64 107 51 108 17 150 2 192 18 192 52 169 65 172 87
Polygon -10899396 true false 85 204 60 233 54 254 72 266 85 252 107 210
Polygon -7500403 true true 119 75 179 75 209 101 224 135 220 225 175 261 128 261 81 224 74 135 88 99
wheel
false
0
Circle -7500403 true true 3 3 294
Circle -16777216 true false 30 30 240
Line -7500403 true 150 285 150 15
Line -7500403 true 15 150 285 150
Circle -7500403 true true 120 120 60
Line -7500403 true 216 40 79 269
Line -7500403 true 40 84 269 221
Line -7500403 true 40 216 269 79
Line -7500403 true 84 40 221 269
x
false
0
Polygon -7500403 true true 270 75 225 30 30 225 75 270
Polygon -7500403 true true 30 75 75 30 270 225 225 270
@#$#@#$#@
NetLogo 6.0
@#$#@#$#@
set population 200
setup
repeat 200 [ go ]
@#$#@#$#@
@#$#@#$#@
@#$#@#$#@
@#$#@#$#@
default
0.0
-0.2 0 0.0 1.0
0.0 1 1.0 0.0
0.2 0 0.0 1.0
link direction
true
0
Line -7500403 true 150 150 90 180
Line -7500403 true 150 150 210 180
@#$#@#$#@
0
@#$#@#$#@

103
demos/octave/FlockingAnalysis/NetLogoExample/loadProperties.m Executable file
View File

@ -0,0 +1,103 @@
%
% This script loads properties for transfer entropy analysis of the data from
% the NetLogo Flocking model.
%
% Author: Joseph T. Lizier, 2019
%
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier
%%
%% This program is free software: you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation, either version 3 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
%%
% This script loads the default properties for the transfer entropy processing:
clear('properties');
%%%%%%%%%
% FILENAMES
%%%%%%%%%
% Input data files:
% properties.files can be:
% a. a cell array of file names, e.g.: {'file1.xlsx', 'file2.xlsx'}
% b. a call to ls or ls with an argument, e.g. ls('*.xlsx')
% c. a space or tab separated character row vector of file names
% d. a character matrix of filenames (each filename on a separate row)
properties.files = 'positions%s-small.txt';
% Function to read in the data files:
% loadScript must point to a function .m file that accepts two arguments
% (the name of a file, and properties object) and returns [x,y,z] (z optional, only when 3D)
% data where each is an array, e.g. x(time, fishIndex) indexed first by time and second by fish index.
% Use the name of the .m file after an "@" character:
properties.loadScript = @loadseparatexy;
% Is the data returned by the loadScript 3D (true) or 2D (false)?
properties.data3d = false;
% Results file - will hold the parsed velocities / relative positions, plus the
% local transfer entropy results
properties.resultsFile = 'results.mat';
%%%%%%%%%
% PARAMETERS
%%%%%%%%%
% Distance within which to consider a pair for the info theoretic analysis (units are as per what is used in the data files)
properties.pairRange = 4; % These ones have a causal range of 3
%%%%%%%%%
% INFORMATION THEORETIC Parameters
% Only lag is used for computing lagged mutual information
% All lag, k and tau are used for transfer entropy
% k - embedding dimension of the past of the destination array.
% tau - embedding delay: time cycles separating each element in the past of the destination.
% lag - time delay between the source and target in cycles
% You can set kRange, tauRange and lagRange to ask that these are optimised by runAnalysis:
% properties.kRange = 1:10;
% properties.tauRange = 1:4;
% properties.lagRange = 1:10;
% You can also set k, tau and lag to values that generateObservations should use
% (although note that if this is called via runAnalysis then it will overwrite them):
properties.k = 1;
properties.tau = 1;
properties.lag = 1;
% Do we include the relative source position in the transfer entropy calculation (true), or
% only the relative source heading (false)
properties.includeSourcePositionInTransfer = false;
% Do we take relative source heading and position with respect to dest heading at that same
% time point (true, this is what we did for Crosato paper) or with respect
% to dest heading just previous to state update (false)?
properties.sourceWrtSameDestTime = true;
% JIDT location:
properties.jidtJarLocation = '../../../../infodynamics.jar';
% Which estimator to use.
% Valid values are 'gaussian' (linear) or 'kraskov' (non-linear)
% properties.estimator = 'gaussian';
properties.estimator = 'kraskov';
% Properties for JIDT estimators:
properties.jidt.kNNs = 4; % Number of nearest neighbours for Kraskov algorithm: just use 4 (default)
properties.jidt.autoDynamicCorrelationExclusion = true; % Exclude nearest neighbours from at least the same target transition from being included in counts for TE
properties.aisNumSurrogates = 0; % Number of surrogate calculations to run for AIS (just to see the noise floor. 0 means skip)
properties.teNumSurrogates = 0; % Number of surrogate calculations to run for TE (just to see the noise floor. 0 means skip)

73
demos/octave/FlockingAnalysis/computeAIS.m Executable file
View File

@ -0,0 +1,73 @@
function ais = computeAIS(D, Dpast, properties)
% Computes active information storage from the pre-processed velocity data
%
% Author: Joseph T. Lizier, 2019
%
% Inputs:
% - D - target samples (may be multivariate as per generateObservations)
% - Dpast - target past samples (multivariate, and embedded up to k previous samples)
% - properties (required) - object with properties for the calculations,
% with sub-members as specificied in the loadProperties.m file. If not supplied
% the properties are loaded from loadProperties.m
%
% Outputs:
% - ais - active information storage value (MI between D and Dpast)
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% This program is free software: you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation, either version 3 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
%%
% -- STEP 0 : create java object
% check out http://lizier.me/joseph/software/jidt/javadocs/v1.3/
% for the description of all classes and methods
javaaddpath(properties.jidtJarLocation); % add JIDT path
% Select correct MI class for this estimator type
if (strcmp('kraskov', properties.estimator))
MI_CLASS = 'infodynamics.measures.continuous.kraskov.MutualInfoCalculatorMultiVariateKraskov1';
else
MI_CLASS = 'infodynamics.measures.continuous.gaussian.MutualInfoCalculatorMultiVariateGaussian';
end
% Compute the information stored in the past of the target.
% Creates java object of the given class
AIScalculator = javaObject(MI_CLASS);
% -- STEP 1 : set properties
AIScalculator.setProperty('k', num2str(properties.jidt.kNNs));
AIScalculator.setProperty('BIAS_CORRECTION', 'true'); % Used for Gaussian only
% -- STEP 2 : initialise
% here the parameters are the dimensionality of the series
% in this case taken directly from the number of columns in each variable
AIScalculator.initialise(size(D,2), size(Dpast,2));
% -- STEP 3 : add in observations
AIScalculator.setObservations(D, Dpast);
% -- STEP 4 : compute the local AIS (will be bias-corrected now for either Gaussian or KSG)
ais = AIScalculator.computeAverageLocalOfObservations(); % global (average) value
fprintf('Mean AIS_%s (k=%d,tau=%d) = %.3f\n', properties.estimator, properties.k, properties.tau, ais);
if (properties.aisNumSurrogates > 0)
% Compute the (statistical significance via) null distribution empirically (e.g. with 100 permutations),
% and use this for empirical bias correction (otherwise we're relying on analytic)
aisMeasDist = AIScalculator.computeSignificance(properties.aisNumSurrogates);
fprintf('Null distribution: %.4f +/- %.4f std dev.; p(surrogate > measured)=%.3f from %d surrogates)\n', ...
aisMeasDist.getMeanOfDistribution(), aisMeasDist.getStdOfDistribution(), ...
aisMeasDist.pValue, properties.aisNumSurrogates);
ais = ais - aisMeasDist.getMeanOfDistribution();
fprintf('Bias corrected Mean AIS_%s (k=%d,tau=%d) = %.3f\n', properties.estimator, properties.k, properties.tau, ais);
end
end

101
demos/octave/FlockingAnalysis/computeTE.m Executable file
View File

@ -0,0 +1,101 @@
function tranEntropy = computeTE(S, D, Dpast, properties)
% Computes transfer entropy from the pre-processed velocity data
%
% Author: Emanuele Crosato, Joseph T. Lizier, 2019
%
% Inputs:
% - S - source samples (may be multivariate as per generate3DObservations)
% - D - target samples (may be multivariate as per generate3DObservations)
% - Dpast - target past samples (multivariate, and embedded up to k previous samples)
% - properties (required) - object with properties for the calculations,
% with sub-members as specificied in the loadProperties.m file. If not supplied
% the properties are loaded from loadProperties.m
%
% Outputs:
% - te - transfer entropy value (conditional MI from S (maybe plus RelSourcePos) to D given Dpast). If not requested, then
% the te (and an array of local values) is saved to properties.resultsFile
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% This program is free software: you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation, either version 3 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
%%
% -- STEP 0 : create java object
% check out http://lizier.me/joseph/software/jidt/javadocs/v1.3/
% for the description of all classes and methods
javaaddpath(properties.jidtJarLocation); % add JIDT path
% the following is the java class for conditional mutual information in JIDT
% transfer entropy is infact mutual information conditioning on the past
% of the destination
if (strcmp('kraskov', properties.estimator))
CMI_CLASS = 'infodynamics.measures.continuous.kraskov.ConditionalMutualInfoCalculatorMultiVariateKraskov1';
else
CMI_CLASS = 'infodynamics.measures.continuous.gaussian.ConditionalMutualInfoCalculatorMultiVariateGaussian';
end
TEcalculator = javaObject(CMI_CLASS); % creates a java object of the given class
% -- STEP 1 : set properties
TEcalculator.setProperty('k', num2str(properties.jidt.kNNs));
TEcalculator.setProperty('BIAS_CORRECTION', 'true'); % Used for Gaussian only
if (isfield(properties.jidt, 'dynamicCorrelationExclusion'))
% We'll ensure samples from the same target transition aren't included in nearest neighbour counts
% (it will exclude some others as well, but this only adds some small noise to the calculation)
TEcalculator.setProperty('DYN_CORR_EXCL', num2str(properties.jidt.dynamicCorrelationExclusion));
end
% -- STEP 2 : initialise
% here the parameters are the dimensionality of the series
% in this case taken directly from the number of columns in each variable
TEcalculator.initialise(size(S,2), size(D,2), size(Dpast,2));
% -- STEP 3 : add in observations
TEcalculator.setObservations(S, D, Dpast);
% -- STEP 4 : compute the local entropies
tranEntropy = TEcalculator.computeAverageLocalOfObservations(); % global (average) value
fprintf('Mean TE_%s (k=%d,tau=%d,lag=%d) = %.3f\n', ...
properties.estimator, properties.k, properties.tau, properties.lag, tranEntropy);
if (properties.teNumSurrogates > 0)
% Compute the (statistical significance via) null distribution empirically (e.g. with 100 permutations):
measDist = TEcalculator.computeSignificance(properties.teNumSurrogates);
fprintf('Null distribution: %.4f +/- %.4f std dev.; p(surrogate > measured)=%.5f from %d surrogates)\n', ...
measDist.getMeanOfDistribution(), measDist.getStdOfDistribution(), ...
measDist.pValue, properties.teNumSurrogates);
pValue = measDist.pValue;
meanOfSurrogates = measDist.getMeanOfDistribution();
stdOfSurrogates = measDist.getStdOfDistribution();
else
pValue = 1;
meanOfSurrogates = 0;
stdOfSurrogates = 0;
end
if (nargout >= 1)
% Supply the samples back to the caller
% (the caller is probably trying to optimise parameters at the moment)
% Nothing to do then actually...
else
% We're going to save the results instead
% First generate the local values to save as well:
localTranEntropy = TEcalculator.computeLocalOfPreviousObservations(); % local values
% save results
save(properties.resultsFile, 'tranEntropy', 'localTranEntropy', 'pValue', 'meanOfSurrogates', 'stdOfSurrogates', '-append');
fprintf('Transfer entropy saved in %s\n', properties.resultsFile);
end
end

378
demos/octave/FlockingAnalysis/generateObservations.m Executable file
View File

@ -0,0 +1,378 @@
function [D, Dpast, S, RelSourcePos, maxSourceSamplesForATarget] = generateObservations(properties)
% This function generates the observations from which
% we can then compute information dynamics with JIDT.
% This will work for either 2D or 3D samples (as specified by the properties)
%
% Author: Emanuele Crosato, Joseph T. Lizier, 2019
%
% Inputs:
% - properties (required) - object with properties for the calculations,
% with sub-members as specificied in the loadProperties.m file.
%
% Outputs:
% - D - target samples (may be multivariate as per below)
% - Dpast - target past samples (multivariate, and embedded up to k previous samples)
% - S - source relative headings samples (may be multivariate as per below)
% - RelSourcePos - relative source position (may be multivariate as below)
% - maxSourceSamplesForATarget - maximum number of sources in range for a specific target sample. Used for dynamic correlation exclusion externally for TE calculations.
% If no outputs are requested, these are saved to properties.resultsFile
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% This program is free software: you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation, either version 3 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
%%
% Call internal utility to put filename lists in a common format
files = processFilenames(properties.files);
% initialize series for storing the samples:
S = []; % initialise source observations
D = []; % initialise destination observations
Dpast = []; % initialise destination past observations
% only S and D are necessary for lagged mutual information
% all S, D and Dpast are necessary for transfer entropy
fileTimeAndPair = []; % to save time and pair
RelSourcePos = []; % initialise relative source positions
sample = 1; % initialise the current sample number
maxSourceSamplesForATarget = 0; % Track the maximum number of in range sources for a given target sample
% Need to loop over fileIndex rather than for file = files (which doesn't work properly for cell array of length 1)
for fileIndex = 1:length(files)
dataFileName = files{fileIndex};
% load preprocessed data using the function specified in properties.loadScript
if (properties.data3d)
[x,y,z] = feval(properties.loadScript, dataFileName, properties);
numMissing = sum(sum(isnan([x,y,z])));
else
[x,y] = feval(properties.loadScript, dataFileName, properties);
numMissing = sum(sum(isnan([x,y])));
end
fprintf('Loading data in %s (%d missing values)\n', dataFileName, numMissing);
samplesBeforeThisFile = sample;
% Translate the raw positions into delta Positions
velX = x(2:end,:) - x(1:end-1,:);
velY = y(2:end,:) - y(1:end-1,:);
if (properties.data3d)
velZ = z(2:end,:) - z(1:end-1,:);
end
% Translate velocities into headings:
if (properties.data3d)
% Spherical polars:
[headingXY,headingZ,speed] = cart2sph(velX,velY,velZ);
else
% Polar coordinates:
[headingXY,speed] = cart2pol(velX, velY);
end
% Manual way is (can verify these are same except for Nans on x,y and z):
% headingXY = atan(velY ./ velX) + (((velX < 0).*(velY>0)) .* pi) + ...
% (((velX < 0).*(velY<=0)) .* (- pi));
% But the above fails if x and y are *both* zero leaving Nan:
% velYOnVelX = velY ./ velX; % Don't replace all nans in headingXY as some mean missing values
% headingXY(isnan(velYOnVelX(:))) = 0;
% xyMagnitude = sqrt(velX.^2 + velY.^2);
% Since xyMagnitude can only be positive, we can take the straight atan to
% compute headingZ
% headingZ = atan(velZ ./ xyMagnitude);
% get number of fish and update time cycles
numCycles = size(velX,1);
numFish = size(velX,2);
% Initialising destPastSample is only important in terms of ensuring it is a row vector.
% If we have 3D data, the vector will get padded out to the appropriate length with the first sample below.
destPastSample = zeros(1, properties.k);
startTime = max(1+properties.lag, (properties.k-1)*properties.tau + 3); % Adding 3: one for target, one for first target past, one for taking differences
for i = startTime : numCycles % cycle over time
timePointForSourceHeading = i-properties.lag; % This is indexed into velX and headingXY, hence no extra +1 !
timePointForSourcePosition = i+1-properties.lag; % This is indexed into x not velX, hence the extra +1 !
% Compute relative position of source (at time timePointForSourcePosition) to
% target either at this same time step or the current time at which it is updating.
% Note this position difference is relative to absolute Cartesian coordinates
% (we'll convert to relative to source heading later):
if (properties.sourceWrtSameDestTime)
destPositionTimePointRef = timePointForSourcePosition;
else
destPositionTimePointRef = i;
end
% cycle over fish pairs
for idxFD = 1 : numFish % Target/Destination
% check destination variable
if isnan(headingXY(i,idxFD)) || isnan(headingXY(i-1,idxFD))
continue;
end
if properties.data3d && (isnan(headingZ(i,idxFD)) || isnan(headingZ(i-1,idxFD)))
continue;
end
% check destination past vector
missingFound = false;
for h = 1 : properties.k
idx = i-1-(h-1)*properties.tau;
if isnan(headingXY(idx,idxFD)) || isnan(headingXY(idx-1,idxFD))
missingFound = true;
break;
end
if properties.data3d && (isnan(headingZ(idx,idxFD)) || isnan(headingZ(idx-1,idxFD)))
missingFound = true;
break;
end
end
if missingFound
continue;
end
% Postcondition: All destination variables are ok
% We will create **source** observation as source headings relative to target headings at appropriate time point:
if (properties.sourceWrtSameDestTime)
% Take reference dest heading at same time as source:
theta_FDXY_ref = headingXY(timePointForSourceHeading,idxFD);
else
% Take reference dest heading at prev time step:
theta_FDXY_ref = headingXY(i-1,idxFD);
end
if properties.data3d
if (properties.sourceWrtSameDestTime)
theta_FDZ_ref = headingZ(timePointForSourceHeading,idxFD);
else
theta_FDZ_ref = headingZ(i-1,idxFD);
end
end
% create **destination** observation as change in headings:
theta_FDXY_curr = headingXY(i,idxFD);
theta_FDXY_prev = headingXY(i-1,idxFD);
if properties.data3d
theta_FDZ_curr = headingZ(i,idxFD);
theta_FDZ_prev = headingZ(i-1,idxFD);
destSample = [angleDifference(theta_FDXY_curr, theta_FDXY_prev), ...
angleDifference(theta_FDZ_curr, theta_FDZ_prev)];
else
destSample = angleDifference(theta_FDXY_curr, theta_FDXY_prev);
end
% create **destination past** observation as changes in headings at each step:
% TODO: we could take differences to previous sample amongst the k rather than only
% one back from each sample: I'm not sure if this would be a more wholistic embedding or not
% (only makes a difference if tau>1)
DpastColIndex = 1;
for h = 1 : properties.k
idx = i-1-(h-1)*properties.tau;
theta_FDXY_curr = headingXY(idx,idxFD);
theta_FDXY_prev = headingXY(idx-1,idxFD);
destPastSample(DpastColIndex) = angleDifference(theta_FDXY_curr, theta_FDXY_prev);
DpastColIndex = DpastColIndex + 1;
if properties.data3d
theta_FDZ_curr = headingZ(idx,idxFD);
theta_FDZ_prev = headingZ(idx-1,idxFD);
destPastSample(DpastColIndex) = angleDifference(theta_FDZ_curr, theta_FDZ_prev);
DpastColIndex = DpastColIndex + 1;
end
end
if (isfield(properties, 'destSamplesOnly'))
if (properties.destSamplesOnly)
% User has asked for [D,Dpast] samples only to be returned,
% so we can do these now (without looping over sources):
% Fill in the destination and destination next samples now from above:
D(sample, :) = destSample;
Dpast(sample, :) = destPastSample;
fileTimeAndPair(sample,:) = [fileIndex i idxFD nan];
% increment sample number
sample = sample + 1;
continue; % skip looping over the sources
end
end
numSourceSamplesForThisTarget = 0;
for idxFS = 1 : numFish % Source
% check not the same fish
if (idxFD == idxFS)
continue;
end
relXOfSource = x(timePointForSourcePosition,idxFS) - x(destPositionTimePointRef,idxFD);
relYOfSource = y(timePointForSourcePosition,idxFS) - y(destPositionTimePointRef,idxFD);
if (properties.data3d)
relZOfSource = z(timePointForSourcePosition,idxFS) - z(destPositionTimePointRef,idxFD);
[xyAbsoluteAngleOfSource,zAbsoluteAngleOfSource,distanceBetween] = ...
cart2sph(relXOfSource,relYOfSource,relZOfSource);
else
[xyAbsoluteAngleOfSource,distanceBetween] = ...
cart2pol(relXOfSource,relYOfSource);
end
% Manually: (verified this matches cart2sph):
% xyAbsoluteAngleOfSource = atan(relYOfSource ./ relXOfSource) + ...
% (((relXOfSource < 0).*(relYOfSource>0)) .* pi) + ...
% (((relXOfSource < 0).*(relYOfSource<=0)) .* (- pi));
% xyRelMagnitude = sqrt(relXOfSource.^2 + relYOfSource.^2);
% zAbsoluteAngleOfSource = atan(relZOfSource ./ xyRelMagnitude);
% distanceBetween = sqrt(relXOfSource.^2 + relYOfSource.^2 + ...
% relZOfSource.^2);
% check in range
if (distanceBetween > properties.pairRange)
continue;
end
% check source variable
if isnan(headingXY(timePointForSourceHeading,idxFD)) || isnan(headingXY(timePointForSourceHeading,idxFS))
continue;
end
if properties.data3d && (isnan(headingZ(timePointForSourceHeading,idxFD)) || isnan(headingZ(timePointForSourceHeading,idxFS)))
continue;
end
% Postcondition: There are no missing headings so we can generate an observation.
numSourceSamplesForThisTarget = numSourceSamplesForThisTarget + 1;
% Now compose the data that will be saved for this sample:
% fileTimeAndPair is [file_index, time_index, target_index, source_index]
fileTimeAndPair(sample,:) = [fileIndex i idxFD idxFS];
% create **source** observation as source headings relative to target headings at appropriate time point:
theta_FSXY_lag = headingXY(timePointForSourceHeading,idxFS);
% And convert the absolute angular positions into
% relative angular positions compared to the target's
% heading.
xyRelativeAngleOfSource = angleDifference(...
xyAbsoluteAngleOfSource, ...
theta_FDXY_ref);
if properties.data3d
theta_FSZ_lag = headingZ(timePointForSourceHeading,idxFS);
sourceSample = [angleDifference(theta_FDXY_ref, theta_FSXY_lag), ...
angleDifference(theta_FDZ_ref, theta_FSZ_lag)];
% Elevation angle differences need to be in -pi/2,pi/2 range
zRelativeAngleOfSource = angleDifferencePiOn2(...
zAbsoluteAngleOfSource, ...
theta_FDZ_ref);
% Store these relative polar coordinates of source at timePointForSourcePosition
RelSourcePos(sample,:) = [distanceBetween, xyRelativeAngleOfSource, zRelativeAngleOfSource];
else
sourceSample = angleDifference(theta_FDXY_ref, theta_FSXY_lag);
% Store these relative polar coordinates of source at timePointForSourcePosition
RelSourcePos(sample,:) = [distanceBetween, xyRelativeAngleOfSource];
end
if properties.includeSourcePositionInTransfer
S(sample,:) = [sourceSample, RelSourcePos];
else
S(sample,:) = sourceSample;
end
% Fill in the destination and destination next samples now from above:
D(sample, :) = destSample;
Dpast(sample, :) = destPastSample;
% increment sample number
sample = sample + 1;
end
if (numSourceSamplesForThisTarget > maxSourceSamplesForATarget)
maxSourceSamplesForATarget = numSourceSamplesForThisTarget;
end
end
% fprintf('Run time step %d\n', i);
end
fprintf(' added %d samples\n', sample - samplesBeforeThisFile);
end
if (sample - 1 == 0)
% We've added no samples
warning('No samples added for the given parameters!');
end
if (nargout > 1)
% Supply the samples back to the caller
% (the caller is probably trying to optimise parameters at the moment)
% Nothing to do then actually...
else
% We're going to save the samples instead
% display to check
fprintf('Displaying first 5 samples for source, target, target past and fileTimeAndPair as a check:');
disp(S(1:5,:));
disp(D(1:5,:));
disp(Dpast(1:5,:));
disp(fileTimeAndPair(1:5,:));
% input(prompt);
% save series and properties
save(properties.resultsFile, 'S', 'D', 'Dpast', 'files', 'fileTimeAndPair', 'RelSourcePos', ...
'maxSourceSamplesForATarget', 'properties');
fprintf('Series saved in %s (%d samples in total)\n', properties.resultsFile, sample - 1);
end
end
% End function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for computing the difference between two angles.
% the differrence must be between pi and -pi
function [diff] = angleDifference(angleA, angleB)
diff = angleA - angleB; % subtract angles
if abs(diff) > pi % if absolute value is larger than pi
% replace with the complementary angle and switch sign
diff = (2*pi - abs(diff) ) * (-sign(diff));
end
end
% function for computing the difference between two angles in [-pi/2,pi/2].
% the differrence must be returned between pi/2 and -pi/2.
% This is used for differences in elevation angles
function [diff] = angleDifferencePiOn2(angleA, angleB)
diff = angleA - angleB; % subtract angles
% Pre-condition: differences between angles which were in range of
% [-pi/2,pi/2] can only be in range [-pi,pi]
if (diff > pi/2)
diff = pi/2 - (diff - pi/2);
elseif (diff < -pi/2)
diff = -pi/2 + (-pi/2 - diff);
end
end
% Turns the fileList into a cell array of file names. The fileList can be either:
% a. a cell array of file names, e.g.: {'file1.xlsx', 'file2.xlsx'}
% b. a call to ls or ls with an argument, e.g. ls('*.xlsx')
% c. a space or tab separated character row vector of file names
% d. a character matrix of filenames (each filename on a separate row)
function fileCellArray = processFilenames(fileList)
if (iscell(fileList))
% We're done already:
fileCellArray = fileList;
elseif (isvector(fileList))
% We have a row vector of space/tab separate filenames:
fileCellArray = strsplit(strtrim(fileList)); % extra strtrim to remove trailing \n's
elseif (ismatrix(fileList))
fileCellArray = {};
for r = 1 : size(fileList, 1)
fileCellArray{r} = strtrim(fileList(r,:));
end
else
error('fileList appears to be of an incorrect format\n');
end
end

66
demos/octave/FlockingAnalysis/loadBasic2d.m Normal file
View File

@ -0,0 +1,66 @@
function [posX,posY] = loadBasic2d(dataFileName, properties)
% This script loads the raw data from a .txt file,
% preprocesses the data and save it as a .mat file.
% The txt data is assumed to have timestamp in column 1,
% then fish 1's X and Y coordinates in columns 2 and 3,
% then fish 2's X and Y coordinates in columns 4 and 5,
% and so on.
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% This program is free software: you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation, either version 3 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
%%
% prompt = 'Press a key to continue';
%%% LOAD RAW DATA %%%
%fprintf('Loading data in %s\n', dataFileName); % print to check
% input(prompt);
% read the .txt as a Matlab table (assuming .txt is tab-separated)
data = load(dataFileName);
%fprintf('Size of data is %d - %d\n', ... % print table's size - %d is for int
% size(data,1), size(data,2)); % see also %f (real) and %s (string)
% disp(data(1:5,:)); % display first 10 rows to check
% input(prompt);
%%% FORMAT DATA IN A MORE CONVENIENT WAY %%%
numFish = (size(data,2)-1) ./ 2; % number of fish (we know it from the raw data)
numCycles = size(data,1); % number of time steps (we also know it)
fprintf('Number of fish %d and cycles %d\n', numFish, numCycles);
% prepare variables for x and y position
% as a table [numFish x numCycles]
posX = nan(numCycles,numFish);
posY = nan(numCycles,numFish);
% fill the position tables
for f = 1 : numFish % cycle over all fish
% copy into new variables
posX(:,f) = data(:,2+(f-1)*2);
posY(:,f) = data(:,3+(f-1)*2);
end
% display to check
% disp(size(posX));
% disp(posX(1:5,:));
% disp(size(posY));
% disp(posY(1:5,:));
% input(prompt);

40
demos/octave/FlockingAnalysis/loadseparatexy.m Executable file
View File

@ -0,0 +1,40 @@
function [x,y] = loadseparatexy(filename, properties)
% loadseparatexy loads 2D fish data from 2 separate txt files (one for x, one for y)
% where in each file time increases down the
% rows and then across the columns we have position columns for each
% fish in turn, i.e. in position x file:
% <fish1x>, <fish2x>, <fish3x>, etc
%
% Inputs:
% - filename - the name template of the file to load, with %s where 'x' and 'y' should be filled in
% - properties (not required) - properties object (may be required for other file loaders)
% Outputs:
% - x - 2D array, each row contains x position for each fish (in columns)
% - y - as per x
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% This program is free software: you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation, either version 3 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
%%
xfilename = sprintf(filename, 'x');
x = load(xfilename);
yfilename = sprintf(filename, 'y');
y = load(yfilename);
end

43
demos/octave/FlockingAnalysis/loadseparatexyz.m Executable file
View File

@ -0,0 +1,43 @@
function [x,y,z] = loadseparatexy(filename, properties)
% loadseparatexy loads 3D fish data from 3 separate txt files (one for x, one for y, one for z)
% where in each file time increases down the
% rows and then across the columns we have position columns for each
% fish in turn, i.e. in position x file:
% <fish1x>, <fish2x>, <fish3x>, etc
%
% Inputs:
% - filename - the name template of the file to load, with %s where 'x', 'y' and 'z' should be filled in
% - properties (not required) - properties object (may be required for other file loaders)
% Outputs:
% - x - 2D array, each row contains x position for each fish (in columns)
% - y - as per x
% - z - as per x
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% This program is free software: you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation, either version 3 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
%%
xfilename = sprintf(filename, 'x');
x = load(xfilename);
yfilename = sprintf(filename, 'y');
y = load(yfilename);
zfilename = sprintf(filename, 'z');
z = load(zfilename);
end

45
demos/octave/FlockingAnalysis/loadxls3d.m Executable file
View File

@ -0,0 +1,45 @@
function [x,y,z] = loadxls3d(filename, properties)
% loadxls3d loads 3D fish data from an xls file where time increases down the
% rows and then across the columns we have 3 x,y,z position columns for each
% fish in turn, i.e.:
% <timestamp1>, <fish1x>, <fish1y>, <fish1z>, <fish2x>, <fish2y>, <fish2z>, etc
% <timestamp2>, <fish1x>, <fish1y>, <fish1z>, <fish2x>, <fish2y>, <fish2z>, etc
%
% Inputs:
% - filename - the name of the file to load
% - properties (not required) - properties object (may be required for other file loaders)
% Outputs:
% - x - 2D array, each row contains x position for each fish (in columns)
% - y - as per x
% - z - as per x
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% This program is free software: you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation, either version 3 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
%%
[data,txt,raw] = xlsread(filename);
% Make sure we preprocess data to have x,y, and z as 2D arrays of x(timeStep, fishID):
% timeSteps = size(data,1); % number of rows; not required
xFishIndex = 2 : 3 : size(data,2); % which indices are the x values for different fish
x = data(:,xFishIndex);
y = data(:,xFishIndex+1);
z = data(:,xFishIndex+2);
end

206
demos/octave/FlockingAnalysis/plotLocalTEs.m Executable file
View File

@ -0,0 +1,206 @@
function plotLocalTEs(properties)
% Plot the local TEs to show where the information transfer hotspots are from target fish relative to each source
%
% Author: Joseph T. Lizier, 2019
%
% Inputs:
% - properties (required) - object with properties for the calculations,
% with sub-members as specificied in the loadProperties.m file. If not supplied
% the properties are loaded from loadProperties.m
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% This program is free software: you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation, either version 3 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
%%
if (nargin < 1)
fprintf('No properties object supplied, attempting to load properties via a loadProperties script ...');
% By default, just try to load properties locally
if (exist('loadProperties') == 2)
% there is a loadProperties script
loadProperties;
else
% there is not a loadProperties script
error('No properties object supplied, and no loadProperties script found.');
end
end
load(properties.resultsFile);
% Loads:
% S -- source samples
% D -- target samples
% Dpast -- target past samples (embedded up to k previous samples)
% files -- cell array of file names that we took samples from
% fileTimeAndPair -- each row holds file index, time index, target index, source index
% RelSourcePos -- each row holds distance between the pair for this sample,
% their xyRelativeAngleOfSource, and zRelativeAngleOfSource
% lag -- source-target lag that is in use
% k -- embedding length for target that is in use
% tau -- embedding delay for target that is in use
% pairRange -- range within which we've pulled source-target interactions
% tranEntropy -- average transfer entropy
% localTranEntropy -- local TE for each sample
fprintf('%d samples in total for %d fish\n', length(S), length(unique(fileTimeAndPair(:,3))));
relDistance = RelSourcePos(:,1);
relTheta = RelSourcePos(:,2);
% I think this is giving us the right conversions:
if (properties.data3d)
relPhi = RelSourcePos(:,3);
distXY = relDistance .* cos(relPhi);
distZ = relDistance .* sin(relPhi);
else
distXY = relDistance;
end
distInFront = distXY .* cos(relTheta);
distToLeft = distXY .* sin(relTheta);
% Plot where all the raw positions are:
% Will need to turn this off when we have too many
figure()
% polar(relTheta, distXY, '.r'); % This is equivalent to below:
scatter(distInFront, distToLeft, 2, localTranEntropy);
title('Relative position of source in XY plane for target heading, coloured for TE');
colorbar;
% Plot the density of samples
makePolarBinnedPlot(relTheta, distXY, ones(length(relTheta), 1), 12, 10, true, false);
title('Density of samples in each bin (r_{XY},\theta)');
% Plot the TE in XY plane
makePolarBinnedPlot(relTheta, distXY, localTranEntropy, 12, 10, true, true);
title('Average TE in each bin (r_{XY},\theta)');
% Plot raw positions in phi-z:
% figure()
% polar(relPhi, relDistance, '.r');
% title('Relative position of source in Z-phi plane for target heading');
if (properties.data3d)
% Plot the density of samples in distance-phi plane
makePolarBinnedPlot(relPhi, relDistance, ones(length(relTheta), 1), 12, 10, true, false);
title('Density of samples in each bin (r, \phi)');
% Plot the TE in distance-phi plane
makePolarBinnedPlot(relPhi, relDistance, localTranEntropy, 12, 10, true, true);
title('Average TE in each bin (r, \phi)');
end
end
% Inputs:
% - thetas - angles for each sample
% - radii - radius for each sample
% - numAngleBins - how many bins to make across 2*pi
% - numRadialBins - how many bins to make up to the maximum radii
% - useMaxEntBinning - whether to make bins with approx same numbers of points (true)
% or same size (false)
% - plotMean - if true (default) plot the mean within each bin, else plot the total (divded by area)
% The latter is used for densities for example
function makePolarBinnedPlot(thetas, radii, valuesToPlot, numAngleBins, numRadialBins, useMaxEntBinning, plotMean)
if (nargin < 6)
useMaxEntBinning = false;
end
if (nargin < 7)
plotMean = true;
end
if ((min(thetas) < -pi/2) || (max(thetas) > pi/2))
% We're using full angular range -pi : pi
minAngle = -pi;
maxAngle = pi;
extraBinForPlotWrap = true;
else
% We're only using -pi/2:pi/2
minAngle = -pi/2;
maxAngle = pi/2;
extraBinForPlotWrap = false;
end
angleStep = (maxAngle - minAngle) / numAngleBins;
radiusStep = max(radii) / numRadialBins;
% Simple way to do the binning for even bins:
% binnedAngles = floor(thetas ./ angleStep); % Gives the discrete bin for the angle
% binnedRadii = floor(radii ./ radiusStep); % Gives the discrete bin for the radius
% More general, and allowing bins to spread with points:
if (useMaxEntBinning)
% Space the bins for roughly same
% numbers of points (when examined marginally):
sortedAngles = sort(thetas);
binAngleEdges = [minAngle; sortedAngles(floor((1:(numAngleBins-1)).*length(sortedAngles)./numAngleBins)); maxAngle]';
sortedRadii = sort(radii);
binRadiusEdges = [0; sortedRadii(floor((1:(numRadialBins-1)).*length(sortedRadii)./numRadialBins)); max(radii)]';
else
% Space the bins equally
binAngleEdges = minAngle:angleStep:maxAngle;
binRadiusEdges = 0:radiusStep:max(radii);
end
[angleHistCounts,binnedAngles] = histc(thetas, binAngleEdges);
[radiiHistCounts,binnedRadii] = histc(radii, binRadiusEdges);
angleBinValues = 1:numAngleBins; % unique(binnedAngles);
radiusBinValues = 1:numRadialBins; % unique(binnedRadii);
if (extraBinForPlotWrap)
valuesForEachBin = zeros(length(angleBinValues) + 1, length(radiusBinValues));
else
valuesForEachBin = zeros(length(angleBinValues), length(radiusBinValues));
end
numberOfSamples = 0;
for aIndex = 1 : length(angleBinValues)
for rIndex = 1 : length(radiusBinValues)
indicesForThisBin = find((binnedAngles == angleBinValues(aIndex)) & (binnedRadii == radiusBinValues(rIndex)));
if (plotMean)
valueForThisBin = mean(valuesToPlot(indicesForThisBin));
else
% Plot a density: compute total then divide by area.
valueForThisBin = sum(valuesToPlot(indicesForThisBin));
areaOfBin = pi .* (binRadiusEdges(rIndex+1).^2 - binRadiusEdges(rIndex).^2) .* ...
mod(abs(binAngleEdges(aIndex+1) - binAngleEdges(aIndex)), 2.*pi) ./ (2.*pi);
valueForThisBin = valueForThisBin ./ areaOfBin;
end
if (length(indicesForThisBin) == 0)
valueForThisBin = 0;
end
% fprintf('Mean value for r=%.1f+,theta=%.3f+ is %.3f (from %d samples)\n', binRadiusEdges(rIndex), ...
% binAngleEdges(aIndex), valueForThisBin, length(indicesForThisBin));
numberOfSamples = numberOfSamples + length(indicesForThisBin);
valuesForEachBin(aIndex, rIndex) = valueForThisBin;
end
end
% Now convert these so we can plot them:
fprintf('Found %d in total in the bins\n', numberOfSamples);
if (extraBinForPlotWrap)
% And add for first angle again to complete the plot
valuesForEachBin(end,:) = valuesForEachBin(1,:);
[THETA,RR] = meshgrid([(binAngleEdges(1:end-1)+binAngleEdges(2:end))./2, (binAngleEdges(1)+binAngleEdges(2))./2], ...
(binRadiusEdges(1:end-1)+binRadiusEdges(2:end))./2);
else
[THETA,RR] = meshgrid([(binAngleEdges(1:end-1)+binAngleEdges(2:end))./2], ...
(binRadiusEdges(1:end-1)+binRadiusEdges(2:end))./2);
end
[A,B] = pol2cart(THETA,RR);
figure();
surf(A,B,valuesForEachBin','edgecolor','none')
xlabel('x [mm]');
ylabel('y [mm]');
colorbar;
view(0,90)
end

174
demos/octave/FlockingAnalysis/runAnalysis.m Normal file
View File

@ -0,0 +1,174 @@
function runAnalysis(properties)
% This high-level function generates the local transfer entropy results, first optimising
% parameters (i.e. embedding length and delay, and source-target lag), then
% storing local transfer entropy values for the optimised parameters.
%
% Author: Joseph T. Lizier, 2019
%
% Inputs:
% - properties (required) - object with properties for the calculations,
% with sub-members as specificied in the loadProperties.m file. If not supplied
% the properties are attempted to be loaded from loadProperties.m
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% This program is free software: you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation, either version 3 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
%%
if (nargin < 1)
fprintf('No properties object supplied, attempting to load properties via a loadProperties script ...');
% By default, just try to load properties locally
if (exist('loadProperties') == 2)
% there is a loadProperties script so attempt to run it to load a properties object
loadProperties;
else
% there is not a loadProperties script
error('No properties object supplied, and no loadProperties script found.');
end
else
% A properties argument was supplied
if (ischar(properties))
% We're assuming it was the name of a properties file
if ((length(properties) > 2) && (strcmp(properties(end-1:end), '.m')))
% Remove the '.m':
properties(end-1:end) = [];
end
if (exist(properties) == 2)
% attempt to run the properties .m file: (after making sure the properties variable is cleared; not necessary but is clean)
propertiesFile = properties;
clear properties;
eval(propertiesFile);
else
error('%s is not an .m file we can find that can be used to load a properties object', properties);
end
% else
% We assume it was the properties object.
end
end
% Step 1: Auto-embed if required:
if (isfield(properties, 'kRange') || isfield(properties, 'tauRange'))
% If any one of these two ranges weren't supplied, set the range variables
% to the value of the corresponding non-range variable:
if (~isfield(properties, 'kRange'))
properties.kRange = properties.k;
end
if (~isfield(properties, 'tauRange'))
properties.tauRange = properties.tau;
end
% Also set the lag to 1 as a dummy if we are optimising over that later as well:
if (~isfield(properties, 'lag'))
properties.lag = 1;
end
% Ask generateObservations to only return the target samples for the AIS calculation
properties.destSamplesOnly = true;
% Optimise k and tau
maxAIS = -inf;
maxAISk = properties.kRange(1);
maxAIStau = properties.tauRange(1);
for k = properties.kRange
properties.k = k;
minTau = min(properties.tauRange);
for tau = properties.tauRange
if ((k == 1) && (tau > minTau))
% We only need compute k=1 for a single tau
continue;
end
properties.tau = tau;
% Generate the observations for k,tau:
[D, Dpast] = generateObservations(properties);
if (isempty(D))
% There were no samples found for the given parameters, presumably k etc are too long
continue;
end
% Compute the AIS:
ais = computeAIS(D, Dpast, properties);
if (ais > maxAIS)
maxAIS = ais;
maxAISk = k;
maxAIStau = tau;
end
end
end
% Optimisation is complete:
properties.k = maxAISk;
properties.tau = maxAIStau;
properties.ais = maxAIS;
properties.destSamplesOnly = false;
fprintf('*** Optmised k=%d and tau=%d (giving AIS=%.4f)\n', properties.k, ...
properties.tau, properties.ais);
else
fprintf('*** Hard-coded values for k=%d and tau=%d to be used\n', properties.k, ...
properties.tau);
end
% Step 2: automatically select the correct lag if required:
teNumSurrogates = properties.teNumSurrogates; % Store this for later, turn it off now
properties.teNumSurrogates = 0; % No need to run any surrogates during parameter fitting
if (isfield(properties, 'lagRange'))
% Caller asks us to maximise the TE over a given range:
maxTE = -inf;
maxTElag = properties.lagRange(1);
for lag = properties.lagRange
properties.lag = lag;
% Generate the observations for k,tau,lag:
[D, Dpast, S, ~, maxSourceSamplesForATarget] = generateObservations(properties);
if (isempty(S))
% There were no samples found for the given parameters, presumably k etc are too long
continue;
end
% Check if we're turning on dynamic correlation exclusion:
if (isfield(properties.jidt, 'autoDynamicCorrelationExclusion'))
properties.jidt.dynamicCorrelationExclusion = maxSourceSamplesForATarget;
end
% Compute the TE:
te = computeTE(S, D, Dpast, properties);
if (te > maxTE)
maxTE = te;
maxTElag = lag;
end
end
% Optimisation is complete:
properties.lag = maxTElag;
properties.tranEntropy = maxTE;
fprintf('*** Optmised lag=%d (giving TE=%.4f)\n', properties.lag, ...
properties.tranEntropy);
else
fprintf('*** Hard-coded value for lag=%d to be used\n', properties.lag);
end
% 3. Compute TE with the correct parameters
% Now, once again pre-process the positional data into velocities, this time
% saving them into the results file (by not requesting [S,D,Dpast] outputs):
generateObservations(properties);
% And load these samples in from the saved file:
load(properties.resultsFile);
properties.teNumSurrogates = teNumSurrogates; % Allow surrogates to be computed for this final run with correct parameters
% And compute the TE again for the optimal parameters, this time
% saving the files:
% Turn on dynamic correlation exclusion if required:
if (isfield(properties.jidt, 'autoDynamicCorrelationExclusion'))
if (properties.jidt.autoDynamicCorrelationExclusion)
properties.jidt.dynamicCorrelationExclusion = maxSourceSamplesForATarget; % maxSourceSamplesForATarget was loaded from the results file
else
properties.jidt.dynamicCorrelationExclusion = 0; % no dynamic correlation exclusion
end
end
% Compute TE with no output arguments so that results are saved
computeTE(S, D, Dpast, properties);
end