qubx.fast.data

1 """Compiled routines for sampled data. 2 3 Copyright 2008-2011 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 from qubx.fast.fast_utils import * 22 23 qubfast.qub_raster_samples.argtypes = (c_int, c_int, c_float_p, c_double, c_double, 24 c_float_p, c_float_p, c_float_p, c_float_p) 25 qubfast.qub_raster_samples.restype = c_int 26 27 #extern "C" QUBFAST_API int qub_raster_samples(int w, int Nsample, float *samples, double samples_per_pixel, double y_res, 28 # float *lows, float *highs, float *gauss_lows, float *gauss_highs) 29

30 -def raster_samples(w, samples, samples_per_pixel, resolution):

31 """ 32 Returns (lows, highs, gauss_lows, gauss_highs); arrays of at most w y-coords; 33 each pixel represents several samples; the arrays give the minimum, maximum, and mean +/- one std. 34 35 @param resolution: min difference between highs[i] and lows[i] 36 """ 37 def cdata(x, typ): 38 if x is None: return x 39 return x.ctypes.data_as(typ)

40 lows = numpy.zeros(shape=(w,), dtype='float32') 41 highs = numpy.zeros(shape=(w,), dtype='float32') 42 gauss_lows = numpy.zeros(shape=(w,), dtype='float32') 43 gauss_highs = numpy.zeros(shape=(w,), dtype='float32') 44 raster_count = qubfast.qub_raster_samples(w, len(samples), cdata(samples, c_float_p), samples_per_pixel, resolution, 45 cdata(lows, c_float_p), cdata(highs, c_float_p), 46 cdata(gauss_lows, c_float_p), cdata(gauss_highs, c_float_p)) 47 return lows[:raster_count], highs[:raster_count], gauss_lows[:raster_count], gauss_highs[:raster_count] 48 49 50 qubfast.qub_doseresponse_equil_init.argtypes = (c_int,) 51 qubfast.qub_doseresponse_equil_init.restype = c_int_p 52 qubfast.qub_doseresponse_equil_free.argtypes = (c_int_p,) 53 qubfast.qub_doseresponse_equil_free.restype = None 54 qubfast.qub_doseresponse_equil_add_data.argtypes = (c_int_p, c_int, c_float_p, c_float_p) 55 qubfast.qub_doseresponse_equil_add_data.restype = None 56 qubfast.qub_doseresponse_equil_result.argtypes = (c_int_p, c_int, c_float_p, c_float_p, c_float_p) 57 qubfast.qub_doseresponse_equil_result.restype = c_int 58

59 -class DoseResponse_Equil(object):

60 """Very important to feed data chunks in reverse order."""

61 - def __init__(self, Ntail):

62 self.dr = qubfast.qub_doseresponse_equil_init(Ntail)

63 - def __del__(self):

64 qubfast.qub_doseresponse_equil_free(self.dr)

65 - def add_data(self, stim, resp):

66 qubfast.qub_doseresponse_equil_add_data(self.dr, min(len(stim), len(resp)), stim.ctypes.data_as(c_float_p), resp.ctypes.data_as(c_float_p))

67 - def result(self, max_results=100):

68 dose = numpy.zeros(shape=(max_results,), dtype='float32') 69 resp = numpy.zeros(shape=(max_results,), dtype='float32') 70 resp_std = numpy.zeros(shape=(max_results,), dtype='float32') 71 Nresult = qubfast.qub_doseresponse_equil_result(self.dr, max_results, dose.ctypes.data_as(c_float_p), 72 resp.ctypes.data_as(c_float_p), resp_std.ctypes.data_as(c_float_p)) 73 return dose[:Nresult], resp[:Nresult], resp_std[:Nresult]

74 75 qubfast.qub_doseresponse_peak_init.argtypes = (c_double, c_double, c_double) 76 qubfast.qub_doseresponse_peak_init.restype = c_int_p 77 qubfast.qub_doseresponse_peak_free.argtypes = (c_int_p,) 78 qubfast.qub_doseresponse_peak_free.restype = None 79 qubfast.qub_doseresponse_peak_add_data.argtypes = (c_int_p, c_int, c_float_p, c_float_p) 80 qubfast.qub_doseresponse_peak_add_data.restype = None 81 qubfast.qub_doseresponse_peak_result.argtypes = (c_int_p, c_int, c_float_p, c_float_p) 82 qubfast.qub_doseresponse_peak_result.restype = c_int 83

84 -class DoseResponse_Peak(object):

85 """Feed whole segments please."""

86 - def __init__(self, sampling_sec, filter_kHz, sign):

87 self.dr = qubfast.qub_doseresponse_peak_init(sampling_sec, filter_kHz, sign)

88 - def __del__(self):

89 qubfast.qub_doseresponse_peak_free(self.dr)

90 - def add_data(self, stim, resp):

91 qubfast.qub_doseresponse_peak_add_data(self.dr, min(len(stim), len(resp)), stim.ctypes.data_as(c_float_p), resp.ctypes.data_as(c_float_p))

92 - def result(self, max_results=100):

93 dose = numpy.zeros(shape=(max_results,), dtype='float32') 94 resp = numpy.zeros(shape=(max_results,), dtype='float32') 95 Nresult = qubfast.qub_doseresponse_peak_result(self.dr, max_results, dose.ctypes.data_as(c_float_p), 96 resp.ctypes.data_as(c_float_p)) 97 return dose[:Nresult], resp[:Nresult]

98 99 100 qubfast.qub_idlstim_create.argtypes = (c_double_p, c_int, c_int, c_double, c_double, c_double) 101 qubfast.qub_idlstim_create.restype = c_int_p 102 qubfast.qub_idlstim_destroy.argtypes = (c_int_p,) 103 qubfast.qub_idlstim_destroy.restype = None 104 qubfast.qub_idlstim_add.argtypes = (c_int_p, c_float_p, c_int) 105 qubfast.qub_idlstim_add.restype = None 106 qubfast.qub_idlstim_done_add.argtypes = (c_int_p,) 107 qubfast.qub_idlstim_done_add.restype = None 108 qubfast.qub_idlstim_get_amp_count.argtypes = (c_int_p,) 109 qubfast.qub_idlstim_get_amp_count.restype = c_int 110 qubfast.qub_idlstim_get_amps.argtypes = (c_int_p, c_double_p) 111 qubfast.qub_idlstim_get_amps.restype = c_int 112 qubfast.qub_idlstim_get_next_dwells.argtypes = (c_int_p, c_int, c_int, c_int_p, c_int_p, c_int_p, c_float_p) 113 qubfast.qub_idlstim_get_next_dwells.restype = c_int 114

115 -class IdlStim(object):

116 """Idealizes a stimulus signal. 117 """

118 - def __init__(self, known_amps, add_deltas, delta, min_dur_ms, sampling_ms):

119 """ 120 @param known_amps: list of known stimulus levels 121 @param add_deltas: True if the idealizer can find new levels 122 @param delta: minimum difference between found levels 123 @param min_dur_ms: dwells shorter than this will be joined to a neighbor 124 @param sampling_ms: sampling interval in milliseconds 125 """ 126 amps = numpy.array(known_amps, dtype='float64') 127 self.idlstim = qubfast.qub_idlstim_create(amps.ctypes.data_as(c_double_p), len(amps), add_deltas, delta, min_dur_ms, sampling_ms)

128 - def __del__(self):

129 qubfast.qub_idlstim_destroy(self.idlstim)

130 - def add(self, samples):

131 """Feeds in a numpy.array(dtype=float32) of samples; you add all samples, then call done_add, then for each array of samples added, you call get_next_dwells.""" 132 qubfast.qub_idlstim_add(self.idlstim, samples.ctypes.data_as(c_float_p), len(samples))

133 - def done_add(self):

134 """Indicates you're done calling add(), so it can prepare the final amps and dwells.""" 135 qubfast.qub_idlstim_done_add(self.idlstim)

