/*

MIT License

Copyright (c) 2020 Moreno Falaschi, Giulia Palma, Linda Brodo, Roberto Bruni

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NON INFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

*/

%==== This initial part allows to define the 'classical' mechanism to compute in a Reaction System (RS) Framework

/*
We have the following lists:
R is a list of reagents (a set of constants);
I is a list of inhibitors;
P is a list of products;
T = list of entities which represent the current computation state
*/


/*
allIn(R,T):-“checks if all the elements of the list R are
present in the list T (R is a subset of T)”.
*/

allIn([],_).
allIn([X|L1],L2):- member(X,L2), allIn(L1,L2).

/*
both_lists(S,S1,SO) computes the list SO as intersection of the input lists S and S1.
empty_inters(S,S1) checks if the input list S and S1 have an empty intersection.
*/

both_lists([],_,[]).
both_lists([X|L1],L2,[X|L3]):- member(X,L2),!, both_lists(L1,L2,L3).
both_lists([_|L1],L2,L3):- both_lists(L1,L2,L3).
empty_inters(I,T):- both_lists(I,T,L), L=[].

/*
A list is non empty if its length is greater or equal to 1.
*/

non_empty(L):- length(L,N), N>=1.

/*
The predicate enable(R,I,T) checks if a reaction can take place (is enabled) in T.
A reaction can take place iff all reactants in R belong to T
and no inhibitor I belongs to T.
*/

enable(R,I,T):- allIn(R,T), empty_inters(I,T).

/*
result(T,R,I,P,P1):- "given the current state T and the reaction (R,I,P), P1 is computed by result/5"
P1 will be P if enable(R,I,T) is true, otherwise it will be the empty list [].
*/

result(T,R,I,P,P):- enable(R,I,T),!.
result(_,_,_,_,[]).

/* The predicate reactionset/1 takes as input the list of triples which defines the set of reactions */

%%%E.g. reactionset([([lac],[ a],[cya]),([lacI],[a ],[lac2]),([lac2],[a ],[lac3]),([cya],[ a],[cya2])]).


/* the predicate 'unlimitedComputation/2' makes a computation in a Reaction System, starting from the initial state,
and returns a list of states, by taking into account the contribution of the context.
The computation can be finite, and stops as soon as an empty state is computed. It can also be unlimited,
and hence it will not return any value, if an empty state is never encountered.

The user can choose to execute the predicate 'computationLimitedToKSteps', and in this case the maximum
number of steps (say K) is requested initially to the user, and the computation will end when
the state becomes empty or when K steps have been executed.
*/

%please notice that a context sequence is a list of (context) lists. Hence the empty initial context sequence is denoted by C = [ [] ]

context([[a1,a2],[a3,a2,a5]]).


/* the predicate computation/2 proposes a choice to the user between an unlimited computation or a computation limited to a maximum of K steps.
*/

computation(InitialState,L):-
   reactionSet(R),preliminaryCheck(R),
   write('Do you want a computation of possibly unlimited length? (say yes or no followed by a dot) '),
   read(Answer),selectComp(Answer,InitialState,L).

preliminaryCheck([]).
preliminaryCheck([(R,_,_)|_]):- R==[], write('Error: the set of reagents in a rule cannot be empty. Please correct'),!, fail.
preliminaryCheck([(_,I,_)|_]):- I==[], write('Error: the set of inhibitors in a rule cannot be empty. Please correct'),!, fail.
preliminaryCheck([(R,I,_)|_]):- both_lists(R,I,Intersection), Intersection\==[], write('Error: the intersection of reagents and inhibitors in a rule is not empty. Please correct'),!, fail.
preliminaryCheck([(_,_,_)|OtherReacts]):-preliminaryCheck(OtherReacts).

selectComp(Answer,InitialState,L):- Answer=='yes',unlimitedComputation(InitialState,L).
selectComp(Answer,InitialState,L):- Answer\=='yes', computationLimitedToKSteps(InitialState,L).

/* the predicate unlimitedComputation/2 will stop only in case an empty state is encountered in the computation */

unlimitedComputation(InitialState,L):-reactionSet(R),context([C0|Cs]),union(InitialState,C0,SC),computeWithContext(SC,Cs,R,L).

