restart:
 with(linalg):

# Reads in the data
 
 read matrix1: 
 N := %:
 read PDF1:
 dist := %:
 nodes := rowdim(N):
 arcs := coldim(N):

# Puts matrix in usable form
 
 insert := 1:
 for row from 1 to nodes do
   for col from 1 to arcs do
     if N[row, col] = 1 then
       N := swapcol(N, col, insert):
       insert := insert + 1:
     fi: 
   od: 
 od:
 top := 1: 
 insert := 1:
 while top <= nodes do
   for row from 1 to nodes do
     for col from 1 to arcs do
       if N[row, col] = -1 and N[top, col] = 1 then
         N := swapcol(N, col, insert):
         insert := insert + 1:
       fi:
     od:
   od:
   top:= top + 1:
 od:

# Initializes the path matrix
 
 Paths := array(1..50, 1..nodes):
 for i from 1 to 50 do
   for j from 1 to nodes do
     Paths[i, j] := 0:
   od: 
 od:

# Finds number of incoming arcs for each node and the source
# node for each edge.
 
 source := array(1..arcs):
 for i from 1 to arcs do 
   source[i] := 0: 
 od:
 incoming := array(1..nodes):
 for i from 1 to nodes do 
   incoming[i] := 0: 
 od: 
 for i from 1 to nodes do
   for j from 1 to arcs do
     if N[i, j] = -1 then
       incoming[i] := incoming[i] + 1 else
       if N[i, j] = 1 then
         source[j] := i: 
       fi:
     fi:
   od:
 od:

# Finds the paths through the network, Step 1 of Algorithm.
 
 GetPaths := proc(node, step, path, Paths)
 local i, j, found, numpaths, total:
 i := 1:
 found := 0:
 numpaths := 0:
 total := 0:
 while (found - incoming[node] < 0) do
   if N[node,i] = -1 then
     numpaths := GetPaths(source[i], step + 1, path + total, Paths):
       for j from 0 to numpaths-1 do
         Paths[path + j + total, step] := i:
       od:
     total := total + numpaths:
     found := found + 1:
   fi:
   i := i + 1:
 od:
 if total = 0 then
   Paths[path, step] := 0:
   total := 1:
 fi:
 RETURN(total):
 end:
 
 numPaths := GetPaths(nodes, 1, 1, Paths):
 Paths:=delrows(Paths, numPaths+1..50):

# The previous procedure finds the paths, but they are in the 
# wrong order, the next part reverses the order.
 
 tempPaths := array(1..numPaths, 1..nodes):
 for i from 1 to numPaths do
   for j from 1 to nodes do
     tempPaths[i, j] := 0:
   od: 
 od:
 for i from 1 to numPaths do
   count := 0:
   for j from 1 to nodes do 
     if Paths[i, j] > 0 then 
       count := count + 1:
     fi:
   od:
   for j from 1 to count do
     tempPaths[i, j] := Paths[i, count + 1 - j]:
   od: 
 od:
 Paths := tempPaths:

# Distinguishes edges  found on multiple paths from those found on 
# only one, Step 2 of Algorithm.

 S := array(1..arcs):
 for i from 1 to arcs do
   S[i] := 0:
 od:
 for i from 1 to numPaths do
   for j from 1 to nodes do
     if Paths[i, j] > 0 then
       S[Paths[i, j]] := S[Paths[i, j]] + 1:
     fi:
   od:
 od:
 for k from 1 to arcs do
   if S[k] > 1 then 
     S[k] := 1:
   else 
     S[k] := 0:
   fi:
 od:

# Calculates k, the number of multiple-use edges, Step  3 of Algorithm

 k := sum(S['n'], 'n' = 1..arcs):

