BeginPackage["Qlib`"]

(************************************************************

The Qlib program library has resulted from the work conducted 
over several years jointly by 

  Laboratory of Telecommunications Technology of 
  Helsinki University of Technology (HUT) and VTT 
  Information Technology, Telecommunications
  
The following individuals have contributed to the library:

  Samuli Aalto, Esa Hyytiä, Jani Lakkakorpi, Ilkka Norros, 
  Anssi Pirhonen, Vesa Timonen and Jorma Virtamo.

The programs can be used and modified freely.

No claims are made about the correctness of the programs
and no liability is taken for any damage caused by the use
of the programs in the library.

Reports of bugs in the program can be sent to qlib@netlab.hut.fi

__
UpDates:

20.11.2001  Functions 'CreateNet', 'Dijkstra' and 'ShowRoutes' 
             added by Vesa Timonen.
08.02.2002  Added KaufmanRoberts -recursion by Esa Hyytiä
19.04.2002  Function 'Berli' updated by Vesa Timonen.
04.05.2002  Function 'Dijkstra' updated, now functioning also
             with networks not compitely connected
23.05.2002  Functions 'CreateNet' and 'ShowRoutes' updated and
             renamed as 'GabrielNet' and 'ShowNetwork'
30.05.2002  Minor modifications to GabrielNet/ShowNetwork:
             Context of optional arguments forced to 'System`'
             Now because context of opt.args in interface.tm is 
             the same, no shadowing appears while simultaneous use
16.01.2003  Added functions 'NormalizationConstant', 'Throughputs' 
	     and 'TreeGraph' by Vesa (functions implemented by JV)
*************************************************************)

(**************** *********************************************)

qlib::usage =
"QLib: Queueing Library\n\n

Markov Processes: CtMarkovChain, DtMarkovChain\n\n

Loss Systems: MVA, KaufmanRoberts\n\n

Erlang B:  Berli, Berld, Xerl, Aerl\n\n

Queueing Systems: Qmd1, Qsdd1, Qndd1, Qmxd1, Mg1, SumMd1, RecMd1, IntFmd1,
 Fmd1, Mdn, IntFmdn, Fmdn, Fekdn, Fend1\n\n

Networks: GabrielNet, Dijkstra, ShowNetwork\n\n

Trees and balanced fairness: NormalizationConstant, Throughputs ,TreeGraph.\n"

(**************** Markov Processes ****************************)

CtMarkovChain::usage =
"CtMarkovChain[Q_] gives the stationary probabilities of a continuous
time Markov chain with the transition rate matrix Q."

DtMarkovChain::usage =
"DtMarkovChain[P_] gives the stationary probabilities of a discrete
time Markov chain with the transition probability matrix P."

(**************** Loss Systems ********************************)

MVA::usage =
"MVA[R_,mu_,K_] is an algorithm for the Mean Value Analysis of a
closed Jackson network. It returns the average queue lengths and the
average sojourn times in the queues. The branching ratios are given by
the matrix R, i.e. the element R_{i,j} gives the probability for a
customer leaving queue i to join queue j next; the vector mu specifies
the service rates of the queues; and K is the number of customers in
the network."

KaufmanRoberts::usage =
"KaufmanRoberts[ linkC, bandwidth, load], returns blocking probabilities of each class
by using the Kaufman and Roberts recursion formula."

Berli::usage = 
"Berli[n,a] gives Erlang loss probability for integer trunk value n >= 0
and Poisson traffic intensity a."

Berld::usage = 
"Berld[x,a] gives Erlang loss probability for real trunk values x >= 0
and Poisson traffic intensity a."

Xerl::usage = "
Xerl[b,a] solves the trunk value x from 'b = Berld[x,a]'."

Aerl::usage = "
Aerl[x, b] solves the traffic intensity a from 'b = Berld[x,a]'."


(**************** M/D/1 etc Queuing Systems *******************)

Qmd1::usage = "
Qmd1[x, rho] is the M/D/1 queue length distribution function with
system load level rho."

Qsdd1::usage = "
Qsdd1[x, {D}] is the Sigma D/D/1 queue length distribution function
for multiple sources with each having it's own period. Table of periods
is given in D."

Qndd1::usage = "
Qndd1[x, N, D] is the N*D/D/1 queue length distribution function
for N sources all with period D."
  
Qmxd1::usage = "
Qmxd1[x, rho, {p1, p2, p3, ...}] is the M^x/D/1 unfinished work tail
distribution function with load level rho and batch size probabilities
{p1, p2, p3, ... pn}. pn = P(batch size = n)."

Mg1::usage = "
Mg1[x, rho, D] is the M/G/1 queue length probability function
with load level rho and service time distribution function D.
It is assumed that the average service time is 1."

SumMd1::usage = "
SumMd1[x, rho] calculates the state probabilities of a M/D/1 queue
with load level rho."

RecMd1::usage = "
RecMd1[x, rho] calculates the state probabilities of a M/D/1 queue
with load level rho. This function uses a recursive algorithm."

