%
% local propagation evaluator
% for nonlinear expressions in clp(r,q)
%

:- use_module(engine(internals), [term_to_meta/2]).

arith_ground(X) :- 			% substitute for number(X)
  type(X,Xt), normalize(Xt, X, _, _, H),
  H = [].

arith_ground(X, Val) :-
  type(X,Xt), normalize(Xt, X, K, I, H),
  H = [],
  arith_eval(K*I, Val).

:- meta_predicate delay_all(_,_,goal).

delay_all( yes, _, _).
delay_all( no, Vs, G) :- delay_all( Vs, G).

:- meta_predicate delay_all(_,goal).

delay_all( [], _).
delay_all( [V|Vs], G) :-
  system_freeze( V, G),
  delay_all( Vs, G).

% a = abs(b)
%
solve_abs( A, B) :-
  solve_abs( _Mutex, A, B, no).
%
solve_abs( Mux, A, B, Retry) :- var(Mux), !,
  ( arith_ground( A, Av) ->
      % arith_eval( A >= 0),
      ( arith_zero( Av) ->
          Mux = solved,
          solve_one( B, 0)
      ; arith_ground( B, Bv) ->   			% +a ?b
          Mux = solved,
          arith_zero( abs(Bv) - Av)
      ;                                       		% +a -b
          delay_all( Retry, [B], solve_abs(Mux,A,B,yes))
      )
  ; arith_ground( B, Bv) ->                   		% -a +b
      Mux = solved,
      arith_eval( abs(Bv), Val), solve_one( A, Val)
  ;                                                 	% -a -b
      delay_all( Retry, [A,B], solve_abs(Mux,A,B,yes))
  ).
solve_abs( _, _, _, _).

normalize_abs( B, 1, Inhom, Hom) :-
  type(B,Bt), normalize( Bt, B, Kb, Ib, Hb),
  ( Hb = [] ->                   			% -a +b
      arith_eval( abs(Kb*Ib), Inhom), Hom = []
  ;                                                 	% -a -b
      Inhom = 0, Hom = [VA*1],
      eqn_var_new( v, VA),
      var_with_def( VB, v, Kb, Ib, Hb),
      delay_all( [VA,VB], solve_abs(_Mux,VA,VB,yes))
  ).

% a = {min,max}(b,c)
%
solve_mix( MIX, A, B, C) :-
  solve_mix( _Mutex, MIX, A, B, C, no).
%
solve_mix( Mux, MIX, A, B, C, Retry) :- var( Mux), !,
  ( arith_ground( A, Av) ->
     ( arith_ground( B, Bv) ->
        ( arith_ground( C, Cv) ->   			% +a +b +c
	    Mux = solved,
 	    ( MIX = min, arith_zero( Av-min(Bv,Cv))
            ; MIX = max, arith_zero( Av-max(Bv,Cv))
            )
        ;                               		% +a +b -c
            delay_all( Retry, [C], solve_mix(Mux,MIX,A,B,C,yes))
        )
     ; arith_ground( C, Cv) ->                  	% +a -b +c
        delay_all( Retry, [B], solve_mix(Mux,MIX,A,B,C,yes))
     ;                                      		% +a -b -c
       delay_all( Retry, [B,C], solve_mix(Mux,MIX,A,B,C,yes))
     )
  ; arith_ground( B, Bv) ->
     ( arith_ground( C, Cv) ->      			% -a +b +c
         Mux = solved,
 	 ( MIX = min, arith_eval( min(Bv,Cv), Av)
         ; MIX = max, arith_eval( max(Bv,Cv), Av)
         ),
         solve_one( A, Av)
     ;                               			% -a +b -c
         delay_all( Retry, [A,C], solve_mix(Mux,MIX,A,B,C,yes))
     )
  ; arith_ground( C, Cv) ->    				% -a -b +c
      delay_all( Retry, [A,B], solve_mix(Mux,MIX,A,B,C,yes))
  ;                                        		% -a -b -c
      delay_all( Retry, [A,B,C], solve_mix(Mux,MIX,A,B,C,yes))
  ).
solve_mix( _, _, _, _, _, _).

