C++ Class and example
Listing of PhaseFunction.h
1 //--------------------------------------------------------------
2 // PhaseFunction.h
3 //
4 // This is C++ header file, and along with PhaseFunction.cc,
5 // provides a class definition and implementation for computing
6 // phase function as reported in
7 // A Phase Function to Quantify Serial Dependence between Discrete Samples.
8 // Ramana Dodla and Charles J. Wilson.
9 // Biophys J. 2010 February 17; 98(4): L5-L7.
10 // http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2820633/
11 //
12 // Author: Ramana Dodla, 2011.
13 //
14 // Typical calling sequence:
15 // InputSequenceFirst(A); // A is a vector of spike times or discrete events.
16 // InputSequenceSecond(B); // B is the same or another vector of spike times or discrete events.
17 // ComputePhaseFunction(lag_start, lag_end, lag_step);
18 // Then the output is available as vectors in:
19 // Lag and Psi
20 //
21 // At the cost by doing a copying operation for each call,
22 // you could repeately call the following function at single lag values
23 // and get phase function as a return value:
24 // ComputePhaseFunction(lag);
25 //
26 // See the example program listing later in this file.
27 //--------------------------------------------------------------
28 #ifndef PHASEFUNCTION_H_
29 #define PHASEFUNCTION_H_
30 #include <iostream> // cout
31 #include <cmath> // sqrt
32 #include <vector>
33 #include <algorithm> // min
34 #include <numeric> // accumulate
35
36 class PhaseFunction
37 {
38 public:
39 PhaseFunction();
40 ~PhaseFunction();
41 void InputSequenceFirst(const std::vector<double>& A); // input: vector A
42 void InputSequenceSecond(const std::vector<double>& B); // input: vector B
43 double ComputePhaseFunction(double lag); // for computing phase function at one lag value.
44 void ComputePhaseFunction(double lag_start, double lag_end, double lag_step); //... at many values.
45 std::vector<double> Lag; //output: Lag[i] is the lag
46 std::vector<double> Psi; //output: Psi[i] is phase function at Lag[i]
47
48 private:
49 std::vector<double> _A, _B, _Bcopy;
50 std::vector<double> _gamma, _delta, _c;
51 void _ComputePhasePairs();
52 void _ComputePhaseDrifts();
53 double _ComputeAvePhaseDrift();
54 void _err_A();
55 void _err_B();
56 void _err_lag();
57 };
58 #endif // PHASEFUNCTION_H_
59
60 /*
61 // ....example program:
62
63 #include <iostream>
64 #include <vector>
65 #include "PhaseFunction.h"
66
67 int main(int argc, char** argv) {
68
69 double X[] = { 0, 0.18595, 0.33955, 0.4997, 0.61495, 0.7413, 0.9024, 1.08125, 1.2001, 1.37285,
70 1.5084, 1.65105, 1.7984, 1.9283, 2.05115, 2.1797, 2.33295, 2.46965, 2.6048, 2.7382,
71 2.86875, 3.04065, 3.1876, 3.327, 3.4601, 3.6072, 3.76645, 3.95745, 4.11405, 4.2476,
72 4.42185, 4.60135, 4.737, 4.85865, 4.9919, 5.14455, 5.26175, 5.3947, 5.53255, 5.67475,
73 5.8135, 5.9587, 6.1078, 6.24605, 6.36345, 6.5113, 6.65755, 6.8065, 6.9528, 7.12085 };
74
75 // ....Get vector A:
76 // std::vector<double> A(X, X+sizeof(X)/sizeof(double)); // Take all of X
77 unsigned int n=7; // or choose a few data points from X
78 std::vector<double> A(X, X+n);
79
80 // ....Get vector B:
81 // ....none
82
83 // ....Give inputs:
84 PhaseFunction PF;
85 PF.InputSequenceFirst(A);
86 PF.InputSequenceSecond(A); // ....seeking auto-phase function.
87
88 // ....Compute Phase Function:
89 PF.ComputePhaseFunction(0.0, 0.5, 0.001); // ....from, to, step
90 for(unsigned int i=0; i<PF.Lag.size(); i++) {
91 std::cout << PF.Lag.at(i) << " " << PF.Psi.at(i) << std::endl;
92 }
93
94 // ....Compute Phase Function (second way):
95 //for(int i=0; i<50; i++) {
96 // std::cout << i*0.01 << " " << PF.ComputePhaseFunction(i*0.01) << std::endl;
97 //}
98
99 return 0;
100 }
101
102 */
Listing of PhaseFunction.cc
1 //--------------------------------------------------------------
2 // PhaseFunction.cc
3 //
4 // This is C++ class implementation file, and along with PhaseFunction.h,
5 // provides a class definition and implementation for computing
6 // phase function as reported in
7 // A Phase Function to Quantify Serial Dependence between Discrete Samples.
8 // Ramana Dodla and Charles J. Wilson.
9 // Biophys J. 2010 February 17; 98(4): L5-L7.
10 // http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2820633/
11 //
12 // Author: Ramana Dodla, 2011.
13 //
14 // Typical calling sequence:
15 // InputSequenceFirst(A); // A is a vector of spike times or discrete events.
16 // InputSequenceSecond(B); // B is the same or another vector of spike times or discrete events.
17 // ComputePhaseFunction(lag_start, lag_end, lag_step);
18 // Then the output is available as vectors in:
19 // Lag and Psi
20 //
21 // At the cost by doing a copying operation for each call,
22 // you could repeately call the following function at single lag values
23 // and get phase function as a return value:
24 // ComputePhaseFunction(lag);
25 //
26 // See the example program listing at the end of PhaseFunction.h file.
27 //--------------------------------------------------------------
28 #include "PhaseFunction.h"
29 //
30 // InputSequenceSecond() is the constructor. No arguments allowed.
31 //
32 PhaseFunction::PhaseFunction(){}
33 //
34 // ~InputSequenceSecond() is the destructor.
35 //
36 PhaseFunction::~PhaseFunction(){}
37 //
38 // InputSequenceSecond(vector<double>&) takes in a double vector A.
39 //
40 void PhaseFunction::InputSequenceFirst(const std::vector<double>& A) {
41 _A = A;
42 _err_A();
43 }
44 //
45 // InputSequenceSecond(vector<double>&) takes in a double vector B.
46 //
47 void PhaseFunction::InputSequenceSecond(const std::vector<double>& B) {
48 _B = B;
49 _Bcopy = B;
50 _err_B();
51 }
52 //
53 // ComputePhaseFunction(lag) computes the phase function at a single value of
54 // the lag (given as argument) and returns it as a double.
55 // The vectors Lag and Psi are also accessible as member variables, but will
56 // be of length 1.
57 //
58 double PhaseFunction::ComputePhaseFunction(double lag) {
59 Lag.clear();
60 Psi.clear();
61 for(size_t i = 0; i < _B.size(); i++){
62 _B[i] = _Bcopy[i] + lag;
63 }
64 Lag.push_back(lag);
65 Psi.push_back(_ComputeAvePhaseDrift());
66 return Psi.at(0);
67 }
68 //
69 // ComputePhaseFunction(lag_start, lag_end, lag_step) computes the phase function
70 // for a range of values starting from lag_start in steps of lag_step up to lag_end,
71 // and stores the lag and phase function values in the vectors: Lag, Psi.
72 // You can access them as member variables of the class instance.
73 //
74 void PhaseFunction::ComputePhaseFunction(double lag_start, double lag_end, double lag_step) {
75 double lag_diff = lag_end-lag_start;
76 size_t nsteps;
77 if(lag_step == 0.0) nsteps = 0;
78 else {
79 double d_nsteps = lag_diff / lag_step;
80 if(d_nsteps < 0) _err_lag();
81 nsteps = size_t (d_nsteps);
82 }
83 Lag.clear();
84 Psi.clear();
85 size_t ix = 0;
86 do {
87 double lag_now = lag_start + (ix++) * lag_step;
88 for(size_t i = 0; i < _B.size(); i++){
89 _B[i] = _Bcopy[i] + lag_now;
90 }
91 Lag.push_back(lag_now);
92 Psi.push_back(_ComputeAvePhaseDrift());
93 } while(ix <= nsteps);
94 }
95 //
96 // _ComputeAvePhaseDrift() computes the phase function
97 // and returns it as a double.
98 // This function makes use of the struct Accum defined below
99 // in computing the sum of absolute values
100 //
101 struct Accum { // This struct is used in _ComputeAvePhaseDrift() below.
102 double operator()(double d1, double d2){return std::fabs(d1) + std::fabs(d2);}
103 };
104 double PhaseFunction::_ComputeAvePhaseDrift() {
105 _ComputePhaseDrifts();
106 double thetadrift = std::accumulate(_c.begin(), _c.end(), 0.0, Accum()) / _c.size();
107 return thetadrift / std::sqrt(2.0);
108 }
109 //
110 // _ComputePhaseDrifts() computes the phase drifts c_i
111 // and stores them in the vector: _c
112 //
113 void PhaseFunction::_ComputePhaseDrifts() {
114 _c.clear();
115 _ComputePhasePairs();
116 double mc = std::accumulate(_gamma.begin(), _gamma.end(), 0.0) / _gamma.size();
117 double md = std::accumulate(_delta.begin(), _delta.end(), 0.0) / _delta.size();
118 double f1 = mc * mc + md * md;
119 double f2 = std::sqrt(f1);
120 for(size_t i = 0; i < _gamma.size(); i++) {
121 _c.push_back((f1-_gamma[i] * mc-_delta[i] * md) / f2);
122 }
123 }
124 //
125 // _ComputePhasePairs() computes the phase pairs (gamma, delta)
126 // and stores them in the two vectors: _gamma, _delta
127 //
128 void PhaseFunction::_ComputePhasePairs() {
129 size_t i = 1, j = 1;
130 double T=std::max(_A[0], _B[0]);
131 _gamma.clear();
132 _delta.clear();
133 while(i < _A.size() && j < _B.size()) {
134 if(_A[i] >= _B[j]) {
135 if(_A[i-1] >= _B[j]) {
136 T = _A[i-1];
137 j = j + 1;
138 }
139 else {
140 _gamma.push_back((_B[j] - T) / (_A[i] - _A[i-1]));
141 _delta.push_back((_B[j] - T) / (_B[j] - _B[j-1]));
142 T = _B[j];
143 j = j + 1;
144 }
145 }
146 else {
147 if(_B[j-1] >= _A[i]) {
148 T = _B[j-1];
149 i = i + 1;
150 }
151 else {
152 _gamma.push_back((_A[i] - T) / (_A[i] - _A[i-1]));
153 _delta.push_back((_A[i] - T) / (_B[j] - _B[j-1]));
154 T = _A[i];
155 i = i + 1;
156 }
157 }
158 }
159 }
160 //
161 // _err_A() is a function that emits a message if vector A's length is < 3
162 //
163 void PhaseFunction::_err_A(){
164 if(_A.size() < 3) {
165 std::cout << "ERROR: size of A vector is less than 3." << std::endl;
166 return;
167 }
168 }
169 //
170 // _err_B() is a function that emits a message if vector B's length is < 3
171 //
172 void PhaseFunction::_err_B(){
173 if(_B.size() < 3) {
174 std::cout << "ERROR: size of B vector is less than 3." << std::endl;
175 return;
176 }
177 }
178 //
179 // _err_lag() is a function that emits a message if you can't reach the final lag.
180 //
181 void PhaseFunction::_err_lag(){
182 std::cout << "ERROR: final lag cannot be reached." << std::endl;
183 return;
184 }
185 //========================= END =========================================
Listing of example_program.cc
1 //
2 // example_program.cc
3 // Computes auto-phase function for spike time sequence of rat globus pallidus in vitro recordings
4 //
5 #include <iostream>
6 #include <vector>
7 #include "PhaseFunction.h"
8
9 int main(int argc, char** argv) {
10
11 double X[] = { 0, 0.18595, 0.33955, 0.4997, 0.61495, 0.7413, 0.9024, 1.08125, 1.2001, 1.37285,
12 1.5084, 1.65105, 1.7984, 1.9283, 2.05115, 2.1797, 2.33295, 2.46965, 2.6048, 2.7382,
13 2.86875, 3.04065, 3.1876, 3.327, 3.4601, 3.6072, 3.76645, 3.95745, 4.11405, 4.2476,
14 4.42185, 4.60135, 4.737, 4.85865, 4.9919, 5.14455, 5.26175, 5.3947, 5.53255, 5.67475,
15 5.8135, 5.9587, 6.1078, 6.24605, 6.36345, 6.5113, 6.65755, 6.8065, 6.9528, 7.12085 };
16
17 // ....Get vector A:
18 // std::vector<double> A(X, X + sizeof(X) / sizeof(double)); // Take all of X
19 unsigned int n = 7; // ....or choose a few data points from X
20 std::vector<double> A(X, X+n);
21
22 // ....Get vector B:
23 // ....none
24
25 // ....Give inputs:
26 PhaseFunction PF;
27 PF.InputSequenceFirst(A);
28 PF.InputSequenceSecond(A); // ....seeking auto-phase function.
29
30 // ....Compute Phase Function:
31 PF.ComputePhaseFunction(0.0, 0.5, 0.001); // ....from, to, step
32 for(unsigned int i = 0; i < PF.Lag.size(); i++) {
33 std::cout << PF.Lag.at(i) << " " << PF.Psi.at(i) << std::endl;
34 }
35
36 // ....Compute Phase Function (second way):
37 //for(int i = 0; i < 50; i++) {
38 // std::cout << i * 0.01 << " " << PF.ComputePhaseFunction(i * 0.01) << std::endl;
39 //}
40
41 return 0;
42 }