%% Redlog file for NAV (Navigation benchmarks) relational abstraction.
%% Results:
%% For continuous mode that goes EAST; i.e.,
%% 4-d system with dynamics:
%% dx/dt = vx,  dvx/dt = -1.2 (vx-vx0) + 0.1 (vy-vy0)
%% dy/dt = vy,  dvy/dt =  0.1 (vx-vx0) - 1.2 (vy-vy0)
%% Assume we are going east; so vx0 = 1, vy0 = 0
%% We get following invariants: R(xinit,yinit,vxinit,vyinit,x,y,vx,vy):
%%  (vxinit-1)+-vyinit >= 0 => 0 <= (vx-1) +- vy <= (vxinit-1) +- vyinit
%%  (vxinit-1)+-vyinit <= 0 => 0 >= (vx-1) +- vy >= (vxinit-1) +- vyinit
%%         +- means we can choose either + or -
%%  (y - yinit) + 10/143*(vx - vxinit) + 120/143*(vy - vyinit) = 0
%%  (x - xinit) + (vx - vxinit) + 2*(vy - vyinit) + 2.3*(y - yinit) >= 0
%% There are infinitely many invariants of the last form.
%% I picked coefficients 1,2,2.3 arbitrarily.
%% We just need to ensure b=12a-10 and 10c = 12b-a
%% where a is coeff of (vx-vxinit), b of (vy-vyinit), c of (y-yinit)

%% Can we synthesize a relational invariant for this.
%% Template:
%%  (vxinit-1)+a*vyinit >= 0 => 0 <= (vx-1) + a*(vy) <= (vxinit-1) + a*vyinit
%%  (vxinit-1)+a*vyinit <= 0 => 0 >= (vx-1) + a*(vy) >= (vxinit-1) + a*vyinit

%% Result: Quantifier elimination gives
%% [ a + 1 >= 0 /\ a - 1 <= 0 /\ [ a - 1 = 0 \/ a + 1 = 0 ] ]
%% which means either a = 1 or a = -1. Both work.  
%% We can easily check that [1;1] and [1;-1] are both eigenvectors.
%% We get the invariants:
%%  (vxinit-1)+-vyinit >= 0 => 0 <= (vx-1) +- vy <= (vxinit-1) +- vyinit
%%  (vxinit-1)+-vyinit <= 0 => 0 >= (vx-1) +- vy >= (vxinit-1) +- vyinit

%% Now we need invariants on x and y.
%% x + a*vx + b*vy = xinit + a*vxinit + b*vyinit ? 


on echo;
load_package redlog;
rlset r;
load_package qepcad;
on rlqefbqepcad;
rlqepcadn(100*10^6);
rlqepcadl(200*10^3);
on rlverbose;
on time;
%off rlqefb;

%% Init is trivially true
phi1 := (
  (vxinit-1)+a*vyinit >= 0 impl  0 <= (vxinit-1) + a*vyinit and
  (vxinit-1)+a*vyinit <= 0 impl  0 >= (vxinit-1) + a*vyinit );

phi2 := ((
  (vxinit-1)+a*vyinit >= 0 and 
  (vx-1) + a*vy = (vxinit-1) + a*vyinit )
  impl
  (-12/10 * (vx-1) + 1/10 * vy + a * ( 0.1 * (vx - 1) - 12/10 * vy ) <= 0));
  
phi3 := ((
  (vxinit-1)+a*vyinit >= 0 and 
  (vx-1) + a*vy = 0 )
  impl
  (-12/10 * (vx-1) + 1/10 * vy + a * ( 0.1 * (vx - 1) - 12/10 * vy ) >= 0));

phi4 := ((
  (vxinit-1)+a*vyinit <= 0 and 
  (vx-1) + a*vy = (vxinit-1) + a*vyinit )
  impl
  (-12/10 * (vx-1) + 1/10 * vy + a * ( 0.1 * (vx - 1) - 12/10 * vy ) >= 0));
  
phi5 := (
 (
  (vxinit-1)+a*vyinit <= 0 and 
  (vx-1) + a*vy = 0 )
  impl
  (-12/10 * (vx-1) + 1/10 * vy + a * ( 0.1 * (vx - 1) - 12/10 * vy ) <= 0) );