IntFmd1::usage = "
IntFmd1[N, rho] is the M/D/1 (virtual) waiting time distribution
function (for integral values of x) with load level rho. Service time
h is assumed to be 1."

Fmd1::usage = "
Fmd1[x, rho] is the M/D/1 (virtual) waiting time distribution
function (it calculates the probability that waiting time W is smaller
than x) with load level rho. Service time h is assumed to be 1."

Mdn::usage = "
Mdn[x, rho, n] calculates the state probabilities of a M/D/n queue
with load level rho and n servers. Accuracy goal is 10^-12."

IntFmdn::usage = "
IntFmdn[N, rho, n] is the M/D/n (virtual) waiting time distribution
function (for integral values of x) with load level rho. Service time
h is assumed to be 1."

Fmdn::usage = "
Fmdn[x, rho, n] is the M/D/n (virtual) waiting time distribution
function (it calculates the probability that waiting time W is smaller
than x) with load level rho. Service time h is assumed to be 1."

Fekdn::usage = "
Fekdn[x, rho, k, n] is the E_k/D/n (virtual) waiting time distribution
function (it calculates the probability that waiting time W is smaller
than x) with load level rho. Service time h is assumed to be 1."

Fend1::usage = "
Fend1[x, n, rho] gives the E_n/D/1 waiting time distribution (it
calculates the probability that waiting time W is smaller than x)
with n arrival phases and system load level rho. Service time is
assumed to be 1."



(********** Routines for Network generation etc ****************)

GabrielNet::usage =
"GabrielNet[NbrOfNodes] generates a random network containing 'NbrOfNodes' 
nodes. The positions of the nodes are drawn from uniform distribution and 
node pairs satisfying the gabriel criterion are connected. Function returns 
a list containing two or three elements: distance/cost matrix, list of the node 
coordinates and, if wanted, a Graphics object containing an image of the network.\n
Optional arguments:\n
 Shape      \tCircle/Box\t the shape of the network (default 'Circle')\n
 Scale      \t          \t the size of the circle radius/box edge (def. '1')\n
 CreateGraph\tFalse/True\t whether to create the graph (def. 'False')\n
 Index      \tFalse/True\t whether to include node indices to the graph (def. 'False')\n
 Graph      \tFalse/True\t whether to show the graph (def. 'False')"

Dijkstra::usage =
"Dijkstra[dist, root, dest] finds out the shortest path from node 'root' to 
node 'dest' using Dijkstras algorithm. If 'dest' is not given, function finds 
out paths to all nodes in the net. Matrix 'dist' is the distance/cost matrix 
between all nodes in the network (if there is no connection, the distance/cost 
must be 'Infinity'). Function returns a list containing three elements: list 
of nodes that belong to the path, list of previous nodes of corresponding node 
in the first list and total distance/cost of corresponding node in the first 
list to the node 'root'. \nIf there is no connection to some nodes, those are left
out from the results. If there is no connection from 'root' to specific 'dest',
result will be list containing only the root node."

ShowNetwork::usage =
"ShowNetwork[netData, pathData] List 'netData' contains information of the network 
in the output form of function 'GabrielNet' and list 'pathData' contains 
information of the route(s) in the output form of function 'Dijkstra'. 
Function returns a list of two Graphics objects that contains an image of the 
network and the path(s).\nIf list 'pathData' is empty i.e. {}, returns the function
only the graph of the network.\n
Optional arguments:\n
 Index \tFalse/True \t whether to include node indices to the graph (default 'False')\n
 Graph \tFalse/True \t whether to show the created graph (def. 'False')"


