3:- module(cplint_util,[
    4  bar/2,
    5  bar/3,
    6  bar1/2,
    7  argbar/2,
    8  histogram/2,
    9  histogram/3,
   10  densities/3,
   11  densities/4,
   12  density/2,
   13  density/3,
   14  density2d/2,
   15  density2d/3,
   16  to_pair/2,
   17  key_pair/2,
   18  value_pair/2,
   19  bin/5,
   20  beta/2,
   21  to_atom/2,
   22  average/2,
   23  variance/2,
   24  variance/3,
   25  std_dev/2,
   26  std_dev/3,
   27  agg_val/3,
   28  swi_builtin/1]).   29
   30
   31:- use_module(library(lists)).   32:- use_module(library(apply)).   33:- use_module(library(clpr)).   34:- use_module(library(matrix)).   35:- use_module(library(clpfd)).   36
   37:- use_module(library(matrix)).   38:- use_module(highlight).

cplint_util

Utility module for cplint

author
- Fabrizio Riguzzi
copyright
- Fabrizio Riguzzi */
license
- Artistic License 2.0 https://opensource.org/licenses/Artistic-2.0



 bar(+Probability:float, -Chart:dict) is det
The predicate returns a dict for rendering with c3 as a bar chart with
a bar for the probability and a bar for one minus the probability.
/


   53bar(P,Chart):-
   54  must_be(float,P),
   55  PF is 1.0-P,
   56  Chart = c3{data:_{x:elem, rows:[elem-prob,'T'-P,'F' -PF], type:bar},
   57          axis:_{x:_{type:category}, rotated: true,
   58                 y:_{min:0.0,max:1.0,padding:_{bottom:0.0,top:0.0},
   59             tick:_{values:[0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]}}},
   60	           size:_{height: 100},
   61	          legend:_{show: false}}.



 bar1(+Probability:float, -Chart:dict) is det
The predicate returns a dict for rendering with c3 as a bar chart with
a bar for the probability
/


   69bar1(P,Chart):-
   70  must_be(float,P),
   71  Chart = c3{data:_{x:elem, rows:[elem-prob,'T'-P], type:bar},
   72          axis:_{x:_{type:category}, rotated: true,
   73                 y:_{min:0.0,max:1.0,padding:_{bottom:0.0,top:0.0},
   74             tick:_{values:[0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]}}},
   75	           size:_{height: 100},
   76	          legend:_{show: false}}.



 bar(+Successes:int, +Failures:int, -Chart:dict) is det
The predicate returns a dict for rendering with c3 as a bar chart with
a bar for the number of successes and a bar for the number of failures
/


   83bar(S,F,Chart):-
   84  must_be(nonneg,S),
   85  must_be(nonneg,F),
   86  Chart = c3{data:_{x:elem, rows:[elem-prob,'T'-S,'F' -F], type:bar},
   87          axis:_{x:_{type:category}, rotated: true,
   88                 y:_{min:0.0,padding:_{bottom:0.0}}},
   89	           size:_{height: 100},
   90	          legend:_{show: false}}.



 argbar(+Values:list, -Chart:dict) is det
Values is a list of pairs V-N where
V is the value and N is the number of samples
returning that value.
The predicate returns a dict for rendering with c3 as a bar chart with
a bar for each value V.
The size of the bar is given by N.
/


  101argbar(ValList,Chart):-
  102  must_be(list,ValList),
  103  maplist(to_atom,ValList,ValList1),
  104  Chart = c3{data:_{x:elem, rows:[elem-prob|ValList1], type:bar},
  105          axis:_{x:_{type:category}, rotated: true,
  106                 y:_{min:0.0,padding:_{bottom:0.0}}},
  107             %  size:_{height: 100},
  108              legend:_{show: false}}.



 to_atom(+In:pair, -Out:pair) is det