# Finds the order in which the integration of each edge must occur, 
# combination of Steps 4 and 5 of Algorithm

 intOrder := array(1..k, []):
 marker := array(1..numPaths): 
 for i from 1 to numPaths do 
   marker[i] := 1
 od: 
 counter := 0: 
 for i from 1 to k do
   intOrder[i] := 0
 od:
 while intOrder[k] = 0 do
   first := {}:
   later := {}:
   for i from 1 to numPaths do
     if Paths[i, marker[i]] > 0 then
       if S[Paths[i, marker[i]]] = 1 then
         first := first union {Paths[i, marker[i]]}:
       fi:
     fi:
   od:
   for i from 1 to numPaths do
     for j from marker[i] + 1 while Paths[i, j] > 0 do
       if S[Paths[i, j]] = 1 then
         later := later union {Paths[i, j]}:
       fi:
     od:
   od:
   first:= first minus (first intersect later):
   for i from 1 to nops(first) do
     counter:= counter + 1:
     intOrder[counter] := first[i]:
   od:
   for i from 1 to numPaths do
     if Paths[i, marker[i]] > 0 then
       if {Paths[i, marker[i]]} intersect first = {Paths[i,marker[i]]} or
          S[Paths[i,marker[i]]] = 0 then
         marker[i] := marker[i] + 1:
       fi: 
     fi:
   od:
 od:

# Sets the lower limits of integration and starts calculating the
# upper limits of integration, Step 6a of Algorithm

 numIntegrals := 1: 
 upper1 := array(1..4, 1..k): 
 lower1 := array(1..3, 1..k):
 for i from 1 to k do
   upper1[1, i] := t:
   for j from 2 to 4 do
     upper1[j, i] := 0:  
     lower1[j-1, i] := 0:
   od:
 od:

# Finishes calculating the upper limits of integration and simplifies 
# the upper limits, Steps 6b and 7 of Algorithm

 for i from 1 to k do 
   first := {}:
   second := {}: 
   for j from 1 to numPaths do
     for m from 2 to nodes do
       if Paths[j, m] = intOrder[i] then
         temp := {}:
         for l from m-1 to 1 by -1 do
           if S[Paths[j, l]] = 1 then
             for n from 1 to k do
               if intOrder[n] = Paths[j, l] then
                 p := n:
               fi:
             od:
             temp := temp union {x[p]}:
           fi:
         od:
         if nops(first) = 0 then
           first := temp:
         else
           second := temp:
         fi:
       fi:
     od:
   od:
   if second = {} then
     upper1[2, i] := first:
     upper1[3, i] := {}:
     upper1[4, i] := {}:
   else
     temp := first intersect second:
     first := first minus temp:
     second := second minus temp:
     upper1[2, i] := temp:
     upper1[3, i] := first:
     upper1[4, i] := second:
   fi:
 od:
 for i from 1 to 4 do
   for j from 1 to k do
     upper1[i, j] := convert(upper1[i, j], `+`):
   od:
 od:

# Eliminates the maximum expressions in the limits of integration
# (upper limits first, then lower limits), Steps 8 and 9 of Algorithm

 finished:= 1:
 while finished > 0 do
   for i from 1 to numIntegrals do
     for j from 1 to k do
       if evalb(upper||i[3, j] = 0) = false or
          evalb(upper||i[4, j] = 0) = false then
         numIntegrals := numIntegrals + 1:
         temp:= convert(convert([upper||i[3, j], upper||i[4, j]],
            `+`), set):
         temp1 := {}:
         for l from 1 to k do
           if evalb(temp intersect {x[l]} = {x[l]}) then
             temp1 := temp1 union {l}: 
           fi: 
         od: 
         temp1:= sort(temp1):
         lim := solve (upper||i[3, j] = upper||i[4, j], x[temp1[-1]]):
         a:=upper||i[3, j]: 
         b:=upper||i[4, j]:
         upper||numIntegrals := copy(upper||i): 
         lower||numIntegrals := copy(lower||i):
         for m from j to k do
           if upper||i[3, m] = a and upper||i[4, m] = b then
             if evalb(upper||i[1, m] = 0) = false then
               upper||i[2, m] := upper||i[2, m] + a: 
               upper||i[3, m] := 0: 
               upper||i[4, m] := 0:
               upper||numIntegrals[2, m] := upper||numIntegrals[2, m] + b: 
               upper||numIntegrals[3, m] := 0: 
               upper||numIntegrals[4, m] := 0:
             else  
               upper||i[1, m] := a: 
               upper||i[3, m] := 0: 
               upper||i[4, m] := 0:
               upper||numIntegrals[1, m] := b:
               upper||numIntegrals[3, m] := 0: 
               upper||numIntegrals[4, m] := 0:
             fi:
           fi: 
         od:
         a1:= convert(convert(a, `+`), set):
         if nops(lim) > 1 then
           if a1 intersect {x[temp1[-1]]} = {x[temp1[-1]]} or a1 intersect
              {-x[temp1[-1]]} = {-x[temp1[-1]]} then
             lower||i[1, temp1[-1]] := 0: 
             lower||i[2, temp1[-1]] := 0: 
             lower||i[3, temp1[-1]] := lim:
             upper||numIntegrals[1, temp1[-1]] := 0:
             upper||numIntegrals[2, temp1[-1]] := 0:
             upper||numIntegrals[3, temp1[-1]] := lim:
             upper||numIntegrals[4, temp1[-1]] := 0:
           else
             lower||numIntegrals[1, temp1[-1]] := 0:
             lower||numIntegrals[2, temp1[-1]] := lim:
             lower||numIntegrals[3, temp1[-1]] := 0:
             upper||i[1, temp1[-1]] := 0: 
             upper||i[2, temp1[-1]] := 0:
             upper||i[3, temp1[-1]] := lim: 
             upper||i[4, temp1[-1]] := 0:
           fi: 
         else 
           if a1 intersect {x[temp1[-1]]} = {x[temp1[-1]]} or a1 intersect
              {-x[temp1[-1]]} = {-x[temp1[-1]]} then
             lower||i[1, temp1[-1]] := lim:
             lower||i[2, temp1[-1]] := 0:
             lower||i[3, temp1[-1]] := 0:
             upper||numIntegrals[1, temp1[-1]] := lim:
             upper||numIntegrals[2, temp1[-1]] := 0:
             upper||numIntegrals[3, temp1[-1]] := 0: 
             upper||numIntegrals[4, temp1[-1]] := 0:
           else
             lower||numIntegrals[1, temp1[-1]] := lim:
             lower||numIntegrals[2, temp1[-1]] := 0:
             lower||numIntegrals[3, temp1[-1]] := 0:
             upper||i[1, temp1[-1]] := lim: 
             upper||i[2, temp1[-1]] := 0: 
             upper||i[3, temp1[-1]] := 0:
             upper||i[4, temp1[-1]] := 0:
           fi:
         fi:
       fi:
     od:
     for j from 1 to k do
       if evalb(lower||i[2, j] = 0) = false or evalb(lower||i[3, j] =
            0) = false then
         numIntegrals := numIntegrals + 1:
         temp:= convert(convert([lower||i[2, j], lower||i[3, j]],
            `+`), set):
         temp1 := {}:
         for l from 1 to k do
           if evalb(temp intersect {x[l]} = {x[l]}) then
             temp1 := temp1 union {l}:
           fi:
         od:
         lim := solve (lower||i[2, j] = lower||i[3, j], x[temp1[-1]]):
         a := lower||i[2, j]:
         b := lower||i[3, j]:
         upper||numIntegrals := copy(upper||i):
         lower||numIntegrals := copy(lower||i):
         for m from j to k do
           if lower||i[2, m] = a and lower||i[3, m] = b then
             lower||i[1, m] := lower||i[1, m] + a:
             lower||i[2, m] := 0:
             lower||i[3, m] := 0:
             lower||numIntegrals[1, m] := lower||numIntegrals[1, m] + b:
             lower||numIntegrals[2, m] := 0:
             lower||numIntegrals[3, m] := 0:
           fi:
         od:
         a1:= convert(convert(a, `+`), set): 
         if nops(lim) > 1 then
           if a1 intersect {x[temp1[-1]]} = {x[temp1[-1]]} or a1 intersect
              {-x[temp1[-1]]} = {-x[temp1[-1]]} then
             lower||i[1, temp1[-1]] := 0:
             lower||i[2, temp1[-1]] := 0: 
             lower||i[3, temp1[-1]] := lim:
             upper||numIntegrals[1, temp1[-1]] := 0:
             upper||numIntegrals[2, temp1[-1]] := 0:
             upper||numIntegrals[3, temp1[-1]] := lim:
             upper||numIntegrals[4, temp1[-1]] := 0:
           else
             lower||numIntegrals[1, temp1[-1]] := 0:
             lower||numIntegrals[2, temp1[-1]] := lim:
             lower||numIntegrals[3, temp1[-1]] := 0:
             upper||i[1, temp1[-1]] := 0: 
             upper||i[2, temp1[-1]] := 0:
             upper||i[3, temp1[-1]] := lim:
             upper||i[4, temp1[-1]] := 0:
           fi: 
         else 
           if a1 intersect {x[temp1[-1]]} = {x[temp1[-1]]} or a1 intersect
              {-x[temp1[-1]]} = {-x[temp1[-1]]} then
             lower||i[1, temp1[-1]] := lim:
             lower||i[2, temp1[-1]] := 0: 
             lower||i[3, temp1[-1]] := 0:
             upper||numIntegrals[1, temp1[-1]] := lim:
             upper||numIntegrals[2, temp1[-1]] := 0:
             upper||numIntegrals[3, temp1[-1]] := 0:
             upper||numIntegrals[4, temp1[-1]] := 0:
           else
             lower||numIntegrals[1, temp1[-1]] := lim:
             lower||numIntegrals[2, temp1[-1]] := 0: 
             lower||numIntegrals[3, temp1[-1]] := 0:
             upper||i[1, temp1[-1]] := lim:
             upper||i[2, temp1[-1]] := 0:
             upper||i[3, temp1[-1]] := 0:
             upper||i[4, temp1[-1]] := 0:
           fi:
         fi:
       fi:
     od:
     finished := 0:
     for g from 1 to numIntegrals do
       for j from 1 to k do
         if evalb(lower||g[2, j] = 0) = false or 
            evalb(lower||g[3, j] = 0) = false or
            evalb(upper||g[3, j] = 0) = false or 
            evalb(upper||g[4, j] = 0) = false then
           finished := 1: 
         fi:
       od:
     od:
   od:
 od:
   