normalize_mix( MIX, B, C, 1, Inhom, Hom) :-
  type(B,Bt), normalize( Bt, B, Kb, Ib, Hb),
  type(C,Ct), normalize( Ct, C, Kc, Ic, Hc),
  ( Hb = [] ->
     ( Hc = [] ->      					% -a +b +c
         ( MIX = min, arith_eval( min(Kb*Ib,Kc*Ic), Inhom)
         ; MIX = max, arith_eval( max(Kb*Ib,Kc*Ic), Inhom)
         ),
         Hom = []
     ;                               			% -a +b -c
         Inhom = 0, Hom = [VA*1],
         eqn_var_new( v, VA),
         arith_eval( Kb*Ib, Vb),
         var_with_def( VC, v, Kc, Ic, Hc),
         delay_all( [VA,VC], solve_mix(Mux,MIX,VA,Vb,VC,yes))
     )
  ; Hc = [] ->    					% -a -b +c
      Inhom = 0, Hom = [VA*1],
      eqn_var_new( v, VA),
      var_with_def( VB, v, Kb, Ib, Hb),
      arith_eval( Kc*Ic, Vc),
      delay_all( [VA,VB], solve_mix(Mux,MIX,VA,VB,Vc,yes))
  ;                                        		% -a -b -c
      Inhom = 0, Hom = [VA*1],
      eqn_var_new( v, VA),
      var_with_def( VB, v, Kb, Ib, Hb),
      var_with_def( VC, v, Kc, Ic, Hc),
      delay_all( [VA,VB,VC], solve_mix(_Mux,MIX,VA,VB,VC,yes))
  ).

% a = b * c
%
solve_mult( A, B, C) :-
  solve_mult( _Mutex, A, B, C, no).
%
solve_mult( Mux, A, B, C, Retry) :-
  var( Mux),
  !,
  ( arith_ground( A, Av) ->
     ( arith_ground( B, Bv) ->       			% +a +b ?c
	Mux = solved,
	( arith_zero( Bv) ->      arith_zero( Av)	% otherwise: zero division
	; arith_eval( Av/Bv, Vc), solve_one( C, Vc)
        )
     ; arith_ground( C, Cv) ->                  	% +a -b +c
	Mux = solved,
        ( arith_zero( Cv)   ->    arith_zero( Av)	% otherwise: zero division
	; arith_eval( Av/Cv, Vb), solve_one( B, Vb)
        )
     ;                                      		% +a -b -c
       delay_all( Retry, [B,C], solve_mult(Mux,A,B,C,yes))
     )
  ; arith_ground( B, Bv) ->
     Mux = solved,
     ( arith_ground( C, Cv) ->      			% -a +b +c
        arith_eval( Bv*Cv, Va), solve_one( A, Va)
     ; arith_zero( Bv) ->   	            		% -a +b -c
 	solve_one( A, 0)
     ;
        solve_two( A, C, Bv)
     )
  ; arith_ground( C, Cv) ->    				% -a -b +c
     Mux = solved,
     ( arith_zero( Cv) ->
        solve_one( A, 0)
     ;
        solve_two( A, B, Cv)
     )
  ;                                        		% -a -b -c
      delay_all( Retry, [A,B,C], solve_mult(Mux,A,B,C,yes))
  ).
solve_mult( _, _, _, _, _).

normalize_mult( B,C, K, Inhom, Hom) :-
  type(B,Bt), normalize( Bt, B, Kb, Ib, Hb),
  type(C,Ct), normalize( Ct, C, Kc, Ic, Hc),
  ( Hb = [] ->
     ( Hc = [] ->      					% -a +b +c
        arith_eval( Kb*Kc, K), arith_eval( Ib*Ic, Inhom), Hom = []
     ; arith_zero( Ib) ->
 	K = 1, Inhom = 0, Hom = []
     ;                               			% -a +b -c
        arith_eval( Kb*Ib*Kc, K), Inhom = Ic, Hom = Hc
     )
  ; Hc = [] ->    					% -a -b +c
     ( arith_zero( Ic) ->
        K = 1, Inhom = 0, Hom = []
     ;
        arith_eval( Kc*Ic*Kb, K), Inhom = Ib, Hom = Hb
     )
  ;                                        		% -a -b -c
      arith_eval( Kb*Kc, K), Inhom = 0, Hom = [VA*1],
      eqn_var_new( v, VA),
      var_with_def( VB, v, 1, Ib, Hb),
      var_with_def( VC, v, 1, Ic, Hc),
      delay_all( [VA,VB,VC], solve_mult(_Mux,VA,VB,VC,yes))
  ).

