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:- init_expansors.   50
   51/*******************************
   52		 *	 TOPLEVEL PRINTING	*
   53		 *******************************/
   54
   55:- multifile
   56	prolog:message/3,
   57        prolog_colour:syntax_message//1.   58
   59prolog_colour:syntax_message(constraint(clpcd(Sub))) -->
   60	[ 'clp(~w) constraint'-[Sub] ].
   61prolog_colour:syntax_message(type_error(constraint(clpcd(Sub)))) -->
   62	[ 'Only clp(~w) constraints may appear inside {}'-[Sub] ].
   63
   64prolog:message(query(YesNo,Bindings)) --> !,
   65	{dump_toplevel_bindings(Bindings,Constraints)},
   66	dump_format(Constraints),
   67        '$messages':prolog_message(query(YesNo,Bindings)).
   68
   69cd_invert(sin(  X),T,X,Y,R) :- compare_d(T, =<, abs(Y), 1), int(I), eval_d(T, asin(Y)+2*pi*I, R).
   70cd_invert(cos(  X),T,X,Y,R) :- compare_d(T, =<, abs(Y), 1), int(I), eval_d(T, acos(Y)+2*pi*I, R).
   71cd_invert(tan(  X),T,X,Y,R) :- int(I), eval_d(T,atan(Y)+pi*I,R).
   72cd_invert(asin( X),T,X,Y,R) :- eval_d(T, sin(  Y), R).
   73cd_invert(acos( X),T,X,Y,R) :- eval_d(T, cos(  Y), R).
   74cd_invert(atan( X),T,X,Y,R) :- eval_d(T, tan(  Y), R).
   75cd_invert(sinh( X),T,X,Y,R) :- eval_d(T, asinh(Y), R).
   76cd_invert(cosh( X),T,X,Y,R) :- compare_d(T, >=, Y, 1), eval_d(T, acosh(Y), R).
   77cd_invert(tanh( X),T,X,Y,R) :- eval_d(T, atanh(Y), R).
   78cd_invert(asinh(X),T,X,Y,R) :- eval_d(T, sinh( Y), R).
   79cd_invert(acosh(X),T,X,Y,R) :- eval_d(T, cosh( Y), R).
   80cd_invert(atanh(X),T,X,Y,R) :- eval_d(T, tanh( Y), R).
   81cd_invert(log(  X),T,X,Y,R) :- eval_d(T, exp(  Y), R).
   82cd_invert(log10(X),T,X,Y,R) :- eval_d(T, exp10(Y), R).
   83cd_invert(log1p(X),T,X,Y,R) :- eval_d(T, expm1(Y), R).
   84cd_invert(log2( X),T,X,Y,R) :- eval_d(T, exp2( Y), R).
   85cd_invert(sqrt( X),T,X,Y,R) :- eval_d(T, Y^2, R).
   86cd_invert(cbrt( X),T,X,Y,R) :- eval_d(T, Y^3, R).
   87cd_invert(exp10(X),T,X,Y,R) :- compare_d(T, >=, Y,  0), eval_d(T, log10(Y), R).
   88cd_invert(exp2( X),T,X,Y,R) :- compare_d(T, >=, Y,  0), eval_d(T, log2( Y), R).
   89cd_invert(exp(  X),T,X,Y,R) :- compare_d(T, >=, Y,  0), eval_d(T, log(  Y), R).
   90cd_invert(expm1(X),T,X,Y,R) :- compare_d(T, >=, Y, -1), eval_d(T, log1p(Y), R).
   91cd_invert(abs(  X),T,X,Y,R) :-
   92    ( compare_d(T, <, 0, Y)
   93    ->( R = Y
   94      ; eval_d(T, -Y, R)
   95      )
   96    ; compare_d(T, =, 0, Y)
   97    ->R = Y
   98    ).
   99cd_invert(B^C,T,X,A,R) :-
  100    ( nf_constant(B, Kb)
  101    ->compare_d(T, <, 0, A),
  102      compare_d(T, <, 0, Kb),
  103      % Kb =\= 1.0
  104      compare_d(T, \=, Kb, 1),
  105      X = C, % note delayed unification
  106      eval_d(T, log(A)/log(Kb), R)
  107    ; nf_constant(C,Kc),
  108      compare_d(T, \=, 0, A),
  109      compare_d(T, <, 0, Kc),
  110      X = B, % note delayed unification
  111      eval_d(T, A**(1/Kc), R)
  112    ).
  113cd_invert(hypot(B,C),T,X,A,R) :-
  114    ( nf_constant(B, Kb)
  115    ->cd_invert_hypot(Kb,C,T,X,A,R)
  116    ; nf_constant(C, Kc)
  117    ->cd_invert_hypot(Kc,B,T,X,A,R)
  118    ).
  119
  120cd_invert_hypot(Kb,C,T,X,A,R) :-
  121    compare_d(T, >=, abs(A), abs(Kb)),
  122    X = C,
  123    eval_d(T, sqrt((A+Kb)*(A-Kb)), R).
  124
  125cd_invertible(Term) :-
  126    clause(cd_invert(Term, _, _, _, _), _),
  127    neck.
  128
  129% crafted to get float operators
  130
  131:- public curr_num_arithmetic_function/1.  132
  133:- meta_predicate curr_num_arithmetic_function(?).  134
  135setnum(acosh(_), N, X) :- !, arithmetic_expression_value(1.1*N, X).
  136setnum(_, N, X) :- arithmetic_expression_value(0.8/N, X).
  137
  138curr_num_arithmetic_function(Expr) :-
  139    current_arithmetic_function(Expr),
  140    \+ \+ ( mapnargs(setnum(Expr), Expr),
  141            catch(arithmetic_expression_value(Expr, Value),
  142                  _,
  143                  fail),
  144            float(Value)
  145          ).
  146% Extra arithmetic functions not implemented in is/2 but in most of float libraries
  147curr_num_arithmetic_function(Expr) :-
  148    member(Expr, [cbrt(_), exp10(_), exp2(_), expm1(_), log1p(_),
  149                  log2(_), tgamma(_), hypot(_, _)]),
  150    neck.
  151
  152num_arithmetic_function(Expr) :-
  153    curr_num_arithmetic_function(Expr),
  154    neck.
  155
  156:- public curr_cd_nonlin_af/4.  157
  158curr_cd_nonlin_af(Term, AL, Skel, SL) :-
  159    num_arithmetic_function(Term),
  160    \+ member(Term, [_^_, _**_, +_, -_, _-_, _+_, _/_, _*_,
  161                     pow(_, _), float(_), eval(_)]),
  162    functor(Term, F, A),
  163    functor(Skel, F, A),
  164    A >= 0,
  165    Term =.. [_|AL],
  166    Skel =.. [_|SL].
  167
  168cd_nonlin(pow(A, B), [A, B], X^Y, [X, Y]).
  169cd_nonlin(A^B,       [A, B], X^Y, [X, Y]).
  170cd_nonlin(A**B,      [A, B], X^Y, [X, Y]).
  171cd_nonlin(Term, AL, Skel, SL) :-
  172    curr_cd_nonlin_af(Term, AL, Skel, SL),
  173    neck.
  174
  175dump_toplevel_bindings(Bindings,Constraints) :-
  176	dump_vars_names(Bindings,[],Vars,Names),
  177	dump(Vars,Names,Constraints).
  178
  179dump_vars_names([],_,[],[]).
  180dump_vars_names([Name=Term|Rest],Seen,Vars,Names) :-
  181	(   var(Term),
  182	    (   get_attr(Term,clpcd_itf,_)
  183	    ;   get_attr(Term,clpcd_geler,_)
  184	    ),
  185	    \+ memberchk_eq(Term,Seen)
  186	->  Vars = [Term|RVars],
  187	    Names = [Name|RNames],
  188	    NSeen = [Term|Seen]
  189	;   Vars = RVars,
  190	    Names = RNames,
  191	    Seen = NSeen
  192	),
  193	dump_vars_names(Rest,NSeen,RVars,RNames).
  194
  195dump_format([]) --> [].
  196dump_format([X|Xs]) -->
  197	['{~w}'-[X], nl],
  198	dump_format(Xs).
  199
  200memberchk_eq(X,[Y|Ys]) :-
  201	(   X == Y
  202	->  true
  203	;   memberchk_eq(X,Ys)
  204	)