136 - def get_amps(self):

137 """After you've called done_add(), returns the amp of each class in the output.""" 138 amps = numpy.zeros(shape=(qubfast.qub_idlstim_get_amp_count(self.idlstim),), dtype='float64') 139 qubfast.qub_idlstim_get_amps(self.idlstim, amps.ctypes.data_as(c_double_p)) 140 return amps

141 - def get_next_dwells(self, offset_in_data, sample_count, get_durations=False):

142 """After you've called done_add(), call this once for each array of samples you added. 143 144 @param offset_in_data: the "first" sample index for the source sample array 145 @param sample_count: length of the source array 146 @param get_durations: whether to return an array of durations 147 @returns: (firsts, lasts, classes[, durations]) 148 """ 149 firsts = numpy.zeros(shape=(sample_count,), dtype='int32') 150 lasts = numpy.zeros(shape=(sample_count,), dtype='int32') 151 classes = numpy.zeros(shape=(sample_count,), dtype='int32') 152 if get_durations: 153 durations = numpy.zeros(shape=(sample_count,), dtype='float32') 154 dd = durations.ctypes.data_as(c_float_p) 155 else: 156 dd = None 157 ndwt = qubfast.qub_idlstim_get_next_dwells(self.idlstim, offset_in_data, sample_count, 158 firsts.ctypes.data_as(c_int_p), 159 lasts.ctypes.data_as(c_int_p), 160 classes.ctypes.data_as(c_int_p), 161 dd) 162 if get_durations: 163 return numpy.array(firsts[:ndwt]), numpy.array(lasts[:ndwt]), numpy.array(classes[:ndwt]), numpy.array(durations[:ndwt]) 164 else: 165 return numpy.array(firsts[:ndwt]), numpy.array(lasts[:ndwt]), numpy.array(classes[:ndwt])