normalize_div( B,C, K, Inhom, Hom) :-			% a = b/c
  type(B,Bt), normalize( Bt, B, Kb, Ib, Hb),
  type(C,Ct), normalize( Ct, C, Kc, Ic, Hc),
  ( Hb = [] ->
     ( Hc = [] ->      					% -a +b +c
        ( arith_zero( Ic) ->
            fail  					% zero divison
        ;
            arith_eval( Kb/Kc, K), arith_eval( Ib/Ic, Inhom), Hom = []
        )
     ; arith_zero( Ib) ->
        var_with_def( VC, nz, 1, Ic, Hc),  		% nonzero( VC)
 	K = 1, Inhom = 0, Hom = []
     ;                               			% -a +b -c
        K = Kb, Inhom = 0, Hom = [VA*1],
        eqn_var_new( v, VA),
        var_with_def( VC, v, Kc, Ic, Hc),
        delay_all( [VA,VC], solve_mult(_Mux,Ib,VA,VC,yes))
     )
  ; Hc = [] ->    					% -a -b +c
     ( arith_zero( Ic) ->
        fail						% zero division
     ;
        arith_eval( Kb/(Kc*Ic), K), Inhom = Ib, Hom = Hb
     )
  ;                                        		% -a -b -c
      arith_eval( Kb/Kc, K), Inhom = 0, Hom = [VA*1],
      eqn_var_new( v, VA),
      var_with_def( VB, v, 1, Ib, Hb),
      var_with_def( VC, nz, 1, Ic, Hc),  		% nonzero( VC)
      delay_all( [VA,VB,VC], solve_mult(_Mux,VB,VA,VC,yes))
  ).


% a = b^c
%
solve_pow( A, B, C) :-
  solve_pow( _Mutex, A, B, C, no).
%
solve_pow( Mux, A, B, C, Retry) :- var(Mux), !,
  ( arith_ground( A, Av) ->      				% +a ?b ?c
     ( arith_zero( Av) ->
        Mux = solved, solve_one( B, 0)
     ; arith_ground( B, Bv) ->                     		% +a +b ?c
	Mux = solved,
	( arith_zero( Bv) ->                fail		% A=0, \+ zero(Va)
        ; arith_zero( 1-Bv) ->              arith_zero( 1-Av)
	; arith_eval( log(Av)/log(Bv), Vc), solve_one( C, Vc)
        )
     ; arith_ground( C, Cv) ->                    		% +a -b +c
	Mux = solved,
        ( arith_zero( Cv)   ->        arith_zero( 1-Av)
        ; arith_zero( 1-Cv) ->        solve_one( B, Av)
	; arith_eval( exp(Av,1/Cv), Vb), solve_one( B, Vb)
        )
     ;                                             		% +a -b -c
       delay_all( Retry, [B,C], solve_pow(Mux,A,B,C,yes))
     )
  ; arith_ground( B, Bv) ->
     ( arith_ground( C, Cv) ->        				% -a +b +c
        Mux = solved, arith_eval( exp(Bv,Cv), Va), solve_one( A, Va)
     ;                                         			% -a +b -c
        delay_all( Retry, [A,C], solve_pow(Mux,A,B,C,yes))
     )
  ; arith_ground( C, Cv) ->        				% -a -b +c
     ( arith_zero( Cv) ->
        Mux = solved,
	% nonzero( B),
        solve_one( A, 1)
     ; arith_zero( 1-Cv) ->
        Mux = solved, solve_two( A, B, 1)
     ;
        delay_all( Retry, [A,B], solve_pow(Mux,A,B,C,yes))
     )
  ;                                                  		% -a -b -c
      delay_all( Retry, [A,B,C], solve_pow(Mux,A,B,C,yes))
  ).
