% =============================================================================
% Automated modeling in process-based qualitative reasoning
%  ** Module for generating a model based on naive dependencies
% 
% April - July 2007
% 
% Hylke Buisman   - hbuisman@science.uva.nl
% =============================================================================

%==============================================================================
% FIND NAIVE DEPENDENCIES
%==============================================================================

% Find naive dependencies
% find_naive_dependencies/0
% find_naive_dependencies.
find_naive_dependencies :-
	retractall(model_dependency(_,_)),
	retractall(dependency_bag(_,_,_)),
	retractall(dependency_group(_,_)),
	retractall(naive_dependency(_,_,_)),

	% Step 1: Find naive dependencies
	generate_quantity_pairs(Pairs),
	findall(S, engine:state(S, _), StateIDs),
	
	find_simple_prop_inf(Pairs, StateIDs, D),
	list_to_set(D, ND),

	assert_model(model(ND)).
	
	
% Find simple influence and proportionality dependencies (naive dependencies)
% find_simple_prop_inf/3
% find_simple_prop_inf(QPairs, States, SimpleModel)
find_simple_prop_inf([], _, []).
find_simple_prop_inf([QPair|RestQPs], States, NDependencies) :-
	find_simple_prop_inf(RestQPs, States, Dependencies),
	writef(' > finding naive dependencies: Pair %q\n', [QPair]),
	flush,
	get_minimal_prop_inf(QPair, States, D),
	append(D, Dependencies, NDependencies).

% Returns the dependencies that are consistent
% in all states
% get_minimal_prop_inf/3
% get_minimal_prop_inf(+QuantityPair, +States, -Dependencies)
get_minimal_prop_inf(QPair, [S], Dependencies) :-
	get_consistent_dependencies(S, QPair, Dependencies).
get_minimal_prop_inf(QPair, [S|Rest], NDependencies) :-
	get_minimal_prop_inf(QPair, Rest, Dependencies),
	get_consistent_dependencies(S, QPair, Dependencies2),
	intersection(Dependencies, Dependencies2, NDependencies).


	
%==============================================================================
% FINDING CALCULI
%==============================================================================
	
% Return triples that possibly have a min relation 
% (based on all naive dependencies)
% min_candidates/2
% min_candidates(+NaiveModel, -MinCandidates)
min_candidates(naive_model(Clusters, _, Dependencies), MinCandidates) :-
	 findall(model_dependency(equal, [C,min(A,B)]), 
		(
		 model_dependency(prop_pos, [C, A]), 
		 model_dependency(prop_neg, [C, B]),
		 % The feedback cannot be from a cluster to the same cluster
		 get_quantity_cluster(A, Clusters, CA),
		 get_quantity_cluster(B, Clusters, CB),
		 get_quantity_cluster(C, Clusters, CC),
		 CA \= CB, CA \= CC, CB \= CC, 
		 
		 % C must be able to serve as a drive
		 % Maybe not C, but the end of the path of the 
		 % cluster of C?
		 is_influence_chk(C),
		 
		 % A and B of same type
		 % TODO: Too constraining??
		 get_quantity_type(A, Type),
		 get_quantity_type(B, Type),
		 
		 % A and B must be end of their causal paths
		 end_of_path(A, Dependencies),
		 end_of_path(B, Dependencies)
		), MinCandidates).


% Return triples that possibly have a plus relation
% Not used yet
% plus_candidates/2
% plus_candidates(+NaiveModel, -MinCandidates)
plus_candidates(naive_model(Clusters, _, Dependencies), PlusCandidates) :-
	 findall(model_dependency(equal, [C,plus(A,B)]), 
		(
		 model_dependency(prop_pos, [C, A]), 
		 model_dependency(prop_pos, [C, B]),
		 % The feedback cannot be from a cluster to the same cluster
		 get_quantity_cluster(A, Clusters, CA),
		 get_quantity_cluster(B, Clusters, CB),
		 get_quantity_cluster(C, Clusters, CC),
		 CA \= CB, CA \= CC, CB \= CC, 
		% No feedback to proportionalities
		 is_influence_chk(C),
		 
		 % A and B of same type
		 % TODO: Too constraining??
		 get_quantity_type(A, Type),
		 get_quantity_type(B, Type),
		 
		 % A and B must be end of their causal paths
		 end_of_path(A, Dependencies),
		 end_of_path(B, Dependencies)
		), PlusCandidates).

		
	 
%==============================================================================
% BUILDING MODELS BASED ON NAIVE DEPENDENCIES
%==============================================================================

% Order of causality must be equal for all quantities under the same 
% entity (assuming no conditions).
% The ordering is important: e.g. clusters are forced to have same
% influence type if they belong to the same entity type
% We first look for a causal path within all clusters that is compatible 
% when combined into a model. For each cluster a 'drive' is sought
% and feedbacks onto this drive within the clusters are sought if necessary.
% Then initial values are set
% generate_naive_model/1
% generate_naive_model(-Model)
generate_naive_model(FinalModel) :-
	get_all_clusters(Cs),
	assert_clusters(Cs, Clusters),
	!,
	cluster_causal_paths(Clusters, PartialModel),
	valid_partial_model(PartialModel),
	
	add_clusterless_quantities(PartialModel, NPartialModel),
	
	add_cluster_drives(NPartialModel, Model),

	find_cluster_feedbacks(Model, Model2),
	
	cluster_connect_list(Model2, Connect),
	
	connect_clusters(Model2, Connect, NModel),
	
	% Check which quantities are not set
	set_unkown_magnitudes(NModel, NModel2),
	
	find_interactions(NModel2, NModel3),
	
	% Call to find feedbacks that can now aris to the newly found new drives
	% Eventually this should be done iteratively: intermix finding new 
	% feedbacks/cluster links/interactions till no changes arise.
	find_cluster_feedbacks(NModel3, naive_model(Cls, NDrives, NDeps)),
	
	append(NDrives, NDeps, FlatModel),	
	
	add_correspondences(FlatModel, Cls, FinalModel).
	
	
