added option "gs" for plotting, adding a gradient of colour for "3dgeo" and...

added option "gs" for plotting, adding a gradient of colour for "3dgeo" and "3d" option in plotting, added arithmetic and geometric smoothing

src/ABMEv.jl

@@ -21,5 +21,5 @@ module ABMEv
     export copy,runWorld_store_WF,runWorld_store_G,clean_world #,runWorld_G!,runWorld_WF!,
     export H_discrete,findclusters,var,covgeo,hamming,get_beta_div, get_alpha_div
     export update_afterbirth_std!,update_afterdeath_std!
-    export generalised_gaussian,gaussian,ma,eth_grad_std
+    export generalised_gaussian,gaussian,ma,geomsmooth,arithsmooth,eth_grad_std
 end

src/ABMEv_plot.jl

@@ -7,7 +7,7 @@ import KernelDensity:kde,pdf
 # ARGS
 - `what = ["x","H"]`: the plots you want to obtain
 - `trait = 1`: the trait that will plotted regarding what you asked. `trait = 0` will plot the geotrait
-- `tplot = false` used when calling xs, as it plots a snapshot of the world at a particular time
+- `tplot = 0` used when calling xs, as it plots a snapshot of the world at a particular time
 It should correspond to an integer, as it indexes the column to plot
 
 # Options available
@@ -89,12 +89,51 @@ It should correspond to an integer, as it indexes the column to plot
             x1_array[:],xt_array[:]
         end
     end
+    if "gs" in what
+        _world = clean_world(world[:, tplot > 0 ? tplot : length(p["tspan"]) ])
+        y = get_x(_world,tend,2)[:]
+        x = get_x(_world,tend,0)[:]
+        X = hcat(x,y)
+        d = kde(X)
+        # by density
+        d_i = diag(pdf(d,X[:,1],X[:,2]))
+        # by value
+        # d_i = y
+        d_i = (d_i .- minimum(d_i)) ./ (maximum(d_i) .- minimum(d_i))
+        # TODO: we stopped here
+        @series begin
+            seriestype := :scatter
+            markercolor := eth_grad_small[d_i]
+            # markercolor := :blue
+            markerstrokewidth := 0
+            # seriesalpha := 1.
+            xaxis := "geotrait"
+            yaxis := "trait value"
+            label := ""
+            grid := false
+            # marker := (:rect,20,1.)
+            x,y
+        end
+    end
     if "3dgeo" in what
+        d_i = []
+        for i in 1:size(world,2)
+            _world = clean_world(world[:,i])
+            x = get_x(_world,p["tspan"][i],2)[:]
+            y = get_x(_world,p["tspan"][i],0)[:]
+            X = hcat(x,y)
+            # d = kde(X)
+            # di_temp = diag(pdf(d,X[:,1],X[:,2]))
+            di_temp = y
+            di_temp = (di_temp .- minimum(di_temp)) ./ (maximum(di_temp) .- minimum(di_temp))
+            # here we normalise with respect to maximum value at each time step
+            append!(d_i,di_temp)
+        end
         @series begin
         xarray = get_geo.(world_sm,tspan_ar)
         yarray = get_x(world_sm,2)
             seriestype := :scatter3d
-            markercolor := "blue"
+            markercolor := eth_grad_std[d_i ./ 1.]
             markerstrokewidth := 0
             seriesalpha :=.1
             xlabel := "time"
@@ -106,11 +145,21 @@ It should correspond to an integer, as it indexes the column to plot
         end
     end
     if "3d" in what
+        d_i = []
+        for i in 1:size(world,2)
+            x = get_x(clean_world(world[:,i]),p["tspan"][i],1)[:]
+            y = get_x(clean_world(world[:,i]),p["tspan"][i],2)[:]
+            X = hcat(x,y)
+            d = kde(X)
+            di_temp = diag(pdf(d,X[:,1],X[:,2]))
+            di_temp = (di_temp .- minimum(di_temp)) ./ (maximum(di_temp) .- minimum(di_temp))
+            append!(d_i,di_temp)
+        end
         @series begin
         xarray = get_x(world_sm,1)
         yarray = get_x(world_sm,2)
             seriestype := :scatter3d
-            markercolor := "blue"
+            markercolor := eth_grad_small[d_i ./ maximum(d_i)]
             markerstrokewidth := 0
             seriesalpha :=.1
             xlabel := "time"

src/ABMEv_utils.jl

@@ -21,6 +21,26 @@ function ma(x::Array{T},f) where T <: Number
     return conv(x,ones(f)./f)[_s:_s+_N-1]
 end
 
+"""
+    function geomsmooth(x,smooth)
+Geometric smoothing, cf `https://en.wikipedia.org/wiki/Exponential_smoothing`
+
+"""
+function geomsmooth(x,smooth)
+    return [prod(x[i-smooth + 1:i])^(1/smooth) for i in smooth:length(x)]
+end
+
+"""
+    function arithsmooth(x,smooth)
+arithmetic smoothing
+
+"""
+function arithsmooth(x,smooth)
+    return [sum(x[i-smooth+1:i])/smooth for i in smooth:length(x)]
+end
+
 import Plots:cgrad
+# asymmetry towards red, blue is only a fifth of the color range
 const eth_grad_small = cgrad([colorant"#1F407A", RGB(0.671,0.851,0.914),RGB(1.0,1.0,0.749), RGB(0.992,0.682,0.38),RGB(0.647,0.0,0.149),],[.0,.1])
+# symmetry between red and blue
 const eth_grad_std = cgrad([colorant"#1F407A", RGB(0.671,0.851,0.914),RGB(1.0,1.0,0.749), RGB(0.992,0.682,0.38),RGB(0.647,0.0,0.149),],[.0,1.])