# Constructs the integrand, Step 10 of Algorithm

 for i from 1 to k do
   dist[intOrder[i]] := subs(x = x[i], dist[intOrder[i]]): 
 od:
 for i from 1 to arcs do
   if S[i] = 0 then
     temp := {}:
     for j from 1 to numPaths do
       for m from 1 to nodes do
         if Paths[j, m] = i then
           for l from 1 to nodes do
             if Paths[j, l] > 0 and S[Paths[j, l]] = 1 then
               for n from 1 to k do
                 if intOrder[n]= Paths[j, l] then
                   temp := temp union {x[n]}:
                 fi:
               od: 
             fi: 
           od:
         fi:
       od:
     od:
     temp := convert(temp, `+`):
     temp1 := t - temp:
     dist[i] := int(dist[i], x = 0..z):
     dist[i] := subs(z = temp1, dist[i]):
   fi:
 od:
 integrand := 1:
 for i from 1 to arcs do
   integrand := integrand * dist[i]:
 od:

# Sets up and solves the integrals to find the solution,
# Step 11 of algorithm

 solution := 0:
 for i from 1 to numIntegrals do
   temp := integrand:
   for j from k to 1 by -1 do
     temp := int(temp, x[j] = lower||i[1, j]..upper||i[1, j] -
        upper||i[2, j]):
   od:
   solution := solution + temp:
 od:
 simplify(solution);
 mean := int(1 - solution, t = 0 .. infinity);