% Figure out which clusters have not yet been connected
% by the previous steps.
% cluster_connect_list/2
% cluster_connect_list(+Model, -ToBeConnectedClusters)
cluster_connect_list(naive_model(Clusters, Drives, _), ConnectClusters) :-
	findall([CTo, CFrom],
		(
		member(model_dependency(_, [To, From]), Drives),
		get_quantity_cluster(To, Clusters, CTo),
		get_quantity_cluster(From, Clusters, CFrom)	
		), List),
	flatten_one_level(List, NotConnect),
	list_to_set(NotConnect, Out),
	subtract(Clusters, Out, ConnectClusters), !.


	
% Add necessary correspondences to the mode
% This includes all correspondences that are alongside
% a proportionality that has already been added
% add_correspondences/3
% add_correspondences(+FlatModel, +Cluster, -FinalModel)
add_correspondences(FlatModel, Clusters, FinalModel) :-
	% Add the correspondences
	% Only adds correspondences that have a proportionality alongside it
	% This prevents that we get cycles of correspondences in clusters
	findall(model_dependency(D, [Q1, Q2]), 
		(
		model_dependency(D, [Q1, Q2]), 
		member(D, [q_correspondence, mirror_q_correspondence]),
		get_quantity_cluster(Q1, Clusters, C),
		get_quantity_cluster(Q2, Clusters, C),
		proportional_connected(Q2, Q1, FlatModel)
		), Corrs),
	append(FlatModel, Corrs, FinalModel), !.
	
	
% Are given quantitities connected using a P?
% proportional_connected/3
% proportional_connected(+QuantityOne, +Quantity2, +Model)
proportional_connected(Q2, Q1, Model) :-
	split_dependency_atom(Prop, prop, _, _),
	member(model_dependency(Prop, [Q1, Q2]), Model), !.
proportional_connected(Q2, Q1, Model) :-
	split_dependency_atom(Prop, prop, _, _),
	member(model_dependency(Prop, [Q2, Q1]), Model), !.

% Connect the remaining clusters with proportionalities 
% Starting from every drive a linking path is formed
% This may be very short, as long as all links combined link all clusters.
% Backtracking let each drive cover other clusters in other orders
% connect_clusters/3
% connect_clusters(+Model, +ToBeConnected, -NewModel)
connect_clusters(naive_model(Clusters, [], Deps), [], naive_model(Clusters, [], Deps)) :- !.
connect_clusters(naive_model(Clusters, [model_dependency(D3, [To, From])|Rest], Deps), ToBeConnected, 
		naive_model(Clusters, [model_dependency(D3, [To, From])|Rest], NDeps2)) :-	
	get_quantity_cluster(To, Clusters, CTo),
	do_connect_clusters(CTo, Clusters, Deps, ToBeConnected, NToBeConnected, ConnectingDeps),
	append(ConnectingDeps, Deps, NDeps),
	connect_clusters(naive_model(Clusters, Rest, NDeps), NToBeConnected, naive_model(Clusters, _, NDeps2)), !.
connect_clusters(Model, _, Model).
	
% Execute the connection finding among the clusters
% Enumerates all possible serial orders of linkage on backtracking
% do_connect_clusters/6
% do_connect_clusters(+CurrentCluster, +AllClusters, +Dependencies, 
%					  +ToBeConnected, -NewToBeConnected, -NewDependencies)
do_connect_clusters(_, _, _, [], [], []).
do_connect_clusters(CurrCluster, AllClusters, Deps, ToBeConnected, NToBeConnectedOut, 
		[model_dependency(Prop, [CBegin, QEnd])|NDeps]) :-
	% Go to next cluster
	member(C, ToBeConnected),
	
	% Get From point
	get_cluster_path(CurrCluster, Deps, QPath),
	end_of_path(QEnd, QPath),
	
	% Get To point
	get_cluster_path(C, Deps, CPath),
	beginning_of_path(CBegin, CPath),
	
	split_dependency_atom(Prop, prop, _, _),
	model_dependency(Prop, [CBegin, QEnd]),
	delete(ToBeConnected, C, NToBeConnected),
	do_connect_clusters(C, AllClusters, Deps, NToBeConnected, NToBeConnectedOut,  NDeps).
do_connect_clusters(_, _, _, C, C, []) :-
	C \= [].
	

% Checks whether this partial model is valid
% This amounts to checking whether causal
% order within equivalent clusters (clusters under the same entity) is the same
% valid_partial_model/1
% valid_partial_model(+Model)
valid_partial_model(naive_model(Clusters, _, Dependencies)) :-
	findall((SC1, SC2),
		(
		member(C1, Clusters),
		member(C2, Clusters),
		C1 \= C2,
		equivalent_clusters(C1, C2),
		sort([C1, C2], [SC1, SC2])
		),	EquivalentClusters),
	list_to_set(EquivalentClusters, SEquivalentClusters),
	have_same_causal_order(SEquivalentClusters, Dependencies).
	
	
% Adds quantities that are not considered as an official cluster 
% to the cluster list as singleton clusters
% add_clusterless_quantities/2
% add_clusterless_quantities(+Model, -NewModel)
add_clusterless_quantities(naive_model(Clusters, Drives, Dependencies), 
		naive_model(NClusters, Drives, Dependencies)) :-
	get_quantity_names(Qs),
	flatten(Clusters, ClusterQs),
	findall([Q], 
		(
		member(Q, Qs),
		\+ member(Q, ClusterQs)
		), NonClusterQs),
	append(NonClusterQs, Clusters, NClusters).
	
	
	
%------------------------------------------------------------------------------
% Predicates for adding cluster drives
%------------------------------------------------------------------------------
	
% Determines which quantities act as drives in this system
% It then adds them at (one of) the appropriate places
% add_cluster_drives/2
% add_cluster_drives(+Model, -NewModel)
add_cluster_drives(Model, Model3) :-
	find_equilibrium_drives(Model, EqQuantities, Model2),
	find_external_drives(Model2, EqQuantities, Model3).
	
	
% Find all drives that form an equilibrium seeking mechanism
% Also return which quanties were set by ESMs
% find_equilibrium_drives/3
% find_equilibrium_drives(+Model, -SetQuantities, -NewModel)
find_equilibrium_drives(Model, SetQuantities, NModel) :-
	min_candidates(Model, Mins),
	get_equilibrium_drives(Model, Mins, Quantities, NModel),
	list_to_set(Quantities, SetQuantities).
	
	
