Commit 7c43f53f authored by FAGES Francois's avatar FAGES Francois
Browse files

merginin

Merge branch 'develop' of gitlab.inria.fr:lifeware/biocham into develop
parents d8c1208f 3707e7e2
%% Cell type:markdown id: tags:
 
# M7 Chemical SAT Solver
 
* Deciding the satisfiability of a Boolean formula in conjunctive normal form is NP-complete
* How can we program a chemical SAT solver ?
 
F. Fages, S. Soliman, Apr. 2020
 
%% Cell type:markdown id: tags:
 
# "Generate and Test" Algorithm with Stochastic CRN
* A stochastic CRN can be used to enumerate random Boolean values for a variable
 
%% Cell type:code id: tags:
 
```
parameter(k=1).
```
 
%% Cell type:code id: tags:
 
```
MA(k) for _/a => a. % reactions with inhibitors cannot fire in presence of the inhibitor, here a
```
 
%% Cell type:code id: tags:
 
```
MA(k) for a => _.
```
 
%% Cell type:code id: tags:
 
```
option(method:spn, stochastic_conversion:1).
```
 
%% Cell type:code id: tags:
 
```
numerical_simulation.
plot.
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%% Cell type:markdown id: tags:
 
## Question 1) Write a random generator of Boolean values for a vector of 3 variables a, b, c
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:markdown id: tags:
 
## Question 2) Write a CRN to find values satisfying the formula (x⋁¬y)⋀(y⋁¬z)⋀(z⋁¬x)
* Use an event to stop the simulation when the formula is satisfied
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:markdown id: tags:
 
# Guided Search Algorithm with Stochastic CRN
## Question 3) Improve your CRN to decrease the probability of moves of the variables that belong to satisfied clauses
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:markdown id: tags:
 
## Question ) Any idea to decide unsatisfiability ? statistical test question ?
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:markdown id: tags:
 
## Question 4) Determine the phase transition threshold in 3-SAT
* the density of a SAT instance is the ratio of the number of clauses divided by the number of variables
* a phase transition phenomenon is an asymptotic result showing the existence of a density threshold
* under the threshold the instances are almost surely satisfiable
* above the threshold the instances are almost surely unsatisfiable
* the hard instances are around the density threshold
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:markdown id: tags:
 
## Question ) Evaluation on 2-SAT and Horn-SAT ?
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:markdown id: tags:
 
# Guided Search Algorithm by Continuous CRN
 
We can see the SAT solving problem as a (continous) global optimization problem and try to solve it by gradient descent.
 
The idea is to find an energy function $E$ that is minimal when all clauses are satisfied, and then to simply enforce that $$\frac{dx}{dt} = -\frac{\partial E}{\partial x}$$
 
For a SAT problem involving variables $x_i\in [0, 1], 1\leq i\leq N$ and clauses $C_j, 1\leq j\leq M$ (with $C_{ji} = 1$ if $x_i$ appears positively in $C_j$, $C_{ji} = -1$ if $x_i$ appears negatively, and $0$ otherwise), we will define our energy function as a sum of squares of sub-energies for each clause.
 
$$E = \sum_{1\leq j\leq M}K_m^2$$
$$E = \sum_{1\leq j\leq M}K_j^2$$
 
 
%% Cell type:markdown id: tags:
 
## Question 5) Write $K_m$
## Question 5) Write $K_j$
 
Define (formally) $K_m$ as a function of the $C_ji$ and of the $x_i$, such that $K_m = 0$ iff clause $m$ is satisfied, and $K_m = 2^N$ if all $N$ variables appear in clause $m$ and are currently at the wrong value.
Define (formally) $K_j$ as a function of the $C_{ji}$ and of the $x_i$, such that $K_j = 0$ iff clause $j$ is satisfied, and $K_j = 2^N$ if all $N$ variables appear in clause $j$ and are currently at the _wrong_ value.
 
One might want to define $s_i\in[-1, 1]$ as a function of $x_i\in[0, 1]$ for ease of writing.
One might want to define $s_i\in[-1, 1]$ as a function of $x_i\in[0, 1]$ and then $K_j$ as a function of the $s_i$ for ease of writing.
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:markdown id: tags:
 
## Question 6) Obtain $dx_i/dt$
 
Obtain the formal expression for $\displaystyle\frac{dx}{dt}$ (if you have used $s_i$ just note that $\frac{\partial E}{\partial x}=\frac{\partial E}{\partial x}\frac{\partial s}{\partial x}$).
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:markdown id: tags:
 
## Question 7) Implement on the above example (x⋁¬y)⋀(y⋁¬z)⋀(z⋁¬x)
 
Using the commands:
`new_ode_system`, `init` (to set $x_1, x_2$ and $x_3$ initial state to 0.5), `ode_parameter` (to set the $c_ji$ corresponding to our 3 clauses), `ode_function` (for the $k_j$ and $s_i$) and `add_ode` (to add the above $dx_i/dt$)
 
Run a numerical simulation and plot the $x_i$ to see the result.
Experiment with different initial states, and variants of the problem (changing only the $c_ji$), what do you observe? [please leave all results and remarks *in* the notebook!]
 
%% Cell type:code id: tags:
 
```
clear_model.
new_ode_system.
init(x_1 = 0.5, x_2 = 0.5, x_3 = 0).
ode_function(
s_1 = 2*x_1 - 1, s_2 = 2*x_2 - 1, s_3 = 2*x_3 - 1,
k_1 = (1-c_11*s_1)*(1-c_12*s_2)*(1-c_13*s_3),
k_2 = (1-c_21*s_1)*(1-c_22*s_2)*(1-c_23*s_3),
k_3 = (1-c_31*s_1)*(1-c_32*s_2)*(1-c_33*s_3)
).
ode_parameter(
c_11 = 1, c_12 = -1, c_13 = 0,
c_21 = 0, c_22 = 1, c_23 = -1,
c_31 = -1, c_32 = 0, c_33 = 1
).
add_ode(
d(x_1)/dt = c_11*k_1*(1-c_12*s_2)*(1-c_13*s_3) + c_21*k_2*(1-c_22*s_2)*(1-c_23*s_3) + c_31*k_3*(1-c_32*s_2)*(1-c_33*s_3),
d(x_2)/dt = c_12*k_1*(1-c_11*s_1)*(1-c_13*s_3) + c_22*k_2*(1-c_21*s_1)*(1-c_23*s_3) + c_32*k_3*(1-c_31*s_1)*(1-c_33*s_3),
d(x_3)/dt = c_13*k_1*(1-c_12*s_2)*(1-c_11*s_1) + c_23*k_2*(1-c_22*s_2)*(1-c_21*s_1) + c_33*k_3*(1-c_32*s_2)*(1-c_31*s_1)
).
```
 
%%%% Output: execute_result
 
 
%% Cell type:code id: tags:
 
```
numerical_simulation.
plot(show: {x_1, x_2, x_3}).
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%% Cell type:code id: tags:
 