computationLimitedToKSteps(InitialState,L):-reactionSet(R),context([C0|Cs]),union(InitialState,C0,SC),
     write('Give me the maximun number of computation steps (a positive integer, followed by a dot) '),read(MaxSteps),
     computeWithContextKSteps(MaxSteps,SC,Cs,R ,L).



/* the predicate union(S,S1,SO) computes the list SO as the set union of the lists S and S1 */
union([],L,L).
union([X|L],L2,L3):-member(X,L2), ! , union(L,L2,L3).
union([X|L],L2,[X|L3]):- union(L,L2,L3).

/* computeWithContext(ComputationState,Context,Reactions,ComputationStateSequence)
takes in input the current ComputationState, the Context sequence, the list of Reactions in the Reaction System,
and returns the computed ComputationStateSequence. */

computeWithContext([],_,_,[]).
computeWithContext([X|L],[],R,[S1|S]):- resultallreactions([X|L],R,S1), computeWithContext(S1,[],R,S).
computeWithContext([X|L],[C|Cs],R,[S1|S]):- resultallreactions([X|L],R,S1),union(S1,C,S2),computeWithContext(S2,Cs,R,S).

/* resultallreactions(CurrentState,ReactionSet,NewState) applies all reactions in ReactionSet to CurrentState and computes the NewState */

resultallreactions(_,[],[]).
resultallreactions(T,[(R,I,P)|OtherReacts],T1):-result(T,R,I,P,P1),resultallreactions(T,OtherReacts,T2),append(P1,T2,T1).

/* computeWithContextKSteps(MaxSteps,State,Context,ReactionSet,FullComputation) it makes the same computation as
computeWithContext(State,Context,ReactionSet,FullComputation). However, the additional parameter 'MaxSteps' allows to
stop the computation after at most MaxSteps steps. The context is used in the 4th clause. When the context is empty
the first three clauses are used.
*/

computeWithContextKSteps(K,[],_,_,[]):-K>0.
computeWithContextKSteps(K,_,_,_,[]):-K==0.
computeWithContextKSteps(K,[X|L],[],R,[S1|S]):- K>0,resultallreactions([X|L],R,S1), K1 is K-1,
                                                computeWithContextKSteps(K1,S1,[],R,S).
computeWithContextKSteps(K,[X|L],[C|Cs],R,[S1|S]):- resultallreactions([X|L],R,S1),union(S1,C,S2),
                                                    K1 is K-1, computeWithContextKSteps(K1,S2,Cs,R,S).


%=================== NEW STUFF FOR PAPER ON SOS RULES =================

/* we recall the syntax for context processes:
nil stop
rec(X) recursive invocation of constant X
pre(C,K) make C available then behaves as K
plus(K1,K2) chooses between K1 and K2
*/

/* find(X,Delta,K) returns the context process K associated with the contant X in the environment Delta
we recall that an environment Delta is a list of (possibly recursive) context declarations def(X,K)
*/
find(_,[],[]).
find(X,[def(X,K)|_],K) :- !.
find(X,[def(_,_)|Ds],K) :- find(X,Ds,K).

/* unfold(Delta,K,Choices) returns the list of choices for the context K given the process definitions Delta
Choices is a list of context moves choice(C,K) where C is a set of entities and K is the continuation
unfold can be applied to single context process K or to lists of context process Ks
*/
% handle single process
unfold(_,nil,[]).
unfold(Delta,rec(X),Choices) :- find(X,Delta,K), unfold(Delta,K,Choices).
unfold(_,pre(C,K),[choice(C,[K])]).
unfold(Delta,plus(K1,K2),Choices) :- !,unfold(Delta,K1,Choices1),unfold(Delta,K2,Choices2),union(Choices1,Choices2,Choices).
% handle lists of processes
unfold(_,[],[choice([],[])]).
unfold(Delta,[K|Ks],Choices) :- unfold(Delta,K,CK),unfold(Delta,Ks,CKs),shuffle(CK,CKs,Choices).

/* shuffle(CKs1,CKs2,Choices) given two lists of choices Cks1 and Cks2,
returns the list Choices of all possible combinations of the choices in CKs1 with those in CKs2
shuffle can be applied to single choices or to lists of choices
*/
% handle lists of choices
shuffle([],_,[]).
shuffle([choice(C,K)|CK],CKs,Choices) :- shuffle(choice(C,K),CKs,Ch1),
                                         shuffle(CK,CKs,Ch2),
                                         append(Ch1,Ch2,Choices).