% Given a list of mins, return the drive dependencies of the ESMs
% get_equilibrium_drives/4
% get_equilibrium_drives(+Model, +Mins, -SetQuantities, -NewModel)
get_equilibrium_drives(Model, [], [], Model).
get_equilibrium_drives(
		naive_model(Clusters, DrivesIn, Dependencies),
		[model_dependency(equal, [C, min(A,B)])|Rest], 
		[C|Qs],
		naive_model(Clusters, 
			[model_dependency(inf_neg_by, [InfA, C]), 
			model_dependency(inf_pos_by, [InfB, C])|Drives], 
			[model_dependency(equal, [C, min(A,B)])|Dependencies])) :-
	get_quantity_cluster(A, Clusters, ClusterA),
	get_quantity_cluster(B, Clusters, ClusterB),
	get_cluster_path(ClusterA, Dependencies, PathA),
	get_cluster_path(ClusterB, Dependencies, PathB),
	beginning_of_path(InfA, PathA),
	beginning_of_path(InfB, PathB),
	model_dependency(inf_neg_by, [InfA, C]),
	model_dependency(inf_pos_by, [InfB, C]), !,
	get_equilibrium_drives(naive_model(Clusters, DrivesIn, Dependencies), 
		Rest, 
		Qs, 
		naive_model(Clusters, Drives, Dependencies)).
get_equilibrium_drives(Model, [_|Rest], Qs, NModel) :-
	get_equilibrium_drives(Model, Rest, Qs, NModel).
	
	
% Determine which quantities can act as a external drive
% find_external_drives/3
% find_external_drives(+Model, +SetQuantities, -NewModel)
find_external_drives(naive_model(Clusters, DrivesIn, Dependencies), 
		EqQuantities, naive_model(Clusters, NDrives, Dependencies)) :-
	findall((C1, C2, model_dependency(Dep, [Q1, Q2])),
		(
		split_dependency_atom(Dep, inf, _, _),
		model_dependency(Dep, [Q1, Q2]),
		\+ member(Q1, EqQuantities),
		\+ member(Q2, EqQuantities),
		get_quantity_cluster(Q1, Clusters, C1),
		get_quantity_cluster(Q2, Clusters, C2),
		C1 \= C2
		), ExtDrives),
		
	findall((C1,C2), 
		(
		member((C1, C2, _), ExtDrives)
		), List),
	list_to_set(List, ClusterLinks),

	get_exteral_drives(ClusterLinks, Dependencies, ExtDrivesOut),

	filter_in_drives(Clusters, ExtDrivesOut, NExtDrivesOut),
	filter_out_drives(Clusters, NExtDrivesOut, NExtDrivesOut2),
	
	append(DrivesIn, NExtDrivesOut2, NDrives).
	
	
% Filter drives in such a way that there are no multiple influences coming into
% one cluster
% filter_in_drives/3
% filter_in_drives(+Clusters, +Drives, -NewDrives)
filter_in_drives(_, [], []).
filter_in_drives(Clusters, [model_dependency(Dep, [To, From])|T], [D|DrivesOut]) :-
	% One inf incoming per cluster 
	findall(model_dependency(Dep2, [To, From2]),
		(
		member(model_dependency(Dep2, [To, From2]), T)
		), SameClusterTo),
	removeall(T, SameClusterTo, NewT),
	
	filter_in_drives(Clusters, NewT, DrivesOut),
	member(D, [model_dependency(Dep, [To, From])|SameClusterTo]).
	
	
% This is the final selection, drives selected here are included in the model
% Filters out multiple infl from a quantitity.
% filter_out_drives/3
% filter_out_drives(+Clusters, +Drives, -NewDrives)
filter_out_drives(_, [], []).
filter_out_drives(Clusters, [model_dependency(Dep, [To, From])|T], [D|DrivesOut]) :-
	% One inf incoming per cluster 
	findall(model_dependency(Dep2, [To2, From]),
		(
		member(model_dependency(Dep2, [To2, From]), T)
		), SameClusterFrom),
	removeall(T, SameClusterFrom, NewT),
	
	member(D, [model_dependency(Dep, [To, From])|SameClusterFrom]),
	
	% Check: is this sufficient to drive the model?
	% if so, we are finished
	% Sum up which clusters can be driven (directly or indirectly) using this drive

	drives_clusters(D, Clusters, DrivenClusters),	
	get_quantity_cluster(From, Clusters, CFrom),
	get_quantity_cluster(To, Clusters, CTo),
	DrivenClusters2 = [CFrom, CTo|DrivenClusters],
	
	% Stop if everything is covered, else continue
	(
		(subtract(Clusters, DrivenClusters2, []), DrivesOut = [], !)
	;	
		filter_out_drives(Clusters, NewT, DrivesOut)
	).
	
	
% Checks which clusters can (directly or indirectly) be driven by D
% Checks this by following Ps. We do not have to check if there are Ps
% between all members of the cluster, since this is always the case 
% (due to the Ps and Qs in the cluster).
% drives_clusters/3
% drives_clusters(+Drive, +Clusters, -ClustersDriven)
drives_clusters(model_dependency(_, [To, From]), Clusters, [C|DrivenClusters]) :-
	member(C, Clusters),
	memberchk(Q, C),
	split_dependency_atom(Dep, prop, _, _),
	model_dependency(Dep, [Q, To]),
	delete(Clusters, C, NClusters),!,
	drives_clusters(model_dependency(Dep, [Q, From]), NClusters, DrivenClusters).
drives_clusters(_, _, []).
	

% Find a drive between clusters where possible
% get_exteral_drives/4
% get_exteral_drives(+ClusterPairs, +Dependencies, -Drives)
get_exteral_drives([], _, []).
get_exteral_drives([(C1, C2)|T], Dependencies,  
		[model_dependency(GoodDep, [To, From])|ExtDrivesOut]) :-
	cluster_drive_relation(Dep, C1, C2),
	get_exteral_drives(T, Dependencies, ExtDrivesOut),
	get_cluster_path(C2, Dependencies, C2Path),

	end_of_path(From, C2Path),
	get_cluster_path(C1, Dependencies, C1Path),
	beginning_of_path(To, C1Path),
	adapt_dependency_perspective(To, Dep, GoodDep),!.
get_exteral_drives([_|T], Dependencies, ExtDrivesOut) :-
	get_exteral_drives(T, Dependencies, ExtDrivesOut).
	

% Based on whether the given quantities Side is 1 or -1 within
% the cluster, change the given dependency. This assumes
% that the given dependency is for Side 1 quantities
% The side indicates if a quantity in the cluster
% inversely corresponds to the other
% adapt_dependency_perspective/3
% adapt_dependency_perspective(+Quantity, +Dependency, -NewDependency)
adapt_dependency_perspective(Quantity, Dependency, Dependency) :-
	cluster_member(Quantity, 1, _).
adapt_dependency_perspective(Quantity, Dependency, DependencyOut) :-
	cluster_member(Quantity, -1, _),
	split_dependency_atom(Dependency, Type, Sign, Dir),
	invert_sign(Sign, InvSign),
	split_dependency_atom(DependencyOut, Type,InvSign, Dir).
	
	
% find possible links (drives) between clusters
% The returned result is based on a drive from any
% quantity in C1 to a quantity of Side = 1.
% To the -1 we obviously have an opposite influence 
% (adaptable using adapt_dependency_perspective/3)
% cluster_drive_relation/3
% cluster_drive_relation(?Drive, +Cluster1, +Cluster2)
cluster_drive_relation(inf_pos_by, C2, C1) :-
	findall(1,
		(
			member(Q1, C1),
			member(Q2, C2),
			((
			cluster_member(Q2, 1, C2), 
			model_dependency(inf_pos_by, [Q2, Q1])
			)
			;
			(
			cluster_member(Q2, -1, C2), 
			model_dependency(inf_neg_by, [Q2, Q1])
			))
		), Links),
	length(Links, Length),
	length(C1, C1Len),
	length(C2, C2Len),
	Length is C1Len * C2Len.
cluster_drive_relation(inf_neg_by, C2, C1) :-
	findall(1,
		(
			member(Q1, C1),
			member(Q2, C2),
			((
			cluster_member(Q2, 1, C2), 
			model_dependency(inf_neg_by, [Q2, Q1])
			)
			;
			(
			cluster_member(Q2, -1, C2), 
			model_dependency(inf_pos_by, [Q2, Q1])
			))
		), Links),
	length(Links, Length),
	length(C1, C1Len),
	length(C2, C2Len),
	Length is C1Len * C2Len.	
	

	
%------------------------------------------------------------------------------
% Finding feedbacks in cluster causal paths
%------------------------------------------------------------------------------
	
% For all clusters add feedback loops to their paths
% If there is one, we always take the one from the end of the path
% if there is a feedback, the last one will also be there.
% find_cluster_feedbacks/2
% find_cluster_feedbacks(+Model, -NewModel
find_cluster_feedbacks(naive_model(Cs, [], Deps), naive_model(Cs, [], Deps)).
find_cluster_feedbacks(
		naive_model(Clusters, [model_dependency(D, [To, Drive])|T], Deps), 
		naive_model(Clusters, [model_dependency(D, [To, Drive])|T], 
			[model_dependency(D2, [Drive, End])|DepsOut])) :-
	get_quantity_cluster(To, Clusters, C),
	get_cluster_path(C, Deps, CPath),
	end_of_path(End, CPath),
	split_dependency_atom(D2, prop, _, _),
	model_dependency(D2, [Drive, End]),	
	\+ member(model_dependency(D2, [Drive, End]),Deps),	!,
	find_cluster_feedbacks(naive_model(Clusters, T, Deps), naive_model(Clusters, T, DepsOut)).
find_cluster_feedbacks(naive_model(Clusters, [H|T], Deps), naive_model(Clusters, [H|T], NDeps)) :-
	find_cluster_feedbacks(naive_model(Clusters, T, Deps), naive_model(Clusters, T, NDeps)).
	

	
%------------------------------------------------------------------------------
% Identifying and assigning initial values
%------------------------------------------------------------------------------
		
% Set the magnitudes that are not set by the scenario
% or by the generated correspondences
% set_unkown_magnitudes/2
% set_unkown_magnitudes(+Model, -NewModel)
set_unkown_magnitudes(Model, NewModel) :-
	get_quantity_names(Qs), 
	quantity_initial_status(Qs, Model, Knowns, Unknowns),
	resolve_unknowns(Unknowns, Knowns, Model, NewModel).
	
	
% Identify which quantities are known (or will become known by correspondences) 
% and which quantities are unkown but should get a value 
% for simulation to succeed
% quantity_initial_status/3
% quantity_initial_status(+Quantities, -KnownQuantities, -UnknownQuantities)
quantity_initial_status([], _, [], []) :- !.
quantity_initial_status([Q|Rest], M, [Q/Reason|Knowns], Unknowns) :- % Given in scenario
	magnitude_known(Q, Reason),
	quantity_initial_status(Rest, M, Knowns, Unknowns), !.
quantity_initial_status([Q|Rest], naive_model(Clusters, Drives, Deps),
		[Q/correspondence|Knowns], Unknowns) :- % Known through correspondence within cluster
	get_quantity_cluster(Q, Clusters, C),
	get_cluster_path(C, Deps, CPath),
	\+ beginning_of_path(Q, CPath),
	quantity_initial_status(Rest, naive_model(Clusters, Drives, Deps), Knowns, Unknowns), !.
quantity_initial_status([Q|Rest], naive_model(Clusters, Drives, Deps),
		[Q/correspondence|Knowns], Unknowns) :- % Known correspondence with other clusters
	get_quantity_cluster(Q, Clusters, C),
	get_cluster_path(C, Deps, CPath),
	beginning_of_path(Q, CPath),
	dependency_bag(Q, Q2, QDeps),
	member(q_correspondence, QDeps),
	get_quantity_cluster(Q2, Clusters, C2),
	C2 \= C,
	quantity_initial_status(Rest, naive_model(Clusters, Drives, Deps), Knowns, Unknowns), !.
quantity_initial_status([Q|Rest], naive_model(Clusters, Drives, Deps),
		Knowns, [Q/calculus|Unknowns]) :- % Calculus but unkown until proven otherwise	
	member(model_dependency(equal, [Q, min(_, _)]), Deps), 
	quantity_initial_status(Rest, naive_model(Clusters, Drives, Deps), Knowns, Unknowns), !.
quantity_initial_status([Q|Rest], naive_model(Clusters, Drives, Deps),
		Knowns, [Q/calculus|Unknowns]) :- % Calculus but unkown until proven otherwise	
	member(model_dependency(equal, [Q, plus(_, _)]), Deps), 
	quantity_initial_status(Rest, naive_model(Clusters, Drives, Deps), Knowns, Unknowns), !.
quantity_initial_status([Q|Rest], naive_model(Clusters, Drives, Deps),
		Knowns, [Q/causal_start|Unknowns]) :- 
	get_quantity_cluster(Q, Clusters, C),
	get_cluster_path(C, Deps, CPath),
	beginning_of_path(Q, CPath), 
	quantity_initial_status(Rest, naive_model(Clusters, Drives, Deps), Knowns, Unknowns), !.
quantity_initial_status([Q|Rest], naive_model(Clusters, Drives, Deps),
		Knowns, [Q/unkown|Unknowns]) :- % Just unkown
	quantity_initial_status(Rest, naive_model(Clusters, Drives, Deps),Knowns, Unknowns), !.
	
	
% Resolve unkown intial magnitudes by doing one of the following:
% - If a calculus element was added, check whether this calculus
%   can initially be evaluated. If not makes sure it can be by adding
%   (in)equalities
% - In case of a causal start, add a value assignment without a condition
% - If that is not possible, find one with a condition
% resolve_unknowns/4
% resolve_unknowns(+UnkownQuantities, +Model, -NewModel)
resolve_unknowns([], _, naive_model(Cs, Drives, Deps), naive_model(Cs, Drives, Deps)).
resolve_unknowns([Q/calculus|Rest], Knowns, naive_model(Cs, Drives, Deps), naive_model(Cs, Drives, NDeps)) :-
	member(model_dependency(equal, [Q, Calc]), Deps),
	Calc =.. [_, _, _],
	calculus_evaluatable(model_dependency(equal, [Q, Calc]), naive_model(Cs, Drives, Deps)),!,
	% Nothing changes, this one is actually known
	resolve_unknowns(Rest, Knowns, naive_model(Cs, Drives, Deps), naive_model(Cs, Drives, NDeps)).
resolve_unknowns([Q/calculus|Rest], Knowns, naive_model(Cs, Drives, Deps), naive_model(Cs, Drives, NDeps2)) :-
	member(model_dependency(equal, [Q, Calc]), Deps),
	Calc =.. [Op, Q2, Q3],
	\+ calculus_evaluatable(model_dependency(equal, [Q, Calc]), naive_model(Cs, Drives, Deps)),
	% Calculus not evaluatable, try to make the calculus evaluatable
	% TODO: Implement the following function
	facilitate_calc_evalutation(Op, Q2, Q3, SolutionDeps),!,
	resolve_unknowns(Rest, Knowns, naive_model(Cs, Drives, Deps), naive_model(Cs, Drives, NDeps)),
	append(SolutionDeps, NDeps, NDeps2).
resolve_unknowns([Q/causal_start|Rest], Knowns, naive_model(Cs, Drives, Deps), 
		naive_model(Cs, Drives, [VA|NDeps])) :-	
	unconditional_value_assignment(Q, VA), !, % TODO: remove cut when DEBUG SKIP is removed
	resolve_unknowns(Rest, Knowns, naive_model(Cs, Drives, Deps), naive_model(Cs, Drives, NDeps)).
resolve_unknowns([Q/causal_start|Rest], Knowns, naive_model(Cs, Drives, Deps), 
		naive_model(Cs, Drives, NDeps2)) :-	
	conditional_value_assignment(Q, Knowns, MFs), !, % TODO: remove cut when DEBUG SKIP is removed
	resolve_unknowns(Rest, Knowns, naive_model(Cs, Drives, Deps), naive_model(Cs, Drives, NDeps)),
	append(NDeps, MFs, NDeps2).
resolve_unknowns([_|Rest], Knowns, Model, Out) :- % DEBUG SKIP (NB SKIP CAUSES FAULTY OUTPUT)
	resolve_unknowns(Rest, Knowns, Model, Out).	  % Used when an unresolvable init val is found


% Determine whether a calculus element is evaluatable with the 
% information provided by the scenario.
% calculus_evaluatable/2
% calculus_evaluatable(+Calculus, +Model)
calculus_evaluatable(model_dependency(equal, [_, Calc]), _) :-
	Calc =.. [_, Q1, Q2],
	inequality_derivable(Q1, Q2).
	
% Checks whether an inequality between Q1 and Q2 can be derived
% directly or inderectly
% TODO: Implement this functions. It should check the initial values
% of the state graph to check whether there are patterns that indicated
% inequalities as conditions in MFs.
inequality_derivable(_Q1, _Q2).

% TODO: Implement the following function
facilitate_calc_evalutation(_Op, _Q2, _Q3, _SolutionDeps) :-
	fail.


% Check if QuantityName is known (either in scenario or via correspondences)
% magnitude_known/2
% magnitude_known(+Quantity, ?KnownType)
magnitude_known(QuantityName, scenario) :-
	scenario_given(QuantityName).
magnitude_known(QuantityName, correspondence) :-
	get_quantity_names(Qs),
	member(QuantityName, Qs),
	member(Q2, Qs),
	Q2 \= QuantityName,
	scenario_given(Q2),
	correspondence_linked(QuantityName, Q2).

% Checks if two quantities are directly linked with a correspondence
% direct_correspondence_linked/2
% direct_correspondence_linked(+Quantity1, +Quantity2)
direct_correspondence_linked(Q1, Q2) :-
	dependency_bag(Q1, Q2, Deps),
	intersection([q_correspondence, mirror_q_correspondence], Deps, Result),
	Result \= [].
	
% Wrapper for correspondence_linked/3
% correspondence_linked/2
% correspondence_linked(+Quantity1, +Quantity2)
correspondence_linked(Q1, Q2) :-
	correspondence_linked(Q1, Q2, [Q1, Q2]).
	
% correspondence_linked/3
% correspondence_linked(+Quantity1, +Quantity2, +AccVisited)
correspondence_linked(Q1, Q2, _) :-
	direct_correspondence_linked(Q1, Q2).
correspondence_linked(Q1, Q2, Visited) :-
	direct_correspondence_linked(Q1, Q3),
	\+ member(Q3, Visited),
	correspondence_linked(Q3, Q2, [Q3|Visited]).
	
% Check if QuantityName is given in the scenario
% scenario_given/1
% scenario_given(+Quantity)
scenario_given(QuantityName) :-
	% get all assigned values from scenario
	get_input_system_name(Name),
	get(@app, findExportedObject, mf, Name,  ScenarioObj),
	get(ScenarioObj?elements, find_all, message(@prolog, ==, @arg1?class_name, 'value'), Values),
	chain_list(Values, ValuesList),
	quantity_value_member(QuantityName, ValuesList).
	
	
% Find a value assignment for Quantity without conditions, 
% in order to set its initial magnitude. Fails if none is found.
% unconditional_value_assignment/2
% unconditional_value_assignment(+Quantity, -ValueAssignment)
unconditional_value_assignment(Quantity, value(Quantity,_,Value,_)) :-
	get_quantity_values(1, magnitude, Quantity, Value),
	findall(S, get_quantity_values(S, magnitude, Quantity, Value), States),
	findall(S, engine:state(S, _), States2),
	sort(States, Sorted),
	sort(States2, Sorted).

	
% Is there a value in ValueList that belongs to the 
% quantity indicated by QuantityDefinition
% quantity_value_member/2
% quantity_value_member(+QuantityName, +ValueList
quantity_value_member(Name, [H|_]) :-
	get_quantity_owner(Name, Owner),
	pretty_atom_to_string(Owner, COwner),
	
	get_quantity_type(Name, Type),
	pretty_atom_to_string(Type, CType),
	
	get(H, quantity, Q), 
	get(Q, name, CType),
	get(Q, garpInstance, D), 
	get(D, name, COwner).
quantity_value_member(Def, [_|T]) :-
	quantity_value_member(Def, T).
	
	
	
%------------------------------------------------------------------------------
% Adding Causal Paths
%------------------------------------------------------------------------------

% For each cluster generate the causal links between the quantities
% within this cluster
% cluster_causal_paths/2
% cluster_causal_paths(+Clusters, -Model)
cluster_causal_paths([], naive_model([], [], [])).
cluster_causal_paths([cluster(_, C)|T], naive_model([C|Clusters], Drives, NPaths)) :-
	cluster_causal_paths(T, naive_model(Clusters, Drives, Paths)),
	causal_path(C, P),
	append(P, Paths, NPaths).


% Checks whether two clusters are equivalent
% Two clusters are equivalent if they contain the same types
% of quantities belonging to the same type of entity.
% equivalent_clusters/2
% equivalent_clusters(+Cluster1, +Cluster2)
equivalent_clusters(C1, C2) :-
	rewrite_cluster(C1, NC1),
	rewrite_cluster(C2, NC2),
	sort(NC1, Sorted),
	sort(NC2, Sorted).

	
% Rewrite cluster for easy cluster comparison
% Every quantity is rewritten as its type and its 
% owners type
% rewrite_cluster/2
% rewrite_cluster(+Cluster, -RewrittenCluster)
rewrite_cluster([], []).
rewrite_cluster([H|T], [OwnerType/Type|Cluster]) :-
       rewrite_cluster(T, Cluster),
       get_quantity_owner(H, Owner),
       get_entity_type(Owner, OwnerType), !,
       get_quantity_type(H, Type).

	   
% Verifies that a list of EquivalentClusters have the 
% same causal orderings
% have_same_causal_order/
% have_same_causal_order(+EquivalentClusters, Depenendencies)
have_same_causal_order([], _).
have_same_causal_order([(C1, C2)|Rest], Dependencies) :-
	get_cluster_path(C1, Dependencies, Path1),
	get_cluster_path(C2, Dependencies, Path2),
	same_causal_order(Path1, Path2),
	have_same_causal_order(Rest, Dependencies).
	
	
% Returns the causal path of a cluster
% get_cluster_path/3
% get_cluster_path(+Cluster, +Dependencies, -CausalPath)
get_cluster_path([Q], _, [null_path/Q]) :- !.
get_cluster_path(Cluster, Dependencies, Path) :-
	findall(model_dependency(D, [Q1, Q2]),
		(
		member(model_dependency(D, [Q1, Q2]), Dependencies),
		member(Q1, Cluster),
		member(Q2, Cluster)
		), Path).

		
% Return the cluster of a given quantity
% get_quantity_cluster/3
% get_quantity_cluster(+Quantity, +Clusters, -Cluster)
get_quantity_cluster(Q, Clusters, C) :-
	member(C, Clusters),
	member(Q, C), !.
	   
	   
% Verfies that Path1 and Path2 have the same causal order
% Prepares paths to right format for comparison
% same_causal_order/2
% same_causal_order(+Path1, +Path2)
same_causal_order(Path1, Path2) :-
	beginning_of_path(B1, Path1),
	beginning_of_path(B2, Path2),
	path_to_list(Path1, B1, [B1], L1),
	path_to_list(Path2, B2, [B2], L2),
	check_same_causal_order(L1, L2), !.

% Verfies that Path1 and Path2 have the same causal order
% check_same_causal_order/2
% check_same_causal_order(+PathList, +PathList)
check_same_causal_order([], []).
check_same_causal_order([H1|T1], [H2|T2]) :-
	get_quantity_type(H1, Type),
	get_quantity_type(H2, Type),
	check_same_causal_order(T1, T2).

	
% Return a dependency form path to a list of 
% quantities in the path in the right order
% path_to_list/4
% path_to_list(+Path, +CurrentNode, +AccPathList, +PathList)
path_to_list(Path, CurrNode, AccList, List) :-
	next_in_path(CurrNode, Path, Next),
	path_to_list(Path, Next, [Next|AccList], List), !.
path_to_list(_Path, _CurrNode, List, List) :- !.

	
% Given a set of Quantities, find a causal path between these quantities
% causal_path/2
% causal_path(+Quantities, -CausalPath)
causal_path(Quantities, Path) :-
	member(Q, Quantities),
	split_dependency_atom(D2, prop, _, _),
	model_dependency(D2, [To, Q]),
	member(To, Quantities),
	build_causal_path(Quantities, [model_dependency(D2, [To, Q])], [To, Q], Path).
	
	
% Find the full path starting at an influence
% build_causal_path/4
% build_causal_path(+Quantities, +AccPath, +AccVisited, -Path)	
build_causal_path(Quantities, [model_dependency(D1, [B, A])|Rest], Visited, Path)	:-
	split_dependency_atom(D2, prop, _, _),
	model_dependency(D2, [C, B]),
	\+ memberchk(C, Visited),
	member(C, Quantities),
	build_causal_path(Quantities, [model_dependency(D2, [C, B]), model_dependency(D1, [B, A])|Rest], [C|Visited], Path).
build_causal_path(Quantities, [model_dependency(D1, [B, A])|Rest], Visited, [model_dependency(D1, [B, A])|Rest]) :-
	length(Visited, L),
	length(Quantities, L).
	
	
% Returns beginning of a path
% beginning_of_path/2
% beginning_of_path(+Path, -BeginningOfPath)
beginning_of_path(Q, [null_path/Q]).
beginning_of_path(A, Path) :-
	split_dependency_atom(D, inf, _, _),
	member(model_dependency(D, [B, A]), Path), 
	in_path_chk(B, Path), !.
beginning_of_path(A, Path) :-
	in_path(A, Path),
	\+ member(model_dependency(_D, [A, _]), Path), !.
	
% Returns end of a path
% end_of_path/2
% end_of_path(+Path, -EndOfPath)
end_of_path(Q, [null_path/Q]).
end_of_path(A, Path) :-
	in_path(A, Path),
	\+ member(model_dependency(_D, [_, A]), Path), !.	
	
		
% Return next quantity from the causal path,
% given the current quantity and the path
% next_in_path/3
% next_in_path(+CurrentQuantity, +Path, -NextQuantity)
next_in_path(Current, Path, Next) :-
	split_dependency_atom(D, inf, _, _),
	member(model_dependency(D, [Next, Current]), Path).	
next_in_path(Current, Path, Next) :-
	split_dependency_atom(D, prop, _, _),
	member(model_dependency(D, [Next, Current]), Path).
	
	
	
% In path check with choice point
% in_path/2
% in_path(?Quantity, +Path)
in_path(Q, [model_dependency(_, Qs)|_T]) :-
	member(Q, Qs).
in_path(Q, [_|T]) :-
	in_path(Q, T).	
	
	
% In path without choice point
% in_path_chk/2
% in_path_chk(+Quantity, +Path)
in_path_chk(Q, [model_dependency(_, Qs)|_T]) :-
	member(Q, Qs), !.
in_path_chk(Q, [_|T]) :-
	in_path(Q, T).	
	

%==============================================================================
% DETERMINING CLUSTERS
%==============================================================================

% Identify all clusters within the set of naive dependencies
% get_all_clusters/1
% get_all_clusters(-Clusters)
get_all_clusters(NClusters) :-
	preprocess_inequalities,
	get_quantity_names(QNames),
	do_get_clusters(QNames, [], [], Clusters), 
	filter_false_clusters(Clusters, NClusters).
	
	
% Execute the cluster finding process.
% do_get_clusters/4
% do_get_clusters(+Quantities, +AccClusters, +AccVisited, -Clusters)
do_get_clusters([], Clusters, _, Clusters).
do_get_clusters([H|T], AccClusters, Visited, Clusters) :-
	\+ member((H,_), Visited),
	get_cluster([(H, 1)], [(H, 1)], cluster(Owner, Members)),
	Members \= [],!,
	append(Visited, Members, NVisited),
	do_get_clusters(T, [cluster(Owner, Members)|AccClusters], NVisited, Clusters).
do_get_clusters([_|T], AccClusters, Visited, Clusters) :-
	do_get_clusters(T, AccClusters, Visited, Clusters).


% For a given quantity expand its cluster by checking which quantities are connected
% to it with correspondences 
% get_cluster/3
% get_cluster(+AccCluster1, +AccCluster2, -Cluster)
get_cluster([], [(Q, Side)|T], cluster(Owner, [(Q, Side)|T])) :-
	T \= [],
	valid_cluster([(Q, Side)|T]),
	get_quantity_owner(Q, Owner).
get_cluster([], _, cluster(_, [])).
get_cluster([(A, Side)|T], AccCluster, Cluster) :-
	get_quantity_names(QNames),
	member(A, QNames),
	get_quantity_owner(A, Entity),

	findall((B, Side2), 
		(
		split_dependency_atom(Dep, prop, _, _),
		dependency_bag(B, A, Ds),
		member(Dep, Ds),
		corresponds(Ds, Side, Side2),
		\+ member((B, Side2), AccCluster),
		valid_cluster([(B, Side2)|AccCluster]),
		get_quantity_owner(B, Entity)
		), NewCenters),
	append(NewCenters, T, NewT),
	append(NewCenters, AccCluster, AccCluster2),
	get_cluster(NewT, AccCluster2, Cluster).
	
	
% Returns wether two quantities with given inter-dependencies correspond
% In addition returns what the side of Q2 is given Q1's Side
% This is based how they are connected with a correspondence
% corresponds/3
% corresponds(+InterDependencies, +Q1Side, -Q2Side)
corresponds(Dependecies, Side, Side2) :-
	member(q_correspondence, Dependecies), 
	member(prop_pos, Dependecies), 
	Side2 = Side.
corresponds(Dependecies, Side, Side2) :-
	member(mirror_q_correspondence, Dependecies), 
	member(prop_neg, Dependecies), 
	Side2 is -Side.

		
% Remove overlapping clusters.
% filter_false_clusters/
% filter_false_clusters(+Clusters, -FilteredClusters)
filter_false_clusters(C, NClusters) :-
	findall(Cls, member(cluster(_, Cls), C), Clusters),
	member(C1, Clusters),
	member(C2, Clusters),
	C1 \= C2,
	intersection(C1, C2, I),
	I \= [],
	delete(Clusters, C1, Cls2),
	delete(Cls2, C2, Cls3), !,
	filter_false_clusters(Cls3, NClusters).
filter_false_clusters(Clusters, Clusters).
	

% Wrapper for valid_cluster/2
% valid_cluster/1
% valid_cluster(
valid_cluster(Cluster) :-
	valid_cluster(Cluster, Cluster).
	
% Check whether Cluster is a valid cluster
% This is based on checking whether all members
% are connected to all the other with correspondences
% valid_cluster/2
% valid_cluster(+Clusters, +Clusters)
valid_cluster([], _).
valid_cluster([(H, _)|T], Cluster) :-
	complies_with(H, Cluster),
	valid_cluster(T, Cluster).

	
% Does Quantity comply with cluster?
% That is, does it contain links to all other cluster members
% complies_with/2
% complies_with(+Quantity, +OtherQuantities)
complies_with(_, []).
complies_with(Q, [(Q, _)|T]) :-
	!,
	complies_with(Q, T).
complies_with(Q, [(H, _)|T]) :-
	dependency_bag(Q, H, Deps),
	split_dependency_atom(D, prop, _, _),
	member(D, Deps),
	(member(q_correspondence, Deps) ; member(mirror_q_correspondence, Deps)),
	complies_with(Q, T).
	
	
% Preprocess (in)equalities to make it easier to check
% whether an inequality holds between two quantities.
% Asserts if in inequality holds between two quantities
% somewhere in the stategraph
% preprocess_inequalities/0
% preprocess_inequalities
preprocess_inequalities :-
	retractall(unequal_somewhere(_)),
	findall(unequal_somewhere(Sorted),
		(
		get_property(_, 6, par_relations(Ineq)),
		member(Inq,Ineq),
		Inq =.. [H, Q1, Q2], 
		member(H, [greater, smaller]),
		% There should be no min(_,_) or plus(_,_) in the (in)equality
		Q1 =.. Q1List,
		Q2 =.. Q2List,
		length(Q1List, LQ1),
		length(Q2List, LQ2),
		LQ1 < 3,
		LQ2 < 3,
		sort([Q1, Q2], Sorted),
		\+ unequal_somewhere(Sorted)
		),Uneq),
	list_to_set(Uneq, SetUneq),
	remove_naive_correspondences(SetUneq).

% Remove Q-correspondences between quantities that were found
% to be unequal somewhere. These quantities can not naively be 
% considered as corresponding
% remove_naive_correspondences/1
% remove_naive_correspondences(+SomewhereUnequalQuantities)
remove_naive_correspondences([]).
remove_naive_correspondences([unequal_somewhere([Q1, Q2])|T]) :-
	dependency_bag(Q1, Q2, Deps),!,
	retractall(model_dependency(q_correspondence, [Q1, Q2])),
	retractall(model_dependency(q_correspondence, [Q2, Q1])),
	retractall(dependency_bag(Q1, Q2, Deps)),
	retractall(dependency_bag(Q2, Q1, Deps)),
	delete(Deps, q_correspondence, NDeps),
	assert(dependency_bag(Q1, Q2, NDeps)),
	assert(dependency_bag(Q2, Q1, NDeps)),
	remove_naive_correspondences(T).
remove_naive_correspondences([_|T]) :-	
	remove_naive_correspondences(T).
	
	
	
%==============================================================================
% AUXIALLIARY PREDICATES
%==============================================================================
	
% Assert the model for better accesibility
% assert_model/1
% assert_model(+Model)
% TODO: add ternary preds
assert_model(model([])).
assert_model(model([H|T])) :-
	H =.. [Dependency, Q1, Q2],
	assert(model_dependency(Dependency, [Q1, Q2])),
	(dependency_bag(Q1, Q2, D1);D1 = []),
	(dependency_bag(Q2, Q1, D2);D2 = []),
	retractall(dependency_bag(Q1, Q2, _)),
	retractall(dependency_bag(Q2, Q1, _)),
	sort([Dependency|D1], SD1),
	sort([Dependency|D2], SD2),
	assert(dependency_bag(Q1, Q2, SD1)),
	assert(dependency_bag(Q2, Q1, SD2)),
	assert_model(model(T)).
	
	
% assert_clusters/2
% assert_clusters(+Clusters, -ReformattedClusters)
assert_clusters([], []) :-
		retractall(cluster_member(_,_,_)).
assert_clusters([cluster(Owner, Members)|T], [cluster(Owner,Qs)|Clusters]) :-
	assert_clusters(T, Clusters),
	findall(Q, 
		(
		member((Q, Side), Members)
		), Qs),
	% Loop hack
	findall(_, 
		(
		member((Q, Side), Members),
		assert(cluster_member(Q, Side, Qs))
		), _).
	
	
% Return all dependencies that are consistent with the current states
% values, derivatives and inequalities
% get_consistent_dependencies/4
% get_consistent_dependencies(+State, +Quantities, -Dependencies)
get_consistent_dependencies(StateId, QuantList, Dependencies) :-
	reverse(QuantList, RevQuants),
	setof(D, 
		(
		dependency_consistency_constraint(StateId, QuantList, D);
		dependency_consistency_constraint(StateId, RevQuants, D)
		), Dependencies), !.
get_consistent_dependencies(_, _, []).
	

% Generate all pairs of quantities
% generate_quantity_pairs/1
% generate_quantity_pairs(+QuantityPairs)
generate_quantity_pairs(Pairs2) :-
	get_state_quantities(1, Quant),
	findall(S, 
		(
		member(quantity(_, Q1, _, _), Quant),
		member(quantity(_, Q2, _, _), Quant),
		Q1 \= Q2,
		sort([Q1, Q2], S)
		), Pairs),
	list_to_set(Pairs, Pairs2).
	
	
% Check whether Quantity serves as an influence
% according to the naive dependencies (no backtracking)
% is_influence_chk/1
% is_influence_chk(+Quantity)
is_influence_chk(Q) :-
	model_dependency(inf_pos_by, [_To, Q]), !.
is_influence_chk(Q) :-
	model_dependency(inf_neg_by, [_To, Q]), !.
	
	
% Flatten list only one level
% flatten_one_level/2
% flatten_one_level(+List, -FlattenedList)
flatten_one_level([], []) :- !.
flatten_one_level([L|T], Out) :-
	is_list(L),
	flatten_one_level(T, TOut),
	append(L, TOut, Out), !.
flatten_one_level([A|T], [A|TOut]) :-
	atom(A),
	flatten_one_level(T, TOut), !.
	
	
% Remove all items from RemoveList from List
% removeall/3
% removeall(+List, +RemoveList, -NewList)
removeall(List, [], List) :- !.
removeall(List, [H|T], NList2) :-
	removeall(List, T, NList),
	delete(NList, H, NList2), !.
	
	
% Invert the sign of a dependency
% invert_sign/1
% invert_sign(?Sign, ?InvSign)
invert_sign(pos, neg).
invert_sign(neg, pos).