166 167 168 169 qubfast.qub_adaptive_resample.argtypes = (c_float_p, c_int, c_double, c_float_p, c_float_p, c_int_p, c_int_p, c_int_p) 170 qubfast.qub_adaptive_resample.restype = c_int 171 172 # QUBFAST_API int qub_adaptive_resample(float *data, int ndata, double MaxStd, 173 # float *out_mean, float *out_std, int *out_first, int *out_last, 174 # int *out_closest); 175

176 -def adaptive_resample(samples, max_std):

177 """ 178 Divides samples into intervals such that all(stds[i] < max_std). Also finds the closest sample to the mean in each interval. 179 180 @param samples: numpy.array(float32) of sampled data 181 @param max_std: desired ceiling for all stds[i] 182 @returns: (means, stds, firsts, lasts, closests) -- array[interval_index] 183 """ 184 def cdata(x, typ): 185 if x is None: return x 186 return x.ctypes.data_as(typ)

187 ndata = samples.shape[0] 188 out_means = numpy.zeros(shape=(ndata,), dtype='float32') 189 out_stds = numpy.zeros(shape=(ndata,), dtype='float32') 190 out_firsts = numpy.zeros(shape=(ndata,), dtype='int32') 191 out_lasts = numpy.zeros(shape=(ndata,), dtype='int32') 192 out_closests = numpy.zeros(shape=(ndata,), dtype='int32') 193 out_count = qubfast.qub_adaptive_resample(cdata(samples, c_float_p), ndata, max_std, 194 cdata(out_means, c_float_p), cdata(out_stds, c_float_p), 195 cdata(out_firsts, c_int_p), cdata(out_lasts, c_int_p), 196 cdata(out_closests, c_int_p)) 197 return (numpy.array(out_means[:out_count]), 198 numpy.array(out_stds[:out_count]), 199 numpy.array(out_firsts[:out_count]), 200 numpy.array(out_lasts[:out_count]), 201 numpy.array(out_closests[:out_count])) 202 203 204 qubfast.qub_adres_create.argtypes = (c_double, c_double, c_double) 205 qubfast.qub_adres_create.restype = c_int_p 206 qubfast.qub_adres_destroy.argtypes = (c_int_p,) 207 qubfast.qub_adres_destroy.restype = None 208 qubfast.qub_adres_add.argtypes = (c_int_p, c_float_p, c_int) 209 qubfast.qub_adres_add.restype = None 210 qubfast.qub_adres_done_add.argtypes = (c_int_p,) 211 qubfast.qub_adres_done_add.restype = None 212 qubfast.qub_adres_get_class_count.argtypes = (c_int_p,) 213 qubfast.qub_adres_get_class_count.restype = c_int 214 qubfast.qub_adres_get_dist.argtypes = (c_int_p, c_double_p, c_double_p) 215 qubfast.qub_adres_get_dist.restype = c_int 216 qubfast.qub_adres_get_next_dwells.argtypes = (c_int_p, c_int, c_int, c_int_p, c_int_p, c_int_p, c_float_p) 217 qubfast.qub_adres_get_next_dwells.restype = c_int 218

219 -class Idl_AdaptiveResample(object):

220 """Idealizes a stimulus signal. 221 """

222 - def __init__(self, delta, min_dur_ms, sampling_ms):

223 """ 224 @param delta: minimum difference between found levels 225 @param min_dur_ms: dwells shorter than this will be joined to a neighbor 226 @param sampling_ms: sampling interval in milliseconds 227 """ 228 self.adres = qubfast.qub_adres_create(delta, min_dur_ms, sampling_ms)