% handle single choices
shuffle(choice(_,_),[],[]).
shuffle(choice(C1,K1),[choice(C2,K2)|CKs],[choice(C,K)|Choices]) :- union(C1,C2,C),
                                                                    append(K1,K2,K),
                                                                    shuffle(choice(C1,K1),CKs,Choices).
/* a system process S is sys(Delta,E,Ks,Rs) consists of
- an environment Delta (a list of constant declarations def(X,K))
- the set E of currently available entities (produced from the previous step)
- the list Ks of context processes
- the list Rs of reaction rules react(R,I,P) (where R,I,P are list of entities)
oneTransition(S,L,S1) holds if the system S has one transition with label L to system S1
when S=sys(Delta,E,Ks,Rs)
- the label L has the form obs(T,R,RI,I,IR,P)
- the target system process has the form sys(Delta,P,Ks',Rs)
i.e. Delta and Rs are not changed and the product set P in the label becomes the set of available entities
*/
oneTransition(sys(Delta,E,Ks,Rs),L,S) :- unfold(Delta,Ks,Choices),
                                         oneTransition(Delta,E,Choices,Rs,[tr(L,S)]).
oneTransition(Delta,E,[choice(C,K)|_],Rs,M) :- union(E,C,T),transition(Delta,T,K,Rs,M).
oneTransition(Delta,E,[_|Cs],Rs,M) :- oneTransition(Delta,E,Cs,Rs,M).


/* allTransitions(S,Moves) returns the list Moves of all possible transitions tr(obs(...),S1) of system S
*/
allTransitions(sys(Delta,E,Ks,Rs),Moves) :- unfold(Delta,Ks,Choices),
                                            allTransitions(Delta,E,Choices,Rs,Moves).
allTransitions(_,_,[],_,[]).
allTransitions(Delta,E,[choice(C,K)|Cs],Rs,Moves) :- union(E,C,T),transition(Delta,T,K,Rs,M),
                                                     allTransitions(Delta,E,Cs,Rs,Ms),
                                                     union(M,Ms,Moves).

/* transition(Delta,T,K,Rs,Move) returns the unique possible transition Move=tr(L,S1) 
given the current Delta, available entities T, context continuation K and reactions Rs
*/                                                  
transition(Delta,T,K,Rs,[tr(obs(T,R,RI,I,IR,P),sys(Delta,P,K,Rs))]) :- result(T,Rs,obs(T,R,RI,I,IR,P)).

/* result(T,Rs,L) returns the label originating from the execution of reactions Rs with available entities T
*/
result(T,[],obs(T,[],[],[],[],[])).
result(T,[react(R,I,P)|Rs],L):- enable(R,I,T),!,
                                result(T,Rs,L2),
                                unionObs(obs(T,R,[],I,[],P),L2,L).
result(T,[react(R,I,_)|Rs],L):- both_lists(I,T,RI),minusSet(R,T,IR),
                                result(T,Rs,L2),
                                unionObs(obs(T,[],RI,[],IR,[]),L2,L).

/* unionObs(L1,L2,L) combines the labels L1 and L2 into L
*/
unionObs(obs(T,R1,RI1,I1,IR1,P1),obs(T,R2,RI2,I2,IR2,P2),obs(T,R,RI,I,IR,P)) :- union(R1,R2,R),
                                                                                union(RI1,RI2,RI),
                                                                                union(I1,I2,I),
                                                                                union(IR1,IR2,IR),
                                                                                append(P1,P2,P).


%==== BioHML

/* checkAssertion(L,F) holds if the label L satisfies the assertion F
*/
checkAssertion(_,true).
checkAssertion(L,sub(S,N)) :- selectSet(N,L,X), allIn(S,X).
checkAssertion(L,nonempty(N)) :- selectSet(N,L,X), non_empty(X).
checkAssertion(L,and(F1,F2)) :- checkAssertion(L,F1),checkAssertion(L,F2).
checkAssertion(L,or(F1,_)) :- checkAssertion(L,F1),!.
checkAssertion(L,or(_,F2)) :- checkAssertion(L,F2).
checkAssertion(L,xor(F1,F2)) :- checkAssertion(L,F1), \+ checkAssertion(L,F2).
checkAssertion(L,xor(F1,F2)) :- checkAssertion(L,F2), \+ checkAssertion(L,F1).
checkAssertion(L,not(F)) :- \+ checkAssertion(L,F).