solve_pow( _, _, _, _, _).

normalize_pow( B, C, K, Inhom, Hom) :-
  type(B,Bt), normalize( Bt, B, Kb, Ib, Hb),
  type(C,Ct), normalize( Ct, C, Kc, Ic, Hc),
  ( Hb = [] ->
     ( Hc = [] ->        					% -a +b +c
        Mux = solved, K = 1, arith_eval( exp(Kb*Ib,Kc*Ic), Inhom), Hom = []
     ;                                         			% -a +b -c
        K = 1, Inhom = 0, Hom = [VA*1],
        eqn_var_new( v, VA),
        arith_eval( Kb*Ib, Vb),
        var_with_def( VC, v, Kc, Ic, Hc),
        delay_all( [VA,VC], solve_pow(_Mux,VA,Vb,VC,yes))
     )
  ; Hc = [] ->        						% -a -b +c
     ( arith_zero( Ic) ->
        Mux = solved,
	var_with_def( VB, nz, 1, Ib, Hb),  			% nonzero( VB)
        K = 1, Inhom = 1, Hom = []
     ; arith_zero( 1-Kc*Ic) ->
        Mux = solved, K = Kb, Inhom = Ib, Hom = Hb
     ;
        K = 1, Inhom = 0, Hom = [VA*1],
        eqn_var_new( v, VA),
        var_with_def( VB, v, Kb, Ib, Hb),
        arith_eval( Kc*Ic, Vc),
        delay_all( [VA,VB], solve_pow(_Mux,VA,VB,Vc,yes))
     )
  ;                                                  		% -a -b -c
      K = 1, Inhom = 0, Hom = [VA*1],
      eqn_var_new( v, VA),
      var_with_def( VB, v, Kb, Ib, Hb),
      var_with_def( VC, v, Kc, Ic, Hc),
      delay_all( [VA,VB,VC], solve_pow(_Mux,VA,VB,VC,yes))
  ).

% a = TRIG(b)
%
solve_trig( Trig, A, B) :-
  solve_trig( _Mutex, Trig, A, B, no).
%
solve_trig( Mux, Trig, A, B, Retry) :- var( Mux), !,
  ( arith_ground( A, Av) ->               		% +a ?b
      Mux = solved,
      trig( Trig, y, Av, Val), solve_one( B, Val)
  ; arith_ground( B, Bv) ->                     	% -a +b
      Mux = solved,
      trig( Trig, x, Bv, Val), solve_one( A, Val)
  ;                                                   	% -a -b
      delay_all( Retry, [A,B], solve_trig( Mux, Trig, A, B, yes))
  ).
solve_trig( _, _, _, _, _).

trig( sin, x, X, Y) :- arith_eval( sin(X),  Y).
trig( sin, y, X, Y) :- arith_eval( asin(X), Y).
trig( cos, x, X, Y) :- arith_eval( cos(X),  Y).
trig( cos, y, X, Y) :- arith_eval( acos(X), Y).
trig( tan, x, X, Y) :- arith_eval( tan(X),  Y).
trig( tan, y, X, Y) :- arith_eval( atan(X), Y).

normalize_trig( Trig, B, 1, Inhom, Hom) :-
  type(B,Bt), normalize( Bt, B, K, Ib, Hb),
  ( Hb = [] ->                     			% -a +b
      arith_eval( K*Ib, Ibk),
      trig( Trig, x, Ibk, Inhom), Hom = []
  ;                                                   	% -a -b
      Inhom = 0, Hom = [VA*1],
      eqn_var_new( v, VA),
      var_with_def( VB, v, K, Ib, Hb),
      delay_all( [VA,VB], solve_trig(_Mux,Trig,VA,VB,yes))
  ).


% See clpr_freeze/2 in clprt.pl

system_freeze(Var, Goal) :-
  % nonvar(Goal),					% assert
  term_to_meta(G,Goal),
  attach_attribute(V, clpr_frozen(V,G,true)),
  Var = V.
