/* different from 2011_01_22:
No no_change: this means there is no enzyme limting --'no_change' reaction state. 
In fact, when no observation about no_change is given, those clauses will be automatically not used.

different from 2011_01_2:
Substrate Limiting is summarized into 3 cases. ?? Add the two Product noChange cases ??
*/


:- dynamic reaction_state/4.

:-['topCla'].
%:-['counting.pl'].
%:-['counting_includeLength'].
:-['b_noCycle'].



%set(depthLimit,s(s(s(s(s(s(s(s(0))))))))). %PianoNumDepthLimit):- numMap(6,PianoNumDepthLimit).
%set(depthLimit,s(0)).
set(depthLimit,s(s(0))).
%set(depthLimit,s(s(s(0)))).
%set(depthLimit,s(s(s(s(0))))).
%set(depthLimit,s(s(s(s(s(s(s(0)))))))). % depth 7
%set(depthLimit,s(s(s(s(s(0)))))).
%set(depthLimit,s(s(s(s(s(s(s(s(s(s(0))))))))))).
set(etoB,yes).








% it is actually minus 2 now
peanoMinus1(s(X),X).

% greater than
peanoGth(X,X).
peanoGth(s(Y),X):-peanoGth(Y,X).


/*********** using other examples ***********/
concentration(DepthLimit,Tr,Tr,VisitedNodes,X,Change,Time):-	
	set(depthLimit,DepthLimit0),
	DepthLimit \== DepthLimit0,
	ex(EI,concentration(X,Change,Time),1),!,
	ex0(EI,concentration(Xname,Change,Time),1),
	write('%***using this example***'), write({EI,Xname,Change}),nl.
	%NewVisitedNodes=[X|VisitedNodes], % since it's the last one, so it doesn't matter in terms of format
	%write('--- one candidate Hypothesis: '),portray_clause(VisitedNodes),pathwayDetails(VisitedNodes).


concentration(DepthLimit,Tr,Tr,VisitedNodes,X,Change,Time):- 
	set(depthLimit,DepthLimit0),
	DepthLimit \== DepthLimit0,
	threeStatesOpposite(Change,Opp1,Opp2),
	(ex(EI,concentration(X,Opp1,Time),1);ex(EI,concentration(X,Opp2,Time),1)),
	%write('!!! this metabolite is measured to be '), ex0(EI,concentration(Xname,ActualChange,Time),1),write({Xname,ActualChange}),nl,
	!,
	fail.


threeStatesOpposite(up,down,no_change).
threeStatesOpposite(down,up,no_change).
threeStatesOpposite(no_change,up,down).





/***************************** Regulation Rules ******************************/
/****************** Direct Effect ******************/

/*---------- Blocked/Slower down ----------*/
concentration(s(DepthLimit),TrSoFar,Tr,VisitedNodes,PID,down,Time):-
	produced_by(PID,ReactionID,'IMPD'),
	reaction_state(TrSoFar,Tr,ReactionID,enzymeLimiting,cataDecreased,Time).

concentration(s(CurrentDepth),TrSoFar,Tr,VisitedNodes,SID,up,Time):-				
	set(depthLimit,DepthLimit),peanoMinus1(DepthLimit,DepthLimit1),peanoGth(CurrentDepth,DepthLimit1),%write('succeed'),
	consumed_by(SID,ReactionID,'IMPD'),	
	reaction_state(TrSoFar,Tr,ReactionID,enzymeLimiting,cataDecreased,Time).


/*---------- Speeded Up ----------*/
concentration(s(DepthLimit),TrSoFar,Tr,VisitedNodes,PID,up,Time):-
	produced_by(PID,ReactionID,'IMPD'),
	reaction_state(TrSoFar,Tr,ReactionID,enzymeLimiting,cataIncreased,Time).

concentration(s(CurrentDepth),TrSoFar,Tr,VisitedNodes,SID,down,Time):-					
	set(depthLimit,DepthLimit),peanoMinus1(DepthLimit,DepthLimit1),peanoGth(CurrentDepth,DepthLimit1),%write('succeed'),
	consumed_by(SID,ReactionID,'IMPD'),	
	reaction_state(TrSoFar,Tr,ReactionID,enzymeLimiting,cataIncreased,Time).



/*---------- Catalytically No change ----------*/
/*concentration(s(DepthLimit),TrSoFar,Tr,VisitedNodes,PID,no_change,Time):-
	produced_by(PID,ReactionID,'IMPD'),
	reaction_state(TrSoFar,Tr,ReactionID,enzymeLimiting,cataNoChange,Time).

concentration(s(CurrentDepth),TrSoFar,Tr,VisitedNodes,SID,no_change,Time):-	
	set(depthLimit,DepthLimit),peanoMinus1(DepthLimit,DepthLimit1),peanoGth(CurrentDepth,DepthLimit1),%write('succeed'),
	consumed_by(SID,ReactionID,'IMPD'),	
	%produced_by(SID,ReactionID_Any,'IMPD'), % complement to the next one, just ensure that there is no redundant %% no matter there is sth to produce or 
	reaction_state(TrSoFar,Tr,ReactionID,enzymeLimiting,cataNoChange,Time).


% This rule doesn't largely increased the candidates, because it requites not produced_by any reactions.
concentration(s(CurrentDepth),TrSoFar,Tr,VisitedNodes,SID,no_change,Time):-					
	set(depthLimit,DepthLimit),peanoMinus1(DepthLimit,DepthLimit1),peanoGth(CurrentDepth,DepthLimit1),%write('succeed'),
	consumed_by(SID,ReactionID,'IMPD'),	
	\+produced_by(SID,ReactionID_Any,'IMPD'),
	reaction_state(TrSoFar,Tr,ReactionID,enzymeLimiting,cataDecreased,Time). */