Given In=A0-N, to_atom/2 returns Out=A-N
where A is an atom representing A0
/


  117to_atom(A0-N,A-N):-
  118  term_to_atom(A0,A).



 histogram(+List:list, -Chart:dict) is det
Equivalent to histogram/3 with an empty option list.
/


  125histogram(L0,Chart):-
  126  histogram(L0,Chart,[]).



 histogram(+List:list, -Chart:dict, +Options:list) is det
Draws a histogram of the samples in List. List must be a list of pairs of the form [V]-W or V-W
where V is a sampled value and W is its weight, or a list of values.



Options is a list of options, the following are recognised by histogram/3:





min(+Min:float)
the minimum value of domain, default value the minimum in List


max(+Max:float)
the maximum value of domain, default value the maximum in List


nbins(+NBins:int)
the number of bins for dividing the domain, default value 40
/







  142histogram(L0,Chart,Options):-
  143  must_be(list,L0),
  144  must_be(list,Options),
  145  (L0=[_-_|_] ->
  146    maplist(to_pair,L0,L1)
  147  ;
  148    maplist([X,Y]>>(Y=(X-1)),L0,L1)
  149  ),
  150  maplist(key_pair,L1,L2),
  151  max_list(L2,DMax0),
  152  min_list(L2,DMin),
  153  DMax is DMax0+(DMax0-DMin)/12,
  154  option(max(Max),Options,DMax),
  155  option(min(Min),Options,DMin),
  156  option(nbins(NBins),Options,40),
  157  histogram(L1,NBins,Min,Max,Chart).



 histogram(+List:list, +NBins:int, +Min:float, +Max:float, -Chart:dict) is det
Draws a histogram of the samples in List dividing the domain in
NBins bins. List must be a list of pairs of the form V-W
where V is a sampled value and W is its weight. The minimum and maximum
values of the domains must be provided.
/


  167histogram(L0,NBins,Min,Max,Chart):-
  168  keysort(L0,L),
  169  D is Max-Min,
  170  BinWidth is D/NBins,
  171  bin(NBins,L,Min,BinWidth,LB),
  172  maplist(key_pair,LB,X),
  173  maplist(value_pair,LB,Y),
  174  Chart = c3{data:_{x:x,
  175    columns:[[x|X],[freq|Y]], type:bar},
  176    axis:_{ x:_{ tick:_{fit:false}}},
  177     bar:_{
  178     width:_{ ratio: 1.0 }},
  179     legend:_{show: false}}.



 densities(+PriorList:list, +PostList:list, -Chart:dict) is det
Equivalent to densities/4 with an empty option list.
/


  186densities(Pri,Post,Chart):-
  187  densities(Pri,Post,Chart,[]).



 densities(+PriorList:list, +PostList:list, -Chart:dict, +Options:list) is det
Draws a line chart of the density of two sets of samples, usually
prior and post observations. The samples from the prior are in PriorList
while the samples from the posterior are in PostList.
PriorList and PostList must be lists of pairs of the form [V]-W or V-W
where V is a sampled value and W is its weight, or lists of values V.
Options is a list of options, the following are recognised by histogram/3:




nbins(+NBins:int)
the number of bins for dividing the domain, default value 40
*/







  201densities(Pri0,Post0,Chart,Options):-
  202  must_be(list,Pri0),
  203  must_be(list,Post0),
  204  must_be(list,Options),
  205  (Pri0=[_-_|_] ->
  206    maplist(to_pair,Pri0,Pri1)
  207  ;
  208    maplist([X,Y]>>(Y=(X-1)),Pri0,Pri1)
  209  ),
  210  (Post0=[_-_|_] ->
  211    maplist(to_pair,Post0,Post1)
  212  ;
  213    maplist([X,Y]>>(Y=(X-1)),Post0,Post1)
  214  ),
  215  option(nbins(NBins),Options,40),
  216  maplist(key_pair,Pri1,Pri),
  217  maplist(key_pair,Post1,Post),
  218  append(Pri,Post,All),
  219  max_list(All,DMax0),
  220  min_list(All,DMin),
  221  DMax is DMax0+(DMax0-DMin)/12,
  222  option(max(Max),Options,DMax),
  223  option(min(Min),Options,DMin),
  224  D is Max-Min,
  225  BinWidth is D/NBins,
  226  keysort(Pri1,Pr),
  227  keysort(Post1,Po),
  228  bin(NBins,Pr,Min,BinWidth,LPr),
  229  bin(NBins,Po,Min,BinWidth,LPo),
  230  maplist(key_pair,LPr,X),
  231  maplist(value_pair,LPr,YPr),
  232  maplist(value_pair,LPo,YPo),
  233  Chart = c3{data:_{x: x,
  234  columns: [[x|X],
  235    [pre|YPr],[post|YPo]]},
  236   axis:_{ x:_{ tick:_{fit:false}}}
  237  }.



 density(+List:list, +NBins:int, +Min:float, +Max:float, -Chart:dict) is det
Draws a line chart of the density of a sets of samples.
The samples are in List
as pairs [V]-W or V-W where V is a value and W its weigth.
The lines are drawn dividing the domain in
NBins bins. The X axis goes from Min to Max.
/


  247density(Post,NBins,Min,Max,Chart):-
  248  D is Max-Min,
  249  BinWidth is D/NBins,
  250  keysort(Post,Po),
  251  bin(NBins,Po,Min,BinWidth,LPo),
  252  maplist(key_pair,LPo,X),
  253  maplist(value_pair,LPo,YPo),
  254  Chart = c3{data:_{x: x,
  255  columns: [[x|X],
  256    [dens|YPo]]},
  257   axis:_{ x:_{ tick:_{fit:false}}}
  258  }.



 density2d(+List:list, +NBins:int, +XMin:float, +XMax:float, +YMin:float, +YMax:float, -Dist:list) is det
Returns the density of a sets of two dimensional samples.
The samples are in List
as pairs [V]-W or V-W where V is a value and W its weigth.
The lines are drawn dividing the domain in
NBins bins. The X axis goes from Min to Max.
/


  269density2d(Post0,NBins,XMin,XMax,YMin,YMax,D):-
  270  maplist(to_pair,Post0,Post),
  271  DX is XMax-XMin,
  272  XBinWidth is DX/NBins,
  273  DY is YMax-YMin,
  274  YBinWidth is DY/NBins,
  275  bin2D(NBins,Post,XMin,YMin,XBinWidth,YBinWidth,D).



 density(+List:list, -Chart:dict) is det
Equivalent to density/3 with an empty option list.
/


  282density(Post,Chart):-
  283  density(Post,Chart,[]).



 density(+List:list, -Chart:dict, +Options:list) is det
Draws a line chart of the density of a sets of samples.
The samples are in List
as pairs [V]-W or V-W where V is a value and W its weigth.



Options is a list of options, the following are recognised by density/3:





min(+Min:float)
the minimum value of domain, default value the minimum in List


max(+Max:float)
the maximum value of domain, default value the maximum in List


nbins(+NBins:int)
the number of bins for dividing the domain, default value 40
/







  300density(Post,Chart,Options):-
  301  must_be(list,Post),
  302  must_be(list,Options),
  303  (Post=[_-_|_] ->
  304    maplist(to_pair,Post,Post0)
  305  ;
  306    maplist([X,Y]>>(Y=(X-1)),Post,Post0)
  307  ),
  308  maplist(key_pair,Post0,PoK),
  309  max_list(PoK,DMax0),
  310  min_list(PoK,DMin),
  311  DMax is DMax0+(DMax0-DMin)/12,
  312  option(max(Max),Options,DMax),
  313  option(min(Min),Options,DMin),
  314  option(nbins(NBins),Options,40),
  315  density(Post0,NBins,Min,Max,Chart).



 density2d(+List:list, -Dens:list) is det