/* select(N,L,E) returns the set of entities extracted from the label L according to the index N
*/
selectSet(1,obs(T,_,_,_,_,_),T).
selectSet(2,obs(_,R,RI,_,_,_),U):-union(R,RI,U).
selectSet(3,obs(_,_,_,I,IR,_),U):-union(I,IR,U).
selectSet(4,obs(_,_,_,_,_,P),P).


/* checkBioHML(S,G,B) returns B=ok if the system S satisfies the BioHML formula G
returns an explanation B of the reasion why G is not satisfied otherwise
*/
checkBioHML(_,true,ok).
checkBioHML(S,false,no(S,false)).
checkBioHML(S,and(G1,_),B) :- checkBioHML(S,G1,B), B\==ok, !.
checkBioHML(S,and(_,G2),B) :- checkBioHML(S,G2,B), B\==ok, !.
checkBioHML(_,and(_,_),ok).
checkBioHML(S,or(G1,_),ok) :- checkBioHML(S,G1,ok),!.
checkBioHML(S,or(_,G2),ok) :- checkBioHML(S,G2,ok),!.
checkBioHML(S,or(G1,G2),no(S,[B1,B2])) :- checkBioHML(S,G1,B1),checkBioHML(S,G2,B2).
checkBioHML(S,diamond(F,G),B) :- allTransitions(S,Moves),filter(Moves,F,Ms),checkOne(S,F,Ms,G,B).
checkBioHML(S,box(F,G),B) :- allTransitions(S,Moves),filter(Moves,F,Ms),checkAll(S,Ms,G,B).

/* filter(Moves,F,Ms) returns the list Ms of all the moves tr(L,S) in Moves such that L satisfies F
*/
filter([],_,[]).
filter([tr(L,S)|Moves],F,[tr(L,S)|Ms]) :- checkAssertion(L,F),!,filter(Moves,F,Ms).
filter([tr(_,_)|Moves],F,Ms) :- filter(Moves,F,Ms).

/* checkOne(S,F,Ms,G,B) checks if there is one transition tr(L,S1) in the list Ms such that S1 satisfies G
if this is the case B=ok, otherwise B explains why the check fails
S and F are needed to build the explanation when there are no moves available
*/
checkOne(S,F,[],_,no(S,miss(F))).
checkOne(_,_,[tr(_,S1)|_],G,ok) :- checkBioHML(S1,G,ok),!.
checkOne(S,F,[tr(_,_)|Ms],G,B) :- checkOne(S,F,Ms,G,B).

/* checkAll(S,Ms,G,B) checks if all transitions tr(L,S1) in the list Ms are such that S1 satisfies G
if this is the case B=ok, otherwise B explains why the check fails
S is needed to build the explanation when there is one transitions that violates the property
+*/
checkAll(_,[],_,ok).
checkAll(S,[tr(L,S1)|_],G,no(S,has(L,B))) :- checkBioHML(S1,G,B), B\==ok, !.
checkAll(S,[tr(_,_)|Ms],G,B) :- checkAll(S,Ms,G,B).


%==== simplified syntax for system processes, labels and contexts



resultallreactionsLabel(_,[],[],[[],[],[],[]]).
resultallreactionsLabel(T,[(R,I,P)|OtherReacts],T1,LabelOut):- resultLabel(T,R,I,P,P1,Label1),
                           resultallreactionsLabel(T,OtherReacts,T2,Label2),append(P1,T2,T1),unionLabel(Label1,Label2,LabelOut).

unionLabel([],[],[]).
unionLabel([L1|L1s],[L2|L2s],[LO|LOs]):- union(L1,L2,LO),unionLabel(L1s,L2s,LOs).

/*
resultLabel(T,R,I,P,P1,Label):- "given the current state T and the reaction (R,I,P), P1 and Label are computed by resultLabel/6, applying reaction (R,I,P)"
P1 will be P if enable(R,I,T) is true, otherwise it will be the empty list [].
*/

