& eQ);
/*
begin_html
As discussed in (Milescu 2005), it can be preferable to start an experiment at equilibrium.
Entry probabilities (P0) are then no longer constant, but a function of rate constants.
We use the "direct" method given in Neher and Sakmann:
Defining \(S = [Q | 1]\) and \(u\) a row vector of ones,
\[P_{eq} = u \cdot (S \cdot S^T)^{-1}\]
end_html
*/
void QtoPe(matrix& Q, matrix& Pe);
/*
begin_html
As discussed in (Qin 1997), the forward probability vectors alpha must be
frequently rescaled to 1 in order to stay within machine precision. Conveniently,
the log likelihood is the sum of these log scaling factors.
end_html
*/
inline double alpha_logscale(vector& alpha)
{
double tot = 0.0;
for (vector::iterator ii = alpha.begin(); ii != alpha.end(); ++ii)
tot += *ii;
double scale = tot ? (1.0 / tot) : 1e10;
for (vector::iterator ii = alpha.begin(); ii != alpha.end(); ++ii)
*ii *= scale;
// std::cerr << scale << std::endl;
return - log(scale);
}
#endif
/*
What are the optimization parameters? (some of) the off-diagonal elements of K0 and K1.
Actually, we use log(K0) to keep them positive.
Breaking with past practice, we include all the unconnected zeros,
and put all the K0s before all the K1s. As we'll see, this is a column vector.
*/
extern "C" MAXILL_API int __stdcall CountK(int Ns);
extern "C" MAXILL_API void __stdcall K01toK(int Ns, double **K0, double **K1, double *K);
// 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.
extern "C" MAXILL_API int __stdcall KtoK01(int Ns, double *K, double **K0, double **K1);
// 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 (11) and derive the transformations (12, 13).
// Ncns_max = upper_bound(Ncns) = Nk [fixed] + 2*Ncns_user
// Ain[Ncns_max][Nk]
// Bin[Ncns_max]
// Ncns = 0 (or how many pre-existing constraints in Ain,Bin)
// Some constraints are more a part of the representation than the model:
// * where K0 == 0 there's no connection; fix k0 and k1
// * where V == 0 fix k1
extern "C" MAXILL_API void __stdcall AddAutoConstraints(int Ns, double **K0, double **K1, int **L, int **V,
int *Ncns, double **Ain, double *Bin);
// The user might have some constraints too:
extern "C" MAXILL_API void __stdcall AddFixRateConstraint(int Ns, double **K0, double **K1, int **L, int **V,
int *Ncns, double **Ain, double *Bin,
int a, int b);
extern "C" MAXILL_API void __stdcall AddFixExpConstraint(int Ns, double **K0, double **K1, int **L, int **V,
int *Ncns, double **Ain, double *Bin,
int a, int b);
extern "C" MAXILL_API void __stdcall AddScaleRateConstraint(int Ns, double **K0, double **K1, int **L, int **V,
int *Ncns, double **Ain, double *Bin,
int a, int b, int c, int d);
extern "C" MAXILL_API void __stdcall AddScaleExpConstraint(int Ns, double **K0, double **K1, int **L, int **V,
int *Ncns, double **Ain, double *Bin,
int a, int b, int c, int d);
extern "C" MAXILL_API void __stdcall AddGeneralizedConstraint(int Ns, int *Ncns, double **Ain, double *Bin,
double *Aentry, double Bentry);
extern "C" MAXILL_API int __stdcall AddLoopBalConstraint(int Ns, double **K0, double **K1, int **L, int **V,
int *Ncns, double **Ain, double *Bin,
int Nloop, int *loopstates,
callbk_reportfun report_cb, void *report_data);
extern "C" MAXILL_API int __stdcall AddLoopImbalConstraint(int Ns, double **K0, double **K1, int **L, int **V,
int *Ncns, double **Ain, double *Bin,
int Nloop, int *loopstates,
callbk_reportfun report_cb, void *report_data);
/*
SetupLinearConstraints requires that all constraints be linearly independent.
First call ReduceConstraints to delete any dependent constraints.
*/
extern "C" MAXILL_API void __stdcall ReduceConstraints(int Ns, int *Ncns, double **Ain, double *Bin);
/*
Now that we have K0, K1, L, V and rate vector K,
we do some magic with the SVD to change your linear constraints
A * K = B
into a transformation between K and the column vector P of free parameters
K = A*P + B
P = Ainv * K
*/
extern "C" MAXILL_API int __stdcall SetupLinearConstraints(int Nc, double **Ain, double *Bin,
int Nk, double *K, double *P_out, // Nk = Nc + Np
double **Acns, double *Bcns, double **Ainv,
callbk_reportfun report_cb, void *report_data);
// Acns[Nk][Np]; Bcns[Nk]; Ainv[Np][Nk]
extern "C" MAXILL_API void __stdcall KtoPars(int Nk, int Np, double *K, double **A, double **Ainv, double *B, double *pars);
extern "C" MAXILL_API void __stdcall ParsToK(int Np, int Nk, double *P, double **A, double **Ainv, double *B, double *K);
/*
And yet 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 = As.I * (P - Bs)
*/
extern "C" MAXILL_API void __stdcall SetupStartAtOnes(int Np, double *pars, double **Ascal, double *Bscal, double **AscalInv);
#endif