Equivalent to density2d/3 with an empty option list.
/


  322density2d(Post0,D):-
  323  density2d(Post0,D,[]).



 density2d(+List:list, -Dens:list, +Options:list) is det
Returns a set of 3-dimensional points representing the plot of the
density of a sets of 2-dimensional samples.
The samples are in List
as pairs [X,Y]-W where (X,Y) is a point and W its weigth.



Options is a list of options, the following are recognised by density2d/3:





xmin(+XMin:float)
the minimum value of the X domain, default value the minimum in List


xmax(-XMax:float)
the maximum value of the X domain, default value the maximum in List


ymin(-YMin:float)
the minimum value of the Y domain, default value the minimum in List


ymax(-YMax:float)
the maximum value of the Y domain, default value the maximum in List


nbins(+NBins:int)
the number of bins for dividing the X and Y domains, default value 40
/







  345density2d(Post0,D,Options):-
  346  must_be(list,Post0),
  347  must_be(list,Options),
  348  maplist(key_x_y,Post0,X,Y),
  349  max_list(X,DxMax),
  350  min_list(X,DxMin),
  351  max_list(Y,DyMax),
  352  min_list(Y,DyMin),
  353  option(xmax(XMax),Options,DxMax),
  354  option(xmin(XMin),Options,DxMin),
  355  option(ymax(YMax),Options,DyMax),
  356  option(ymin(YMin),Options,DyMin),
  357  option(nbins(NBins),Options,40),
  358  density2d(Post0,NBins,XMin,XMax,YMin,YMax,D).



 to_pair(+Pair:pair, -FlattenedPair:pair) is det
Given a pair E-W, returns a pair Ep-W where
Ep=EE if E=[EE], otherwise Ep=E
/


  366to_pair([E]-W,E-W):- !.
  367to_pair(E-W,E-W).



 key_pair(+Pair:pair, -Key:term) is det
Given a pair Key-Vaule, returns its first element Key
/


  374key_pair(K-_,K).



 value_pair(+Pair:pair, -Value:term) is det
Given a pair Key-Vaule, returns its second element Value
/


  381value_pair(_ - Y,Y).
  382
  383key_x_y([X,Y]-_,X,Y).
  384
  385
  386bin2D(NBins,Post,XMin,YMin,XBinWidth,YBinWidth,D):-
  387  binX(NBins,NBins,Post,XMin,YMin,XBinWidth,YBinWidth,D).
  388
  389binX(0,_NBins,_Post,_XMin,_YMin,_XBinWidth,_YBinWidth,[]):-!.
  390
  391binX(N,NBins,L,XLower,YMin,XBW,YBW,[R|D]):-
  392  V is XLower+XBW/2,
  393  XUpper is XLower+XBW,
  394  binY(NBins,L,V,XLower,XUpper,YMin,YBW,R),
  395  N1 is N-1,
  396  binX(N1,NBins,L,XUpper,YMin,XBW,YBW,D).
  397
  398binY(0,_Post,_XV,_XMin,_YMin,_XBinWidth,_YBinWidth,[]):-!.
  399
  400binY(N,L,XV,XLower,XUpper,YLower,YBW,[[XV,YV]-Freq|D]):-
  401    YV is YLower+YBW/2,
  402    YUpper is YLower+YBW,
  403    count_bin2d(L,XLower,XUpper,YLower,YUpper,0,Freq),
  404    N1 is N-1,
  405    binY(N1,L,XV,XLower,XUpper,YUpper,YBW,D).
  406
  407count_bin2d([],_XL,_XU,_YL,_YU,F,F).
  408
  409count_bin2d([[X,Y]-W|T0],XL,XU,YL,YU,F0,F):-
  410  ((X>=XL,X<XU,Y>=YL,Y<YU)->
  411    F1 is F0+W
  412  ;
  413    F1 = F0
  414  ),
  415  count_bin2d(T0,XL,XU,YL,YU,F1,F).



 bin(+N:int, +Values:list, +Lower:number, +BinWidth:number, -Couples:list) is det
