This function computes the log-likelihood of data given a model.
See also:
Up: IndexOur Markov model is a graph with colored vertices. A vertex is called a "state," and its color is a nonnegative integer called its "class." States in the same class are indistinguishable (same amp). To describe the vertices, provide the array
int clazz[Ns] = [class of each state]
int Ns = number of states
Each edge is labeled with its transition rate (probability per second). These form the matrix
Q, a Ns x Ns matrix with
\(Q_{a,b}\) = rate from state a to state b
\(Q_{a,a} = - \sum_i Q_{a,i}\) where \(i \neq a\)
Each \(Q_{a,b} = K0_{a,b} * Ligand_{a,b} * e^{K1_{a,b} * Voltage_{a,b}}\). You provide the Ns x Ns matrices
double **K0, **K1 of kinetic parameters
int **Ligand, **Voltage index of the ligand or voltage signal influencing each rate, or 0
with the diagonals undefined.
A pair of states is either connected in both directions or neither. To indicate un-connectedness, set
\(K0_{a,b} = K0_{b,a} = 0\)
For constant state entry probabilities, as described in (Qin...1996), provide the array
int P0[Ns]
For equilibrium entry probabilities, as described in (Milescu...2005), provide
P0 = NULLTODO: talk about amp/cond model and channel count
This is wrong -- leftover from MSL -- The data consist of one or more parallel signals. The first (index 0) is the Markovian one to be analyzed. Each additional signal describes a ligand or voltage variable. Signals are in model order; i.e., if (v = Voltage[a][b]) != 0, then the signal at index v is the voltage controlling a->b. For each signal you provide
int DwellCount number of events int Classes[DwellCount] class index of each event float Durations[DwellCount] length of each dwell in milliseconds double Amps[ClassCount] amplitude (or ligand/voltage value) of each class
and you give the signals together as
int Nsignal number of signals int *dwellCounts int **classeses float **durationses double **ampses
You provide one or more segments of data, each with the same number of signals. All together they are given as
int Nseg int Nsignal int **dwellCountses [segment][signal] int ***classeseses float ***durationseses double ***ampseses
We multiplex the signals to create a single idealized signal which changes class whenever any source signal changes class. Each plexi-class denotes a Markov class with a specific set of experimental constants (the other signals' idealized amplitudes).
Then, as in MIL, we subtract tdead. tdead is the longest duration of events that can't be reliably detected. MSL deletes any such events, by merging them with their prior. Then for computational reasons, tdead is subtracted from each event, and they're converted to seconds.Sadly, someone has to allocate memory to store this processed signal. It has the form
int Nseg int Nsignal int dwellCounts[Nseg] int **classeses float **durationses int Nplex int plexiclasses[2*Nplex] alternating markov-class[i/2], stimclass[i/2] int Nstim double stimclasses[Nstim*Nsignal] [stimcls * Nsignal + signal_ix], signal 0 undefined
First, we can scan your data to compute upper bounds for dwellCounts, Nplex and Nstim. You allocate dwellCounts[Nseg] and call:
end_html */ // MAXMLL_API void __stdcall mac_multiplex_bounds(int Nseg, int Nsignal, // int **dwellCountses, int ***classeseses, // float ***durationseses, double ***ampseses, // int *out_dwellCounts, int *out_Nplex, int *out_Nstim); // /* begin_htmlThen you allocate classeses[Nseg][dwellCounts[i]] and durationses, plexiclasses and stimclasses, and call this function to multiplex and process the signals:
end_html */ // MAXMLL_API void __stdcall mac_multiplex(int Nseg, int Nsignal, // int **dwellCountses, int ***classeseses, // float ***durationseses, double ***ampseses, // double tdead_ms, int *out_dwellCounts, // int **out_classeses, float **out_durationses, // int *out_Nplex, int *out_plexiclasses, // int *out_Nstim, double *out_stimclasses); typedef struct { float time, timeRel, cur;//x, var; } qubx_mac_sample; typedef struct { float fit, fit_var, ll, scaled_err; } qubx_mac_sample_out; // opaque type qubx_mac_data = (void *): // stores data (uneven sampling possible, or no data) and stimulus idealization QUBFAST_API void* __stdcall qubx_mac_data_create(int Nstimsig, int data_cap); QUBFAST_API void __stdcall qubx_mac_data_reset(void *data, int Nstimsig, int data_cap); QUBFAST_API void __stdcall qubx_mac_data_free(void *data); QUBFAST_API int __stdcall qubx_mac_data_append_segment(void *data); QUBFAST_API void __stdcall qubx_mac_data_append_chunk(void *data, int seg_ix, double dur_sec, double t0, int Nd, float *T, float *I,/* float *V,*/ int *dwellCounts, int **classeses, float **durationses, double **ampses); QUBFAST_API double* __stdcall qubx_mac_data_get_stimclasses(void *data, int *Nstimclass); // returns double[Nstimclass * Nsignal], Nsignal = Nstimsig + 1; QUBFAST_API double* __stdcall qubx_mac_data_get_stimclass_frac(void *data); // after get_stimclasses // returns double[Nstimclass]; weight of each stimclass; sums to 1.0; valid after calling get_stimclasses() QUBFAST_API qubx_mac_sample_out* __stdcall qubx_mac_data_get_buffer(void *data); typedef struct { void *storage; // you set: void *data; int useVar; int optNchannel; int accel; // and add_model()s // ready after setup_constraints(): int Npar, NfPar; double *pars, *fPars; // to gray out the GPU checkbox int gpu_fail; // output after qubx_mac_ll() double ll; // this is what's actually used (initialized from fastmodel nchannel) double float_Nchannel; } qubx_mac_param; QUBFAST_API qubx_mac_param* __stdcall qubx_mac_param_create(int Nmodel); QUBFAST_API void __stdcall qubx_mac_param_free(qubx_mac_param *p); QUBFAST_API void __stdcall qubx_mac_param_add_model(qubx_mac_param *p, qubx_model *model, qubx_stim_amps *ampm, qubx_stim_rates *ratesm); QUBFAST_API int __stdcall qubx_mac_setup_constraints(qubx_mac_param *p); QUBFAST_API void __stdcall qubx_mac_setup_nchannel(qubx_mac_param *p); QUBFAST_API void __stdcall qubx_mac_do_fPars_to_model(qubx_mac_param *p); QUBFAST_API int __stdcall qubx_mac_ll(qubx_mac_param *p); QUBFAST_API void __stdcall qubx_mac_calc_std(qubx_mac_param *p, double *InvHessian); #ifdef __cplusplus } #endif #endif