```
```
 
%%%% Output: display_data
 
\begin{align*}
{x_1}_0 &= 0.5\\
{x_2}_0 &= 0.5\\
{x_3}_0 &= 0.5\\
c_11 &= 1\\
c_12 &= -1\\
c_13 &= 0\\
c_21 &= 0\\
c_22 &= 1\\
c_23 &= -1\\
c_31 &= -1\\
c_32 &= 0\\
c_33 &= 1\\
s_1 &= 2*x_1-1\\
s_2 &= 2*x_2-1\\
s_3 &= 2*x_3-1\\
k_1 &= (1-c_11*s_1)*(1-c_12*s_2)*(1-c_13*s_3)\\
k_2 &= (1-c_21*s_1)*(1-c_22*s_2)*(1-c_23*s_3)\\
k_3 &= (1-c_31*s_1)*(1-c_32*s_2)*(1-c_33*s_3)\\
\frac{dx_1}{dt} &= c_11*k_1*(1-c_12*s_2)*(1-c_13*s_3)+c_21*k_2*(1-c_22*s_2)*(1-c_23*s_3)+c_31*k_3*(1-c_32*s_2)*(1-c_33*s_3)\\
\frac{dx_2}{dt} &= c_12*k_1*(1-c_11*s_1)*(1-c_13*s_3)+c_22*k_2*(1-c_21*s_1)*(1-c_23*s_3)+c_32*k_3*(1-c_31*s_1)*(1-c_33*s_3)\\
\frac{dx_3}{dt} &= c_13*k_1*(1-c_12*s_2)*(1-c_11*s_1)+c_23*k_2*(1-c_22*s_2)*(1-c_21*s_1)+c_33*k_3*(1-c_32*s_2)*(1-c_31*s_1)\\
\end{align*}
 
