// This program illustrates the principle of interval constraint propagation
// Consider an electronic circuit made with one battery E and deux resistors R1 and R2
// The current is denoted by I, the voltages at the bounds of the resistors by  
// R1 and R2, and the power delivered by the battery by P
// These quantities are related by the following (primitive) constraints
// 	E=R*I,R=R1+R2,E*I=P,U2=R2*I,U1=R1*I,U1+U2=E 
// For information, the true values for the variables are : 
//   	R1=2,R3=3,R=5,I=5,E=25,U1=10,U2=15,P=125.

//--------------------------------------------------------
/////////////// INTERVAL ARITHMETIC
//--------------------------------------------------------
function z=inter(x,y)
	z=[max(x(1),y(1)),min(x(2),y(2))]
	if (z(2)<z(1)) then write ('intervalle vide'); end;
endfunction
//-----------------------------------------
function z=add(x,y); 	
	z=x+y; 
endfunction
//-----------------------------------------
function z=dif(x,y)
	z=[x(1)-y(2),x(2)-y(1)];
endfunction
//-----------------------------------------
function z=mult(x,y)
	v=[x(1)*y(1),x(1)*y(2),x(2)*y(1),x(2)*y(2)];
	z=[min(v),max(v)];
endfunction
//-----------------------------------------
function y=invert(x)
	if and([x(1)<=0,x(2)>=0]) then y=[-%inf,%inf];
     				  else y=[1/x(2),1/x(1)];
	end 
endfunction
//-----------------------------------------
function z=div(x,y)
	a=invert(y); 
	z=mult(x,a);
endfunction
//-----------------------------------------------------------------------
/////////////// PROJECTION OPERATORS
//-----------------------------------------------------------------------
function [z1,x1,y1]=padd(z,x,y)
	z1=inter(z,add(x,y));
	x1=inter(x,dif(z,y));
	y1=inter(y,dif(z,x));
endfunction
//-----------------------------------------
function [z1,x1,y1]=pmult(z,x,y)
	z1=inter(z,mult(x,y))
	x1=inter(x,div(z,y))
	y1=inter(y,div(z,x))
endfunction
//------------------------------------------------------------------------
/////////////// 	MAIN PROGRAM
//------------------------------------------------------------------------

// prior domains for the variables : 
R=[-%inf,%inf]; R1=[0,%inf]; R2=[0,%inf]; E=[23,26]; I=[4,8];
U1=[10,11]; U2=[14,17]; P=[124,130];

// Propagation
for i=1:3,
     [R,R1,R2]=padd(R,R1,R2);
     [P,E,I]=pmult(P,E,I);
     [E,R,I]=pmult(E,R,I);
     [U2,R2,I]=pmult(U2,R2,I);
     [U1,R1,I]=pmult(U1,R1,I);
     [E,U1,U2]=padd(E,U1,U2);
end; 

// display the posterior domains
R1,R2,I,U1,U2,E,P