1/*  Constraint logic programming over continuous domains
    2
    3    Author:        Edison Mera
    4    E-mail:        efmera@gmail.com
    5    WWW:           https://github.com/edisonm/assertions
    6    Copyright (C): 2020, Process Design Center, Breda, The Netherlands.
    7    All rights reserved.
    8
    9    Redistribution and use in source and binary forms, with or without
   10    modification, are permitted provided that the following conditions
   11    are met:
   12
   13    1. Redistributions of source code must retain the above copyright
   14       notice, this list of conditions and the following disclaimer.
   15
   16    2. Redistributions in binary form must reproduce the above copyright
   17       notice, this list of conditions and the following disclaimer in
   18       the documentation and/or other materials provided with the
   19       distribution.
   20
   21    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   22    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   23    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
   24    FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
   25    COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
   26    INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
   27    BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
   28    LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
   29    CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
   30    LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
   31    ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
   32    POSSIBILITY OF SUCH DAMAGE.
   33*/
   34
   35:- module(clpcd_inv,
   36          [cd_invertible/1,
   37           cd_invert/5,
   38           cd_nonlin/4,
   39           num_arithmetic_function/1]).   40
   41:- use_module(library(arithmetic)).   42:- use_module(library(lists)).   43:- use_module(library(neck)).   44:- use_module(library(mapnargs)).   45:- use_module(library(typeprops)).   46:- use_module(library(clpcd/domain_ops)).   47:- use_module(library(clpcd/nf)).   48:- use_module(library(clpcd/dump)).   49
   50		 /*******************************
   51		 *	 TOPLEVEL PRINTING	*
   52		 *******************************/
   53
   54:- multifile
   55	prolog:message/3,
   56        prolog_colour:syntax_message//1.   57
   58prolog_colour:syntax_message(constraint(clpcd(Sub))) -->
   59	[ 'clp(~w) constraint'-[Sub] ].
   60prolog_colour:syntax_message(type_error(constraint(clpcd(Sub)))) -->
   61	[ 'Only clp(~w) constraints may appear inside {}'-[Sub] ].
   62
   63prolog:message(query(YesNo,Bindings)) --> !,
   64	{dump_toplevel_bindings(Bindings,Constraints)},
   65	dump_format(Constraints),
   66        '$messages':prolog_message(query(YesNo,Bindings)).
   67
   68cd_invert(sin(  X),T,X,Y,R) :- compare_d(T, =<, abs(Y), 1), int(I), eval_d(T, asin(Y)+2*pi*I, R).
   69cd_invert(cos(  X),T,X,Y,R) :- compare_d(T, =<, abs(Y), 1), int(I), eval_d(T, acos(Y)+2*pi*I, R).
   70cd_invert(tan(  X),T,X,Y,R) :- int(I), eval_d(T,atan(Y)+pi*I,R).
   71cd_invert(asin( X),T,X,Y,R) :- eval_d(T, sin(  Y), R).
   72cd_invert(acos( X),T,X,Y,R) :- eval_d(T, cos(  Y), R).
   73cd_invert(atan( X),T,X,Y,R) :- eval_d(T, tan(  Y), R).
   74cd_invert(sinh( X),T,X,Y,R) :- eval_d(T, asinh(Y), R).
   75cd_invert(cosh( X),T,X,Y,R) :- compare_d(T, >=, Y, 1), eval_d(T, acosh(Y), R).
   76cd_invert(tanh( X),T,X,Y,R) :- eval_d(T, atanh(Y), R).
   77cd_invert(asinh(X),T,X,Y,R) :- eval_d(T, sinh( Y), R).
   78cd_invert(acosh(X),T,X,Y,R) :- eval_d(T, cosh( Y), R).
   79cd_invert(atanh(X),T,X,Y,R) :- eval_d(T, tanh( Y), R).
   80cd_invert(log(  X),T,X,Y,R) :- eval_d(T, exp(  Y), R).
   81cd_invert(log10(X),T,X,Y,R) :- eval_d(T, exp10(Y), R).
   82cd_invert(log1p(X),T,X,Y,R) :- eval_d(T, expm1(Y), R).
   83cd_invert(log2( X),T,X,Y,R) :- eval_d(T, exp2( Y), R).
   84cd_invert(sqrt( X),T,X,Y,R) :- eval_d(T, Y^2, R).
   85cd_invert(cbrt( X),T,X,Y,R) :- eval_d(T, Y^3, R).
   86cd_invert(exp10(X),T,X,Y,R) :- compare_d(T, >=, Y,  0), eval_d(T, log10(Y), R).
   87cd_invert(exp2( X),T,X,Y,R) :- compare_d(T, >=, Y,  0), eval_d(T, log2( Y), R).
   88cd_invert(exp(  X),T,X,Y,R) :- compare_d(T, >=, Y,  0), eval_d(T, log(  Y), R).
   89cd_invert(expm1(X),T,X,Y,R) :- compare_d(T, >=, Y, -1), eval_d(T, log1p(Y), R).
   90cd_invert(abs(  X),T,X,Y,R) :-
   91    ( compare_d(T, <, 0, Y)
   92    ->( R = Y
   93      ; eval_d(T, -Y, R)
   94      )
   95    ; compare_d(T, =, 0, Y)
   96    ->R = Y
   97    ).
   98cd_invert(B^C,T,X,A,R) :-
   99    ( nf_constant(B, Kb)
  100    ->compare_d(T, <, 0, A),
  101      compare_d(T, <, 0, Kb),
  102      % Kb =\= 1.0
  103      compare_d(T, \=, Kb, 1),
  104      X = C, % note delayed unification
  105      eval_d(T, log(A)/log(Kb), R)
  106    ; nf_constant(C,Kc),
  107      compare_d(T, \=, 0, A),
  108      compare_d(T, <, 0, Kc),
  109      X = B, % note delayed unification
  110      eval_d(T, A**(1/Kc), R)
  111    ).
  112cd_invert(hypot(B,C),T,X,A,R) :-
  113    ( nf_constant(B, Kb)
  114    ->cd_invert_hypot(Kb,C,T,X,A,R)
  115    ; nf_constant(C, Kc)
  116    ->cd_invert_hypot(Kc,B,T,X,A,R)
  117    ).
  118
  119cd_invert_hypot(Kb,C,T,X,A,R) :-
  120    compare_d(T, >=, abs(A), abs(Kb)),
  121    X = C,
  122    eval_d(T, sqrt((A+Kb)*(A-Kb)), R).
  123
  124cd_invertible(Term) :-
  125    clause(cd_invert(Term, _, _, _, _), _),
  126    neck.
  127
  128% crafted to get float operators
  129
  130:- public curr_num_arithmetic_function/1.  131
  132:- meta_predicate curr_num_arithmetic_function(?).  133
  134setnum(acosh(_), N, X) :- !, arithmetic_expression_value(1.1*N, X).
  135setnum(_, N, X) :- arithmetic_expression_value(0.8/N, X).
  136
  137curr_num_arithmetic_function(Expr) :-
  138    current_arithmetic_function(Expr),
  139    \+ \+ ( mapnargs(setnum(Expr), Expr),
  140            catch(arithmetic_expression_value(Expr, Value),
  141                  _,
  142                  fail),
  143            float(Value)
  144          ).
  145% Extra arithmetic functions not implemented in is/2 but in most of float libraries
  146curr_num_arithmetic_function(Expr) :-
  147    member(Expr, [cbrt(_), exp10(_), exp2(_), expm1(_), log1p(_),
  148                  log2(_), tgamma(_), hypot(_, _)]),
  149    neck.
  150
  151num_arithmetic_function(Expr) :-
  152    curr_num_arithmetic_function(Expr),
  153    neck.
  154
  155:- public curr_cd_nonlin_af/4.  156
  157curr_cd_nonlin_af(Term, AL, Skel, SL) :-
  158    num_arithmetic_function(Term),
  159    \+ member(Term, [_^_, _**_, +_, -_, _-_, _+_, _/_, _*_,
  160                     pow(_, _), float(_), eval(_)]),
  161    functor(Term, F, A),
  162    functor(Skel, F, A),
  163    A >= 0,
  164    Term =.. [_|AL],
  165    Skel =.. [_|SL].
  166
  167cd_nonlin(pow(A, B), [A, B], X^Y, [X, Y]).
  168cd_nonlin(A^B,       [A, B], X^Y, [X, Y]).
  169cd_nonlin(A**B,      [A, B], X^Y, [X, Y]).
  170cd_nonlin(Term, AL, Skel, SL) :-
  171    curr_cd_nonlin_af(Term, AL, Skel, SL),
  172    neck.
  173
  174dump_toplevel_bindings(Bindings,Constraints) :-
  175	dump_vars_names(Bindings,[],Vars,Names),
  176	dump(Vars,Names,Constraints).
  177
  178dump_vars_names([],_,[],[]).
  179dump_vars_names([Name=Term|Rest],Seen,Vars,Names) :-
  180	(   var(Term),
  181	    (   get_attr(Term,clpcd_itf,_)
  182	    ;   get_attr(Term,clpcd_geler,_)
  183	    ),
  184	    \+ memberchk_eq(Term,Seen)
  185	->  Vars = [Term|RVars],
  186	    Names = [Name|RNames],
  187	    NSeen = [Term|Seen]
  188	;   Vars = RVars,
  189	    Names = RNames,
  190	    Seen = NSeen
  191	),
  192	dump_vars_names(Rest,NSeen,RVars,RNames).
  193
  194dump_format([]) --> [].
  195dump_format([X|Xs]) -->
  196	['{~w}'-[X], nl],
  197	dump_format(Xs).
  198
  199memberchk_eq(X,[Y|Ys]) :-
  200	(   X == Y
  201	->  true
  202	;   memberchk_eq(X,Ys)
  203	)