%% y + a*(vx-1) + b*vy = yinit + a*(vxinit-1) + b*vyinit
phi6 := ( 
 ( vy + a*( -12/10*(vx-1) + 1/10 * vy) + b*( 1/10 * (vx-1) - 12/10 * vy) = 0));
%% We get a=10/143, b=120/143, ... y + 10/143*vx + 120/143*vy is a constant.

%% x + a*vx + b*vy = xinit + a*vxinit + b*vyinit ? 
phi6 := ( vx + a*( -12/10*(vx-1) + 1/10 * vy) + b*( 1/10 * (vx-1) - 12/10 * vy) = 0);
% No equality invariant of this form.

%% x + c*y + a*vx + b*vy = xinit + c*yinit + a*vxinit + b*vyinit ? 
phi6 := ( vx + c*vy + a*( -12/10*(vx-1) + 1/10 * vy) + b*( 1/10 * (vx-1) - 12/10 * vy) = 0);
% No equality invariant of this form too.

%% x + c*y + a*(vx-1) + b*vy <= xinit + c*yinit + a*(vxinit-1) + b*vyinit ? 
phi6 := ( 
 xinit + c*yinit + a*(vxinit-1) + b*vyinit >= 0 and
 x + c*y + a*vx + b*vy = xinit + c*yinit + a*vxinit + b*vyinit impl 
 ( vx + c*vy + a*( -12/10*(vx-1) + 1/10 * vy) + b*( 1/10 * (vx-1) - 12/10 * vy) <= 0));
%% We get FALSE... so no such invariant with a in there....

%% x + c*y + a*(vx-1) + b*vy >= xinit + c*yinit + a*(vxinit-1) + b*vyinit
phi6 := ( 
 ( vx + c*vy + a*( -12/10*(vx-1) + 1/10 * vy) + b*( 1/10 * (vx-1) - 12/10 * vy) >= 0));
%% We get b=12a-10 and 10c = 12b-a;  c = 14.3a-12; 
%% x + vx + 2*vy + 2.3*y >= its initial value...

%% x + c*y + a*(vx-1) + b*vy <= xinit + c*yinit + a*(vxinit-1) + b*vyinit
phi6 := ( 
 ( vx + c*vy + a*( -12/10*(vx-1) + 1/10 * vy) + b*( 1/10 * (vx-1) - 12/10 * vy) <= 0));
%% THIS IS FALSE!!!

%% For going south....
%% x + a*(vx-1) + b*vy = xinit + a*(vxinit-1) + b*vyinit
phi6 := ( 
 ( vx + a*( -12/10*vx + 1/10*(vy+1)) + b*( 1/10*vx - 12/10*(vy+1)) = 0));

%% y + c*x + a*vx + b*vy >= yinit + c*xinit + a*vxinit + b*vyinit
phi6 := ( 
 ( vy + c*vx + a*( -12/10*vx + 1/10 *(vy+1)) + b*( 1/10 * vx - 12/10 *(vy+1)) <= 0));

%% For going south-east....
%% x + a*(vx-1) + b*vy + c*y = xinit + a*(vxinit-1) + b*vyinit + c*yinit
phi6 := ( 
 ( vx + a*( -12/10*(vx-7/10) + 1/10*(vy+7/10)) + b*( 1/10*(vx-7/10) - 12/10*(vy+7/10)) + c*vy = 0));

%% For going south-east....
%% a*vx + b*vy + y <= 0
phi6 := ( 
 ( a*( -12/10*(vx-7/10) + 1/10*(vy+7/10)) + b*( 1/10*(vx-7/10) - 12/10*(vy+7/10)) + vy <= 0));

%% For going south-east....
%% a*vx + b*vy + x >= 0
phi6 := ( 
 ( a*( -12/10*(vx-7/10) + 1/10*(vy+7/10)) + b*( 1/10*(vx-7/10) - 12/10*(vy+7/10)) + vx >= 0));


%% We get a=10/143, b=120/143, ... y + 10/143*vx + 120/143*vy is a constant.
phi := phi1 and phi2 and phi3 and phi4 and phi5;
phi := phi6;

phi := all({x,y,vx,vy,vxinit,vyinit,xinit,yinit}, phi);

psi := rlqe(phi);

rlqepcad( rlex(psi), "nav.qepcad");

end;