max_subi_ll.h.html mathcode2html   
 Source file:   max_subi_ll.h
 Converted:   Sat May 9 2015 at 14:46:13
 This documentation file will not reflect any later changes in the source file.

$$\phantom{******** If you see this on the webpage then the browser could not locate *********}$$
$$\phantom{******** jsMath/easy/load.js or the variable root is set wrong in this file *********}$$
$$\newcommand{\vector}[1]{\left[\begin{array}{c} #1 \end{array}\right]}$$ $$\newenvironment{matrix}{\left[\begin{array}{cccccccccc}} {\end{array}\right]}$$ $$\newcommand{\A}{{\cal A}}$$ $$\newcommand{\W}{{\cal W}}$$ 

/* Copyright 2008-2013 Research Foundation State University of New York   */

/* This file is part of QUB Express.                                      */

/* QUB Express is free software; you can redistribute it and/or modify    */
/* it under the terms of the GNU General Public License as published by   */
/* the Free Software Foundation, either version 3 of the License, or      */
/* (at your option) any later version.                                    */

/* QUB Express is distributed in the hope that it will be useful,         */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of         */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          */
/* GNU General Public License for more details.                           */

/* You should have received a copy of the GNU General Public License,     */
/* named LICENSE.txt, in the QUB Express program directory.  If not, see  */
/* <http://www.gnu.org/licenses/>.                                        */

#ifndef MAX_SUBI_LL_H
#define MAX_SUBI_LL_H

#include "maxill.h"
#include "callbk_reportfun.h"
#include "../qubfast/qubx_model.h"

/*

Maximum Sub-Interval Likelihood

[adaptation of (Qin et al, 1996)] by Chris Nicolai 2009

This function computes the forward log-likelihood of data given a model.

See also:

Up: Index

The Model

Our 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 = NULL

The Data

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: 

 */

#ifdef __cplusplus
extern "C" {
#endif

extern "C" MAXILL_API void msl_multiplex_bounds(int Nseg, int Nsignal,
							  int **dwellCountses, int ***classeseses,
							  float ***durationseses, double ***ampseses,
							  int *out_dwellCounts, int *out_Nplex, int *out_Nstim);

/*

Then you allocate classeses[Nseg][dwellCounts[i]] and durationses, plexiclasses and stimclasses, and call this function to multiplex and process the signals: 

*/

extern "C" MAXILL_API void msl_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);

/*

Now you can call the subi_ll function below.

double *LL: upon return, contains log(prob. that model generated this data)

Returns: 0 on success, error codes to be defined.

Limitations

  • no gradients
  • no maximization
  • no support for multiple identical models (channel count > 1)

*/

MAXILL_API int subi_ll_arr(
	int Ns, int *clazz, double *P0_in,  double *K0_in, double *K1_in, double *K2_in,
	int *Ligand_in, int *Voltage_in, int *Pressure_in,
	double tdead_sec, int Nseg, int *dwellCounts, int **classeses, float **durationses,
	int Nsig, int Nplex, int *plexiclasses, int Nstim, double *stimclasses,
	double *ll, eigen_func custom_eigen);

MAXILL_API int subi_ll(
	int Ns, int *clazz, double *P0,  double **K0, double **K1, double **K2,
	int **Ligand, int **Voltage, int **Pressure,
	double tdead_sec, int Nseg, int *dwellCounts, int **classeses, float **durationses,
	int Nsig, int Nplex, int *plexiclasses, int Nstim, double *stimclasses,
	double *LL, eigen_func custom_eigen,
	double *segLL, callbk_reportfun output_cb, void *output_data, int *stop_flag);


#define MSL_ACCEL_NONE 0
#define MSL_ACCEL_OPENMP 1
#define MSL_ACCEL_OPENCL 2

MAXILL_API void* msl_accel_context_init(int *accel);
MAXILL_API void msl_accel_context_free(void *ctx_);
MAXILL_API void* msl_accel_data_init(void *ctx_, int Nseg, int *dwellCounts, int **classeses, float **durationses, int Nsig, int Nplex, int *plexiclass, int Nstim, double *stimclasses);
MAXILL_API void msl_accel_data_free(void *data_);
MAXILL_API double* msl_accel_data_get_ll(void *data_);
MAXILL_API void* msl_accel_models_init(void *ctx_, int dim, int Nmodel, int Nsc);
MAXILL_API void msl_accel_models_free(void *models_);
MAXILL_API void msl_accel_models_reset(void *models_);
MAXILL_API int msl_accel_models_setup_model(void *models_, void *data_, int Ns, int *clazz, double *P0_in, double **K0, double **K1, double **K2,
					    int **Ligand, int **Voltage, int **Pressure, double tdead_sec,
					    eigen_func custom_eigen, callbk_reportfun report_cb, void *report_data);
MAXILL_API int msl_accel_models_setup_model_arr(void *models_, void *data_, int Ns, int *clazz, double *P0_in, double *K0, double *K1, double *K2,
						int *Ligand, int *Voltage, int *Pressure, double tdead_sec,
						eigen_func custom_eigen, callbk_reportfun report_cb, void *report_data);
MAXILL_API double* msl_accel_models_get_ll(void *models_);
MAXILL_API int msl_accel_ll(void *ctx_, void *data_, void *models_, int accel);
MAXILL_API int mil_accel_ll(void *ctx_, void *data_, void *models_, int accel);

#ifdef __cplusplus
}
#endif

#endif