1 """Compiled routines for HMM simulation.
2
3 Copyright 2008-2014 Research Foundation State University of New York
4 This file is part of QUB Express.
5
6 QUB Express is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 QUB Express is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License,
17 named LICENSE.txt, in the QUB Express program directory. If not, see
18 <http://www.gnu.org/licenses/>.
19
20 """
21
22 from qubx.fast.fast_utils import *
23 from qubx.fast.model import *
24 from itertools import izip
25 import random
26
27
28 qubfast.qubx_simulate.argtypes = (PModel, c_int_p, c_int, c_int, c_int, c_double, c_int, c_int, c_int_pp, c_float_ppp,
29 c_int_p, c_float_pp, c_int_p, c_int_pp, c_int_pp, c_double_pp, c_double_pp,
30 c_int, StatusCallback, c_int_p)
31 qubfast.qubx_simulate.restype = c_int
32
33
34
35
36
37
38 -def simulate(fastmodel, start_states, add_baseline, use_Peq, v_signal, sampling, stimulus_of_seg, Nsample_of_seg, samples_of_seg,
39 on_pct = lambda frac: 0, force_Q=False):
40 """Returns independent segments of simulated HMM, as sampled data and as state-idealization.
41 @param fastmodel: qubx.fast.Model
42 @param start_states: numpy.array[Nchannel], type int32, of entry state indices, or -1 to pick fresh. On output, contains exit states for sequential simulations.
43 @param add_baseline: False to skip base amp and std, for adding multi channels
44 @param v_signal: if nonzero, multiply exc_amp and exc_std by (stimulus[v] - vRev) * 1e-3 (in volts)
45 @param sampling: interval between samples, in seconds
46 @param stimulus_of_seg: list[Nseg][Nsignal] of stimulus values or sample arrays. Nsignal > all x in union(L, V). Since 0 indicates 'no stimulus,' stimulus[0] is ignored. Entries should be either a constant (float) or a list or array[nsample] containing a stimulus waveform.
47 @param nsample: list[Nseg] of number of data points to return
48 @param samples_of_seg: list[Nseg] of array[Nsample_of_seg[i]] of float32, will be filled with output
49 @param on_pct: callback function on_pct(frac), 0 <= frac <= 1. Return nonzero to stop the simulation.
50 @param force_Q: uses Q matrix method (instead of A matrix) even if it would be slower
51 @return: (states as numpy.array[Nevent] of int32, counts as numpy.array[Nevent] of int32, amp, std as numpy.array[Nevent] of float64)
52 """
53 def cdata(x, typ):
54 if x is None: return x
55 return x.ctypes.data_as(typ)
56 temps = []
57 def temparray(shape, dtype, init):
58 a = numpy.zeros(shape=shape, dtype=dtype)
59 a += init
60 temps.append(a)
61 return a
62 Nseg = len(stimulus_of_seg)
63 result = []
64 states = []
65 counts = []
66 amps = []
67 stds = []
68 stim_counts = []
69 stimuli = []
70 samples = []
71 for stimulus, samples1, Nsample in izip(stimulus_of_seg, samples_of_seg, Nsample_of_seg):
72 stim_counts_ = temparray(shape=(len(stimulus),), dtype='int32', init=0)
73 stimuli1 = []
74 for i,stim in enumerate(stimulus):
75 if hasattr(stim,'__iter__'):
76 stim_counts_[i] = len(stim)
77 stimuli1.append(cdata(stim,c_float_p))
78 else:
79 stimuli1.append(cdata(temparray(shape=(1,),dtype='float32', init=stim),c_float_p))
80 stim_counts.append(cdata(stim_counts_, c_int_p))
81 stimuli.append((c_float_p*len(stimuli1))(*stimuli1))
82 samples.append(cdata(samples1, c_float_p))
83 tt = temparray(shape=(Nsample,), dtype='int32', init=0)
84 cc = temparray(shape=(Nsample,), dtype='int32', init=0)
85 aa = temparray(shape=(fastmodel.Nstate,), dtype='float64', init=0.0)
86 ss = temparray(shape=(fastmodel.Nstate,), dtype='float64', init=0.0)
87 states.append(cdata(tt, c_int_p))
88 counts.append(cdata(cc, c_int_p))
89 amps.append(cdata(aa, c_double_p))
90 stds.append(cdata(ss, c_double_p))
91 result.append([tt, cc, aa, ss])
92 stim_counts_ = (c_int_p*len(stim_counts))(*stim_counts)
93 stimuli_ = (c_float_pp*len(stimuli))(*stimuli)
94 samples_ = (c_float_p*len(samples))(*samples)
95 states_ = (c_int_p*len(states))(*states)
96 counts_ = (c_int_p*len(counts))(*counts)
97 amps_ = (c_double_p*len(amps))(*amps)
98 stds_ = (c_double_p*len(stds))(*stds)
99 ndwt = temparray(shape=(Nseg,), dtype='int32', init=0)
100 Nsample_ = numpy.array(Nsample_of_seg, dtype='int32', copy=True)
101
102 on_status = StatusCallback(lambda ptr, frac: on_pct(frac))
103
104 if stimulus_of_seg:
105 qubfast.qubx_simulate(fastmodel.obj, cdata(start_states, c_int_p), int(add_baseline), int(use_Peq), int(v_signal),
106 sampling, len(stimulus_of_seg[0]), len(stimulus_of_seg), stim_counts_, stimuli_,
107 cdata(Nsample_, c_int_p), samples_, cdata(ndwt, c_int_p), states_, counts_, amps_, stds_,
108 force_Q, on_status, None)
109
110 for rr, nd in izip(result, ndwt):
111 tt, cc = rr[:2]
112 rr[:2] = tt[:nd], cc[:nd]
113 return result
114
115
116 qubfast.qubx_simulate_make_idealized_flc.argtypes = (c_int_p, c_int_p, c_int_p, c_int, c_int_p, c_int_p, c_int_p, c_int)
117 qubfast.qubx_simulate_make_idealized_flc.restype = c_int
118
120 return qubfast.qubx_simulate_make_idealized_flc(cdata(ff, c_int_p), cdata(ll, c_int_p), cdata(cc, c_int_p),
121 len(states), cdata(clazz, c_int_p), cdata(states, c_int_p), cdata(st_counts, c_int_p), f0)
122
123
124
125
126
127
128
129
130
131
132
133 qubfast.qubx_simulate_baselines.argtypes = (c_int, c_int_p, c_float_pp, c_double,
134 c_int, c_int, c_double, c_double,
135 c_int, c_double, c_double,
136 c_int, c_double, c_double, c_double,
137 c_int, c_double,
138 c_double, c_double, c_double, c_double,
139 c_int, c_int, c_double,
140 c_double, c_double)
141 qubfast.qubx_simulate_baselines.restype = None
142
143 -def simulate_baselines(sample_arrays, sampling,
144 baselineFluct, baselineFluctN, baselineFluctLifetime, baselineFluctMaxAmp,
145 baselineDrift, baselineDriftSlope, baselineDriftIntercept,
146 baselinePeriodic, baselinePeriodicFreq, baselinePeriodicPhase, baselinePeriodicAmp,
147 baselineJumps, baselineJumpsAmpStd,
148 baselineJumpsFastW, baselineJumpsFast, baselineJumpsSlow, baselineJumpsLifetime,
149 baselineLGM, baselineLGM2State, baselineLGMProcStd, baselineLGMX0, baselineLGMX0p):
150
151 def cdata(x, typ):
152 if x is None: return x
153 return x.ctypes.data_as(typ)
154 temps = []
155 def temparray(shape, dtype, init):
156 a = numpy.zeros(shape=shape, dtype=dtype)
157 a += init
158 temps.append(a)
159 return a
160
161 Nsamples = numpy.array([len(samples) for samples in sample_arrays], dtype='int32')
162 arrs = (c_float_p*len(sample_arrays))(*[cdata(samples, c_float_p) for samples in sample_arrays])
163
164 qubfast.qubx_simulate_baselines(len(sample_arrays), cdata(Nsamples, c_int_p), arrs, sampling,
165 baselineFluct, baselineFluctN, baselineFluctLifetime, baselineFluctMaxAmp,
166 baselineDrift, baselineDriftSlope, baselineDriftIntercept,
167 baselinePeriodic, baselinePeriodicFreq, baselinePeriodicPhase, baselinePeriodicAmp,
168 baselineJumps, baselineJumpsAmpStd,
169 baselineJumpsFastW, baselineJumpsFast, baselineJumpsSlow, baselineJumpsLifetime,
170 baselineLGM, baselineLGM2State, baselineLGMProcStd, baselineLGMX0, baselineLGMX0p)
171