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
   47r_hook(true).
   48r_hook(false).
   49
   50%
   51% Call R 
   52%
   53int_hook(r, r1(atomic), _, [evaluate(false)]).
   54r1(atomic(A), Res, _Flags) :-
   55    eval_hook(r(A), Res1),
   56    !,
   57    clean(Res1, Res).
   58
   59int_hook(r, r2(_), _, [evaluate(false)]).
   60r2(A, Res, Flags) :-
   61    compound(A),
   62    compound_name_arguments(A, Name, Args1),
   63    maplist(interval__(Flags), Args1, Args2),
   64    compound_name_arguments(A1, Name, Args2),
   65    unwrap_r(A1, A2),
   66    !,
   67    eval_hook(r(A2), Res1),
   68    clean(Res1, Res).
   69
   70r2(A, Res, Flags) :-
   71    interval_(A, Res, Flags).
   72
   73%
   74% Binomial distribution
   75%
   76int_hook(pbinom, pbinom(atomic, atomic, ..., atomic), ..., []).
   77
   78% lower tail
   79pbinom(atomic(X), atomic(N), P, atomic(true), Res, Flags) :-
   80    !,
   81    interval_(pbinom0(atomic(X), atomic(N), P), Res, Flags).
   82
   83% upper tail
   84pbinom(atomic(X), atomic(N), P, atomic(false), Res, Flags) :-
   85    interval_(pbinom1(atomic(X), atomic(N), P), Res, Flags).
   86
   87r_hook(pbinom0/3).
   88mono(pbinom0/3, [+, -, -]).
   89
   90r_hook(pbinom1/3).
   91mono(pbinom1/3, [-, +, +]).
   92
   93%
   94% Quantile function - hier weiter
   95%
   96int_hook(qbinom, qbinom(..., ..., ..., atomic), ..., []).
   97
   98% lower tail
   99qbinom(Alpha, N, P, atomic(true), Res, Flags) :-
  100    !,
  101    interval_(qbinom0(Alpha, N, P), Res, Flags).
  102
  103% upper tail
  104qbinom(Alpha, N, P, atomic(false), Res, Flags) :-
  105    interval_(qbinom1(Alpha, N, P), Res, Flags).
  106
  107r_hook(qbinom0/3).
  108mono(qbinom0/3, [+, +, +]).
  109
  110r_hook(qbinom1/3).
  111mono(qbinom1/3, [-, +, +]).
  112
  113%
  114% Density
  115%
  116int_hook(dbinom, dbinom(..., ..., ...), ..., []).
  117
  118% left to X / N
  119dbinom(X1...X2, N1...N2, P1...P2, Res, Flags) :-
  120    X2 < N1 * P1,
  121    !,
  122    interval_(dbinom0(X1...X2, N1...N2, P1...P2), Res, Flags).
  123
  124% right to X / N
  125dbinom(X1...X2, N1...N2, P1...P2, Res, Flags) :-
  126    X1 > N2 * P2,
  127    !,
  128    interval_(dbinom1(X1...X2, N1...N2, P1...P2), Res, Flags).
  129
  130% otherwise
  131dbinom(X1...X2, N1...N2, P1...P2, Res, _Flags) :-
  132    r(dbinom2(X1, X2, N1, N2, P1, P2), ##(L, U)),
  133    Res = L...U.
  134
  135r_hook(dbinom0/3).
  136mono(dbinom0/3, [+, -, -]).
  137
  138r_hook(dbinom1/3).
  139mono(dbinom1/3, [-, +, +]).
  140
  141%
  142% Normal distribution
  143%
  144r_hook(pnorm0/1).
  145mono(pnorm0/1, [+]).
  146
  147int_hook(pnorm, pnorm(..., ..., ...), ..., []).
  148pnorm(X, Mu, Sigma, Res, Flags) :-
  149     interval_((X - Mu)/Sigma, Z, Flags),
  150     interval_(pnorm0(Z), Res, Flags).
  151
  152int_hook(pnorm, pnorm1(...), ..., []).
  153pnorm1(Z, Res, Flags) :-
  154     interval_(pnorm0(Z), Res, Flags).
  155
  156int_hook(pnorm, pnorm2(atomic), atomic, []).
  157pnorm2(atomic(Z), atomic(Res), Flags) :-
  158     eval(pnorm0(Z), Res, Flags).
  159
  160%
  161% Quantile function
  162%
  163r_hook(qnorm0/1).
  164mono(qnorm0/1, [+]).
  165
  166int_hook(qnorm, qnorm(..., ..., ...), ..., []).
  167qnorm(P, Mu, Sigma, Res, Flags) :-
  168     interval_(qnorm0(P), Z, Flags),
  169     interval_(Mu + Z * Sigma, Res, Flags).
  170
  171int_hook(qnorm, qnorm1(...), ..., []).
  172qnorm1(P, Res, Flags) :-
  173     interval_(qnorm0(P), Res, Flags).
  174
  175int_hook(qnorm, qnorm2(atomic), atomic, []).
  176qnorm2(atomic(P), atomic(Res), Flags) :-
  177     eval(qnorm0(P), Res, Flags).
  178
  179%
  180% Density
  181%
  182r_hook(dnorm1/1).
  183mono(dnorm1/1, [+]).
  184
  185r_hook(dnorm2/1).
  186mono(dnorm2/1, [-]).
  187
  188int_hook(dnorm, dnorm(..., ..., ...), ..., []).
  189dnorm(X, Mu, Sigma, Res, Flags) :-
  190    interval_((X - Mu)/Sigma, Z, Flags),
  191    interval_(atomic(1)/Sigma * dnorm0(Z), Res, Flags).
  192
  193int_hook(dnorm0, dnorm0(...), ..., []).
  194dnorm0(A...B, Res, Flags) :-
  195    B =< 0,
  196    !,
  197    interval_(dnorm1(A...B), Res, Flags).
  198
  199dnorm0(A...B, Res, Flags) :-
  200    A >= 0,
  201    !,
  202    interval_(dnorm2(A...B), Res, Flags).
  203
  204% mixed
  205dnorm0(A...B, Res, Flags) :-
  206    Max is max(abs(A), B),
  207    interval_(dnorm2(0...Max), Res, Flags).
  208
  209%
  210% t distribution
  211%
  212int_hook(pt, pt(..., ..., atomic), ..., []).
  213
  214r_hook(pt0/2).
  215mono(pt0/2, [+,-]).
  216
  217r_hook(pt1/2).
  218mono(pt1/2, [+,+]).
  219
  220r_hook(pt2/2).
  221mono(pt2/2, [-,+]).
  222
  223r_hook(pt3/2).
  224mono(pt3/2, [-,-]).
  225
  226% lower tail
  227pt(L...U, Df, atomic(true), Res, Flags) :-
  228    U =< 0,
  229    !,
  230    interval_(pt0(L...U, Df), Res, Flags).
  231
  232pt(L...U, Df, atomic(true), Res, Flags) :-
  233    L >= 0,
  234    !,
  235    interval_(pt1(L...U, Df), Res, Flags).
  236
  237pt(L...U, Df, atomic(true), Res, Flags) :-
  238    Max is max(abs(L), U), 
  239    interval_(pt1(0...Max, Df), Res, Flags).
  240
  241% upper tail
  242pt(L...U, Df, atomic(false), Res, Flags) :-
  243    U =< 0,
  244    !, 
  245    interval_(pt2(L...U, Df), Res, Flags).
  246
  247pt(L...U, Df, atomic(false), Res, Flags) :-
  248    L >= 0,
  249    !, 
  250    interval_(pt3(L...U, Df), Res, Flags).
  251
  252pt(L...U, Df, atomic(false), Res, Flags) :-
  253    Max is max(abs(L), U), 
  254    interval_(pt3(0...Max, Df), Res, Flags).
  255
  256%
  257% Quantile function
  258%
  259r_hook(qt0/2).
  260mono(qt0/2, [+,-]).
  261
  262int_hook(qt, qt(..., ...), ..., []).
  263qt(P, Df, Res, Flags) :-
  264    interval_(qt0(P, Df), Res, Flags).
  265
  266%
  267% Density
  268%
  269r_hook(dt0/2).
  270mono(dt0/2, [+, +]).
  271
  272r_hook(dt1/2).
  273mono(dt1/2, [-, +]).
  274
  275int_hook(dt, dt(..., ...), ..., []).
  276dt(L...U, Df, Res, Flags) :-
  277    U =< 0,
  278    !,
  279    interval_(dt0(L...U, Df), Res, Flags).
  280
  281dt(L...U, Df, Res, Flags) :-
  282    L >= 0,
  283    !,
  284    interval_(dt1(L...U, Df), Res, Flags).
  285
  286% mixed
  287dt(L...U, Df, Res, Flags) :-
  288    Max is max(abs(L), U),
  289    interval_(dt1(0...Max, Df), Res, Flags). 
  290
  291%
  292% chisq
  293%
  294
  295int_hook(pchisq, pchisq(..., atomic, atomic), ..., []).
  296
  297r_hook(pchisq0/2).
  298mono(pchisq0/2, [+,-]).
  299
  300r_hook(pchisq1/2).
  301mono(pchisq1/2, [-,+]).
  302
  303% lower tail
  304pchisq(L...U, Df, atomic(true), Res, Flags):-
  305    !,
  306    interval_(pchisq0(L...U, Df), Res, Flags).
  307
  308% upper tail
  309pchisq(L...U, Df, atomic(false), Res, Flags):-
  310    !,
  311    interval_(pchisq1(L...U, Df), Res, Flags).
  312
  313%
  314% quantile function
  315%
  316int_hook(qchisq, qchisq(..., atomic, atomic), ..., []).
  317
  318r_hook(qchisq0/2).
  319mono(qchisq0/2, [+,+]).
  320
  321r_hook(qchisq1/2).
  322mono(qchisq1/2, [-,+]).
  323
  324qchisq(L...U, Df, atomic(true), Res, Flags):-
  325    !,
  326    interval_(qchisq0(L...U, Df), Res, Flags).
  327
  328qchisq(L...U, Df, atomic(false), Res, Flags):-
  329    interval_(qchisq1(L...U, Df), Res, Flags).
  330
  331%
  332% density
  333%
  334int_hook(dchisq, dchisq(..., atomic), ..., []).
  335
  336r_hook(dchisq0/2).
  337mono(dchisq0/2, [-,/]).
  338
  339r_hook(dchisq1/2).
  340mono(dchisq1/2, [+,/]).
  341
  342% for df<=2
  343dchisq(L...U, atomic(Df), Res, Flags):-
  344    Df =< 2,
  345    !,
  346    interval_(dchisq0(L...U, atomic(Df)), Res, Flags).
  347
  348% for df>2
  349dchisq(L...U, atomic(Df), Res, Flags):-
  350    dchisq_A(L...U, atomic(Df), Res, Flags).
  351
  352% for x < mode
  353dchisq_A(L...U, atomic(Df), Res, Flags) :-
  354    Mode is Df - 2,
  355    U =< Mode,
  356    !,
  357    interval_(dchisq1(L...U, atomic(Df)), Res, Flags).
  358
  359% for x > mode
  360dchisq_A(L...U, atomic(Df), Res, Flags) :-
  361    Mode is Df - 2,
  362    L >= Mode,
  363    !,
  364    interval_(dchisq0(L...U, atomic(Df)), Res, Flags).
  365
  366% for L < mode, U > mode
  367dchisq_A(L...U, atomic(Df), Res, Flags) :-
  368    interval_(dchisq(atomic(L), atomic(Df)), X1..._, Flags),
  369    interval_(dchisq(atomic(U), atomic(Df)), X3..._, Flags),
  370    L1 is min(X1, X3),
  371    Mode is Df - 2,
  372    interval_(dchisq(atomic(Mode), atomic(Df)), U1..._, Flags),
  373    Res = L1...U1.
  374
  375%
  376% Assignment
  377%
  378r_hook('<-'/2).
  379int_hook('<-', assign(_, _), _, [evaluate(false)]).
  380assign(Var, A, Res, Flags) :-
  381    interval_(A, A1, Flags),
  382    assign_(Var, A1, Res).
  383
  384assign_(atomic(Var), L...U, Res) :-
  385    eval_hook(Var <- call("...", L, U), Res),
  386    !.
  387
  388assign_(atomic(Var), atomic(A), Res) :-
  389    eval_hook(Var