Given a list of numeric Values, a Lower value and BinWidth, returns in Couples
a list of N pairs V-Freq where V is the midpoint of a bin and Freq is the number
of values that are inside the bin interval [V-BinWidth/2,V+BinWidth/2)
starting with the bin where V-BinWidth/2=Lower
/


  425bin(0,_L,_Min,_BW,[]):-!.
  426
  427bin(N,L,Lower,BW,[V-Freq|T]):-
  428  V is Lower+BW/2,
  429  Upper is Lower+BW,
  430  count_bin(L,Lower,Upper,0,Freq,L1),
  431  N1 is N-1,
  432  bin(N1,L1,Upper,BW,T).
  433
  434count_bin([],_L,_U,F,F,[]).
  435
  436count_bin([H-_W|T0],L,U,F0,F,T):-
  437  H<L,!,
  438  count_bin(T0,L,U,F0,F,T).
  439
  440count_bin([H-W|T0],L,U,F0,F,T):-
  441  (H>=U->
  442    F=F0,
  443    T=[H-W|T0]
  444  ;
  445    F1 is F0+W,
  446    count_bin(T0,L,U,F1,F,T)
  447  ).



 beta(+Alphas:list, -Beta:float) is det
Computes the value of the multivariate beta function for vector Alphas
https://en.wikipedia.org/wiki/Beta_function#Multivariate_beta_function
Alphas is a list of floats
/


  456beta(Par,B):-
  457  sumlist(Par,Sum),
  458  maplist(comp_lgamma,Par,LnGPar),
  459  LnG is lgamma(Sum),
  460  sumlist([-LnG|LnGPar],Exp),
  461  B is exp(Exp).
  462
  463comp_lgamma(X,LnG):-
  464  LnG is lgamma(X).



 average(+Values:list, -Average:float) is det
Computes the average of Values.
Values can be




  481average([H|T],Av):-
  482  number(H),!,
  483  sum_list([H|T],Sum),
  484  length([H|T],N),
  485  Av is Sum/N.
  486
  487average([H|T],E):-
  488  is_list(H),!,
  489  length(H,N),
  490  list0(N,L0),
  491  foldl(vector_sum,[H|T],L0,Sum),
  492  length([H|T],NV),
  493  matrix_div_scal([Sum],NV,[E]).
  494
  495average([H-W|T],E):-
  496  is_list(H),!,
  497  length(H,N),
  498  list0(N,L0),
  499  foldl(single_value_vect,[H-W|T],L0,Sum),
  500  foldl(agg_val,[H-W|T],0,SW),
  501  matrix_div_scal([Sum],SW,[E]).
  502
  503average(ValList,E):-
  504  foldl(single_value_cont,ValList,0,Sum),
  505  foldl(agg_val,ValList,0,SW),
  506  E is Sum/SW.
  507
  508
  509single_value_cont(H-N,S,S+N*H).
  510
  511single_value_vect(H-N,S0,S):-
  512  matrix_mult_scal([H],N,[H1]),
  513  matrix_sum([H1],[S0],[S]).



 agg_val(+Couple:atom, +PartialSum:number, -Sum:number) is det
Aggregate values by summation. The first argument is a couple _-N with
N the new value to sum to PartialSum
/


  521agg_val(_ -N,S,S+N).
  522
  523vector_sum(A,B,C):-
  524  matrix_sum([A],[B],[C]).



 variance(+Values:list, -Variance:float) is det
Computes the variance of Values.
Values can be




  539variance(L,Var):-
  540  variance(L,_Av,Var).



 variance(+Values:list, -Average:float, -Variance:float) is det
Computes the variance the average of Values.
Values can be




  555variance(L,Av,Var):-
  556  average(L,Av),
  557  maplist(sq_diff(Av),L,LS),
  558  average(LS,Var).



 std_dev(+Values:list, -Dev:float) is det
Computes the standard deviation of Values.
Values can be




  573std_dev(L,Dev):-
  574  std_dev(L,_Av,Dev).



 std_dev(+Values:list, -Average:float, -Dev:float) is det
Computes the standard deviation and the average of Values.
Values can be




  589std_dev(L,Av,Dev):-
  590  variance(L,Av,Var),
  591  root(Var,Dev).
  592
  593root(Var,Dev):-
  594  number(Var),!,
  595  Dev is sqrt(Var).
  596
  597root(Var,Dev):-
  598  maplist(sqroot,Var,Dev).
  599
  600sqroot(A,B):-
  601  B is sqrt(A).
  602
  603sq_diff(Av,A,S):-
  604  number(A),!,
  605  S is (A-Av)^2.
  606
  607sq_diff(Av,A-W,S):-
  608  number(A),!,
  609  S is W*(A-Av)^2.
  610
  611sq_diff(Av,A-W,S):-
  612  maplist(sq_diff,Av,A,S0),
  613  matrix_mult_scal([S0],W,[S]).



 swi_builtin(+Goal:atom) is det
Succeeds if Goal is an atom whose predicate is defined in Prolog
(either builtin or defined in a standard library).
/


  621:- dynamic
  622  builtin_cache/2.  623
  624swi_builtin(G):-
  625  builtin_cache(G, IsBuiltin),
  626  !,
  627  IsBuiltin == true.
  628swi_builtin(G):-
  629  functor(G, Name, Arity),
  630  functor(Gen, Name, Arity),
  631  (   builtin_int(Gen)
  632  ->  asserta(builtin_cache(Gen, true))
  633  ;   asserta(builtin_cache(Gen, false)),
  634      fail
  635  ).
  636
  637builtin_int(G):-
  638  predicate_property(G,built_in).
  639builtin_int(G):-
  640  predicate_property(G,imported_from(lists)).
  641builtin_int(G):-
  642  predicate_property(G,imported_from(apply)).
  643builtin_int(G):-
  644  predicate_property(G,imported_from(nf_r)).
  645builtin_int(G):-
  646  predicate_property(G,imported_from(matrix)).
  647builtin_int(G):-
  648  predicate_property(G,imported_from(clpfd)).
  649
  650
  651:- multifile sandbox:safe_primitive/1.  652sandbox:safe_primitive(cplint_util:bar(_,_)).
  653sandbox:safe_primitive(cplint_util:bar(_,_,_)).
  654sandbox:safe_primitive(cplint_util:bar1(_,_)).
  655sandbox:safe_primitive(cplint_util:argbar(_,_)).
  656sandbox:safe_primitive(cplint_util:histogram(_,_)).
  657sandbox:safe_primitive(cplint_util:histogram(_,_,_)).
  658sandbox:safe_primitive(cplint_util:densities(_,_,_)).
  659sandbox:safe_primitive(cplint_util:densities(_,_,_,_)).
  660sandbox:safe_primitive(cplint_util:density(_,_)).
  661sandbox:safe_primitive(cplint_util:density(_,_,_)).
  662sandbox:safe_primitive(cplint_util:density2d(_,_)).
  663sandbox:safe_primitive(cplint_util:density2d(_,_,_)).
  664sandbox:safe_primitive(cplint_util:to_pair(_,_)).
  665sandbox:safe_primitive(cplint_util:key_pair(_,_)).
  666sandbox:safe_primitive(cplint_util:value_pair(_,_)).
  667sandbox:safe_primitive(cplint_util:bin(_,_,_,_,_)).
  668sandbox:safe_primitive(cplint_util:beta(_,_)).
  669
  670
  671
  672:- multifile license:license/3.  673
  674license:license(artisticv2, permissive,
  675                [ comment('Artistic License 2.0'),
  676                  url('https://opensource.org/licenses/Artistic-2.0')
  677                ]).
  678
  679:- license(artisticv2).