229 - def __del__(self):

230 qubfast.qub_adres_destroy(self.adres)

231 - def add(self, samples):

232 """Feeds in a numpy.array(dtype=float32) of samples; you add all samples, then call done_add, then for each array of samples added, you call get_next_dwells.""" 233 qubfast.qub_adres_add(self.adres, samples.ctypes.data_as(c_float_p), len(samples))

234 - def done_add(self):

235 """Indicates you're done calling add(), so it can prepare the final amps and dwells.""" 236 qubfast.qub_adres_done_add(self.adres)

237 - def get_dist(self):

238 """After you've called done_add(), returns the (amps, stds).""" 239 n = qubfast.qub_adres_get_class_count(self.adres) 240 amps = numpy.zeros(shape=(n,), dtype='float64') 241 stds = numpy.zeros(shape=(n,), dtype='float64') 242 qubfast.qub_adres_get_dist(self.adres, amps.ctypes.data_as(c_double_p), stds.ctypes.data_as(c_double_p)) 243 return amps, stds

244 - def get_next_dwells(self, offset_in_data, sample_count, get_durations=False):

245 """After you've called done_add(), call this once for each array of samples you added. 246 247 @param offset_in_data: the "first" sample index for the source sample array 248 @param sample_count: length of the source array 249 @param get_durations: whether to return an array of durations 250 @returns: (firsts, lasts, classes[, durations]) 251 """ 252 firsts = numpy.zeros(shape=(sample_count,), dtype='int32') 253 lasts = numpy.zeros(shape=(sample_count,), dtype='int32') 254 classes = numpy.zeros(shape=(sample_count,), dtype='int32') 255 if get_durations: 256 durations = numpy.zeros(shape=(sample_count,), dtype='float32') 257 dd = durations.ctypes.data_as(c_float_p) 258 else: 259 dd = None 260 ndwt = qubfast.qub_adres_get_next_dwells(self.adres, offset_in_data, sample_count, 261 firsts.ctypes.data_as(c_int_p), 262 lasts.ctypes.data_as(c_int_p), 263 classes.ctypes.data_as(c_int_p), 264 dd) 265 if get_durations: 266 return numpy.array(firsts[:ndwt]), numpy.array(lasts[:ndwt]), numpy.array(classes[:ndwt]), numpy.array(durations[:ndwt]) 267 else: 268 return numpy.array(firsts[:ndwt]), numpy.array(lasts[:ndwt]), numpy.array(classes[:ndwt])

269 270 271 272 273 qubfast.decode_acquirefile_samples.argtypes = (c_char_p, c_short_p, c_int, c_int) 274 qubfast.decode_acquirefile_samples.restype = None 275

276 -def decode_acquirefile_samples(src_8, src_16, src_16s, dest, size, count):

277 qubfast.decode_acquirefile_samples(cdata(src_8, c_char_p), cdata(dest, c_short_p), size, count)

278 279 280 281 qubfast.qub_idlstats_create.argtypes = (c_double,) 282 qubfast.qub_idlstats_create.restype = c_void_p 283 qubfast.qub_idlstats_free.argtypes = (c_void_p,) 284 qubfast.qub_idlstats_free.restype = None 285 qubfast.qub_idlstats_reset.argtypes = (c_void_p,) 286 qubfast.qub_idlstats_reset.restype = None 287 qubfast.qub_idlstats_add.argtypes = (c_void_p, c_float_p, c_int, c_int, c_int_p, c_int_p, c_int_p) 288 qubfast.qub_idlstats_add.restype = None 289 qubfast.qub_idlstats_add_up_means.argtypes = (c_void_p,) 290 qubfast.qub_idlstats_add_up_means.restype = c_int 291 qubfast.qub_idlstats_get_means.argtypes = (c_void_p,) 292 qubfast.qub_idlstats_get_means.restype = c_double_p 293 qubfast.qub_idlstats_get_lifetimes.argtypes = (c_void_p,) 294 qubfast.qub_idlstats_get_lifetimes.restype = c_double_p 295 qubfast.qub_idlstats_get_median_lifetimes.argtypes = (c_void_p,) 296 qubfast.qub_idlstats_get_median_lifetimes.restype = c_double_p 297 qubfast.qub_idlstats_get_occupancies.argtypes = (c_void_p,) 298 qubfast.qub_idlstats_get_occupancies.restype = c_double_p 299 qubfast.qub_idlstats_get_nevents.argtypes = (c_void_p,) 300 qubfast.qub_idlstats_get_nevents.restype = c_int_p 301 qubfast.qub_idlstats_add_up_stds.argtypes = (c_void_p,) 302 qubfast.qub_idlstats_add_up_stds.restype = c_int 303 qubfast.qub_idlstats_get_stds.argtypes = (c_void_p,) 304 qubfast.qub_idlstats_get_stds.restype = c_double_p 305