%% Cell type:markdown id: tags:
 
## Question 8) Improving the search
 
To avoid getting stuck in some local minima, we can add *Lagrange multipliers* $a_j, 1\leq j\leq M$.
To avoid getting stuck in some local minima, we can add *Lagrange multipliers* $a_j, 1\leq j\leq M$ so that the energy becomes $$E = \sum_{1\leq j\leq M}a_j K_j^2$$
These are new variables that will have an exponential increase proportional to $K_j$.
 
Add the 3 new variables and their ODEs with `add_ode`.
Do you notice any difference?
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:code id: tags:
 
```
```
 
%%%% Output: display_data
 
\begin{align*}
{x_1}_0 &= 0.5\\
{x_2}_0 &= 0.5\\
{x_3}_0 &= 0.5\\
{a_1}_0 &= 1\\
{a_2}_0 &= 1\\
{a_3}_0 &= 1\\
c_1_1 &= 0\\
c_1_2 &= 0\\
c_1_3 &= 0\\
c_2_1 &= 0\\
c_2_2 &= 0\\
c_2_3 &= 0\\
c_3_1 &= 0\\
c_3_2 &= 0\\
c_3_3 &= 0\\
s_1 &= 2*x_1-1\\
s_2 &= 2*x_2-1\\
s_3 &= 2*x_3-1\\
k_1 &= (1-c_1_1*s_1)*(1-c_1_2*s_2)*(1-c_1_3*s_3)\\
k_2 &= (1-c_2_1*s_1)*(1-c_2_2*s_2)*(1-c_2_3*s_3)\\
k_3 &= (1-c_3_1*s_1)*(1-c_3_2*s_2)*(1-c_3_3*s_3)\\
\frac{dx_1}{dt} &= a_1*c_1_1*k_1*((1-c_1_2*s_2)*(1-c_1_3*s_3))+a_2*c_2_1*k_2*((1-c_2_2*s_2)*(1-c_2_3*s_3))+a_3*c_3_1*k_3*((1-c_3_2*s_2)*(1-c_3_3*s_3))\\
\frac{dx_2}{dt} &= a_1*c_1_2*k_1*((1-c_1_1*s_1)*(1-c_1_3*s_3))+a_2*c_2_2*k_2*((1-c_2_1*s_1)*(1-c_2_3*s_3))+a_3*c_3_2*k_3*((1-c_3_1*s_1)*(1-c_3_3*s_3))\\
\frac{dx_3}{dt} &= a_1*c_1_3*k_1*((1-c_1_1*s_1)*(1-c_1_2*s_2))+a_2*c_2_3*k_2*((1-c_2_1*s_1)*(1-c_2_2*s_2))+a_3*c_3_3*k_3*((1-c_3_1*s_1)*(1-c_3_2*s_2))\\
\frac{da_1}{dt} &= a_1*k_1\\
\frac{da_2}{dt} &= a_2*k_2\\
\frac{da_3}{dt} &= a_3*k_3\\
\end{align*}
 
