#Use:> read "c:/cygwin/home/cr569/vTool.mpl";

with(LinearAlgebra):

Vols := module()
     export Dag,pa,ch,fillts,Tobin,BlanketOf,b2,b10,posit,vpa,Tprod,ProductProb,Prob,merge,Wi,JP,VolMCMC,VolExact,nparams,DofI,Winaive,FInt,showVolsTool,B,Explode,Collapse,Line,LowVol,HighVol,SCformula,MeanCurvature,LS,fill_g,IntRicci;

Dag := proc(n,ls)
   description "returs an n by n adjacency matrix with links given in the list ls=[i1,j1,i2,j2,...]";
   local i,k,TS;
   global dag,T,S;
   k := nops(ls); if (is(k,odd)) then print("WRONG list of connectors!"); end if;
   with(LinearAlgebra):
   dag := Matrix(n,n);
   for i from 1 to k by 2 do
       dag[ls[i],ls[i+1]] := 1;
   end do;
   TS := fillts(dag);
   T := TS[1];
   S := TS[2];
   return(dag);
end:

pa := proc(j)
 description "returns the set of parents of node j";
 local s,i;
 s := {};
 for i from 1 to op(1,dag)[1] do
   if dag[i,j]=1 then s := s union {i} end if;
 end do;
return s
end:

ch := proc(i)
  description "returns the set of children of node i";
  local s,j;
  s := {};
  for j from 1 to op(1,dag)[1] do
    if dag[i,j]=1 then s := s union {j} end if;
  end do;
return s
end:

fillts := proc(dag)
 description "returns two vectors: tt=first index of the parameter for each node. st= 2^nops(pa(i))";
 local i,j,tt,st,nnodes,s,np;
 nnodes := op(1,dag)[1]: s := 0: 
 tt := Vector[row](nnodes):
 st := Vector[row](nnodes):
 for i from 1 to nnodes do
    np := nops(pa(i));
    st[i] := 2^np;
    tt[i] := s + 1;
    s := s + st[i];
 end do;
 return tt,st
end:

Tobin := proc(i,k)
  description "transform to binary parents of i =k";
  local j,np,rs,q,r; 
  np := nops(pa(i));
  rs := Vector[row](np); q := k;
  for j from 1 to np do
   q := iquo(q,2,'r');
   rs[np-j+1] := r;
  end do;
  return rs
  end:

BlanketOf := proc(s)
  description "returns the Markov Blanket of the set s";
 local i,j,b;
 b := {};
 for i in s do
  b := b union pa(i)
 end do;
 return b minus s
end:

b2 := proc(k,j)
  description "returns the k least significant binary digits of j";
  local i,q,r,rs;
  rs := Vector[row](k);
  q := j;
  for i from 1 to k do
    q := iquo(q,2,'r');
    rs[k-i+1] := r;
  end do;
  return rs
end:

b10 := proc(v)
  description "returns the base 10 equiv. of the vector of binary digits in v";
  local i,n,x,u; u := convert(v,list);
  x := 0; 
  n := nops(u); 
  for i from 1 to n do
    if (u[i]=1) then x := x+2^(n-i) end if;
  end do;
  return x;
end:

posit := proc(a,s)
  description "returns the position of a in the list s";
  local k;
  member(a,s,'k');
  return k;
end:

vpa := proc(i,s,v)
   description "returns the value of the parents of node i when the set of nodes s has values specified by the binary vector v";
   local par,j,k,r;   
   par := pa(i);
   k := nops(par);
   r := Vector[row](k);
   for j from 1 to k do
     r[j] := v[posit(par[j],s)];
   end do;
   return r
end:

Tprod := proc(i,s,v,T)
  description "returns (1-tn) if ith node has value 0 and tn otherwise. The value of n is computed assuming the set of nodes s has given vector of values v which must include all the parents of i. T is the output of fillts()";
  local n,u,xi; u := vpa(i,s,v); xi := v[posit(i,s)];
  n := T[i] + b10(u);
  if (xi=0) then return (1-t||n) else return t||n end if
end:

ProductProb := proc(s1,s2,v)
  description "returns the product of the probabilities of nodes in the set s1 conditional of their parents when all the values of s1 union s2 are given by the binary vector v. It is assumed that s2 contains the BlanketOf(s1)";
  local i,prod,s; prod := 1;
  s := s1 union s2;
  for i from 1 to nops(s1) do prod := prod*Tprod(s1[i],s,v,T); end do;
  return prod
end:

Prob := proc(s,j)
  description "returns the probability that the nodes in the set s have the values j when written in binary";
  local m,tot,bs,v,vj,vm,k,n; n := nops(s); if (n=0) then return 1 end if;
  bs := BlanketOf(s);
  tot := 0;
  k := nops(bs);
  vj := b2(n,j);
  for m from 0 to 2^k-1 do
    vm := b2(k,m);
    v := merge(s,vj,bs,vm);
    tot := tot + ProductProb(s,bs,v)*Prob(bs,m);
  end do;
  return tot;
end:

merge := proc(s1,v1,s2,v2)
  description "returns the values of s1 union s2 when the nodes in set s1 have values v1 and those in s2 have values v2";
  local i,ns,si,s,v;
  s := s1 union s2;
  ns := nops(s);
  v := Vector[row](ns);
  for i from 1 to ns do
    si := s[i];
    if (si in s1) then v[i] := v1[posit(si,s1)];
      else v[i] := v2[posit(si,s2)];
    end if;
  end do;
  return v;
end:

Wi := proc(i)
  description "returns the contribution of ith node to the det of Fisher info";
  local np,j,tn,s,uprod,dprod;
  s := pa(i);
  np := nops(s); uprod := 1; dprod := 1;
  for j from 0 to 2^np-1 do 
    uprod := uprod*Prob(s,j);
    tn := t || (T[i]+j);
    dprod := dprod*tn*(1-tn);
    end do;
  return uprod/dprod;
end:

JP := proc(dag)
  description "Jeffreys Prior for dag: sqrt(det(I))"; 
  local T,J,i,n,prod;
  T := fillts(dag);
  n := op(1,dag)[1]; prod := 1;
  for i from 1 to n do
    prod := prod*Wi(i);
  end do;
  return sqrt(prod);
end:

VolMCMC := proc(dag,eps)
  description "returns the volume of the dag computed by MCMC with epsilon = eps";
  local i,j,n,s,T,N,prod,nc,np,npis,uprod,tts,childless;
  T := fillts(dag); childless := {};
  n := op(1,dag)[1]; N:= nparams(dag); prod := 1; npis := 0; 
  for i from 1 to n do
    nc := nops(ch(i));
    s := pa(i);   
    np := nops(s);
    if (nc=0) then
        uprod := 1; childless := childless union {i};
        for j from 0 to 2^np-1 do 
        uprod := uprod*Prob(s,j);
        end do;
        prod := prod*uprod;
        npis := npis+2^np;
    else prod := prod*Wi(i);
    end if;
  end do;
  tts := DofI({seq(i,i=1..n)} minus childless);
  return Pi^npis*Int(sqrt(prod),tts,method=_MonteCarlo, epsilon=eps);
end:

VolExact := proc(dag)
  description "returns the volume of the dag as an exact multiple integral";
  local i,j,n,s,T,N,prod,nc,np,npis,uprod,tts,childless;
  T := fillts(dag); childless := {};
  n := op(1,dag)[1]; N:= nparams(dag); prod := 1; npis := 0; 
  for i from 1 to n do
    nc := nops(ch(i));
    s := pa(i);   
    np := nops(s);
    if (nc=0) then
        uprod := 1; childless := childless union {i};
        for j from 0 to 2^np-1 do 
        uprod := uprod*Prob(s,j);
        end do;
        prod := prod*uprod;
        npis := npis+2^np;
    else prod := prod*Wi(i);
    end if;
  end do;
  tts := {seq(i,i=1..n)} minus childless;
  return Pi^npis*FInt(sqrt(prod),tts);
end:


nparams := proc(dag)
  description "returns the total number of parameters of the dag";
  local N,n,i;
  n := op(1,dag)[1]; N := 0;
  for i from 1 to n do 
     N := N + 2^nops(pa(i));     
  end do;
return N;
end:

DofI := proc(s)
   description "returns the domain of integration for all the nodes in the set s";
   local i,j,r;
   r := {};
   for i in s do
      for j from T[i] to T[i]+2^nops(pa(i))-1 do
         r := r union{t||j=0..1};
      end do;
   end do;