resultLabel(T,R,I,P,P,[R,[],I,[]]):- enable(R,I,T),!.
resultLabel(T,R,I,_,[], [[],RI,[],IR] ):- both_lists(I,R,RI),minusSet(R,T,IR).

/* minusSet(S,S1,SO) :- "SO is the list of the elements in the set S-S1 (notice that it may contain repeated elements if S1 contains repeated elements)"
*/

minusSet([],_,[]).
minusSet([X|R],T,Rs):- member(X,T),!,minusSet(R,T,Rs).
minusSet([X|R],T,[X|Rs]):- minusSet(R,T,Rs).

/* verifyLabel(E,[R,RI,I,IR],P,Assertion):-"it verifies if the label E,[R,RI,I,IR],P |= Assertion. "
verifyLabel/4 follows the structure of the Assertion Language in the paper on SOS rules for Reaction Systems (see Definition 10) */


verifyLabel(E,Label,P,sub(S,N)):- selectLabel(N,E,Label,P,Selected), allIn(S,Selected).
verifyLabel(E,Label,P,nonempty(N)):- selectLabel(N,E,Label,P,Selected), non_empty(Selected).
verifyLabel(E,Label,P,and(F1,F2)):- verifyLabel(E,Label,P,F1), verifyLabel(E,Label,P,F2).
verifyLabel(E,Label,P,or(F1,_)):- verifyLabel(E,Label,P,F1).
verifyLabel(E,Label,P,or(_,F2)):-  verifyLabel(E,Label,P,F2).
verifyLabel(E,Label,P,xor(F1,F2)):-  verifyLabel(E,Label,P,F1), \+ verifyLabel(E,Label,P,F2).
verifyLabel(E,Label,P,xor(F1,F2)):-  \+ verifyLabel(E,Label,P,F1), verifyLabel(E,Label,P,F2).
verifyLabel(E,Label,P,not(F)):-  \+ verifyLabel(E,Label,P,F).

/* selectLabel(X,E, [R,RI,I,IR], P, L) takes as input X an integer between 1 and 4
and returns L=E if X==1, L=R \cup RI if X=2, L=I \cup IR if X=3, L=P if X=4
*/
selectLabel(1,E,[_,_,_,_],_,E).
selectLabel(2,_,[R,RI,_,_],_,U):-union(R,RI,U).
selectLabel(3,_,[_,_,I,IR],_,U):-union(I,IR,U).
selectLabel(4,_,[_,_,_,_],P,P).


/*  verifyBioHML(CurrentState,ReactionSet,bioHMLFormula,ContextSequence,VerifiedFormulaSequence)
checks if the bioHML formula (represented as a list of subformulas -- see definition of bioHML formulas)
is verified in the current Reaction System framework.
In each computation step one element of the 'bioHMLFormula' list is checked w.r.t. the current computation step */

verifyBioHML(_,_,[t|_],_,[],[]).
verifyBioHML([],_,[f|_],_,[],[]).
verifyBioHML(S,R,[box(F)|_],plus([A|_],[B|_]),[box(F)],[S1]):-union(S,A,S1),computeNewStateAndLabel(S1,R,SO,Label1),
                                             \+ verifyLabel(S1,Label1,SO,F),
                                             union(S,B,S2),computeNewStateAndLabel(S2,R,SO1,Label2),
                                             \+ verifyLabel(S2,Label2,SO1,F),!,true.
verifyBioHML(S,R,[box(F)|As],plus([A|A1],[B|_]),[box(F)|VF],[S1|Ss]):-union(S,A,S1),computeNewStateAndLabel(S1,R,SO1,Label1),
                     verifyLabel(S1,Label1,SO1,F),
                     union(S,B,S2),computeNewStateAndLabel(S2,R,SO2,Label2),
                     \+ verifyLabel(S2,Label2,SO2,F),!,
                     verifyBioHML(SO1,R,As,A1,VF,Ss).
