| qub_constraints.cpp.html | mathcode2html |
| Source file: qub_constraints.cpp | |
| Converted: Tue May 19 2015 at 15:35:14 | |
| 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-2014 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/>. */ /*
|
|
See also: Up: Index |
*/
#ifdef _WIN32
#include <windows.h>
#else
#include <stdlib.h>
#define BOOL int
#define TRUE 1
#define FALSE 0
#endif
#include <iostream>
#include <fstream>
using std::cout;
using std::cerr;
using std::endl;
#include <set>
using std::set;
#include <map>
using std::map;
#include <string.h>
#include "qub_constraints.h"
#include "ublas_matrixutil.h"
#include "callbk_reportstream.h"
//std::ofstream logg;
/*
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.
UPDATE: we omit any fixed zero from the final K vector
*/
inline int k_ix_of(int Ns, int a, int b)
{
return (a * (Ns - 1)) + b - ((b > a) ? 1 : 0);
}
extern "C" QUBFAST_API int CountK(int Ns)
{
return 3*(Ns*Ns - Ns); // the off-diagonal elements, twice
}
extern "C" QUBFAST_API void K01toK(int Ns, double **K0, double **K1, double *K)
{
int Noffd = CountK(Ns) / 3;
int i=0;
for (int a=0; a<Ns; ++a)
for (int b=0; b<Ns; ++b)
if ( a != b ) {
double x = K0[a][b];
if ( x > 0.0 ) // tricky; 1->0 and 0->0; see KtoK01
x = log(x);
else
x = 0.0;
K[i] = x;
K[Noffd+i] = K1[a][b];
++i;
}
}
extern "C" QUBFAST_API void K012toK(int Ns, double **K0, double **K1, double **K2, double *K)
{
int Noffd = CountK(Ns) / 3;
int i=0;
for (int a=0; a<Ns; ++a)
for (int b=0; b<Ns; ++b)
if ( a != b ) {
double x = K0[a][b];
if ( x > 0.0 ) // tricky; 1->0 and 0->0; see KtoK01
x = log(x);
else
x = 0.0;
K[i] = x;
K[Noffd+i] = K1[a][b];
K[2*Noffd+i] = K2[a][b];
++i;
}
}
// K alone isn't enough to reconstruct K0 and K1, since 0 could mean log(1), or else unconnected.
// We disambiguate by checking the original K0 for zeros, which mean unconnected.
#define MIN_LOG_K0 -30.0
#define MAX_LOG_K0 30.0
extern "C" QUBFAST_API int KtoK01(int Ns, double *K, double **K0, double **K1)
{
int Noffd = CountK(Ns) / 3;
int i = 0;
for (int a=0; a<Ns; ++a)
for (int b=0; b<Ns; ++b)
if ( a != b ) {
if ( K0[a][b] ) {
if ( (K[i] > MIN_LOG_K0) && (K[i] < MAX_LOG_K0) )
K0[a][b] = exp(K[i]);
else
return -(i+1);
}
K1[a][b] = K[Noffd+i];
++i;
}
return 0;
}
extern "C" QUBFAST_API int KtoK012(int Ns, double *K, double **K0, double **K1, double **K2)
{
int Noffd = CountK(Ns) / 3;
int i = 0;
for (int a=0; a<Ns; ++a)
for (int b=0; b<Ns; ++b)
if ( a != b ) {
if ( K0[a][b] ) {
if ( (K[i] > MIN_LOG_K0) && (K[i] < MAX_LOG_K0) )
K0[a][b] = exp(K[i]);
else
return -(i+1);
}
K1[a][b] = K[Noffd+i];
K2[a][b] = K[2*Noffd+i];
++i;
}
return 0;
}
extern "C" QUBFAST_API int dK_to_dK01(int Ns, double *K, double *dK, double **K0, double **d_K0, double **d_K1)
{
int Noffd = CountK(Ns) / 3;
int i = 0;
for (int a=0; a<Ns; ++a)
for (int b=0; b<Ns; ++b)
if ( a != b ) {
if ( K0[a][b] ) {
d_K0[a][b] = dK[i];
if ( (K[i] > MIN_LOG_K0) && (K[i] < MAX_LOG_K0) )
d_K0[a][b] *= exp(K[i]);
}
d_K1[a][b] = dK[Noffd+i];
++i;
}
return 0;
}
extern "C" QUBFAST_API int dK_to_dK012(int Ns, double *K, double *dK, double **K0, double **d_K0, double **d_K1, double **d_K2)
{
int Noffd = CountK(Ns) / 3;
int i = 0;
for (int a=0; a<Ns; ++a)
for (int b=0; b<Ns; ++b)
if ( a != b ) {
if ( K0[a][b] ) {
d_K0[a][b] = dK[i];
if ( (K[i] > MIN_LOG_K0) && (K[i] < MAX_LOG_K0) )
d_K0[a][b] *= exp(K[i]);
}
d_K1[a][b] = dK[Noffd+i];
d_K2[a][b] = dK[2*Noffd+i];
++i;
}
return 0;
}
// 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).
// 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
// Ncns_max = upper_bound(Ncns) = Nk [fixed] + 3*Ncns_user
// Ain[Ncns_max][Nk]
// Bin[Ncns_max]
// Ncns = 0 (or how many pre-existing constraints in Ain,Bin)
extern "C" QUBFAST_API void AddAutoConstraints(int Ns, double **K0, QUBFAST_VAR_NOT_USED double ** K1, QUBFAST_VAR_NOT_USED double ** K2,
QUBFAST_VAR_NOT_USED int ** L, int **V, int **P,
int *Ncns, double **Ain, double *Bin)
// """Returns the constraints (Ain, Bin) [11] which fix unconnected k0 and unsensitive k1."""
{
int Nk = CountK(Ns);
int Noffd = Nk/3;
int i=0;
for (int a=0; a<Ns; ++a)
for (int b=0; b<Ns; ++b)
if (a != b) {
if ( ! K0[a][b] ) {
for (int j=0; j<Nk; ++j)
Ain[*Ncns][j] = 0.0;
Ain[*Ncns][i] = 1.0;
Bin[*Ncns] = 0.0;
V[a][b] = 0;
++*Ncns;
}
if ( ! V[a][b] ) {
for (int j=0; j<Nk; ++j)
Ain[*Ncns][j] = 0.0;
Ain[*Ncns][Noffd+i] = 1.0;
Bin[*Ncns] = 0.0;
++*Ncns;
}
if ( ! P[a][b] ) {
for (int j=0; j<Nk; ++j)
Ain[*Ncns][j] = 0.0;
Ain[*Ncns][2*Noffd+i] = 1.0;
Bin[*Ncns] = 0.0;
++*Ncns;
}
++i;
}
}
// The user might have some constraints too:
extern "C" QUBFAST_API void AddFixRateConstraint(int Ns, double **K0, QUBFAST_VAR_NOT_USED double ** K1, QUBFAST_VAR_NOT_USED double ** K2,
QUBFAST_VAR_NOT_USED int ** L, QUBFAST_VAR_NOT_USED int ** V, QUBFAST_VAR_NOT_USED int ** P,
int *Ncns, double **Ain, double *Bin,
int a, int b)
{
int Nk = CountK(Ns);
int k_ix = k_ix_of(Ns, a, b);
for (int j=0; j<Nk; ++j)
Ain[*Ncns][j] = 0.0;
Ain[*Ncns][k_ix] = 1.0;
Bin[*Ncns] = log(K0[a][b]);
++*Ncns;
}
extern "C" QUBFAST_API void AddFixExpConstraint(int Ns, QUBFAST_VAR_NOT_USED double ** K0, double **K1, QUBFAST_VAR_NOT_USED double ** K2,
QUBFAST_VAR_NOT_USED int ** L, QUBFAST_VAR_NOT_USED int ** V, QUBFAST_VAR_NOT_USED int ** P,
int *Ncns, double **Ain, double *Bin,
int a, int b)
{
int Nk = CountK(Ns);
int k_ix = k_ix_of(Ns, a, b);
for (int j=0; j<Nk; ++j)
Ain[*Ncns][j] = 0.0;
Ain[*Ncns][k_ix + Nk/3] = 1.0;
Bin[*Ncns] = K1[a][b];
++*Ncns;
}
extern "C" QUBFAST_API void AddFixPresConstraint(int Ns, QUBFAST_VAR_NOT_USED double ** K0, QUBFAST_VAR_NOT_USED double **K1, double ** K2,
QUBFAST_VAR_NOT_USED int ** L, QUBFAST_VAR_NOT_USED int ** V, QUBFAST_VAR_NOT_USED int ** P,
int *Ncns, double **Ain, double *Bin,
int a, int b)
{
int Nk = CountK(Ns);
int k_ix = k_ix_of(Ns, a, b);
for (int j=0; j<Nk; ++j)
Ain[*Ncns][j] = 0.0;
Ain[*Ncns][k_ix + 2*Nk/3] = 1.0;
Bin[*Ncns] = K2[a][b];
++*Ncns;
}
// we optimize log(k0); this constraint actually means log(k0ab) - log(k0cd) = const
extern "C" QUBFAST_API void AddScaleRateConstraint(int Ns, double **K0, QUBFAST_VAR_NOT_USED double ** K1, QUBFAST_VAR_NOT_USED double ** K2,
QUBFAST_VAR_NOT_USED int ** L, QUBFAST_VAR_NOT_USED int ** V, QUBFAST_VAR_NOT_USED int ** P,
int *Ncns, double **Ain, double *Bin,
int a, int b, int c, int d)
{
int Nk = CountK(Ns);
int k_ix1 = k_ix_of(Ns, a, b);
int k_ix2 = k_ix_of(Ns, c, d);
for (int j=0; j<Nk; ++j)
Ain[*Ncns][j] = 0.0;
Ain[*Ncns][k_ix1] = 1.0;
Ain[*Ncns][k_ix2] = - 1.0;
Bin[*Ncns] = log(K0[a][b]) - log(K0[c][d]);
++*Ncns;
}
// ratio = (k1ab / k1cd)
// k1ab - ratio * k1cd = 0
extern "C" QUBFAST_API void AddScaleExpConstraint(int Ns, QUBFAST_VAR_NOT_USED double ** K0, double **K1, QUBFAST_VAR_NOT_USED double ** K2,
QUBFAST_VAR_NOT_USED int ** L, QUBFAST_VAR_NOT_USED int ** V, QUBFAST_VAR_NOT_USED int ** P,
int *Ncns, double **Ain, double *Bin,
int a, int b, int c, int d)
{
int Nk = CountK(Ns);
int k_ix1 = k_ix_of(Ns, a, b);
int k_ix2 = k_ix_of(Ns, c, d);
for (int j=0; j<Nk; ++j)
Ain[*Ncns][j] = 0.0;
Ain[*Ncns][k_ix1 + Nk/3] = 1.0;
Ain[*Ncns][k_ix2 + Nk/3] = - K1[a][b] / K1[c][d];
Bin[*Ncns] = 0.0;
++*Ncns;
}
extern "C" QUBFAST_API void AddScalePresConstraint(int Ns, QUBFAST_VAR_NOT_USED double ** K0, QUBFAST_VAR_NOT_USED double **K1, double ** K2,
QUBFAST_VAR_NOT_USED int ** L, QUBFAST_VAR_NOT_USED int ** V, QUBFAST_VAR_NOT_USED int ** P,
int *Ncns, double **Ain, double *Bin,
int a, int b, int c, int d)
{
int Nk = CountK(Ns);
int k_ix1 = k_ix_of(Ns, a, b);
int k_ix2 = k_ix_of(Ns, c, d);
for (int j=0; j<Nk; ++j)
Ain[*Ncns][j] = 0.0;
Ain[*Ncns][k_ix1 + 2*Nk/3] = 1.0;
Ain[*Ncns][k_ix2 + 2*Nk/3] = - K2[a][b] / K2[c][d];
Bin[*Ncns] = 0.0;
++*Ncns;
}
// nonzero if there's a problem, reported through report_cb
extern "C" QUBFAST_API int AddLoopBalConstraint(int Ns, QUBFAST_VAR_NOT_USED double ** K0, double **K1, double **K2, int **L, int **V, int **P,
int *Ncns, double **Ain, double *Bin,
int Nloop, int *loopstates,
callbk_reportfun report_cb, void *report_data)
{
callbk_reportstream output(report_cb, report_data);
set<int> ligs, volts, press;
map<int,int> ligsF, ligsB, voltsPlus, voltsMinus, pressPlus, pressMinus;
int Nk = CountK(Ns);
if ( loopstates[Nloop-1] == loopstates[0] )
--Nloop; // don't repeat first as last
for (int j=0; j<Nk; ++j)
Ain[*Ncns][j] = 0.0;
for (int i=0; i<Nloop; ++i) {
int a = loopstates[i];
int b = loopstates[(i+1)%Nloop];
int k_ix1 = k_ix_of(Ns, a, b);
int k_ix2 = k_ix_of(Ns, b, a);
Ain[*Ncns][k_ix1] = 1.0;
Ain[*Ncns][k_ix2] = - 1.0;
if ( L[a][b] ) {
ligs.insert(L[a][b]);
++ligsF[L[a][b]];
}
if ( L[b][a] ) {
ligs.insert(L[b][a]);
++ligsB[L[b][a]];
}
if ( V[a][b] ) {
volts.insert(V[a][b]);
if ( K1[a][b] < 0 )
++voltsMinus[V[a][b]];
else
++voltsPlus[V[a][b]];
}
if ( V[b][a] ) {
volts.insert(V[b][a]);
if ( K1[b][a] < 0 )
++voltsPlus[V[b][a]];
else
++voltsMinus[V[b][a]];
}
if ( P[a][b] ) {
press.insert(P[a][b]);
if ( K2[a][b] < 0 )
++pressMinus[P[a][b]];
else
++pressPlus[P[a][b]];
}
if ( P[b][a] ) {
press.insert(P[b][a]);
if ( K2[b][a] < 0 )
++pressPlus[P[b][a]];
else
++pressMinus[P[b][a]];
}
}
Bin[*Ncns] = 0.0;
for (set<int>::iterator ligi = ligs.begin(); ligi != ligs.end(); ++ligi) {
if ( ligsF[*ligi] != ligsB[*ligi] ) {
output << "In loop";
for (int i=0; i<Nloop; ++i)
output << " " << (loopstates[i]+1);
output << ": Ligand (signal " << *ligi
<< ") must influence the same number of rates forward and back." << endl;
return -1;
}
}
for (set<int>::iterator voli = volts.begin(); voli != volts.end(); ++voli) {
if ( (! voltsPlus[*voli]) || (! voltsMinus[*voli]) ) {
output << "In loop";
for (int i=0; i<Nloop; ++i)
output << " " << (loopstates[i]+1);
output << ": Voltage (signal " << *voli
<< ") must influence at least two rates,"
<< " and at least one must have opposite sign or direction." << endl;
return -2;
}
}
for (set<int>::iterator presi = press.begin(); presi != press.end(); ++presi) {
if ( (! pressPlus[*presi]) || (! pressMinus[*presi]) ) {
output << "In loop";
for (int i=0; i<Nloop; ++i)
output << " " << (loopstates[i]+1);
output << ": Pressure (signal " << *presi
<< ") must influence at least two rates,"
<< " and at least one must have opposite sign or direction." << endl;
return -2;
}
}
++*Ncns;
if ( volts.size() ) {
for (int j=0; j<Nk; ++j)
Ain[*Ncns][j] = 0.0;
for (int i=0; i<Nloop; ++i) {
int a = loopstates[i];
int b = loopstates[(i+1)%Nloop];
int k_ix1 = k_ix_of(Ns, a, b);
int k_ix2 = k_ix_of(Ns, b, a);
Ain[*Ncns][k_ix1 + Nk/3] = 1.0;
Ain[*Ncns][k_ix2 + Nk/3] = - 1.0;
}
Bin[*Ncns] = 0.0;
++*Ncns;
}
if ( press.size() ) {
for (int j=0; j<Nk; ++j)
Ain[*Ncns][j] = 0.0;
for (int i=0; i<Nloop; ++i) {
int a = loopstates[i];
int b = loopstates[(i+1)%Nloop];
int k_ix1 = k_ix_of(Ns, a, b);
int k_ix2 = k_ix_of(Ns, b, a);
Ain[*Ncns][k_ix1 + 2*Nk/3] = 1.0;
Ain[*Ncns][k_ix2 + 2*Nk/3] = - 1.0;
}
Bin[*Ncns] = 0.0;
++*Ncns;
}
return 0;
}
extern "C" QUBFAST_API int AddLoopImbalConstraint(int Ns, double **K0, double **K1, double **K2, int **L, int **V, int **P,
int *Ncns, double **Ain, double *Bin,
int Nloop, int *loopstates,
callbk_reportfun report_cb, void *report_data)
{
callbk_reportstream output(report_cb, report_data);
set<int> ligs, volts, press;
map<int,int> ligsF, ligsB, voltsPlus, voltsMinus, pressMinus, pressPlus;
int Nk = CountK(Ns);
if ( loopstates[Nloop-1] == loopstates[0] )
--Nloop; // don't repeat first as last
for (int j=0; j<Nk; ++j)
Ain[*Ncns][j] = 0.0;
double imbal = 0.0;
for (int i=0; i<Nloop; ++i) {
int a = loopstates[i];
int b = loopstates[(i+1)%Nloop];
int k_ix1 = k_ix_of(Ns, a, b);
int k_ix2 = k_ix_of(Ns, b, a);
imbal += log(K0[a][b]) - log(K0[b][a]);
Ain[*Ncns][k_ix1] = 1.0;
Ain[*Ncns][k_ix2] = - 1.0;
if ( L[a][b] ) {
ligs.insert(L[a][b]);
++ligsF[L[a][b]];
}
if ( L[b][a] ) {
ligs.insert(L[b][a]);
++ligsB[L[b][a]];
}
if ( V[a][b] ) {
volts.insert(V[a][b]);
if ( K1[a][b] < 0 )
++voltsMinus[V[a][b]];
else
++voltsPlus[V[a][b]];
}
if ( V[b][a] ) {
volts.insert(V[b][a]);
if ( K1[b][a] < 0 )
++voltsPlus[V[b][a]];
else
++voltsMinus[V[b][a]];
}
if ( P[a][b] ) {
press.insert(P[a][b]);
if ( K2[a][b] < 0 )
++pressMinus[P[a][b]];
else
++pressPlus[P[a][b]];
}
if ( P[b][a] ) {
press.insert(P[b][a]);
if ( K2[b][a] < 0 )
++pressPlus[P[b][a]];
else
++pressMinus[P[b][a]];
}
}
Bin[*Ncns] = imbal;
for (set<int>::iterator ligi = ligs.begin(); ligi != ligs.end(); ++ligi) {
if ( ligsF[*ligi] != ligsB[*ligi] ) {
output << "In loop";
for (int i=0; i<Nloop; ++i)
output << " " << (loopstates[i]+1);
output << ": Ligand (signal " << *ligi
<< ") must influence the same number of rates forward and back." << endl;
return -1;
}
}
for (set<int>::iterator voli = volts.begin(); voli != volts.end(); ++voli) {
if ( (! voltsPlus[*voli]) || (! voltsMinus[*voli]) ) {
output << "In loop";
for (int i=0; i<Nloop; ++i)
output << " " << (loopstates[i]+1);
output << ": Voltage (signal " << *voli
<< ") must influence at least two rates,"
<< " and at least one must have opposite sign or direction." << endl;
return -2;
}
}
for (set<int>::iterator presi = press.begin(); presi != press.end(); ++presi) {
if ( (! pressPlus[*presi]) || (! pressMinus[*presi]) ) {
output << "In loop";
for (int i=0; i<Nloop; ++i)
output << " " << (loopstates[i]+1);
output << ": Pressure (signal " << *presi
<< ") must influence at least two rates,"
<< " and at least one must have opposite sign or direction." << endl;
return -2;
}
}
++*Ncns;
if ( volts.size() ) {
for (int j=0; j<Nk; ++j)
Ain[*Ncns][j] = 0.0;
for (int i=0; i<Nloop; ++i) {
int a = loopstates[i];
int b = loopstates[(i+1)%Nloop];
int k_ix1 = k_ix_of(Ns, a, b);
int k_ix2 = k_ix_of(Ns, b, a);
Ain[*Ncns][k_ix1 + Nk/3] = 1.0;
Ain[*Ncns][k_ix2 + Nk/3] = - 1.0;
}
Bin[*Ncns] = 0.0;
++*Ncns;
}
if ( press.size() ) {
for (int j=0; j<Nk; ++j)
Ain[*Ncns][j] = 0.0;
for (int i=0; i<Nloop; ++i) {
int a = loopstates[i];
int b = loopstates[(i+1)%Nloop];
int k_ix1 = k_ix_of(Ns, a, b);
int k_ix2 = k_ix_of(Ns, b, a);
Ain[*Ncns][k_ix1 + 2*Nk/3] = 1.0;
Ain[*Ncns][k_ix2 + 2*Nk/3] = - 1.0;
}
Bin[*Ncns] = 0.0;
++*Ncns;
}
return 0;
}
/*
|
|
To derive the transformation we will use singular value decomposition (svd) which is really well described at http://www.uwlax.edu/faculty/will/svd/index.html And derive some useful extras:
|
*/
void DeleteConstraint(int Nk, int *Ncns, double **Ain, double *Bin, int ic)
{
for (int jc=ic+1; jc<*Ncns; ++jc) {
memcpy(Ain[jc-1], Ain[jc], Nk*sizeof(double));
Bin[jc-1] = Bin[jc];
}
--*Ncns;
}
/*
|
|
Now we have K with 2*Nstate*Nstate elements, and a constraint system in which most of which are fixed.
We remove the fixed elements to speed up the constraint engine.
K[Kr_indices] = Kr kR_indices: length Nk constraint_indices: length Ncns */
extern "C" QUBFAST_API void IdentifyFixedParams_p(int Nk, int Ncns, double *K, double *Ain, double *Bin,
/*output: */ int *NkR, int *kR_indices, int *NcnsR, int *constraint_indices)
{
double **Ain_pp = new double*[Ncns];
for (int i=0; i
extern "C" QUBFAST_API void Kr_to_K(double *Kr, int NkR, double *K, int *kR_indices)
{
for (int i=0; i
#define PRETTY_SMALL 1e-5
/*
|