%% Cell type:markdown id: tags:
 
link to M4 chemical_logic
 
%% Cell type:markdown id: tags:
 
## Question ) on the link to the first Modal session below ?
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:code id: tags:
 
```
clear_model.
```
 
%% Cell type:code id: tags:
 
```
present(A,a). present(B, b).
```
 
%% Cell type:code id: tags:
 
```
parameter(a=1, b=1).
```
 
%% Cell type:markdown id: tags:
 
# And
$A\wedge B=A*B$
 
%% Cell type:code id: tags:
 
```
A+B=>A+B+AandB.
```
 
%% Cell type:code id: tags:
 
```
AandB=>_.
```
 
%% Cell type:code id: tags:
 
```
list_ode.
```
 
%%%% Output: display_data
 
\begin{align*}
{AandB}_0 &= 0\\
{A}_0 &= 1\\
{B}_0 &= 1\\
a &= 1\\
b &= 1\\
\frac{dAandB}{dt} &= A*B-AandB\\
\frac{dA}{dt} &= 0\\
\frac{dB}{dt} &= 0\\
\end{align*}
 
%% Cell type:code id: tags:
 
```
%slider a b
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%% Cell type:markdown id: tags:
 
# Or
$A\vee B = A+B-A*B$
 
%% Cell type:code id: tags:
 
```
A=>A+AorB.
```
 
%% Cell type:code id: tags:
 
```
B=>B+AorB.
```
 
%% Cell type:code id: tags:
 
```
A+B+AorB=>A+B.
```
 
%% Cell type:code id: tags:
 
```
AorB=>_.
```
 
%% Cell type:code id: tags:
 
```
%slider a b
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%% Cell type:markdown id: tags:
 
# Not
$\neg A=1-A$
 
%% Cell type:code id: tags:
 
```
_ => notA.
```
 
%% Cell type:code id: tags:
 
```
notA => _.
```
 
%% Cell type:code id: tags:
 
```
A+ notA => A.
```
 
%% Cell type:code id: tags:
 
```
list_ode.
```
 
%%%% Output: display_data
 
\begin{align*}
{notA}_0 &= 0\\
{A}_0 &= 0.0\\
a &= 0.0\\
\frac{dnotA}{dt} &= 1-notA*A-notA\\
\frac{dA}{dt} &= 0\\
\end{align*}
 
%% Cell type:code id: tags:
 
```
%slider a
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%% Cell type:code id: tags:
 
```
```
 
%% Cell type:code id: tags:
 
```
clear_model.
```
 
%% Cell type:code id: tags:
 
```
present(A,a). parameter(a=1).
```
 
%% Cell type:code id: tags:
 
```
present(B, b). parameter(b=1).
```
 
%% Cell type:code id: tags:
 
```
A=>D.
```
 
%% Cell type:code id: tags:
 
```
B=>D.
```
 
%% Cell type:code id: tags:
 
```
MA(k) for A+B=>_. parameter(k=2.5).
```
 
%% Cell type:code id: tags:
 
```
list_ode.
```
 
%%%% Output: display_data
 
\begin{align*}
{B}_0 &= 1\\
{A}_0 &= 1\\
{D}_0 &= 0\\
a &= 1\\
b &= 1\\
k &= 2.5\\
\frac{dB}{dt} &= -B-k*A*B\\
\frac{dA}{dt} &= -A-k*A*B\\
\frac{dD}{dt} &= A+B\\
\end{align*}
 
%% Cell type:code id: tags:
 
```
%slider a b k
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%% Cell type:code id: tags:
 
```
validity_domain(F(G(D=v))).
```
 
%%%% Output: execute_result
 
v=1.00221
 
%% Cell type:code id: tags:
 
```
```
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment