%% bouncingBall.pl %%
%% Using Constraint Logic Programming to solve a planning problem for bouncing balls.
%% This example adapts a DFS planner without heuristics and reduces the planning
%% problem to solving a Linear Programming problem with an external EPLEX library.
%% This library is designed to interface with Linear Programming solvers only. 
%% For this reason, the trajectory of the bouncing ball is approximated as a
%% linear function of time, since otherwise a different interface to an external
%%  non-linear programming solver would be required. 
%% 
%% Author: Mikhail Soutchanski. 
%% Contact Email:   soutchanski (at) protonmail (dot) com
%% Copyright (c) 2020, 2023 Mikhail Soutchanski
%%
%% This software is available with the MIT Licence  https://opensource.org/license/MIT

%% See details about physics at
%% https://www.physicsclassroom.com/class/vectors/Lesson-2/Characteristics-of-a-Projectile-s-Trajectory 
%% This example works with ECLiPSe Prolog 6.1 and more recent releases.

/* Eplex: Harnessing Mathematical Programming Solvers for Constraint Logic Programming
    by Kish Shen and Joachim Schimpf. In: Proceedings of the 11th International 
    Conference on Principles and Practice of Constraint Programming, pp 622–636.
    CP-2005, LNCS 3709,  Springer Verlag, 2005. 
    https://link.springer.com/chapter/10.1007/11564751_46
    http://eclipseclp.org/reports/eplex_cp2005.pdf
*/
    
:- lib(eplex).  /* EPLEX library with an external linear constraints solver */

/* The following simple iterative deepening DFS planner is taken from the
   Hector Levesque's course notes prepared for the 1st year SCI-199 course.
   A similar program is in Chapter 9 of his textbook "Thinking as Computation", 
   MIT Press, 2017. It turns out that his program can be adapted to solve 
   planning problems for the mixed discrete-continuous (hybrid) systems.
*/

solve_problem(N,Plan) :- C0 is cputime,
                    max_length(L,N),
                    reachable(S,L), goal_state(S),
					myReverse(L,Plan),
                    Cf is cputime, Duration is Cf - C0, nl,
                    write('Elapsed time (sec): '), write(Duration), nl.

max_length([],N) :- N >= 0.
max_length([_|L],N1) :- N1 > 0, N is N1 -1, max_length(L,N).

myReverse(L1,L2) :- myRev(L1, [], L2).

myRev([], L1, L1).
myRev([H|L1], L2, L3) :- myRev(L1, [H|L2], L3).

reachable(S, []) :- initial_state(S).
reachable(S2, [M|ListOfActions]) :- 
		reachable(S1, ListOfActions), legal_move(S2, M, S1).

	%using legal_move predicate with natural actions
legal_move([A|S], A, S) :- natural(Nat), 
	( poss(Nat,S) -> A=Nat ; (agent(A), poss(A,S) )
	).

/* 
% The above rule can be simplified as follows.
legal_move([A|S], A, S) :- natural(A), poss(A,S).
legal_move([A|S], A, S) :- agent(A), poss(A,S).
*/
ball(b1).	ball(b2).
epsilon(E) :- E is (1/10). %approximation

natural(bounce(Ball,Time)).
natural(atPeak(Ball,Time)).
agent(drop(Ball,Time)).
agent(catch(Ball, Time)).

		%--------- Precondition axioms ---------------%

poss(drop(Ball,Time),S) :- ball(Ball), not falling(Ball,S), not flying(Ball,S),
		start(S, T), T $=< Time.

poss(catch(Ball, Time), S) :- ball(Ball), 
		( falling(Ball, S); flying(Ball,S) ), start(S,T), T $=< Time.

poss(bounce(Ball, Time), S) :- ball(Ball), falling(Ball, S), epsilon(E),
	dist(Ball, Distance, Time, S), Distance $= 0, 
	vel(Ball, Velocity, Time, S), (Velocity $>= E ),
		start(S,T), T $=< Time.

poss(atPeak(Ball, Time), S) :- ball(Ball), flying(Ball, S), epsilon(E),
	dist(Ball, Distance, Time, S), Distance $>= E, 
	vel(Ball, Velocity, Time, S), Velocity $=0,
			start(S,T), T $=< Time.

/* falling(B,S): ball B is falling down and accelerating under the Earth gravity.
flying(B,S): ball B bounced and it is flying up and decelerating due to gravity.
Assume the vertical axis is oriented downwards, i.e., if a ball is falling
down, then its speed is positive and increases. But when it bounces, its
speed is negative and decreases. For simplicity: elastic ball, no damping.
*/
			%-----Successor state axioms ------------%

falling(Ball, [A|S]) :- A = drop(Ball, Time).
falling(Ball, [A|S]) :- A = atPeak(Ball, Time).
falling(Ball, [A|S]) :- not A = catch(Ball,Time), not A=bounce(Ball, Time),
		falling(Ball, S).

flying(Ball, [A|S]) :- A = bounce(Ball, Time).
flying(Ball, [A|S]) :- not A = catch(Ball,Time), not A=atPeak(Ball, Time),
		flying(Ball, S).

numberFlights(0,B, []).
numberFlights(N,B, [A | S]) :- A = atPeak(B,_), numberFlights(M,B,S), N is M+1.
numberFlights(N,B, [A | S]) :- not A = atPeak(B,_), numberFlights(N,B,S).

/* The height of the ball from the floor  dist(B,Dist,T,S)   changes 
	continuously, no matter what actions or events happen.			 */
init_dist(Ball,Dist,[Action|S]) :- time(Action, Time), 
    dist(Ball,Dist,Time,S).

init_vel(b1,0,[]).
init_vel(b2,0,[]).

init_vel(Ball,Vel, [Action|S]) :- Action = catch(Ball, Time), Vel $= 0.

init_vel(Ball,Vel, [Action|S]) :- Action = bounce(Ball, Time), 
	time(Action, Time),
    vel(Ball, Old_Vel, Time, S), Vel $= -1*Old_Vel.

init_vel(Ball,Vel,[Action|S]) :- time(Action,T),
	not Action = catch(Ball, _), not Action = bounce(Ball, _), 
	vel(Ball, Vel, T, S).

		/* ------------ State evolution axioms --------------
		See details at  https://arxiv.org/abs/1807.04861 
		Hybrid Temporal Situation Calculus
        Vitaliy Batusov, Giuseppe De Giacomo, Mikhail Soutchanski   */

vel(b1, 0, 0, []).
vel(b2, 0, 0, []).

vel(Ball, Vel, Time, S) :- not falling(Ball, S), not flying(Ball, S), Vel $= 0.

vel(Ball, Vel, Time, S) :- falling(Ball, S), 
	init_vel(Ball, Old_Vel, S), start(S, T1), Time $>= T1,
    Vel $= Old_Vel + 9.81*(Time - T1). /* since Old_vel >0, Ball accelerates */

vel(Ball, Vel, Time, S) :- flying(Ball, S), 
	init_vel(Ball, Old_Vel, S), start(S, T1), Time $>= T1,
    Vel $= Old_Vel + 9.81*(Time - T1). /* since Old_vel < 0, Ball decelerates */

dist(b1, 100, 0, []).
dist(b2, 150, 0, []).
dist(Ball, Dist, Time, S) :- not falling(Ball, S), not flying(Ball, S),
		 init_dist(Ball, Dist, S).

/* Since eplex solves linear constraints only, approximate the law for
   distance with a linear function of time, instead of a quadratic function.
*/
dist(Ball, Dist, Time, S) :- falling(Ball, S), 
	init_dist(Ball, Old_Dist, S), start(S,T1), Time $>= T1,
    Dist $= Old_Dist - 0.5*9.81*(Time - T1). /* Dist decreases if Ball falls down */

dist(Ball, Dist, Time, S) :- flying(Ball, S), 
	init_dist(Ball, Old_Dist, S), start(S,T1), Time $>= T1,
    Dist $= Old_Dist + 0.5*9.81*(Time - T1). /* Dist increases if Ball moves up */

        %--------------Goal conditions and time minimization-------%

firstAction([A], [A]).
/* The following rule makes sure there is no other solution upon backtracking */
firstAction([A1, A2 |S], Sit) :- firstAction([A2 | S], Sit).

duration([A|S], Duration) :- firstAction([A|S], [FirstAct]), 
	time(FirstAct, Tinit), 
	time(A, Tlast), 
    Duration $= (Tlast - Tinit), Duration $>= 0.

/* Linear Programming: schedule actions to execute as soon as possible */
minimizeDuration(TotalDur, Cost) :- 
    eplex_solver_setup(min(TotalDur)), eplex_solve(Cost).

goal_state(S) :- satisfy(S,TotalTime), duration(S,Time), 
	/* Let Delta be extra time that passed since the start of the last action*/
	Delta $>= 0,
	TotalTime $= Time + Delta,
    minimizeDuration(TotalTime, Cost), chooseTimes(S), 
    eplex_get(cost,EndTime), nl, write("EndTime= "), writeln(EndTime).

satisfy(S,Time) :- start(S,Tinit), Time $>= Tinit, 
			falling(b1,S), falling(b2,S), 
	numberFlights(N1,b1,S), N1 >= 1, numberFlights(N2,b2,S), N2 >= 1, 
	vel(b1, V1, Time, S), vel(b2, V2, Time, S), V1 - V2 $= 0,
	dist(b1, D1, Time, S), dist(b2, D2, Time, S), D1 - D2 $= 0 .

		%--------------Auxiliary axioms ---------------------------%

time( drop(Ball,Time), Time). 
time( catch(Ball,Time), Time).
time( bounce(Ball, Time), Time).
time( atPeak(Ball, Time), Time).

initial_state([]).

%Reiter domain independent axioms about start time of situations.
start([], 0).
start([A|S], T) :- time(A,T).

chooseTimes([]).
chooseTimes([A|S]) :- chooseTimes(S), time(A,T), eplex_var_get(T, typed_solution, T).