verifyBioHML(S,R,[box(F)|As],plus([A|_],[B|B1]),[box(F)|VF],[S2|Ss]):-
                     union(S,A,S1),computeNewStateAndLabel(S1,R,SO1,Label1),\+ verifyLabel(S1,Label1,SO1,F),
                     union(S,B,S2),computeNewStateAndLabel(S2,R,SO2,Label2),verifyLabel(S2,Label2,SO2,F),!,
                     verifyBioHML(SO2,R,As,B1,VF,Ss).
verifyBioHML(S,R,[box(F)|As],plus([A|A1],[B|B1]),[box(F)|VF],SsO):-
                     union(S,A,S1),computeNewStateAndLabel(S1,R,SO1,_), verifyBioHML(SO1,R,As,A1,VF,Ss1),
                     union(S,B,S2),computeNewStateAndLabel(S2,R,SO2,_),
                     verifyBioHML(SO2,R,As,B1,VF,Ss2),append([S1|Ss1],[S2|Ss2],SsO).
verifyBioHML(S,R,[box(F)|_],[C|_],[box(F)],[S1]):- C\==plus(_,_),union(S,C,S1),
                     computeNewStateAndLabel(S1,R,SO,Label),\+ verifyLabel(S1,Label,SO,F),!.
verifyBioHML(S,R,[box(F)|As],[C|C1],[box(F)|VF],[S1|Ss]):- C\==plus(_,_),union(S,C,S1),computeNewStateAndLabel(S1,R,SO,Label),
                     verifyLabel(S1,Label,SO,F),verifyBioHML(SO,R,As,C1,VF,Ss).
verifyBioHML(S,R,[diamond(F)|_],plus([A|_],[B|_]),[diamond(F)],[S1]):-union(S,A,S1),computeNewStateAndLabel(S1,R,SO1,Label1),
                     \+ verifyLabel(S1,Label1,SO1,F), union(S,B,S2),
                     computeNewStateAndLabel(S2,R,SO2,Label2), \+ verifyLabel(S2,Label2,SO2,F),!,
                     fail.

verifyBioHML(S,R,[diamond(F)|As],plus([A|A1],[_|_]),[diamond(F)|VF],[S1|Ss]):-
          union(S,A,S1),computeNewStateAndLabel(S1,R,SO,Label),verifyLabel(S1,Label,SO,F),
          verifyBioHML(SO,R,As,A1,VF,Ss).

verifyBioHML(S,R,[diamond(F)|As],plus([_|_],[B|B1]),[diamond(F)|VF],[S1|Ss]):-
          union(S,B,S1),computeNewStateAndLabel(S1,R,SO,Label),verifyLabel(S1,Label,SO,F),verifyBioHML(SO,R,As,B1,VF,Ss).

verifyBioHML(S,R,[diamond(F)|As],[C|C1],[diamond(F)|VF],[S1|Ss]):- %C\==plus(A,B),
          union(S,C,S1),computeNewStateAndLabel(S1,R,SO,Label),
          verifyLabel(S1,Label,SO,F),
          verifyBioHML(SO,R,As,C1,VF,Ss).

verifyBioHML(S,R,[or(A,_)],C,[or(A,_)|VF],[S|Ss]):-verifyBioHML(S,R,A,C,VF,Ss).
verifyBioHML(S,R,[or(_,B)],C,[or(_,B)|VF],[S|Ss]):-verifyBioHML(S,R,B,C,VF,Ss).
verifyBioHML(S,R,[and(A,B)],C,[and(A,B)|VF],[S|SsO]):-verifyBioHML(S,R,A,C,VF1,Ss1),verifyBioHML(S,R,B,C,VF2,Ss2),append(VF1,VF2,VF),append(Ss1,Ss2,SsO).


computeNewStateAndLabel(S1,R,SO,Label):- resultallreactionsLabel(S1,R,SO,Label).





main(Answer) :- myenvironment(Delta),
                myentities(E),
                mycontext(Ks),
                myreactions(Rs),
                mybhml(G),
                checkBioHML(sys(Delta,E,Ks,Rs),G,Answer).


/* examples from paper */
/* default example */
myenvironment([]).
myentities([]).
myreactions([react([a,b],[c],[b])]).
mycontext([plus(pre([a, b], pre([a], pre([a, c], nil))), pre([a, b], pre([a], pre([a], nil))))]) .

mybhml(diamond(not(sub([c], 1)), box(not(sub([c], 1)), diamond(not(sub([c], 1)), true)))).