1:- reexport('../r'), r_initialize.    2
    3:- use_module(cleaning).    4
    5%
    6% Skip R vectors
    7%
    8int_hook(:, colon(_, _), _, []).
    9colon(A, A).
   10
   11%
   12% Obtain atoms or functions from R
   13%
   14eval_hook(r(Expr), Res) :-
   15    eval_hook(Expr, Res).
   16
   17eval_hook(r(Expr), Res) :-
   18    !,
   19    r(Expr, Res).
   20
   21eval_hook(Atom, Res) :-
   22    atomic(Atom),
   23    r_hook(R, Atom),
   24    !,
   25    call(R, Atom, Res).
   26
   27eval_hook(Atom, Res) :-
   28    atomic(Atom),
   29    r_hook(Atom),
   30    !,
   31    r(Atom, Res).
   32
   33eval_hook(Expr, Res) :-
   34    compound(Expr),
   35    compound_name_arity(Expr, Name, Arity),
   36    r_hook(R, Name/Arity),
   37    !,
   38    call(R, Expr, Res).
   39
   40eval_hook(Expr, Res) :-
   41    compound(Expr),
   42    compound_name_arity(Expr, Name, Arity),
   43    r_hook(Name/Arity),
   44    !,
   45    r(Expr, Res).
   46
   47%
   48% Call R 
   49%
   50int_hook(r, r1(atomic), _, [evaluate(false)]).
   51r1(atomic(A), Res, _Flags) :-
   52    eval(r(A), Res1),
   53    !,
   54    clean(Res1, Res).
   55
   56int_hook(r, r2(_), _, [evaluate(false)]).
   57r2(A, Res, Flags) :-
   58    compound(A),
   59    compound_name_arguments(A, Name, Args1),
   60    maplist(interval__(Flags), Args1, Args2),
   61    compound_name_arguments(A1, Name, Args2),
   62    unwrap_r(A1, A2),
   63    !,
   64    eval(r(A2), Res1),
   65    clean(Res1, Res).
   66
   67r2(A, Res, Flags) :-
   68    interval_(A, Res, Flags).
   69
   70%
   71% Binomial distribution
   72%
   73int_hook(pbinom, pbinom_(atomic, atomic, atomic), atomic, []).
   74pbinom_(atomic(X), atomic(N), atomic(P),atomic(Res), _Flags) :-
   75    eval(r(pbinom(X, N, P)), Res).
   76
   77int_hook(pbinom, pbinom_(atomic, atomic, atomic, atomic), atomic, []).
   78pbinom_(atomic(X), atomic(N), atomic(P), atomic(Tail), atomic(Res), _Flags) :-
   79    eval(r(pbinom(X, N, P, Tail)), Res).
   80
   81int_hook(pbinom, pbinom(atomic, atomic, ...), ..., []).
   82pbinom(X, N, P, Res, Flags) :-
   83    pbinom(X, N, P, atomic(true), Res, Flags).
   84
   85int_hook(pbinom, pbinom(atomic, atomic, ..., atomic), ..., []).
   86
   87% lower tail
   88pbinom(X, N, P, atomic(true), Res, Flags) :-
   89    !,
   90    interval_(pbinom0(X, N, P), Res, Flags).
   91
   92% upper tail
   93pbinom(X, N, P, atomic(false), Res, Flags) :-
   94    interval_(pbinom1(X, N, P), Res, Flags).
   95
   96r_hook(pbinom0/3).
   97mono(pbinom0/3, [+, -, -]).
   98
   99r_hook(pbinom1/3).
  100mono(pbinom1/3, [-, +, +]).
  101
  102%
  103% Quantile function
  104%
  105int_hook(qbinom, qbinom_(atomic, atomic, atomic), atomic, []).
  106qbinom_(atomic(Alpha), atomic(N), atomic(P), atomic(Res), _Flags) :-
  107    eval(r(qbinom(Alpha, N, P)), Res).
  108
  109int_hook(qbinom, qbinom_(atomic, atomic, atomic, atomic), atomic, []).
  110qbinom_(atomic(Alpha), atomic(N), atomic(P), atomic(Tail), atomic(Res), _Flags) :-
  111    eval(r(qbinom(Alpha, N, P, Tail)), Res).
  112
  113int_hook(qbinom, qbinom2(..., _, ...), ..., []).
  114qbinom2(Alpha, N, P, Res, Flags) :-
  115    qbinom(Alpha, N, P, atomic(true), Res, Flags).
  116
  117int_hook(qbinom, qbinom(..., _, ..., atomic), ..., []).
  118
  119% lower tail
  120qbinom(Alpha, N, P, atomic(true), Res, Flags) :-
  121    !,
  122    interval_(qbinom0(Alpha, N, P), Res, Flags).
  123
  124% upper tail
  125qbinom(Alpha, N, P, atomic(false), Res, Flags) :-
  126    interval_(qbinom1(Alpha, N, P), Res, Flags).
  127
  128r_hook(qbinom0/3).
  129mono(qbinom0/3, [+, +, +]).
  130
  131r_hook(qbinom1/3).
  132mono(qbinom1/3, [-, +, +]).
  133
  134%
  135% Density
  136%
  137int_hook(dbinom, dbinom_(atomic, atomic, atomic), atomic, []).
  138dbinom_(atomic(X), atomic(N), atomic(P), atomic(Res), _Flags) :-
  139    eval(r(dbinom(X, N, P)), Res).
  140
  141int_hook(dbinom, dbinom3(..., atomic, ...), ..., []).
  142dbinom3(X1...X2, atomic(N), P1...P2, Res, Flags) :-
  143    dbinom(X1...X2, N...N, P1...P2, Res, Flags).
  144
  145int_hook(dbinom, dbinom4(..., atomic, atomic), ..., []).
  146dbinom4(X1...X2, atomic(N), atomic(P), Res, Flags) :-
  147    dbinom(X1...X2, N...N, P...P, Res, Flags).
  148
  149int_hook(dbinom, dbinom(..., ..., ...), ..., []).
  150% left to X / N
  151dbinom(X1...X2, N1...N2, P1...P2, Res, Flags) :-
  152    X2 < N1 * P1,
  153    !,
  154    interval_(dbinom0(X1...X2, N1...N2, P1...P2), Res, Flags).
  155
  156% right to X / N
  157dbinom(X1...X2, N1...N2, P1...P2, Res, Flags) :-
  158    X1 > N2 * P2,
  159    !,
  160    interval_(dbinom1(X1...X2, N1...N2, P1...P2), Res, Flags).
  161
  162% otherwise
  163dbinom(X1...X2, N1...N2, P1...P2, Res, _Flags) :-
  164    eval(r(dbinom2(X1, X2, N1, N2, P1, P2)), ##(L, U)),
  165    Res = L...U.
  166
  167r_hook(dbinom0/3).
  168mono(dbinom0/3, [+, -, -]).
  169
  170r_hook(dbinom1/3).
  171mono(dbinom1/3, [-, +, +]).
  172
  173%
  174% Normal distribution
  175%
  176r_hook(pnorm0/1).
  177mono(pnorm0/1, [+]).
  178
  179r_hook(pnorm1/1).
  180mono(pnorm1/1, [-]).
  181
  182% Atomic, Mu = 0, Sd = 1, lower tail
  183int_hook(pnorm, pnorm_(atomic), atomic, []).
  184pnorm_(atomic(A), atomic(Res), _Flags) :-
  185    eval(r(pnorm(A)), Res).
  186
  187% Atomic, lower tail
  188int_hook(pnorm, pnorm_(atomic, atomic, atomic), atomic, []).
  189pnorm_(atomic(A), atomic(Mu), atomic(Sigma), atomic(Res),_Flags) :-
  190    eval(r(pnorm(A, Mu, Sigma)), Res).
  191
  192% Atomic
  193int_hook(pnorm, pnorm_(atomic, atomic, atomic, atomic), atomic, []).
  194pnorm_(atomic(A), atomic(Mu), atomic(Sigma), atomic(Tail), atomic(Res), _Flags) :-
  195    eval(r(pnorm(A, Mu, Sigma, Tail)), Res).
  196
  197% Interval, Mu = 0, Sd = 1, lower tail
  198int_hook(pnorm, pnorm2(...), ..., []).
  199pnorm2(A, Res, Flags) :-
  200    pnorm6(A, atomic(true), Res, Flags).
  201
  202% Interval, lower tail
  203int_hook(pnorm, pnorm4(..., ..., ...), ..., []).
  204pnorm4(A, Mu, Sigma, Res, Flags) :-
  205    pnorm5(A, Mu, Sigma, atomic(true), Res, Flags).
  206
  207% Interval
  208int_hook(pnorm, pnorm5(..., ..., ..., atomic), ..., []).
  209pnorm5(A, Mu, Sigma, Tail, Res, Flags) :-
  210     interval_((A - Mu)/Sigma, Z, Flags),
  211     pnorm6(Z, Tail, Res, Flags).
  212
  213pnorm6(Z, atomic(true), Res, Flags) :-
  214     interval_(pnorm0(Z), Res, Flags).
  215
  216pnorm6(Z, atomic(false), Res, Flags) :-
  217     interval_(pnorm1(Z), Res, Flags).    
  218
  219%
  220% Quantile function
  221%
  222r_hook(qnorm0/1).
  223mono(qnorm0/1, [+]).
  224
  225r_hook(qnorm1/1).
  226mono(qnorm1/1, [-]).
  227
  228% Atomic, Mu = 0, Sd = 1, lower tail
  229int_hook(qnorm, qnorm_(atomic), atomic, []).
  230qnorm_(atomic(P), atomic(Res), _Flags) :-
  231    eval(r(qnorm(P)), Res).
  232
  233% Atomic, lower tail
  234int_hook(qnorm, qnorm_(atomic, atomic, atomic), atomic, []).
  235qnorm_(atomic(P), atomic(Mu), atomic(Sigma), atomic(Res), _Flags) :-
  236    eval(r(qnorm(P, Mu, Sigma)), Res).
  237
  238% Atomic
  239int_hook(qnorm, qnorm_(atomic, atomic, atomic, atomic), atomic, []).
  240qnorm_(atomic(P), atomic(Mu), atomic(Sigma), atomic(Tail), atomic(Res), _Flags) :-
  241    eval(r(qnorm(P, Mu, Sigma, Tail)), Res).
  242
  243% Interval, Mu = 0, Sd = 1, lower tail
  244int_hook(qnorm, qnorm2(...), ..., []).
  245qnorm2(P, Res, Flags) :-
  246     interval_(qnorm0(P), Res, Flags).
  247
  248% Interval, lower tail  
  249int_hook(qnorm, qnorm3(..., ..., ...), ..., []).
  250qnorm3(P, Mu, Sigma, Res, Flags) :-
  251    qnorm5(P, Mu, Sigma, atomic(true), Res, Flags).
  252
  253% Interval
  254int_hook(qnorm, qnorm4(..., ..., ..., atomic), ..., []).
  255qnorm4(P, Mu, Sigma, Tail, Res, Flags) :-
  256    qnorm5(P, Mu, Sigma, Tail, Res, Flags).
  257
  258qnorm5(P, Mu, Sigma, atomic(true), Res, Flags) :-
  259    interval_(qnorm0(P), Z, Flags),
  260    interval_(Mu + Z * Sigma, Res, Flags).
  261
  262qnorm5(P, Mu, Sigma, atomic(false), Res, Flags) :-
  263    interval_(qnorm1(P), Z, Flags),
  264    interval_(Mu + Z * Sigma, Res, Flags).
  265%
  266% Density
  267%
  268r_hook(dnorm0/1).
  269mono(dnorm0/1, [+]).
  270
  271r_hook(dnorm1/1).
  272mono(dnorm1/1, [-]).
  273
  274% Atomic, Mu = 0, Sd = 1
  275int_hook(dnorm, dnorm_(atomic), atomic, []).
  276dnorm_(atomic(A), atomic(Res), _Flags) :-
  277    eval(r(dnorm(A)), Res).
  278
  279% Atomic
  280int_hook(dnorm, dnorm_(atomic, atomic, atomic), atomic, []).
  281dnorm_(atomic(A), atomic(Mu), atomic(Sigma), atomic(Res), _Flags) :-
  282    eval(r(dnorm(A, Mu, Sigma)), Res).
  283
  284% Interval, Mu = 0, Sd = 1 
  285int_hook(dnorm, dnorm2(...), ..., []).
  286
  287% Interval 
  288int_hook(dnorm, dnorm3(..., ..., ...), ..., []).
  289dnorm3(X, Mu, Sigma, Res, Flags) :-
  290    interval_((X - Mu)/Sigma, Z, Flags),
  291    dnorm2(Z, Res0, Flags),
  292    interval_(atomic(1)/Sigma * Res0, Res, Flags).
  293
  294dnorm2(A...B, Res, Flags) :-
  295    B =< 0,
  296    !,
  297    interval_(dnorm0(A...B), Res, Flags).
  298
  299dnorm2(A...B, Res, Flags) :-
  300    A >= 0,
  301    !,
  302    interval_(dnorm1(A...B), Res, Flags).
  303
  304% mixed
  305dnorm2(A...B, Res, Flags) :-
  306    Max is max(abs(A), B),
  307    interval_(dnorm1(0...Max), Res, Flags).
  308
  309%
  310% t distribution
  311%
  312r_hook(pt0/2).
  313mono(pt0/2, [+,+]).
  314
  315r_hook(pt1/2).
  316mono(pt1/2, [-,-]).
  317
  318% Atomic, lower tail
  319int_hook(pt, pt_(atomic, atomic), atomic, []).
  320pt_(atomic(A), atomic(Df), atomic(Res), _Flags) :-
  321    eval(r(pt(A, Df)), Res).
  322
  323% Atomic
  324int_hook(pt, pt_(atomic, atomic, atomic), atomic, []).
  325
  326pt_(atomic(A), atomic(Df), atomic("lower"), atomic(Res), _Flags) :-
  327    eval(r(pt(A, Df, 'lower.tail'=true)), Res).
  328pt_(atomic(A), atomic(Df), atomic("upper"), atomic(Res), _Flags) :-
  329    eval(r(pt(A, Df, 'lower.tail'=false)), Res).
  330pt_(atomic(A), atomic(Df), atomic("two.sided"), atomic(Res), _Flags) :-
  331    eval(2 * r(pt(abs(A), Df, 'lower.tail'=false)), Res).
  332pt_(atomic(A), atomic(Df), atomic("density"), atomic(Res), _Flags) :-
  333    eval(r(dt(A, Df)), Res).
  334    
  335pt_(atomic(A), atomic(Df), atomic(Tail), atomic(Res), _Flags) :-
  336    eval(r(pt(A, Df, 'lower.tail'=Tail)), Res).
  337
  338% Interval, lower tail
  339int_hook(pt, pt2(..., _), ..., []).
  340pt2(A, Df, Res, Flags) :-
  341    pt(A, Df, atomic(true), Res, Flags).
  342
  343% Interval
  344int_hook(pt, pt(..., _, atomic), ..., []).
  345
  346pt(A, Df, atomic("lower"), Res, Flags) :-
  347    pt(A, Df, atomic(true), Res, Flags).
  348pt(A, Df, atomic("upper"), Res, Flags) :-
  349    pt(A, Df, atomic(false), Res, Flags).
  350pt(A, Df, atomic("two.sided"), Res, Flags) :-
  351    interval_(atomic(2) * pt(abs(A), Df, atomic("upper")), Res, Flags).
  352pt(A, Df, atomic("density"), Res, Flags) :-
  353    interval_(dt(A, Df), Res, Flags).
  354
  355% lower tail
  356pt(L...U, Df, atomic(true), Res, Flags) :-
  357    !,
  358    interval_(pt0(L...U, Df), Res, Flags).
  359
  360% upper tail
  361pt(L...U, Df, atomic(false), Res, Flags) :-
  362    !, 
  363    interval_(pt1(L...U, Df), Res, Flags).
  364
  365%
  366% Quantile function
  367%
  368r_hook(qt0/2).
  369mono(qt0/2, [+,-]).
  370
  371r_hook(qt1/2).
  372mono(qt1/2, [-,+]).
  373
  374% Atomic, lower tail
  375int_hook(qt, qt_(atomic, atomic), atomic, []).
  376qt_(atomic(P), atomic(Df), atomic(Res), _Flags) :-
  377    eval(r(qt(P, Df)), Res).
  378
  379% Atomic
  380int_hook(qt, qt_(atomic, atomic, atomic), atomic, []).
  381qt_(atomic(P), atomic(Df), atomic(Tail), atomic(Res), _Flags) :-
  382    eval(r(qt(P, Df, 'lower.tail'=Tail)), Res).
  383
  384% Interval, lower tail
  385int_hook(qt, qt2(..., _), ..., []).
  386qt2(P, Df, Res, Flags) :-
  387    qt(P, Df, atomic(true), Res, Flags).
  388
  389% Interval
  390int_hook(qt, qt(..., _, atomic), ..., []).
  391
  392% lower tail
  393qt(P, Df, atomic(true), Res, Flags) :-
  394    interval_(qt0(P, Df), Res, Flags).
  395
  396% upper tail
  397qt(P, Df, atomic(false), Res, Flags) :-
  398    interval_(qt1(P, Df), Res, Flags).
  399
  400%
  401% Density
  402%
  403r_hook(dt0/2).
  404mono(dt0/2, [+, +]).
  405
  406r_hook(dt1/2).
  407mono(dt1/2, [-, +]).
  408
  409int_hook(dt, dt_(atomic, atomic), atomic, []).
  410dt_(atomic(A), atomic(Df), atomic(Res), _Flags) :-
  411    eval(r(dt(A, Df)), Res).
  412
  413int_hook(dt, dt(..., _), ..., []).
  414dt(L...U, Df, Res, Flags) :-
  415    U =< 0,
  416    !,
  417    interval_(dt0(L...U, Df), Res, Flags).
  418
  419dt(L...U, Df, Res, Flags) :-
  420    L >= 0,
  421    !,
  422    interval_(dt1(L...U, Df), Res, Flags).
  423
  424% mixed
  425dt(L...U, Df, Res, Flags) :-
  426    Max is max(abs(L), U),
  427    interval_(dt1(0...Max, Df), Res, Flags). 
  428
  429%
  430% chisq
  431%
  432r_hook(pchisq0/2).
  433mono(pchisq0/2, [+,-]).
  434
  435r_hook(pchisq1/2).
  436mono(pchisq1/2, [-,+]).
  437
  438% Atomic, lower tail
  439int_hook(pchisq, pchisq_(atomic, atomic), atomic, []).
  440pchisq_(atomic(A), atomic(Df), atomic(Res), _Flags) :-
  441    eval(r(pchisq(A, Df)), Res).
  442
  443% Atomic
  444int_hook(pchisq, pchisq_(atomic, atomic, atomic), atomic, []).
  445pchisq_(atomic(A), atomic(Df), atomic(Tail), atomic(Res), _Flags) :-
  446    eval(r(pchisq(A, Df, 'lower.tail'=Tail)), Res).
  447
  448% Interval, lower tail
  449int_hook(pchisq, pchisq2(..., atomic), ..., []).
  450pchisq2(A, Df, Res, Flags) :-
  451    pchisq(A, Df, atomic(true), Res, Flags).
  452
  453% Interval
  454int_hook(pchisq, pchisq(..., atomic, atomic), ..., []).
  455
  456% lower tail
  457pchisq(L...U, Df, atomic(true), Res, Flags):-
  458    !,
  459    interval_(pchisq0(L...U, Df), Res, Flags).
  460
  461% upper tail
  462pchisq(L...U, Df, atomic(false), Res, Flags):-
  463    !,
  464    interval_(pchisq1(L...U, Df), Res, Flags).
  465
  466%
  467% quantile function
  468%
  469r_hook(qchisq0/2).
  470mono(qchisq0/2, [+,+]).
  471
  472r_hook(qchisq1/2).
  473mono(qchisq1/2, [-,+]).
  474
  475% Atomic, lower tail
  476int_hook(qchisq, qchisq_(atomic, atomic), atomic, []).
  477qchisq_(atomic(P), atomic(Df), atomic(Res), _Flags) :-
  478    eval(r(qchisq(P, Df)), Res).
  479
  480% Atomic
  481int_hook(qchisq, qchisq_(atomic, atomic, atomic), atomic, []).
  482qchisq_(atomic(P), atomic(Df), atomic(Tail), atomic(Res), _Flags) :-
  483    eval(r(qchisq(P, Df, 'lower.tail'=Tail)), Res).
  484
  485% Interval, lower tail
  486int_hook(qchisq, qchisq2(..., atomic), ..., []).
  487qchisq2(P, Df, Res, Flags) :-
  488    qchisq(P, Df, atomic(true), Res, Flags).
  489
  490% Interval
  491int_hook(qchisq, qchisq(..., atomic, atomic), ..., []).
  492
  493% lower tail
  494qchisq(L...U, Df, atomic(true), Res, Flags):-
  495    !,
  496    interval_(qchisq0(L...U, Df), Res, Flags).
  497
  498% upper tail
  499qchisq(L...U, Df, atomic(false), Res, Flags):-
  500    interval_(qchisq1(L...U, Df), Res, Flags).
  501
  502%
  503% density
  504%
  505r_hook(dchisq0/2).
  506mono(dchisq0/2, [-,/]).
  507
  508r_hook(dchisq1/2).
  509mono(dchisq1/2, [+,/]).
  510
  511int_hook(dchisq, dchisq_(atomic, atomic), atomic, []).
  512dchisq_(atomic(A), atomic(Df), atomic(Res), _Flags) :-
  513    eval(r(dchisq(A, Df)), Res).
  514
  515int_hook(dchisq, dchisq(..., atomic), ..., []).
  516% for df<=2
  517dchisq(L...U, atomic(Df), Res, Flags):-
  518    Df =< 2,
  519    !,
  520    interval_(dchisq0(L...U, atomic(Df)), Res, Flags).
  521
  522% for df>2
  523dchisq(L...U, atomic(Df), Res, Flags):-
  524    dchisq_A(L...U, atomic(Df), Res, Flags).
  525
  526% for x < mode
  527dchisq_A(L...U, atomic(Df), Res, Flags) :-
  528    Mode is Df - 2,
  529    U =< Mode,
  530    !,
  531    interval_(dchisq1(L...U, atomic(Df)), Res, Flags).
  532
  533% for x > mode
  534dchisq_A(L...U, atomic(Df), Res, Flags) :-
  535    Mode is Df - 2,
  536    L >= Mode,
  537    !,
  538    interval_(dchisq0(L...U, atomic(Df)), Res, Flags).
  539
  540% for L < mode, U > mode
  541dchisq_A(L...U, atomic(Df), Res, Flags) :-
  542    interval_(dchisq(atomic(L), atomic(Df)), atomic(X1), Flags),
  543    interval_(dchisq(atomic(U), atomic(Df)), atomic(X2), Flags),
  544    L1 is min(X1, X2),
  545    Mode is Df - 2,
  546    interval_(dchisq(atomic(Mode), atomic(Df)), atomic(U1), Flags),
  547    Res = L1...U1.
  548
  549%
  550% Assignment
  551%
  552r_hook('<-'/2).
  553int_hook('<-', assign0(_, _), _, [evaluate(false)]).
  554assign0(Var, A, Res, Flags) :-
  555    interval_(A, A1, Flags),
  556    interval2_(assign(Var, A1), Res, Flags).
  557
  558int_hook(assign, assign1(atomic, atomic), atomic, []).
  559assign1(atomic(Var), atomic(A), Res, _Flags) :-
  560    eval(Var <- A, Res1),
  561    clean(Res1, Res).
  562
  563int_hook(assign, assign2(atomic, ...), ..., []).
  564assign2(atomic(Var), L...U, Res, _Flags) :-
  565    eval(Var