/****************** Indirect Effect (Substrate Limiting) ******************/
/*------------------------- Substrate Up------------------------*/
concentration(s(DepthLimit),TrSoFar,Tr,VisitedNodes,PID,up,Time):-				
	produced_by(PID,ReactionID,'IMPD'),
	reaction_state(TrSoFar,Tr1,ReactionID,substrateLimiting,sUp,Time),
	consumed_by(SID,ReactionID,'IMPD'), \+partMember(SID,VisitedNodes),
	concentration(DepthLimit,Tr1,Tr,[{ReactionID,*}-PID|VisitedNodes],SID,up,Time).


/*------------------------- Substrate Down ------------------------*/
concentration(s(DepthLimit),TrSoFar,Tr,VisitedNodes,PID,down,Time):-				
	produced_by(PID,ReactionID,'IMPD'),
	reaction_state(TrSoFar,Tr1,ReactionID,substrateLimiting,sDown,Time),
	consumed_by(SID,ReactionID,'IMPD'), \+partMember(SID,VisitedNodes),
	concentration(DepthLimit,Tr1,Tr,[{ReactionID,*}-PID|VisitedNodes],SID,down,Time).



/*
concentration(s(DepthLimit),TrSoFar,Tr,VisitedNodes,PID,no_change,Time):-				
	produced_by(PID,ReactionID,'IMPD'),
	reaction_state(TrSoFar,Tr1,ReactionID,substrateLimiting,cataDecreased,Time),
	consumed_by(SID,ReactionID,'IMPD'),\+partMember(SID,VisitedNodes),
	concentration(DepthLimit,Tr1,Tr,[{ReactionID,cataDecreased}-PID|VisitedNodes],SID,up,Time).

concentration(s(DepthLimit),TrSoFar,Tr,VisitedNodes,PID,no_change,Time):-				
	produced_by(PID,ReactionID,'IMPD'),	
	reaction_state(TrSoFar,Tr1,ReactionID,substrateLimiting,cataIncreased,Time),
	consumed_by(SID,ReactionID,'IMPD'),\+partMember(SID,VisitedNodes),
	concentration(DepthLimit,Tr1,Tr,[{ReactionID,cataIncreased}-PID|VisitedNodes],SID,down,Time).
*/
% ??? What about substrate no_change?


partMember(M,Record):-
	maplist(wholeRecord_M,Record,Ms),
	member(M,Ms).

wholeRecord_M({ReactionID,Change}-PID,PID).
wholeRecord_M(SID-{ReactionID,Change},SID).






/*************** Abduce reaction_state	****************/
/* Simple version -- no incorporation of integrity constraint
reaction_state(TrSoFar,Tr,EC_Number,State):-
	Tr=[[reaction_state-{EC_Number,State}]|TrSoFar].
*/


% incorporate integrity constraint
% 1. one enzyme not to be hypothesized to several states: (1) the existing B; (2) the same 
% 2. the gene expression data. 

propotional(cataIncreased,up).
propotional(cataDecreased,down).
propotional(cataNoChange,noChange).

% This is not suitable when there are alternative similar rules 


% I don't think I really need Sign to be recorded, but to ensure to be used together with deriving hypotheses


reaction_state(TrSoFar,Tr,ReactionID,LimitingType,State,Time):-
	(appeared(ReactionID,TrSoFar,LimitingType2,State2)-> 
		State2==State, LimitingType2==LimitingType, % write('not contradiction'),nl,	% this will rarely happen -- unless it is going through a cycle, which is not allowed
		Tr=TrSoFar;
		Tr=[[rs-{ReactionID,LimitingType,State,Time}]|TrSoFar]	% Sign here will be left as a variable
	).


appeared(ReactionID,TrSoFar,LimitingType,State):-	
	member([rs-{ReactionID,LimitingType,State,Time}],TrSoFar).








p([reaction_state-{Sign,ReactionID,EC_Number,State,Time}],P):-
	pFact([reaction_state-{Sign,ReactionID,EC_Number,State,Time}],P),!.
	

p([reaction_state-{Sign,ReactionID,EC_Number,State,Time}],P):-
	concentration_enzClass(EC_Number,Change,Time),!,
	(propotional(State,Change)->
				P=0.8;  % whatever up/down/no_change, as long as they agree with each other
				P=0.1
			).

p([reaction_state-{Sign,ReactionID,EC_Number,State,Time}],1/3).




	

claLogPrior([reaction_state-{Sign,ReactionID,EC_Number,State,Time}],LogP):-
	concentration_enzClass(EC_Number,Change,Time),!,
	(propotional(State,Change)->
				LogP is log(0.8);  % whatever up/down/no_change, as long as they agree with each other
				LogP is log(0.1)
			).

claLogPrior([reaction_state-{Sign,ReactionID,EC_Number,State,Time}],LogP):- LogP is -log(3).