unit maxill_iface; { /* Copyright 1998-2011 Research Foundation State University of New York */ /* This file is part of QuB. */ /* QuB 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 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 program directory. If not, see */ /* . */ } (* begin_html See also: Up: Index end_html *) interface uses Alloc, AllocTypes, QUB_Python_Eigen; type TMsgOutFunc = function(Msg :PChar; Data :Pointer): Integer; stdcall; (* begin_html

Maximum Interval Likelihood (Qin et al, 1996)

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

See also: max_inter_ll.cpp

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

For data, you provide one or more idealized segments. A segment is described by 

 int   DwellCount            number of events
 int   Classes[DwellCount]   class index of each event
 float Durations[DwellCount] length of each dwell in seconds, less tdead

and you give the segments together as 

 int   Nseg                  number of segments
 int   **classeses
 float **durationses

tdead is the longest duration of events that can't be reliably detected. MIL assumes you have deleted any such events, by merging them with their prior. Also, for computational reasons, tdead should be subtracted from each event. We provide a utility which merges and shortens events, in-place:

end_html *) procedure mil_prep_events(DwellCount :Integer; Classes :TPAInteger; Durations :TPASingle; TDead_ms :Double); stdcall; external 'maxill.dll' name '_mil_prep_events@20'; procedure mil_prep_segments(Nseg :Integer; DwellCount :TPAInteger; Classeses :TPMInteger; Durationses :TPMSingle; TDead_ms :Double); stdcall; external 'maxill.dll' name '_mil_prep_segments@24'; (* begin_html

If any \(Ligand_{a,b}\ \neq 0\) or \(Voltage_{a,b}\ \neq 0\), you must provide an array 

 double *Constants

where e.g. if \(Ligand_{a,b} = 2\) then Constants[2] holds the ligand concentration. We assume all segments were recorded at the same constant conditions. For global fitting, call max_inter_ll separately for each dataset and sum the LL.

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

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

Limitations

end_html *) function inter_ll_arr(Ns :Integer; Clazz :TPAInteger; P0_or_nil :TPADouble; K0, K1 :TPADouble; Ligand, Voltage :TPAInteger; Constants :TPADouble; TDead_sec :Double; NSeg :Integer; DwellCounts :TPAInteger; Classeses :TPMInteger; Durationses :TPMSingle; out LL :Double): Integer; stdcall; external 'maxill.dll' name '_inter_ll_arr@60'; function inter_ll(Ns :Integer; Clazz :TPAInteger; P0_or_nil :TPADouble; K0, K1 :TPMDouble; Ligand, Voltage :TPMInteger; Constants :TPADouble; TDead_sec :Double; NSeg :Integer; DwellCounts :TPAInteger; Classeses :TPMInteger; Durationses :TPMSingle; out LL :Double): Integer; stdcall; external 'maxill.dll' name '_inter_ll@60'; // ********************************************************************** (* begin_html

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: max_inter_ll.cpp

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. 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. You can provide additional constants which will be treated as constant signals index N, N+1, ... All together they are given as 

 int    Nseg
 int    Nsignal
 int    **dwellCountses          [segment][signal]
 int    ***classeseses
 float  ***durationseses
 double ***ampseses
 int     Nconst
 double  *constants

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 *) procedure msl_multiplex_bounds(Nseg, Nsignal :Integer; DwellCountses :TPMInteger; Classeseses, Durationseses, Ampseses :TPAPointer; Out_DwellCounts :TPAInteger; out Nplex, Nstim :Integer); stdcall; external 'maxill.dll' name '_msl_multiplex_bounds@36'; (* begin_html

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

end_html *) procedure msl_multiplex(Nseg, Nsignal :Integer; DwellCountses :TPMInteger; Classeseses, Durationseses, Ampseses :TPAPointer; TDead_ms :Double; DwellCounts :TPAInteger; Out_Classeses :TPMInteger; Out_Durationses :TPMSingle; out Nplex :Integer; Out_Plexiclasses :TPAInteger; out NStim :Integer; Out_Stimclasses :TPADouble); stdcall; external 'maxill.dll' name '_msl_multiplex@60'; (* begin_html

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

end_html *) function subi_ll_arr(Ns :Integer; Clazz :TPAInteger; P0 :TPADouble; K0, K1 :TPADouble; Ligand, Voltage :TPAInteger; TDead_sec :Double; Nseg :Integer; DwellCounts :TPAInteger; Classeses :TPMInteger; Durationses :TPMSingle; Nsignal, Nplex :Integer; Plexiclasses :TPAInteger; Nstim :Integer; Stimclasses :TPADouble; out LL :Double; CustomEigen :BypassEigenFunc): Integer; stdcall; external 'maxill.dll' name '_subi_ll_arr@80'; function subi_ll (Ns :Integer; Clazz :TPAInteger; P0 :TPADouble; K0, K1 :TPMDouble; Ligand, Voltage :TPMInteger; TDead_sec :Double; Nseg :Integer; DwellCounts :TPAInteger; Classeses :TPMInteger; Durationses :TPMSingle; Nsignal, Nplex :Integer; Plexiclasses :TPAInteger; Nstim :Integer; Stimclasses :TPADouble; out LL :Double; CustomEigen :BypassEigenFunc; SegLL :TPADouble; MsgOut :TMsgOutFunc; MsgOutData :Pointer; StopFlag :TPAInteger): Integer; stdcall; external 'maxill.dll' name '_subi_ll@96'; // We build the full parameters K as // K = [k0_0->1, k0_0->2, k0_0->3, ..., K0_1->0, k0_1->2, k0_1->3, ..., // k1_0->1, k1_0->2, k1_0->3, ..., K1_1->0, k1_1->2, k1_1->3, ... ] // Actually, for the k0s we take log(k0) to make negative values impossible. // Breaking with past practice, we include all k_i->j, i<>j whether connected or not. // Also new: all the k0s are before all the k1s. function Maxill_CountK(Ns :Integer): Integer; stdcall; external 'maxill.dll' name '_CountK@4'; procedure Maxill_K01toK(Ns :Integer; K0, K1 :TPMDouble; K :TPADouble); stdcall; external 'maxill.dll' name '_K01toK@16'; // K alone isn't enough to reconstruct K0 and K1, since 0 could mean log(1) or else unconnected. // To disambiguate, you provide a valid K0 to be updated -- its zeros mean unconnected. function Maxill_KtoK01(Ns :Integer; K :TPADouble; K0, K1 :TPMDouble): Integer; stdcall; external 'maxill.dll' name '_KtoK01@16'; // returns 0 on success, or -(i+1) if exp(K[i]) is too tiny or huge // But some of those K are zero and have to stay zero, and the user may have other constraints. // We represent the constraints by the system of linear equations ... and derive the transformations... // You'll have to allocate memory for these constraints. Figure one row per constraint, except // that loop constraints involving V-sensitivity take 2 rows. Here's a safe upper bound: // Ncns_max = upper_bound(Ncns) = Nk [max fixed] + 2*Ncns_user // Ain[Ncns_max][Nk] // Bin[Ncns_max] // Ncns = 0 to start with -- no constraints // Some constraints are more a part of the representation than the specific model: // * where K0 == 0 there's no connection; fix k0 and k1 // * where V == 0 fix k1 procedure Maxill_AddAutoConstraints(Ns :Integer; K0, K1 :TPMDouble; Ligand, Voltage :TPMInteger; out Ncns :Integer; Ain :TPMDouble; Bin :TPADouble); stdcall; external 'maxill.dll' name '_AddAutoConstraints@32'; // the user might have some constraints: procedure Maxill_AddFixRateConstraint(Ns :Integer; K0, K1 :TPMDouble; Ligand, Voltage :TPMInteger; out Ncns :Integer; Ain :TPMDouble; Bin :TPADouble; A, B :Integer); stdcall; external 'maxill.dll' name '_AddFixRateConstraint@40'; procedure Maxill_AddFixExpConstraint(Ns :Integer; K0, K1 :TPMDouble; Ligand, Voltage :TPMInteger; out Ncns :Integer; Ain :TPMDouble; Bin :TPADouble; A, B :Integer); stdcall; external 'maxill.dll' name '_AddFixExpConstraint@40'; procedure Maxill_AddScaleRateConstraint(Ns :Integer; K0, K1 :TPMDouble; Ligand, Voltage :TPMInteger; out Ncns :Integer; Ain :TPMDouble; Bin :TPADouble; A, B, C, D :Integer); stdcall; external 'maxill.dll' name '_AddScaleRateConstraint@48'; procedure Maxill_AddScaleExpConstraint(Ns :Integer; K0, K1 :TPMDouble; Ligand, Voltage :TPMInteger; out Ncns :Integer; Ain :TPMDouble; Bin :TPADouble; A, B, C, D :Integer); stdcall; external 'maxill.dll' name '_AddScaleExpConstraint@48'; procedure Maxill_AddGeneralizedConstraint(Ns :Integer; out Ncns :Integer; Ain :TPMDouble; Bin :TPADouble; Aentry :TPADouble; Bentry :Double); stdcall; external 'maxill.dll' name '_AddGeneralizedConstraint@28'; // these fail nonzero if the constraint can't be satisfied for all stimuli // details are called back to MsgOut function Maxill_AddLoopBalConstraint(Ns :Integer; K0, K1 :TPMDouble; Ligand, Voltage :TPMInteger; out Ncns :Integer; Ain :TPMDouble; Bin :TPADouble; Nloop :Integer; LoopStates :TPAInteger; MsgOut :TMsgOutFunc; MsgOutData :Pointer): Integer; stdcall; external 'maxill.dll' name '_AddLoopBalConstraint@48'; function Maxill_AddLoopImbalConstraint(Ns :Integer; K0, K1 :TPMDouble; Ligand, Voltage :TPMInteger; out Ncns :Integer; Ain :TPMDouble; Bin :TPADouble; Nloop :Integer; LoopStates :TPAInteger; MsgOut :TMsgOutFunc; MsgOutData :Pointer): Integer; stdcall; external 'maxill.dll' name '_AddLoopImbalConstraint@48'; // presto change-o we take // Ain * K = Bin // to get // K = Acns * P + B // P = Ainv * K // Np = Nk - Nc // You alloc // Acns[Nk][Np] // Bcns[Nk] // Ainv[Np][Nk] procedure Maxill_ReduceConstraints(Ns :Integer; out Ncns :Integer; Ain :TPMDouble; Bin :TPADouble); stdcall; external 'maxill.dll' name '_ReduceConstraints@16'; function Maxill_SetupLinearConstraints(Nc :Integer; Ain :TPMDouble; Bin :TPADouble; Nk :Integer; K :TPADouble; P_out :TPADouble; Acns :TPMDouble; Bcns :TPADouble; Ainv :TPMDouble; MsgOut :TMsgOutFunc; MsgOutData :Pointer): Integer; stdcall; external 'maxill.dll' name '_SetupLinearConstraints@44'; // here's another linear transformation, this time to scale initial parameter values to 1.0. // (in the case where p[i] == 0, we scale by 1 and translate by -1) // P = As*S + Bs // S = Asinv * (P - Bs) procedure Maxill_SetupStartAtOnes(Np :Integer; P :TPADouble; Ascal :TPMDouble; Bscal :TPADouble; AscalInv :TPMDouble); stdcall; external 'maxill.dll' name '_SetupStartAtOnes@20'; // these actually work for either transform procedure Maxill_KtoPars(Nk, Np :Integer; K :TPADouble; A, Ainv :TPMDouble; B, P :TPADouble); stdcall; external 'maxill.dll' name '_KtoPars@28'; procedure Maxill_ParsToK(Np, Nk :Integer; P :TPADouble; A, Ainv :TPMDouble; B, K :TPADouble); stdcall; external 'maxill.dll' name '_ParsToK@28'; function Maxill_k_ix_of(Ns, a, b :Integer; IsExp :Boolean=False): Integer; procedure Maxill_FirstLatencies(A :TPMDouble; Peq :TPADouble; Sampling :Double; N :Integer; FL :TPMDouble); stdcall; external 'maxill.dll' name '_FirstLatencies@24'; type TQUB_Kalman_Rec = record Storage :Pointer; NProcess :Integer; NMeasure :Integer; X :TPADouble; Z :TPADouble; P :TPMDouble; A :TPMDouble; Q :TPMDouble; H :TPMDouble; R :TPMDouble; end; TQUB_Kalman_PRec = ^TQUB_Kalman_Rec; function qub_kalman_filter_create(Nprocess, Nmeasure :Integer): TQUB_Kalman_PRec; cdecl; external 'maxill.dll' name 'qub_kalman_filter_create'; procedure qub_kalman_filter_free(F :TQUB_Kalman_PRec); cdecl; external 'maxill.dll' name 'qub_kalman_filter_free'; procedure qub_kalman_filter_next(F :TQUB_Kalman_PRec; Z :TPADouble); cdecl; external 'maxill.dll' name 'qub_kalman_filter_next'; procedure qub_kalman_filter_reset(F :TQUB_Kalman_PRec); cdecl; external 'maxill.dll' name 'qub_kalman_filter_reset'; type TQUB_Kalman_Filter = class public Rec :TQUB_Kalman_PRec; constructor Create(NProcess, NMeasure :Integer); destructor Destroy; override; procedure Next(Z :TPADouble); procedure Reset; end; implementation function Maxill_k_ix_of(Ns, a, b: Integer; IsExp :Boolean): Integer; begin Result := a*(Ns-1) + b; if b > a then Dec(Result); if IsExp then Inc(Result, Ns*(Ns-1)); end; constructor TQUB_Kalman_Filter.Create(NProcess, NMeasure :Integer); begin Rec := qub_kalman_filter_create(NProcess, NMeasure); end; destructor TQUB_Kalman_Filter.Destroy; begin try qub_kalman_filter_free(Rec); finally inherited; end; end; procedure TQUB_Kalman_Filter.Next(Z :TPADouble); begin qub_kalman_filter_next(Rec, Z); end; procedure TQUB_Kalman_Filter.Reset; begin qub_kalman_filter_reset(Rec); end; end.