return convert(r,'list');
end:

Winaive := proc(i)
  description "returns the contribution of ith node to the det of Fisher info for the naive approximation";
  local np,j,tn,dprod,s,k,sdprod;
  s := pa(i);
  np := nops(s); dprod := 1; k := 2^np;
  for j from 0 to k-1 do 
    tn := t || (T[i]+j);
    dprod := dprod*tn*(1-tn);
    end do; 
  return (1/k)^k/dprod;
end:

FInt := proc(e,s)    
  description "FInt(e,s): returns the multiple integral of the expression e over the set of variables of nodes in s";
     local i,j,tn,r; 
     tn := t||(s[1]);
     r := int(e,tn=0..1);
     for j from T[s[1]] to T[s[1]]+2^nops(pa(s[1]))-1 do
         tn := t||j;
	 r := int(r,tn=0..1);
     end do;
     for i in (s minus {s[1]}) do
         for j from T[i] to T[i]+2^nops(pa(i))-1 do 
               tn := t||j;
	       r := int(r,tn=0..1);
          end do;
      end do;
 return r;
end:

showVolsTool := proc()
   description "Prints all the descriptions in the VolsTool";
   local i;
   print ("Dag(n,ls): returs an n by n adjacency matrix with links given in the list ls=[i1,j1,i2,j2,etc]");

 print ("pa(j): returns the set of parents of node j");

  print ("ch(i): returns the set of children of node i");

 print ("fillts(dag): returns the first index of the parameter for each node");

  print ("Tobin(i,k): transform to binary parents of i =k");

  print ("BlanketOf(s): returns the Markov Blanket of the set s");

  print ("b2(k,j): returns the k least significant binary digits of j");

  print ("b10(v): returns the base 10 equiv. of the vector of binary digits in v");

  print ("posit(a,s): returns the position of a in the list s");

   print ("vpa(i,s,v): returns the value of the parents of node i when the set of nodes s has values specified by the binary vector v");

  print ("Tprod(i,s,v,T): returns (1-tn) if ith node has value 0 and tn otherwise. The value of n is computed assuming the set of nodes s has given vector of values v which must include all the parents of i. T is the output of fillts()");

  print ("ProductProb(s1,s2,v): returns the product of the probabilities of nodes in the set s1 conditional of their parents when all the values of s1 union s2 are given by the binary vector v. It is assumed that s2 contains the BlanketOf(s1)");

  print ("Prob(s,j): returns the probability that the nodes in the set s have the values j when written in binary");

  print ("merge(s1,v1,s2,v2): returns the values of s1 union s2 when the nodes in set s1 have values v1 and those in s2 have values v2");

  print ("Wi(i): returns the contribution of ith node to the det of Fisher info");

  print ("JP(dag): Jeffreys Prior for dag: sqrt(det(I))"); 

  print ("VolMCMC(dag,eps): returns the volume of the dag from MCMC. eps is the desired precission");

  print ("VolExact(dag): returns the volume of the dag as an exact multiple integral");

  print ("LowVol(dag); returns a Lower Bound on the Volume of the dag");

  print ("HighVol(dag); returns an Upper Bound on the Volume of the dag");

  print ("nparams(dag): returns the total number of parameters of the dag");

   print ("DofI(s): returns the domain of integration for all the nodes in the set s");

  print ("Winaive(i): returns the contribution of ith node to the det of Fisher info for the naive approximation");

 print ("pa(j): returns the set of parents of node j");


end:


alias(g=GAMMA):
#B(n) = Beta(n/2+1,n/2+1)=int(t^(n/2)*(1-t)^(n/2),t=0..1);
B := n -> GAMMA(n/2+1)^2/GAMMA(n+2):

# Explode(n) is the volume of the exploding star: one node with n children;
Explode := n -> Pi^(2*n)*B(n-1):

# Collapse(n) is the volume of the collapsing star: one node with n parents;
Collapse := n -> Pi^(2^n)*B(2^(n-1)-1)^n:

# Line(n) is the volume fo the straight line of n nodes;
Line := n -> 4*(Pi/2)^(3*n-4):

