#---------------------------------------------------------
# Code for Macaulay2

-- Define a ring of polynomials in 9 variables.

R = ZZ/2[x1,x2,x3,x4,x5,x6,x7,x8,x9]; 

-- Define a quotient ring, where each x_i^2 = x_i.

I = ideal(x1^2-x1, x2^2-x2, x3^2-x3, x4^2-x4, x5^2-x5, x6^2-x6, x7^2-x7, x8^2-x8, x9^2-x9);
Q = R / I;

-- Shortcut for AND and OR functions.

RingElement | RingElement :=(x,y)->x+y+x*y;
RingElement & RingElement :=(x,y)->x*y;

-- Set the parameters (constants).

Ge = 0_Q
Le = 1_Q

-- This is the 9-variable lac operon model.

f1 = (1+x5) & x4;
f2 = x1;
f3 = x1;
f4 = 1+Ge;
f5 = (1+x6) & (1+x7);
f6 = x8 & x3;
f7 = x6 | x8 | x9;
f8 = (1+Ge) & x2 & Le;
f9 = (1+Ge) & (x8 | Le);

(f1,f2,f3,f4,f5,f6,f7,f8,f9)

-- Compute the ideal to find the fixed point(s).

I = ideal(f1+x1, f2+x2, f3+x3, f4+x4, f5+x5, f6+x6, f7+x7, f8+x8, f9+x9)

-- Compute a Groebner basis.

G = gens gb I

#---------------------------------------------------------
# Code for Cyclone in AlgoRun (if Ge=0, Le=1)

NUMBER OF VARIABLES: 9
NUMBER OF STATES :   2

x1 = x4*x5+x4
x2 = x1
x3 = x1
x4 = 1
x5 = x6*x7+x6+x7+1
x6 = x3*x8
x7 = x6*x8*x9+x6*x8+x6*x9+x8*x9+x6+x8+x9
x8 = x2
x9 = 1 

#---------------------------------------------------------
# Code for BoolNet in RStudio

targets, factors
M, !R & C
P, M
B, M
C, !Ge
R, !A & !Am
A, L & B
Am, A | L | Lm
L, !Ge & P & Le
Lm, !Ge & (L | Le)


library(Boolnet)
lacModel <- loadNetwork("lac-operon")
lacAttractors <- getAttractors(lacModel)
plotStateGraph(lacAttractors)