add intermittent_correlation

Project.toml

@@ -11,6 +11,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 EasyFit = "fde71243-0cda-4261-b7c7-4845bd106b21"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 PDBTools = "e29189f1-7114-4dbd-93d0-c5673a921a58"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
@@ -39,6 +40,7 @@ Documenter = "1.3"
 EasyFit = "0.6"
 LaTeXStrings = "1.3"
 LinearAlgebra = "1.9"
+OffsetArrays = "1.13"
 PDBTools = "1.1"
 Plots = "1.39"
 Printf = "1.9"

docs/src/functions.md

@@ -6,7 +6,11 @@ CollapsedDocStrings = true
 
 ```@autodocs
 Modules = [ MolSimToolkit ]
-Pages = [ "distances.jl", "center_of_mass.jl" ]
+Pages = [ 
+    "intermittent_correlation.jl",
+    "distances.jl", 
+    "center_of_mass.jl" 
+]
 ```
 
 

src/MolSimToolkit.jl

@@ -6,6 +6,7 @@ using AtomsBase
 using StaticArrays
 import Chemfiles
 import PDBTools
+import OffsetArrays
 using LinearAlgebra: norm
 using Reexport: @reexport
 using ProgressMeter: Progress, next!
@@ -15,6 +16,7 @@ import LaTeXStrings # only because Aqua complains: used in the Plotting extensio
 export wrap, wrap_to_first
 export distances
 export align, align!, rmsd, rmsd_matrix
+export intermittent_correlation
 
 # Testing module
 include("./Testing.jl")
@@ -25,9 +27,10 @@ include("./datastructures/Positions.jl")
 
 # Basic functions
 include("./wrap.jl")
-include("./distances.jl")
-include("./center_of_mass.jl")
+include("./simple_functions/distances.jl")
+include("./simple_functions/center_of_mass.jl")
 include("./procrustes.jl")
+include("./simple_functions/intermittent_correlation.jl")
 
 # Analysis functions and modules
 include("./BlockAverages.jl")

src/simple_functions/intermittent_correlation.jl

@@ -0,0 +1,59 @@
+"""
+    intermittent_correlation(data::AbstractVector; maxdelta = length(data) ÷ 10)
+
+Calculate the intermittent correlation function of a time series. That is,
+computes the probability of finding a value of the same type at a distance
+`delta` in the time series, given that it was present in step `i`.
+
+Returns an `OffsetArray` with indices `0:maxdelta`, where the value at position
+`0` is `1.0`, corresponding to the normalized count of events. 
+
+# Arguments
+
+- `data::AbstractVector`: The time series to be analyzed. 
+- `maxdelta::Int`: The maximum delta-step to be considered. Defaults to 
+  `length(data) ÷ 10`.
+
+# Examples
+
+Here we produce a time-series of 10,000 elements, as a sequence of 
+1's and 0's (`[1, 0, 1, 0, ...]`), and calculate the intermittent correlation function.
+The probability of finding the same number (0 or 1) after odd steps is 0, and
+the probability of finding the same number after even steps is 1.
+
+```jldoctest
+julia> using MolSimToolkit
+
+julia> data = [ mod(i,2) for i in 1:10^4 ];
+
+julia> intermittent_correlation(data; maxdelta=4)
+5-element OffsetArray(::Vector{Float64}, 0:4) with eltype Float64 with indices 0:4:
+ 1.0
+ 0.0
+ 1.0
+ 0.0
+ 1.0
+```
+
+"""
+function intermittent_correlation(
+    data::AbstractVector; 
+    maxdelta::Integer = min(1, length(data) ÷ 10)
+)
+    types = unique(data)
+    counts = OffsetArrays.OffsetArray(zeros(maxdelta+1), 0:maxdelta)
+    for type in types
+        positions = findall(x -> isequal(x, type), data)
+        np = length(positions)
+        for i in 1:np, j in i:np
+            delta = positions[j] - positions[i]
+            if delta <= maxdelta
+                counts[delta] += 1
+            end
+        end
+    end
+    for i in 0:maxdelta
+        counts[i] /= length(data) - i
+    end
+    return counts
+end