LowVol := proc(dag)   
    description  "LowVol(dag); returns a Lower Bound on the Volume of the dag";
    local p,i,j,ai,nnodes;
    nnodes := op(1,dag)[1];
    p := 1;
    for i to nnodes do
        ai := 0;
	for j in ch(i) do
	    ai := ai + (S[j]/S[i]);
        end do; 
        p := p*B(ai/2-1)^S[i];
    end do;
    return evalf(p);
end:

HighVol := proc(dag)    
     description "Upper bound for volume of dag";
     local i,j,bi,nnodes,p,hits,th,C;
     nnodes := op(1,dag)[1];
     hits := Vector[row](nnodes);
     p := 1;
     for i to nnodes do
         bi := 0;
         for j in ch(i) do
             if (pa(i) subset pa(j)) then bi := bi + S[j]/S[i];
             else hits[j] := hits[j]+1;
             end if;
         end do;
         p := p*B(bi/2-1)^S[i];
     end do;
     C := 1; i:=1;
     for i to nnodes do
         if (hits[i]>0) then
	      th := 2^hits[i];
              C := C*(1/th)^(S[i]/2);
         end if;
     end do;
     return evalf(p*C);
end:

fill_g := proc(dag)
  local j,Ti,tk,k;
  global S,T,dim,nnodes,g_compts; 
  for k from 1 to nnodes do
    Ti := T[k];
    for j from 0 to S[k]-1 do
      tk := t||(Ti+j); 
      g_compts[Ti+j,Ti+j] := Prob(pa(k),j)/(tk*(1-tk));
    end do;
  end do;
  eval(g_compts);
end:

SCformula := proc(n,ls)
description "SCformula(n,ls): returns the formula for the Ricci scalar for the dag of n nodes  and with links given in the list ls=[i1,j1,i2,j2,...]";
 local coord,g,ginv,D1g,D2g,Cf1,RMN,RICCI,RS,TS;
 global S,T,dim,nnodes,g_compts,dag;

 nnodes := n;
 dag := Dag(n,ls);
 TS := fillts(dag);
 T := TS[1];
 S := TS[2];
 with(tensor):
 dim := nparams(dag);
 coord := [seq(t||i,i=1..dim)]:
 g_compts := array(symmetric,sparse, 1..dim, 1..dim):
 
 fill_g(dag):
 alias(g = 'g'): 
 g := create( [-1,-1], eval(g_compts));
 ginv := invert( g, 'detg' ):
 D1g := d1metric ( g, coord ):  D2g := d2metric ( D1g, coord ):
 Cf1 := Christoffel1 ( D1g ):
 RMN := Riemann( ginv, D2g, Cf1 ):
 RICCI := Ricci( ginv, RMN ):
 RS := Ricciscalar( ginv, RICCI ):

 RETURN(-RS[compts]);
end:

IntRicci := proc(n,ls)
description "returns the integral of the Ricci scalar over the manifold represented by the dag of n nodes and links in ls=[i1,j1,i2,j2,...]";
  local i,j,s,T,N,prod,nc,np,npis,uprod,tts,childless;
  global dag;

  dag := Dag(n,ls);
  T := fillts(dag); childless := {};
  N:= nparams(dag); prod := 1; npis := 0; 
  for i from 1 to n do
    nc := nops(ch(i));
    s := pa(i);   
    np := nops(s);
    if (nc=0) then
        uprod := 1; childless := childless union {i};
        for j from 0 to 2^np-1 do 
        uprod := uprod*Prob(s,j);
        end do;
        prod := prod*uprod;
        npis := npis+2^np;
    else prod := prod*Wi(i);
    end if;
  end do;
  tts := {seq(i,i=1..n)} minus childless;
  return Pi^npis*FInt(SCformula(n,ls)*sqrt(prod),tts);
end:

MeanCurvature := proc(dag)
local n,ls,i,j;
n := op(1,dag)[1];
ls := LS(dag);
return (IntRicci(n,ls)/VolExact(dag));
end;

LS := proc(dag)
local n,ls,i,j;

ls := [];

n := op(1,dag)[1];
for i to n do
    for j from i+1 to n do
	if (dag[i,j] > 0) then 
             ls := [op(ls),i,j] end if;
    end do;
end do;
return ls;
end;

end module;

with(Vols):