项目作者: andreabonvini

项目描述 :
A C++ Library for Point Process Analysis
高级语言: C++
项目地址: git://github.com/andreabonvini/pointprocess.git
创建时间: 2021-05-10T16:21:46Z
项目社区:https://github.com/andreabonvini/pointprocess

开源协议:MIT License

下载


Documentation Status
Build Status
codecov.io

ppbig

This repository contains a C++ implementation of the MATLAB software provided by Riccardo Barbieri and Luca Citi here.

Refer to the following papers for more details:

Requirements

The project can be built with CMake, the only dependencies are Eigen
and Boost.

Check here for information on how to install these two packages.

Documentation

The technical and scientific documentation for this repository is available here. (WIP)

In order to build the project you should follow the instructions in the build.sh script.

Quick Tour

This brief tutorial shows how to fit a history-dependent Inverse Gaussian point process model to a temporal series of RR events.

You can find this example as a Jupyter Notebook in the examples/ directory.

Note: In order to call the Python bindings you need to first build the library.

  1. # Basic imports
  2. import sys
  3. import os
  4. import numpy as np
  5. import matplotlib.pyplot as plt
  6. # Path to the generated .so library (MacOs or Linux) or .dll (Windows)
  7. # Should be "build/src/" if you built the project with build.sh
  8. PATH_TO_LIB = "path/to/library/"
  9. sys.path.append(PATH_TO_LIB)
  10. # Now we can import the Python byndings from the pointprocess library
  11. from pointprocess import (
  12. compute_single_regression,
  13. compute_full_regression,
  14. compute_spectral_analysis,
  15. Distributions,
  16. Result,
  17. RegressionResult,
  18. get_ks_coords
  19. )
  1. # Example: load and plot a series of RR events
  2. rr = np.load("events.npy")
  3. events = np.array(rr[75:301])
  4. plt.plot(events[1:], 1000*np.diff(events),"b")
  5. plt.xlabel('time [s]')
  6. plt.ylabel('RR [ms]')

Single Regression

  1. # Compute single regression on a subset of the RR intervals
  2. result_single_regression = compute_single_regression(
  3. events = rr[75:301],
  4. ar_order = 9,
  5. has_theta0 = True,
  6. right_censoring = False,
  7. alpha = 0.02,
  8. distribution=Distributions.InverseGaussian,
  9. max_iter=10000
  10. )
  11. # You can now access to elements as:
  12. # result_single_regression.thetap
  13. # result_single_regression.theta0
  14. # result_single_regression.kappa
  15. # result_single_regression.likelihood
  16. # result_single_regression.mean_interval

Get Spectral info

  1. thetap = result_single_regression.thetap
  2. mean_interval = result_single_regression.mean_interval
  3. kappa = result_single_regression.kappa
  4. variance = result_single_regression.sigma**2
  5. # Retrieve spectral info
  6. analysis = compute_spectral_analysis(thetap, mean_interval, variance, aggregate=False)
  7. # Plot stuff
  8. plt.figure(figsize=(8,6),dpi=100)
  9. colors = []
  10. for comp in analysis.comps:
  11. p = plt.plot(analysis.frequencies, np.real(comp),linewidth=0.7)
  12. colors.append(p[-1].get_color())
  13. for i in range(len(analysis.poles)):
  14. plt.vlines(analysis.poles[i].frequency,0,10000,colors[i],"--")
  15. plt.plot(analysis.frequencies, analysis.powers, "k-",linewidth=0.8,label="Total PSD")
  16. plt.xlabel("f [Hz]")
  17. plt.ylabel("PSD [$ms^2$/Hz]")
  18. plt.legend()

Full Regression

  1. # Fit Inverse Gaussian distribution to RR by moving a 60.0 seconds windows and shifting it by 0.005 s at each step.
  2. # The mean of the ditribution will bi given by a 9th order AR model
  3. result = compute_full_regression(
  4. events=rr,
  5. window_length=60.0,
  6. delta=0.005,
  7. ar_order=9,
  8. has_theta0=True,
  9. right_censoring=True,
  10. alpha = 0.02,
  11. distribution=Distributions.InverseGaussian,
  12. max_iter = 1000
  13. )
  14. # Compute spectral info
  15. result.compute_hrv_indices()
  16. # Convert result to dictionary...
  17. d = result.to_dict()
  1. # Plot first moment of the distribution in time, along with the discrete RR intervals
  2. plt.figure(figsize=(6,4),dpi=120)
  3. plt.plot(d["Time"],d["Mu"],"b",linewidth=0.5,label="First moment of IG regression")
  4. plt.xlabel("Time [s]")
  5. plt.ylabel("$\mu$ [s]")
  6. plt.plot(rr[1:],np.diff(rr),"r*",mew=0.01,label="RR interval")
  7. plt.legend()

  1. # Plot hazard rate
  2. plt.figure(figsize=(6,4),dpi=120)
  3. plt.plot(d["Time"],d["lambda"],"firebrick",linewidth=0.15,label="Hazard rate ($\lambda$)")
  4. plt.xlabel("Time [s]")
  5. plt.ylabel("Hazard rate ($\lambda$)")
  6. plt.legend()
  7. # Same but just for a smaller number of samples...
  8. plt.figure(figsize=(6,4),dpi=120)
  9. s,e = 30200,31000
  10. plt.plot(d["Time"][s:e],d["lambda"][s:e],"firebrick",linewidth=0.5,label="Hazard rate ($\lambda$)")
  11. plt.xlabel("Time [s]")
  12. plt.ylabel("Hazard rate ($\lambda$)")
  13. plt.legend()

  1. # Plot LF/HF ration
  2. plt.figure(figsize=(8,8))
  3. bal = d["powLF"] / d["powHF"]
  4. plt.plot(d["Time"],bal,"blue",linewidth=0.5)
  5. plt.ylim(-1,20)
  6. plt.ylabel("LF/HF")
  7. plt.xlabel("Time [s]")

  1. # Build KS plot to assess goodness of fit
  2. coords = get_ks_coords(result.taus)
  3. plt.figure(figsize=(5,5),dpi=100)
  4. plt.plot(coords.z, coords.inner, "k", linewidth=0.8)
  5. plt.plot(coords.lower, coords.inner, "b", linewidth=0.5)
  6. plt.plot(coords.upper, coords.inner, "b", linewidth=0.5)
  7. plt.plot(coords.inner, coords.inner, "r", linewidth=0.5)