(* Defaults for optitonal arguments of functions GabrielNet, ShowNetwork *)
(*************************************************************************)
Options[GabrielNet] = {System`Scale -> 1,
                       System`Shape -> Circle, 
                       System`CreateGraph -> False, 
                       System`Index -> False, 
                       System`Graph -> False};
Options[ShowNetwork] = {System`Index -> False, System`Graph -> False};
(*************************************************************************)

(********** Routines for Trees, Balanced Fairness etc ****************)

NormalizationConstant::usage = 
"NormalizationConstant[List] Calculates the normalization constant for a tree 
network with balanced fair capacity allocation. Input is a list describing the 
tree and offered flows."

Throughputs::usage =
"Throughputs[List] Calculates the normalization constant and flow throughputs
for a tree network with balanced fair capacity allocation. Input is a list 
describing the tree and offered flows. Output is a list: 
{normalization constant, {throughputs}, {flows}}."

TreeGraph::usage =
"TreeGaph[List] Returns a graphics object for drawing a tree network.
Input is a list describing the tree and offered flows."


(**************** *********************************************)


Begin["`Private`"]

CtMarkovChain[Q_] := Module[{n=Length[Q],EE,ee},
  EE=Table[1,{n},{n}];
  ee=Table[1,{n}];
  Return[ee.Inverse[Q+EE]//Chop]]


DtMarkovChain[P_] := Module[{n=Length[P],EE,ee},
  EE=Table[1,{n},{n}]-IdentityMatrix[n];
  ee=Table[1,{n}];
  Return[ee.Inverse[P+EE]//Chop]]



(**************** *********************************************)

MVA[R_,mu_,K_] := Module[{n=Length[R],EE,ee,la,T,nn},
  EE=Table[1,{n},{n}]-IdentityMatrix[n];
  ee=Table[1,{n}];
  nn=Table[0,{n}];
  la=ee.Inverse[R+EE];
  Do[T=(1+nn)/mu; nn=k*la*T/(la.T), {k,1,K}]; 
  {nn,T}]

(* Kaufman Roberts recursion *)

KaufmanRoberts[ linkC_, bandwidth_, load_]:=Module[{G,Gnorm,c},
    G[x_ /; x < 0] = 0; G[0] = 1;
    (* Laskentarekursio :*)
    Do[G[c] = 1/c Plus @@ (load*(#*G[c - #] & /@ bandwidth)), {c, 1, linkC}];
    Gnorm = Plus @@ (G[#] &) /@ Range[0, linkC];
    (Plus @@ (G[#] & /@ (linkC + 1 - Range[#])))/Gnorm & /@ bandwidth  ];


Berli[n_, a_] := 
  Module[{e = 1}, Do[e = 1 + i*e/a, {i, 1, n}]; 1/e ]

Berld[x_, a_] := N[a^x Exp[-a]/Gamma[x+1,a]] /; (x >= 0) && (a > 0)


Xerl[b_, a_] :=
Module[{b2, s, n, x},
  b2 = 1;
  s = 1/a;
  For[x = s, b2 > b, x+=s, 
      b2 /= x + b2];
  x -= s; 
  n = x*a - 1;
  x /. FindRoot[Berld[x, a] == b, {x, {n, n+1}}]
] /; (0 < b < 1) 
  
  
Aerl[x_, b_] :=
Module[{c = N[Sqrt[2/(Pi x)]], start, a},
  start = If[b >= c,
    x + (b-c)/(c^2),
    ((Sqrt[-2 Log[(2-b/c) b/c]] - 2 Sqrt[x])^2+1)/4
  ];
  a /. FindRoot[a (b - Berld[x,a]) == 0, {a, {start, start + 1}}]
] /; (0 < b < 1)


(**************** M/D/1 etc Queuing Systems *******************)


Getzval[n_, N_, x_, D_, d_, p_] :=
Module[{df0 = 0, fz = 0, z},
  df0 = Sum[p[[i]] / (1 - p[[i]]), {i,N}] ;
  z = N / df0 * (n-d) / (N - (n - d)) ;
  fz = Sum[p[[i]] z / (1 - p[[i]] + p[[i]] z), {i,N}] ;
  
  While[Abs[fz - n + d] > 0.001,
      z = (n - d) / (df0 - (n-d)*(df0/fz - 1/z));
      fz = Sum[p[[i]] z / (1 - p[[i]] + p[[i]] z), {i,N}] ;
  ];
  z ];
  
  
Qsdd1[x_, D_List] := 
Module[{d, z, p, th, s, t, ans = 0, N = Length[D], n, k, rho, Frac, x2 },
  Frac[x2_] := x2 - Floor[x2];
  rho = Sum[1 / D[[i]], {i, N}];
  For[n = x + 1, n <= (N - rho x)/(1 - rho), n++,
      d = Sum[Floor[(n-x) / D[[i]]], {i,N}] ;
      If [(n - d) >= N,
	  Continue[]];
      p = Table[Frac[(n-x)/D[[i]]], {i,N}] ;
      z = Getzval[n, N, x, D, d, p];
      th = Product[1 - p[[i]] + p[[i]] z, {i,N}];
      s = Sum[p[[i]]z(1-p[[i]])/(1-p[[i]]+p[[i]]z)^2, {i,N}];
      t = Sum[1 / D[[i]] / (1-p[[i]]+p[[i]]z),{i,N}];   
      ans += th / z^(n-d) / Sqrt[2 Pi s]  (1 - t);
  ];
  ans
]
  
  
Qmd1[x_, rho_] :=
Module[{n,t},
  If [rho < .3,
      If [x < 9 + 15 Log[10, .3 / rho],
	  Qmd1series[x, rho],
	  Qmd1asymp[x, rho]
      ],
      If [x > 8,
	  Qmd1asymp[x, rho],
	  Qmd1poly[x, rho]
      ]
  ]
] 
      
      
Qmd1series[x_, rho_] :=
Module[{sum = 0, logrho, prev = 0, lf, fa = 1, n},
  For[n = Floor[x] + 1, n < x + 50, n++, 
      sum += (rho(n-x))^n / n! * Exp[-rho n];
  ];
  sum  Exp[rho x] (1 - rho)
]


Qmd1poly[x_, rho_] :=
Module[{P = {(1-rho) Exp[x*rho]}, a = rho Exp[-rho], i = Floor[x], 
        y = x - Floor[x], v = {}, d = {1}, k = 1},
	Do [
	  AppendTo[v, -a/j];
	  P = Prepend[v P, Plus @@ P];
	  AppendTo[d, k *= y],
	{j, i}];
        1 - P . d
]
      
      
Qmd1asymp[x_,r_] := ((1-r)/(r E^z-1)*Exp[-z*x])/. FindRoot[r(E^z-1)-z,{z,(3-r)(1-r)Log[10,10/r]},
  AccuracyGoal->15,MaxIterations->100]


Qndd1[x_, N_, D_] :=
Module[{sum = 0, bin},
  bin = 1;
  For[n = N, n > x, n--,  
      sum += bin ((n-x)/D)^n (1-(n-x)/D)^(N-n) (D-N+x)/(D-n+x);
      bin *= n/(N-n+1);
  ];
  sum
]


Qmxd1[x_, rho_, p_] :=
Module[{n = Floor[x], b, lambda, f, beta, i, j, k, z},
  b = Sum[i*p[[i]], {i,Length[p]}];
  lambda = rho/b;
  beta = p*lambda*Exp[-lambda*Range[Length[p]]];
  f[1,1] = (1 - rho);
  For[j = 2, j <= n+1, j++,
      f[1,j] = 0;
   ];
  For[i = 2, i <= n+1, i++,
      f[i,1] = Sum[f[(i-1),k], {k,n+1}];
      For[j = 2, j <= n+1, j++,
	  If [Length[p] >= (i-1),
	      f[i,j] = -1/(j-1) Sum[beta[[z]]*f[(i-z),(j-1)],{z,i-1}],
	      f[i,j] = -1/(j-1) Sum[beta[[z]]*f[(i-z),(j-1)],{z,Length[p]}]
	  ];  
      ];
  ];
  1 - Exp[lambda x] (f[n+1,1] + Sum[f[n+1,k] (x-n)^(k-1), {k,2,n+1}])
]
      
      
Mg1[x_Integer, rho_, D_] :=
Module[{lambda = rho, n, j, t, k, a, p},
  For[n = 0, n <= x, n++,
      a[n] = NIntegrate[(Exp[-lambda t] (lambda t)^n/n! (1-D[t])), {t, 0, Infinity}]
  ];
  p[0] = (1-rho);
  For[j = 1, j <= x, j++,
      p[j] = 1/(1 - lambda a[0]) (lambda a[j-1] p[0] + lambda Sum[a[j-k] p[k], {k,(j-1)}])
  ]; 
  Table[p[j], {j,0,x}]          
]
      

SumMd1[x_, rho_] :=
Module[{n = 1},
  If [(x == 0),
      1-rho,
      If [(x == 1),
	  ( (1-rho) (E^rho - 1) ),
	  If [(x > 1),
	      ( (1-rho) Sum[(-1)^(x-n) E^(n rho)
	      ( (n rho)^(x-n)/(x-n)! + (n rho)^(x-n-1)/(x-n-1)! ),
	      {n,1,x}] ),
	      0
	  ]
      ]
  ]
]


RecMd1[x_, rho_] :=
Module[{n=1, i=1, P},
  P[0] = 1-rho;
  P[1] = (1-rho) (E^rho-1);
  For[i=1, i<x, i++,
      P[i+1] = 1/(E^-rho) ( P[i] - (P[0]+P[1]) rho^i/i! E^-rho
               - Sum[P[n] rho^(i-n+1)/(i-n+1)! E^-rho,{n,2,i}] )];
  Table[P[i], {i,0,x}]
]


IntFmd1[N_, rho_] := (Plus @@ RecMd1[N,rho])


Fmd1[x_, rho_] :=
Module[{n=0},
  T=Floor[x];
  t=x-T;
  v=RecMd1[T,rho];
  If[ (t == 0),
      Plus @@ v,
      E^(rho t) Sum[(-rho t)^n/n! (Plus @@ Take[v,T-n+1]), {n,0,T}]
  ]
]


Fekdn[x_, rho_, k_, n_] := Fmdn[x, rho, k n]


IntFmdn[x_, rho_, n_] := (Plus @@ Mdn[n(x+1)-1, rho, n])


Fmdn[x_, rho_, n_] :=
Module[{i=0, j=0, v=0, V},
  V=Mdn[n-1, rho, n];
  If[x == 0.,
      Plus @@ V,
      Sum[(Plus @@ Take[V,i+1]) * 
	  (Sum[(((rho n) (v-x))^(v n + n - 1 - i))/(v n + n - 1 - i)! E^(-(rho n) (v-x)),
	  {v,0,Floor[x]}]), {i,0,n-1}]
  ]
]


Mdn[x_, rho_, n_] :=
Module[{i = 0, I = 0, GuessP, Sum1, P, S, Difference, Max = 0},
  
  GuessP[0] = ( Sum[(rho n)^i/i!,{i,0,n-1}] + (rho n)^n/n! 1/(1-(rho n)/n) )^-1;
  For[i=1, i<n, i++,
      GuessP[i] = GuessP[i-1] (rho n)/i;
  ];
  i = n - 1; 
  While[ (GuessP[i] > 10^-12),
      i=i+1;
      GuessP[i] = GuessP[i-1] (rho n)/n;
  ];
  I = i; 
  For[i=(I+1), i<=(I+n), i++,
      GuessP[i] = GuessP[i-1] (rho n)/n;
  ];

  While[ ((Max > 10^-12) || (Max == 0)),
      Max=0;
      Sum1 = Sum[GuessP[k], {k,0,n}];
      For[i=0, i<=I, i++,
	  P[i] = Sum1 (rho n)^i/i! E^-(rho n)
	  + Sum[GuessP[k] (rho n)^(n+i-k)/(n+i-k)! E^-(rho n), {k,n+1,n+i}];
      ];
      For[i=(I+1), i<=(I+n), i++,
	  P[i] = P[i-1] (rho n)/n;  
      ];
      S = Sum[P[k], {k,0,I+n}];
      For[i=0, i<=(I+n), i++,
	  P[i] = P[i]/S;
      ];
      For[i=0, i<=(I+n), i++,
	  Difference = Abs[P[i] - GuessP[i]];
	  If [Difference > Max,
	      Max = Difference
	  ];
	  GuessP[i] = P[i]; 
      ];
  ];
  P[i_]:=0;
  Table[P[i],{i,0,x}]
]


Fend1[x_, n_, rho_] :=
Module[{m = Floor[x], lambda, beta, h, i, j, k, a, P0, f, F, pi, ans},
  
  lambda = n rho;
  beta = lambda Exp[-lambda];
  
  (* Initial values *)  
  
  For[i = 0, i <= (n-1), i++,
      pi[i] = 1/n;
  ];

  For[k = 0, k <= (n-1), k++,
      h=0;
      For[i = 0, i <= (n-1), i++,
	  If[i == k,
	      f[i,0] = 1,
	      f[i,0] = 0;
	  ];
	  For[j = 1, j <= i, j++,
	      f[i,j] = -lambda/j f[(i-1),(j-1)];
	  ];
      ];
      While[( h<10 ),
	  h = h+1;
	  For[i = (h n), i <= (h n + n - 1), i++,
	      f[i,0] = Sum[f[(i-n),j], {j,0,(i-n)}];
	      For[j = 1, j <= i, j++,
		  If [Mod[i, n] == 0,
		      f[i,j] = -beta/j f[(i-1),(j-1)],
		      f[i,j] = -lambda/j f[(i-1), (j-1)]
		  ];
	      ];
	      (* Print[f[i,0]," ", Exp[lambda h] f[i,0]]; *)
	  ];
      ];     
      For[i = 1, i <= n, i++,
	  F[i-1,k] = Exp[lambda h] f[(h n + i - 1), 0];    
      ];
  ]; (* End of k-loop *)

  P0 = Inverse[Table[F[i,j], {i,0,n-1}, {j,0,n-1}]] . Table[pi[i], {i,0,n-1}];
  (* Print[P0]; *)      
  
  (* End of initial values *)
  
  For[i = 0, i <= (m n + n - 1), i++,
      For[j = 0, j <= i, j++,
	  a[i,j] = 0;
      ];
  ];

  a[0,0] = Take[P0, 1];
  For[i = 1, i <= (n-1), i++,
      a[i,0] = Take[P0, {(i+1), (i+1)}];
      For[j = 1, j <= i, j++,
	  a[i,j] = -lambda/j a[(i-1),(j-1)];
      ];
  ];     

  For[i = n, i <= (m n + n - 1), i++,
      a[i,0] = Sum[a[(i-n),j], {j,0,(i-n)}];
      For[j = 1, j <= i, j++,
	  If [Mod[i, n] == 0,
	      a[i,j] = -beta/j a[(i-1),(j-1)],
	      a[i,j] = -lambda/j a[(i-1), (j-1)]
	  ];
      ];
  ];

  If[(x - m) == 0, 
      ans = n Exp[lambda x] a[(m n + n - 1), 0],
      ans = n Exp[lambda x] Sum[a[(m n + n - 1), j] (x-m)^j, {j, 0, (m n + n - 1)}];
  ]; 
  ans
]


(********** Routines for Network generation etc ****************)

(* Function: GabrielNet *)

GabrielNet[nbrNds_, options___] := 
  CNetCore[nbrNds, {Scale, Shape, CreateGraph, Index, Graph} /. {options} /. Options[GabrielNet]];

CNetCore[nbrNds_, {Scale_, Shape_, CreateGraph_, Index_, Graph_}] :=
  Module[{rr := Random[], p, p1, p2, d, points, P, dtemp, xsort, 
      ysort, indx, indy, xlim, ylim},
    
    If[TrueQ[Shape == Circle],
      p := (#1*Scale*{Cos[#2], Sin[#2]} & @@ {Sqrt[rr], 2N[Pi]*rr}), (* Shape : Circle *)
      p := Scale*{rr, rr}];      (* Shape : Box *)
    
    points = Table[p, {nbrNds}];
    d = Table[0, {nbrNds}, {nbrNds}];
    xsort = Ordering[points]; (* Index list of ordered points *)
    ysort = Ordering[points, All, #1[[2]] < #2[[2]] &];
    
    Do[p1 = points[[i]]; p2 = points[[j]];
      indx = ({Position[xsort, i][[1, 1]], Position[xsort, j][[1, 1]]} // Sort);
      indy = ({Position[ysort, i][[1, 1]], Position[ysort, j][[1, 1]]} // Sort);
      
      If[Intersection[Take[xsort, indx], Take[ysort, indy]] != {i, j},
        d[[i, j]] = d[[j, i]] = Infinity (* box defined by p1 and p2 is not empty *),
        
        P = (p1 + p2)/2; dtemp = Sqrt[(p1 - p2).(p1 - p2)];
        xlim = P[[1]] + {-dtemp, dtemp}/2;
        ylim = P[[2]] + {-dtemp, dtemp}/2;
 
 	(* Points that fit in the x - limits *)
        While[indx[[1]] > 1 && points[[xsort[[indx[[1]]]], 1]] > xlim[[1]], 
          indx[[1]] -= 1];
        While[indx[[2]] < nbrNds && points[[xsort[[indx[[2]]]], 1]] < xlim[[2]], 
          indx[[2]] += 1];

        (* Points that fit in the y - limits *)
        While[indy[[1]] > 1 && points[[ysort[[indy[[1]]]], 2]] > ylim[[1]], 
          indy[[1]] -= 1];
        While[indy[[2]] < nbrNds && points[[ysort[[indy[[2]]]], 2]] < ylim[[2]], 
          indy[[2]] += 1];

        (* (indices of) points that fit in the both limits *)
        P = Intersection[Take[xsort, indx], Take[ysort, indy]];
        P = Complement[P, {i, j}]; (* Remove points i and j *)

        (* Check the Gabriel - criteria for points P *)
        d[[i, j]] = d[[j, i]] = If[And @@ ((# - p1).(# - p2) >= 0 & /@ points[[P]]),
              (* connected points *) dtemp,
              (* not connected *) Infinity]], (* end if *)
      {i, 1, nbrNds}, {j, i + 1, nbrNds}];
    
    If[CreateGraph || Graph, (* Create the graph *)
    	{d, points, RCore[{d,points},{},{Index,Graph}]},
       (* do not create or return the graph *)
       	{d, points}]
    ]
    
    (*
    graph = {GrayLevel[.3], Line /@ graph};
    graph = {graph, GrayLevel[0], AbsolutePointSize[5], Point /@ points};

    If[Index, graph = {graph, RGBColor[0, 0.4, 0.2],
          MapIndexed[Text[First[#2], #1, {-2, 0}, Background -> GrayLevel[.95]] &, points]}];
    graph = Graphics[graph];
    If[Graph, Show[graph, AspectRatio -> Automatic, PlotRange -> All];];
    {d, points, graph}]*)

(* End: GabrielNet *)

(* Function: Dijkstra *)

Dijkstra[d_, root_, dstn_:False] :=
  Module[{core, prev, dist, cond, temp, new, x, y, interface, closest, mindist, newnode, from, pos},

    (* set the given root node to the interface *)
    core = interface = {root};(* (indices of) core nodes still having possible neighbourgs *)
    prev = {root};
    peri = Delete[Range[pos = (d // Length)], root]; (* NOTA BENE! here pos used as a temp variable *)
    dist = {0};
    
    (* cases in which the distance/cost is Infinity are excluded *)
    temp = Flatten[Position[#, Infinity]] & /@ d;
    pos = Range[pos]; (* NB! here pos used as a temp variable *);
    closest = (Complement[pos, #] & /@ temp);
    
    (* sort the (indices of) closest nodes for each node *)
    closest = MapIndexed[#1[[Rest[Ordering[d[[#2, #1]] // Flatten]]]] &, closest];

    (* remove cases in which the root is the closest neighbour *)
    closest = (Delete[#, Position[#, root]] & /@ closest);
    If[closest[[root]]=={},interface={}]; (* case that no connection from the root *)

    If[! IntegerQ[dstn],
      (* no destination *)cond := (interface =!= {}),
      (* destination dstn *)cond := ((dstn != Last[core])&&(interface =!= {}))];
    
    While[cond, (* loop is repeated as long as cond is TRUE *)
      (* find the closest node and the corresponding distance *)
      mindist = Infinity;
      MapIndexed[If[#1 < mindist, mindist = #1; from = #2] &,  (* NB! here from used as an index *)
        Flatten[dist[[Position[core, #] // First]] + d[[#, closest[[#, 1]]]] & /@ interface]];
      
      (* the new node of core is now known *)
      newnode = closest[[from = interface[[from // First]], 1]]; (* NB! here from denotes the previous node *)

      (* update the lists *)
      core = {core, newnode} // Flatten;
      interface = {interface, newnode} // Flatten;
      prev = {prev, from} // Flatten;
      dist = {dist, mindist} // Flatten;
      peri = Delete[peri, Position[peri, newnode] // Flatten];

      (* remove newnode from lists in closest and update it *)
      closest = (Delete[#, Position[#, newnode] // Flatten] & /@ closest);

      (* remove nodes that do not have neighbourghs from list interface *)
      interface = Part[interface, Complement[Range[Length[interface]],
            Flatten[Position[interface, #] & /@ Flatten[Position[closest, {}]]]]]
      ]; (* end While loop *)

    If[IntegerQ[dstn], (* one destination dstn *)
    
      If[(y = Position[core, dstn] // Flatten)=!={},

      (* destination is found *)
      x = Position[core, prev // Last] // First;
      While[First[y] != 1, y = {x, y} // Flatten;
          x = Position[core, prev[[x]] // First] // First];
      {core, prev, dist} = {core[[y]], prev[[y]], dist[[y]]};,
      {core, prev, dist} = {{root},{root},{0}};(* destination not found *)];
      ];
    (* returns following lists *)
    {core, prev, dist}
    ]

(* End: Dijkstra *)

(* Function: ShowNetwork *)

ShowNetwork[GraphData_, RouteData_, options___] := 
  RCore[GraphData, RouteData, {Index, Graph} /. {options} /. Options[ShowNetwork]];

RCore[GraphData_, RouteData_, {Index_, Graph_}] := Module[{net={}, graph, root},
    If[Length[GraphData] == 2, (* no graph created *)
       Do[If[GraphData[[1]][[i,j]] != Infinity, (* points i and j connected *)
             net = net~Join~{{GraphData[[2,i]], GraphData[[2,j]]}}],
          {i,1,Length[GraphData[[1]]]},{j,i+1,Length[GraphData[[1]]]}];
       net = {GrayLevel[.3], Line /@ net};
       net = {net, GrayLevel[0], AbsolutePointSize[5], Point /@ GraphData[[2]]};
       If[Index, (* add indices *)
          net = {net, RGBColor[0, 0.4, 0.2],
          MapIndexed[Text[First[#2], #1, {-2, 0}, Background -> GrayLevel[.95]]&, GraphData[[2]]]}];
       net = Graphics[net],
       net = GraphData[[3]] (* else *)];
    
    If[RouteData=!={}, (* path info is known *)
      graph = Line[GraphData[[2, #]]] & /@ Transpose[Take[RouteData, 2]];
      root = {PointSize[.035], RGBColor[1, .5, .1], Point[GraphData[[2, RouteData[[1, 1]]]]]};
      graph = {root, {AbsoluteThickness[2], graph}};
      If[Index, graph = {graph, RGBColor[.3, .1, 0], 
         Text[#, GraphData[[2, #]], {-2, 0}, Background -> GrayLevel[.95]] & /@ RouteData[[1]]}];
      If[Length[GraphData]==3,graph = {GraphData[[3]], Graphics[graph]},graph = {net, Graphics[graph]}],
      graph=net (* else *)];
      
    If[Graph, Show[graph, AspectRatio -> Automatic, PlotRange -> All];];
    graph];

(* End: ShowNetwork *)

(********** Routines for Trees, Balanced Fairness etc ****************)

(* Implementation of KSubsets taken from packet DiscreteMath`Combinatorica` *)

KSubsets[l_List, 0] := {{}}
KSubsets[l_List, 1] := Partition[l, 1]
KSubsets[l_List, k_Integer?Positive] := {l} /; (k == Length[l])
KSubsets[l_List, k_Integer?Positive] := {} /; (k > Length[l])
KSubsets[l_List, k_Integer?Positive] := Join[Map[(Prepend[#, First[l]]) &, KSubsets[Rest[l], k - 1]],
       KSubsets[Rest[l], k]]

(* Function: NormalizationConstant *)

NormalizationConstant[TreeList_] := 
  Module[{tree, uses, offered, saturated, recursion, flows, nfl, isets, rt, tt, g, G, a}, 
    (* Definitions *) 
  tree[n_] := {Flatten[{n, (tree[#][[1]] & /@ Flatten[ Position[TreeList, {n, __}]])}], TreeList[[n]][[2]],
        tree /@ Flatten[ Position[TreeList, {n, __}]]} ;
  offered[t_, ind_] :=  Min[t[[2]], 
        If[MemberQ[ind, t[[1, 1]] ], Infinity , 0] + 
          Plus @@ (offered[#, ind] & /@ Last[t])];
  saturated[{}, ind_] := {Null, Null, False};
  saturated[t_, ind_] :=  If[offered[t, ind] == t[[2]], {Intersection[ind, t[[1]]], t[[2]], True},
         If[(i =  Position[Last[saturated[#, ind]] & /@ Last[t],  True]) == {}, {Null, Null, False},
  saturated[Last[t][[i // First // First]], ind] ]];
  recursion[ind_] := Module[{sat, ii, iii, r, gg, c},
        {ii, c} = Drop[saturated[tt, ind], -1];
        iii = Complement[ind, {#}] & /@ ii;
        gg = g /@ iii; 
        hh = h /@ iii;
        r = Last /@ TreeList[[ii]] ;
        c -= Plus @@ r;
        gg = r.gg/c ];
    (* Initializations *);
  flows = Position[TreeList, {_, _, _}] // Flatten; nfl = Length[flows];
  isets = Table[KSubsets[flows, k], {k, nfl}] // Flatten[#, 1] &;
  rt = Position[TreeList, {Null, __}] // First // First;
    tt = tree[rt];
    (* Actual recursion *)
    G = g[{}] = 1;
    G += (g[#] = recursion[#]) & /@ isets;
  G // Simplify]

(* End: NormalizationConstant *)

(* Function: Throughputs *)

Throughputs[TreeList_] := 
  Module[{tree, uses, offered, saturated, recursion, flows, nfl, isets, rt, tt, g, h, G, H, i}, 
    (* Definitions *) 
  tree[n_] := {Flatten[{n, (tree[#][[1]] & /@ 
                Flatten[ Position[TreeList, {n, __}]])}], TreeList[[n]][[2]],
        tree /@ Flatten[ Position[TreeList, {n, __}]]} ;
  offered[t_, ind_] :=  
      Min[t[[2]], 
        If[MemberQ[ind, t[[1, 1]] ], Infinity , 0] + 
          Plus @@ (offered[#, ind] & /@ Last[t])];
  saturated[{}, ind_] := {Null, Null, False};
  saturated[t_, ind_] := 
      If[offered[t, ind] == t[[2]], {Intersection[ind, t[[1]]], t[[2]], True},
         If[(i = 
                Position[Last[saturated[#, ind]] & /@ Last[t], 
                  True]) == {}, {Null, Null, False},
  saturated[Last[t][[i // First // First]], ind] ]];
  recursion[ind_] := Module[{sat, ii, iii, r, gg, c},
        {ii, c} = Drop[saturated[tt, ind], -1];
        iii = Complement[ind, {#}] & /@ ii;
        gg = g /@ iii; 
        hh = h /@ iii;
        r = Last /@ TreeList[[ii]] ;
        c -= Plus @@ r;
        gg = r.gg/c;
        hh = (If[MemberQ[ii, #], gg + g[Complement[ind, {#}]], 0] & /@ flows +
                 r.hh)/c;
        {gg, hh} ];
    (* Initializations *);
  flows = Position[TreeList, {_, _, _}] // Flatten; nfl = Length[flows];
  isets = Table[KSubsets[flows, k], {k, nfl}] // Flatten[#, 1] &;
  rt = Position[TreeList, {Null, __}] // First // First;
    tt = tree[rt];
    (* Actual recursion *)
    {G, H} = {g[{}], h[{}]} = {1, Table[0, {i, nfl}]};
    {G, H} += ({g[#], h[#]} = recursion[#]) & /@ isets;
  {G, (G/H), flows} // Simplify]

(* End: Throughputs *)

(* Function: TreeGraph *)

TreeGraph[TreeList_] := Module[{link, tree, tg, rt, tt, d = 0.6},
  link[b_, t_] := Module[{a = 0.125, s, p} ,
        s = {-1, -a/2} + b;
        p = {{0, 0}, {1, 0}, {1, a}, {0, a}, {0, 0}} + {s, s, s, s, s};
        {{Hue[0.5], Polygon[p]}, Line[p], 
          Text[t, {1/2, 2a} + s], {GrayLevel[1], Point[{0, 3a} + s]}}];
  tree[n_] := {{Max[d, 
              If[Length[TreeList[[n]]] == 3, d, 0] + 
                Plus @@ (tree[#][[1, 1]] & /@ 
                      Flatten[ Position[TreeList, {n, __}]])], n, 
  TreeList[[n]]} // Flatten,
        tree /@ Flatten[ Position[TreeList, {n, __}]]} ;
  tg[t_, b_] := 
      Module[{bb = b, bbb, a = t[[1, 1]], aa, tt = Last[t], lst = {}},
        {link[bb, 
  StringJoin["C", First[t][[2]] // ToString, "=", 
  First[t][[4]] // ToString]],
          bb -= {2, a/2}; aa = If[tt == {}, 0, tt[[1, 1, 1]] ]/2;
          tg[#, aa = #[[1, 1]]/2; bbb = bb += {0, aa}; 
                lst = {lst, Line[{b - {1, 0}, bb}]}; bb += {0, aa}; bbb] & /@ 
            tt, lst, If[Length[t[[1]]] == 5,
            
            bb += {0, d/2}; {Line[{b - {1, 0}, bb}], 
  Text[StringJoin["x", First[t][[2]] // ToString], 
                bb - {.2, 0}]}, {}]}];
  rt = Position[TreeList, {Null, __}] // First // First;
    tt = tree[rt];
    {tg[tt, {0, 0}], Line[{{0, 0}, {1, 0}}]}]

(* End: TreeGraph *)

(**************** *********************************************)

End[]

EndPackage[]