/* 
ECLiPSe Constraint Logic Programming System [kernel threads]
Kernel and basic libraries copyright Cisco Systems, Inc.
and subject to the Cisco-style Mozilla Public Licence 1.1
(see legal/cmpl.txt or http://eclipseclp.org/licence)
Source available at www.sourceforge.org/projects/eclipse-clp
GMP library copyright Free Software Foundation, see legal/lgpl.txt
For other libraries see their individual copyright notices
Version 7.0 #63 (x86_64_linux), Sun Apr 24 18:46 2022
[eclipse 1]: solve_problem(1, Plan).

No (0.00s cpu)
[eclipse 2]: solve_problem(7, Plan).

No (0.08s cpu)
[eclipse 3]: solve_problem(8, Plan).

EndTime= 91.743119266055

Elapsed time (sec): 0.180136583

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), drop(b1, 50.9683995922528), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), bounce(b1, 71.35575942915392), drop(b2, 91.743119266055047), atPeak(b1, 91.743119266055047)]
Yes (0.18s cpu, solution 1, maybe more) ? 
[eclipse 4]: solve_problem(8, Plan).

EndTime= 91.743119266055

Elapsed time (sec): 0.181632581

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), drop(b1, 50.9683995922528), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), bounce(b1, 71.35575942915392), drop(b2, 91.743119266055047), atPeak(b1, 91.743119266055047)]
Yes (0.18s cpu, solution 1, maybe more) ? ;

EndTime= 91.743119266055

Elapsed time (sec): 0.185122109

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), drop(b1, 50.968399592252787), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), bounce(b1, 71.3557594291539), atPeak(b1, 91.743119266055032), drop(b2, 91.743119266055032)]
Yes (0.19s cpu, solution 2, maybe more) ? ;

EndTime= 101.936799184506

Elapsed time (sec): 0.213021238

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), drop(b1, 61.162079510703364), catch(b2, 71.35575942915392), bounce(b1, 81.549439347604476), drop(b2, 101.9367991845056), atPeak(b1, 101.9367991845056)]
Yes (0.21s cpu, solution 3, maybe more) ? ;

EndTime= 101.936799184506

Elapsed time (sec): 0.216116812

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), drop(b1, 61.162079510703364), catch(b2, 71.355759429153949), bounce(b1, 81.549439347604476), atPeak(b1, 101.9367991845056), drop(b2, 101.9367991845056)]
Yes (0.22s cpu, solution 4, maybe more) ? ;

EndTime= 112.130479102956

Elapsed time (sec): 0.222204706

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), drop(b2, 112.13047910295616), atPeak(b1, 112.13047910295616)]
Yes (0.22s cpu, solution 5, maybe more) ? ;

EndTime= 112.130479102956

Elapsed time (sec): 0.22502153

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), atPeak(b1, 112.13047910295616), drop(b2, 112.13047910295616)]
Yes (0.23s cpu, solution 6, maybe more) ? ;

EndTime= 112.130479102956

Elapsed time (sec): 0.231816381

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), drop(b2, 112.13047910295616), atPeak(b1, 112.13047910295616)]
Yes (0.23s cpu, solution 7, maybe more) ? ;

EndTime= 112.130479102956

Elapsed time (sec): 0.234644073

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), atPeak(b1, 112.13047910295616), drop(b2, 112.13047910295616)]
Yes (0.23s cpu, solution 8, maybe more) ? ;

EndTime= 50.9683995922528

Elapsed time (sec): 0.271247942

Plan = [drop(b2, 0), drop(b1, 10.19367991845056), catch(b2, 10.19367991845056), drop(b2, 10.19367991845056), bounce(b1, 30.581039755351682), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.9683995922528)]
Yes (0.27s cpu, solution 9, maybe more) ? ;

EndTime= 50.9683995922528

Elapsed time (sec): 0.296041525

Plan = [drop(b2, 0), catch(b2, 10.193679918450554), drop(b1, 10.193679918450554), drop(b2, 10.193679918450554), bounce(b1, 30.581039755351679), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.968399592252808)]
Yes (0.30s cpu, solution 10, maybe more) ? ;

EndTime= 50.9683995922528

Elapsed time (sec): 0.322303164

Plan = [drop(b2, 0), catch(b2, 10.19367991845056), drop(b2, 10.19367991845056), drop(b1, 10.19367991845056), bounce(b1, 30.581039755351682), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.9683995922528)]
Yes (0.32s cpu, solution 11, maybe more) ? ;

EndTime= 50.9683995922528

Elapsed time (sec): 0.331785495

Plan = [drop(b2, 0), catch(b2, 10.193679918450554), drop(b1, 10.193679918450554), drop(b2, 10.193679918450554), bounce(b1, 30.581039755351679), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.968399592252808)]
Yes (0.33s cpu, solution 12, maybe more) ? ;

EndTime= 50.9683995922528

Elapsed time (sec): 0.358241188

Plan = [drop(b2, 0), catch(b2, 10.19367991845056), drop(b2, 10.19367991845056), drop(b1, 10.19367991845056), bounce(b1, 30.581039755351682), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.9683995922528)]
Yes (0.36s cpu, solution 13, maybe more) ? ;

EndTime= 91.743119266055

Elapsed time (sec): 0.485687106

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), drop(b1, 50.9683995922528), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), bounce(b1, 71.35575942915392), drop(b2, 91.743119266055047), atPeak(b1, 91.743119266055047)]
Yes (0.49s cpu, solution 14, maybe more) ? ;

EndTime= 91.743119266055

Elapsed time (sec): 0.489053124

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), drop(b1, 50.968399592252787), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), bounce(b1, 71.3557594291539), atPeak(b1, 91.743119266055032), drop(b2, 91.743119266055032)]
Yes (0.49s cpu, solution 15, maybe more) ? ;

EndTime= 101.936799184506

Elapsed time (sec): 0.517981219

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), drop(b1, 61.162079510703364), catch(b2, 71.35575942915392), bounce(b1, 81.549439347604476), drop(b2, 101.9367991845056), atPeak(b1, 101.9367991845056)]
Yes (0.52s cpu, solution 16, maybe more) ? ;

EndTime= 101.936799184506

Elapsed time (sec): 0.521125555

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), drop(b1, 61.162079510703364), catch(b2, 71.355759429153949), bounce(b1, 81.549439347604476), atPeak(b1, 101.9367991845056), drop(b2, 101.9367991845056)]
Yes (0.52s cpu, solution 17, maybe more) ? ;

EndTime= 112.130479102956

Elapsed time (sec): 0.52731797

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), drop(b2, 112.13047910295616), atPeak(b1, 112.13047910295616)]
Yes (0.53s cpu, solution 18, maybe more) ? ;

EndTime= 112.130479102956

Elapsed time (sec): 0.53014014

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), atPeak(b1, 112.13047910295616), drop(b2, 112.13047910295616)]
Yes (0.53s cpu, solution 19, maybe more) ? ;

EndTime= 112.130479102956

Elapsed time (sec): 0.536990292

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), drop(b2, 112.13047910295616), atPeak(b1, 112.13047910295616)]
Yes (0.54s cpu, solution 20, maybe more) ? ;

EndTime= 112.130479102956

Elapsed time (sec): 0.539800036

Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), atPeak(b1, 112.13047910295616), drop(b2, 112.13047910295616)]
Yes (0.54s cpu, solution 21, maybe more) ? ;

EndTime= 50.9683995922528

Elapsed time (sec): 0.576289241

Plan = [drop(b2, 0), drop(b1, 10.19367991845056), catch(b2, 10.19367991845056), drop(b2, 10.19367991845056), bounce(b1, 30.581039755351682), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.9683995922528)]
Yes (0.58s cpu, solution 22, maybe more) ? ;

EndTime= 50.9683995922528

Elapsed time (sec): 0.586033084

Plan = [drop(b2, 0), catch(b2, 10.193679918450554), drop(b1, 10.193679918450554), drop(b2, 10.193679918450554), bounce(b1, 30.581039755351679), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.968399592252808)]
Yes (0.59s cpu, solution 23, maybe more) ? ;

EndTime= 50.9683995922528

Elapsed time (sec): 0.596641336

Plan = [drop(b2, 0), catch(b2, 10.19367991845056), drop(b2, 10.19367991845056), drop(b1, 10.19367991845056), bounce(b1, 30.581039755351682), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.9683995922528)]
Yes (0.60s cpu, solution 24, maybe more) ? ;

EndTime= 50.9683995922528

Elapsed time (sec): 0.606203237

Plan = [drop(b2, 0), catch(b2, 10.193679918450554), drop(b1, 10.193679918450554), drop(b2, 10.193679918450554), bounce(b1, 30.581039755351679), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.968399592252808)]
Yes (0.61s cpu, solution 25, maybe more) ? ;

EndTime= 50.9683995922528

Elapsed time (sec): 0.616902587

Plan = [drop(b2, 0), catch(b2, 10.19367991845056), drop(b2, 10.19367991845056), drop(b1, 10.19367991845056), bounce(b1, 30.581039755351682), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.9683995922528)]
Yes (0.62s cpu, solution 26, maybe more) ? ;

No (0.62s cpu)
[eclipse 5]: halt.
*/