306 -class IdlStats(object):

307 - def __init__(self, sampling):

308 self.sampling = sampling 309 self.obj = qubfast.qub_idlstats_create(sampling)

310 - def __del__(self):

311 qubfast.qub_idlstats_free(self.obj)

312 - def reset(self):

313 qubfast.qub_idlstats_reset(self.obj)

314 - def add(self, data, offset, firsts, lasts, classes):

315 # requires numpy arrays of float, int, int, int 316 if len(firsts): 317 qubfast.qub_idlstats_add(self.obj, cdata(data, c_float_p), len(firsts), offset, 318 cdata(firsts, c_int_p), cdata(lasts, c_int_p), cdata(classes, c_int_p))

319 - def add_up_means(self):

320 self.Nclass = qubfast.qub_idlstats_add_up_means(self.obj) 321 return self.Nclass

322 - def get_means(self):

323 return ptr_to_array(qubfast.qub_idlstats_get_means(self.obj), self.Nclass, c_double, 'float64')

324 - def get_lifetimes(self):

325 return ptr_to_array(qubfast.qub_idlstats_get_lifetimes(self.obj), self.Nclass, c_double, 'float64')

326 - def get_median_lifetimes(self):

327 return ptr_to_array(qubfast.qub_idlstats_get_median_lifetimes(self.obj), self.Nclass, c_double, 'float64')

328 - def get_occupancies(self):

329 return ptr_to_array(qubfast.qub_idlstats_get_occupancies(self.obj), self.Nclass, c_double, 'float64')

330 - def get_nevents(self):

331 return ptr_to_array(qubfast.qub_idlstats_get_nevents(self.obj), self.Nclass, c_int, 'int32')

332 - def add_up_stds(self):

333 return qubfast.qub_idlstats_add_up_stds(self.obj)

334 - def get_stds(self):

335 return ptr_to_array(qubfast.qub_idlstats_get_stds(self.obj), self.Nclass, c_double, 'float64')

336 337 338 qubfast.qub_mean_rss.argtypes = (c_float_p, c_int, c_double_p, c_double_p) 339 qubfast.qub_mean_rss.restype = None 340 qubfast.qub_mean_rss_init.argtypes = (c_float_p, c_int, c_double_p, c_double_p) 341 qubfast.qub_mean_rss_init.restype = c_void_p 342 qubfast.qub_mean_rss_free.argtypes = (c_void_p,) 343 qubfast.qub_mean_rss_free.restype = None 344 qubfast.qub_mean_rss_add.argtypes = (c_void_p, c_float_p, c_int) 345 qubfast.qub_mean_rss_add.restype = None 346 qubfast.qub_mean_rss_read.argtypes = (c_void_p, c_double_p, c_double_p) 347 qubfast.qub_mean_rss_read.restype = c_int 348

349 -def mean_rss(data):

350 if (not hasattr(data, 'dtype')) or (data.dtype != numpy.dtype('float32')): 351 d = numpy.array(data, dtype='float32') 352 else: 353 d = data 354 mean = c_double() 355 rss = c_double() 356 qubfast.qub_mean_rss(cdata(d, c_float_p), len(d), byref(mean), byref(rss)) 357 return mean.value, rss.value

358

359 -class MeanRSS(object):

360 - def __init__(self):

361 self.obj = qubfast.qub_mean_rss_init()

362 - def __del__(self):

363 qubfast.qub_mean_rss_free(self.obj)

364 - def add(self, data):

365 if (not hasattr(data, 'dtype')) or (data.dtype != numpy.dtype('float32')): 366 d = numpy.array(data, dtype='float32') 367 else: 368 d = data 369 qubfast.qub_mean_rss_add(self.obj, cdata(d, c_float_p), len(d))

370 - def read(self):

371 mean = c_double() 372 rss = c_double() 373 N = qubfast.qub_mean_rss_read(self.obj, byref(mean), byref(rss)) 374 return mean.value, rss.value, N

375

Source Code for Module qubx.fast.data