# Define the method for Psi- evaluation,
Psi_minus:=proc(OpSet)
           local i, solution;
           solution:={};
           for i from 1 to nops(OpSet) do
             solution:=solution union op(1,op(i,OpSet))
           od;
           solution;
end:

# Define the method for Psi+ evaluation,
Psi_plus:=proc(OpSet)
          local i, solution;
          solution:={};
          for i from 1 to nops(OpSet) do
            solution:=solution union op(2,op(i,OpSet))
          od;
          solution;
end:
# Define the method for phi- evaluation,
phi_minus:=proc(mSet,OpSet)
           local i, solution, item;
           solution:={};
           for i from 1 to nops(OpSet) do
               item:={op(i,OpSet)};
               if (op(2,op(1,item)) intersect mSet) <> {} then
                  solution:=solution union item  fi;
           od;
           solution;
end:

# Define the method for phi+ evaluation,
phi_plus:=proc(mSet,OpSet)
          local i, solution, item;
          solution:={};
          for i from 1 to nops(OpSet) do
            item:={op(i,OpSet)};
            if (op(1,op(1,item)) intersect mSet) <> {} then
               solution:=solution union item  fi;
          od;
          solution;
end:


# Program MSG algorithm,
MSG:=proc(OpSet,RM,Pr)
     local i, item, opset, Mat, o, r, p, m, ox;
     opset:=OpSet minus phi_minus(RM,OpSet);
     Mat:=Psi_minus(opset) union Psi_plus(opset);
     r:=Psi_minus(opset) minus (Psi_plus(opset) union RM);
     if nops(r) > 0 then 
        while r <> {}  do
              item:={op(1,r)};
              Mat:=Mat minus item;
              o:=phi_plus(item,opset);
              opset:=opset minus o;
              r:=(r union (Psi_plus(o) minus Psi_plus(opset))) minus item;
        od;
     fi;
     if (Pr union Mat) <> Mat then
        lprint(`Error:  No maximal structure exists!`); 
        RETURN(`No maximal structure exists!`);
     fi;
     p:=Pr;
     m:={};
     o:={};
     while p <> {}  do
        item:={op(1,p)};
        m:=m union item;
        ox:=phi_minus(item,opset);
        o:=o union ox;
        p:=(p union Psi_minus(ox)) minus (RM union m);
     od;
     m:=Psi_minus(o) union Psi_plus(o);
     [m,o];
end:

# Convert the Pgraph structure into a network,
graph_network:=proc(Structure,naming)
   local i, j, rm, pr, mid, m, mat_order, listing, opset, Gr,
         mat_ins, mat_outs, el, unitnum, inset, outset, inset2, outset2,
         unit_order, unit_order2, unitname, location;
   global  Raw_Materials, Products, OpData, LabelsList;
   m:=op(1,Structure);
   opset:=op(2,Structure);
   # check for material inlets and outlets of the structure
   mat_ins:=Psi_minus(opset);
   mat_outs:=Psi_plus(opset);
   # find all raw materials (they are only inlets to unit ops)
   rm:={};
   # find all products (they are only outlets of unit ops)
   pr:={};
   for el in m  do
       if member(el,mat_ins) and member(el,mat_outs) = false then
          rm:=rm union {el}
       fi;
       if member(el,mat_ins) = false and member(el,mat_outs) then
          pr:=pr union {el}
       fi;
   od;
   mid:=m minus rm minus pr;
   mat_order:=[op(pr)],[op(mid)],[op(rm)];
   # get labels from OpData information
   j:=0;
   for el in opset do
       j:=j+1;
       inset:=op(1,el);
       outset:=op(2,el);
       for i from 1 to nops(OpData) do
           inset2:={op(op(2,op(i,OpData)))};
           outset2:={op(op(3,op(i,OpData)))};
           if inset = inset2 and outset = outset2 then break fi;
       od;
       unitname:=op(1,op(i,OpData));
       member(unitname,LabelsList,'location');
       unitnum[j]:=location;
   od;
   unit_order:={};
   if (naming = 'numbers') then
      unit_order:=[seq(unitnum[j],j=1..nops(opset))];
   else
      unit_order:=[seq(cat(`:`,op(unitnum[j],LabelsList)),
                   j=1..nops(opset))]
   fi;
   unit_order2:=convert(convert(unit_order,set),list);
   listing:=mat_order,unit_order2;
   new(Gr);
   addvertex(op(1,Structure),Gr);
   addvertex(unit_order,Gr);
   for i from 1 to nops(opset) do
       connect(op(1,op(i,opset)),[op(i,unit_order)],directed,Gr);
       connect([op(i,unit_order)],op(2,op(i,opset)),directed,Gr);
   od;
   Linear(listing),Gr;
end:

# Program SSG algorithm,
SSG:=proc(by_prod)
	local mats, e, maps;
	global ox_vals, decision_map, by_prods, Op, Products;
        by_prods:=by_prod;
        mats:=Psi_minus(Op) union Psi_plus(Op);
	ox_vals:=table();
        for e in mats do
            ox_vals[{e}]:=phi_minus({e},Op);
        od;
	SSGmethod(Products,{},{});
        maps:=nops(decision_map);
	`# of decision maps found = `.maps
end:
SSGmethod:=proc(p,m,delta_m)
   local i, j, k, item_x, set_c, testing, item_c, item_y, ox, oy, delta_y, 
         delta_m2, p2, m2, op_y, mat_outs, mat_ins, opmat, pr, el;
   global Raw_Materials, Products, Op, decision_map, ox_vals, by_prods;
   if nops(p) = 0 then
     if by_prods = true then
       # by-products are allowed; accept the solution
       decision_map:=decision_map union {delta_m}
     elif by_prods = false then
       # accept a solution if no by-products are present
       # check for material inlets and outlets of the decision map
       opmat:={seq(op(op(2,op(k,delta_m))),k=1..nops(delta_m))};
       mat_ins:=Psi_minus(opmat);
       mat_outs:=Psi_plus(opmat);
       # subtract the raw materials and products
       mat_ins:=mat_ins minus Raw_Materials minus Products;
       mat_outs:=mat_outs minus Raw_Materials minus Products;
       # for a valid mapping, the intermediate materials must be both
       # inlet and outlet materials
       if ( (mat_ins minus mat_outs) = {} and
	    (mat_outs minus mat_ins) = {} ) then 
          decision_map:=decision_map union {delta_m};
       fi;
     else
       # accept the solution if the by-products are within the input set
       # of allowed by-products.
       # check for material inlets and outlets of the structure
       opmat:={seq(op(op(2,op(k,delta_m))),k=1..nops(delta_m))};
       mat_ins:=Psi_minus(opmat);
       mat_outs:=Psi_plus(opmat);
       # find all products (they are only outlets of unit ops)
       pr:={};
       for el in (mat_ins union mat_outs) do
          if member(el,mat_ins) = false and member(el,mat_outs) then
             pr:=pr union {el}
          fi;
       od;
       # test that the by-products are contained in the by_prods set.
       # accept the solution if true.
       if ((pr minus Products) minus by_prods) = {} then
          decision_map:=decision_map union {delta_m} fi;
     fi;
     if (nops(decision_map) mod 100) = 0 and 
         nops(decision_map) > 0              then
        lprint(nops(decision_map), ` decision maps found so far`)
     fi;
     RETURN(0);
   fi;
   item_x:={op(1,p)};
   ox:=ox_vals[item_x]; 
   set_c:=choose(ox) minus {{}}; 
   for i from 1 to nops(set_c) do
       item_c:=op(i,set_c);
       testing:=0;
       for j from 1 to nops(m) do
           item_y:={op(j,m)};
           oy:=ox_vals[item_y];
           delta_y:={};
           if nops(delta_m) > 0 then 
              op_y:={seq(op(op(2,op(k,delta_m))),k=1..nops(delta_m))};
              delta_y:=phi_minus(item_y,op_y);
           fi;
           if (item_c intersect (oy minus delta_y)) <> {}  or   
              ((ox minus item_c) intersect delta_y) <> {}  then
                testing:=1;
                break;
           fi;
       od;
       if testing = 1  then  next  fi; 
       delta_m2:=delta_m union {[item_x,item_c]};
       p2:=p;
       for j from 1 to nops(item_c) do
           p2:=p2 union op(1,op(j,item_c))
       od;
       p2:=p2 minus (Raw_Materials union m union item_x);
       m2:=m union item_x;
       SSGmethod(p2, m2, delta_m2);
   od;
end:

# Convert a decision mapping into a P graph,
Pgraph:=proc(delta_m)
        local k, m, o;
        o:={seq(op(op(2,op(k,delta_m))),k=1..nops(delta_m))};
        m:=Psi_minus(o) union Psi_plus(o);
        [m,o];
end:

# updating the draw package to show different colors on the graph
`draw/Linear`:=
proc(partitions::specfunc(list,Linear),G::GRAPH,
    origin::list(numeric),xrng::name,yrng::name)
local e,max_x,middle,n,nx,part,shift,t1,text,pos,x,y,
    t,vset,lines,points,i,j,v,center, 
    a,b,c,d,a1,a2,a3,angle,dist; #** add new local vars
options `Copyright 1993 by Waterloo Maple Software`;
    center := [origin[1 .. 2]];
    vset := {};
    n := 0;
    part := table();
    shift := table();
    for t in partitions do
        n := n+1;
        if not type(t,list) then
            ERROR(`partition should be a list`)
        fi;
        shift[n] := 0;
        t1 := select(
            hastype,t,'identical'('offset') = numeric)
            ;
        if t1 <> [] then
            shift[n] := subs({op(t1)},'offset');
            t :=
              select(proc(x) not hastype(x,`=`) end,t)
        else
            if {op(t)} intersect vset <> {} then
                ERROR(
                 `intersecting partitions involving`,t
                 )
            fi
        fi;
        if not type(t,'list'('VERTEX'(G))) then
            ERROR(`not a list of vertices`,t)
        fi;
        part[n] := t;
        vset := vset union {op(t)}
    od;
    userinfo(
        5,'networks',`allocated vertices:`,print(vset)
        );
    if networks['vertices'](G) minus vset <> {} then
        n := n+1;
        t :=
        sort([op(networks['vertices'](G) minus vset)])
        ;
        part[n] := t;
        shift[n] := 0;
        userinfo(
            1,'networks',`adding a new partition`,t)
    fi;
    max_x := origin[1];
    for j to n do
        if max_x < nops(part[j]) then
            max_x := nops(part[j])
        fi
    od;
    middle := 1/2*max_x-1/2+origin[1];
    points := table();
    text := table();
    pos := table();
    x := origin[1];
    for j to n do
        vset := part[j];
        nx := nops(vset);
        x := x+1;
        #** update to shift the alignment of columns
        y := middle-1/2*nx+1/2+shift[j]+origin[2]+0.2*((j mod 3) - 1);
        for i to nops(vset) do
            v := vset[i];
            pos[v] := [x,y];
            points[v] := POINTS(pos[v]);
            if j = n then
                text[v] := 'TEXT'(
                    [pos[v][1]+1.6*10^(-1),pos[v][2]],
                    convert(v,string))
            else
                text[v] := 'TEXT'(
                    [pos[v][1]-.1 ,pos[v][2]],
                    convert(v,string))
            fi;
            y := y+1
        od
    od;
    lines := table();
    for e in edges(G) do
        x := networks['ends'](e,G);
        if 1 < nops(x) then y := x[2]; x := x[1]
        else x := x[1]
        fi;
	if ( type(networks['head'](e,G),integer) or
             substring(convert(networks['head'](e,G),string),1..1) = `:` )
        then
        lines[e] :=
        #** update the line colour and linestyle
        CURVES([pos[x],pos[y]],COLOUR('RGB',1.0,0,0),LINESTYLE(0))
        else
        lines[e] :=
        #** update the line colour and linestyle
        CURVES([pos[x],pos[y]],COLOUR('RGB',0,0,1.0),LINESTYLE(3))
        fi;
    od;
    t := map(
     op,{entries(points),entries(lines),entries(text)}
     );
    xrng := origin[1]+1/2 .. origin[1]+n+1/2;
    yrng := origin[2]-1/2 .. origin[2]+max_x-1/2;
    RETURN(t)
end:

# Calculate statistics for a process network,
netstats:=proc(first,Gr)
          local mat, mats, unitop, unitops, nodes,
                streams, systems, subsystems, loops;
          mat:=[op(op(1,op(1,[first]))),
               op(op(2,op(1,[first]))),
               op(op(3,op(1,[first])))];
          mat;
          mats:=nops(mat);
          unitop:=[op(op(4,op(1,[first])))];
          unitops:=nops(unitop);
          nodes:=mats+unitops;
          streams:=nops(edges(Gr));
          systems:=components(Gr);
          subsystems:=nops(systems);
          loops:=streams-nodes+subsystems;
          `# of materials = `.mats, 
          `# of unit ops = `.unitops,
          `# of network nodes (materials+unit ops) = `.nodes,
          `# of connections = `.streams,
          `# of subsystems = `.subsystems,
          `# of loops = `.loops; 
end:

# Add routines related to generating a graphic flow diagram
# of the process

# Generate Op data structure on available unit operations from
# the full unit operations set data.
OpGen:=proc(OpData)
          local i, Op, count, Labels;
          global LabelsList;
          count:=nops(OpData);
          # setup unique list for the unit operation names
          Labels:={};  
          for i from 1 to count do
              Labels:=Labels union {op(1,op(i,OpData))}
          od;
          LabelsList:=convert(Labels,list);
          # check for the proper number of inlet ports
          for i from 1 to count do
              if nops(op(2,op(i,OpData))) <> 
                 nops(op(6,op(i,OpData)))  then
                 RETURN(`Error: Wrong number of inlet stream port values.  Line=`,i)
              fi;
          od;
          # check for the proper number of outlet ports
          for i from 1 to count do
              if nops(op(3,op(i,OpData))) <> 
                 nops(op(7,op(i,OpData)))  then
                 RETURN(`Error: Wrong number of outlet stream port values.  Line=`,i)
              fi;
          od;
          Op:={};
          for i from 1 to nops(OpData)  do
              Op := Op union {[ {op(op(2,op(i,OpData)))},
                               {op(op(3,op(i,OpData)))} ]};
	  od;
          Op;
end:

# Convert the Pgraph structure into data for a Visio flowsheet plot.
Visio_data:=proc(PGraph,filename)
   local i, j, k, m, opset, rm, pr, mid, el, el2,
         mat_ins, mat_outs, unitnum, place,  
         inset, outset, inset2, outset2, nextunit,
         unitname, matname, port1, port2,  vals,
         max_els, unitprint, routes, route_count, Tree, 
         AP, route_pair, max_route, table_data, route_value,
         route, Gr, Gr2, boxes, location, Streams, nextunitname,
         fileID;
   global  OpData, LabelsList, none, numbers;
   boxes:={};
   m:=op(1,PGraph);
   opset:=op(2,PGraph);
   # generate the network from the PGraph
   Gr:=op(2,[graph_network(PGraph,numbers)]);
   # check for material inlets and outlets of the structure
   mat_ins:=Psi_minus(opset);
   mat_outs:=Psi_plus(opset);
   # find all raw materials (they are only inlets to unit ops)
   rm:={};
   # find all products (they are only outlets of unit ops)
   pr:={};
   for el in m  do
       if member(el,mat_ins) and member(el,mat_outs) = false then
          rm:=rm union {el}
       fi;
       if member(el,mat_ins) = false and member(el,mat_outs) then
          pr:=pr union {el}
       fi;
   od;
   mid:=m minus rm minus pr;
   
   # Print info for all boxes on the diagram,
   fileID:=fopen(`visio.dat`,WRITE,TEXT);
   fprintf(fileID,`%s %s\n`,`File:`, filename);
   # raw materials have the Feed icon
   for el in rm do
       vals[el]:=1;
       fprintf(fileID,`"%s" "towers" "Flowsheet Feed"\n`,el);
       boxes:=boxes union {el};
   od;
   # mid materials use a box icon if they have multiple feeds 
   for el in mid do
       Tree:=shortpathtree(Gr,el);
       if nops(daughter(el,Tree)) > 1 or nops(ancestor(el,Tree)) > 1 then
          vals[el]:=1;
          fprintf(fileID,`"%s" "towers" "Diamond"\n`,el);
          boxes:=boxes union {el};
       else
          vals[el]:=0;
          fprintf(fileID,`"%s" "none" "none"\n`,el);
       fi;
   od;
   # products have the Product icon
   for el in pr do
       vals[el]:=1;
       fprintf(fileID,`"%s" "towers" "Flowsheet Feed"\n`,el);
       boxes:=boxes union {el};
   od;
   # the unit ops have their special icons from OpData
   j:=0;
   for el in opset do
       j:=j+1;
       inset:=op(1,el);
       outset:=op(2,el);
       for i from 1 to nops(OpData) do
           inset2:={op(op(2,op(i,OpData)))};
           outset2:={op(op(3,op(i,OpData)))};
           if inset = inset2 and outset = outset2 then break fi;
       od;
       fprintf(fileID,`"%s" "%s" "%s"\n`,op(1,op(i,OpData)),op(4,op(i,OpData)),
                 op(5,op(i,OpData)) ) ;
       member(op(1,op(i,OpData)),LabelsList,'location');
       unitnum[j]:=i;
       boxes:=boxes union {location};
   od;
   # print the stream connectivity data
   # source node, destination node, source port, dest. port, label
   fprintf(fileID,`Stream Data:\n`);
   # cycle through all units in the process
   Streams:={};
   for i from 1 to nops(opset) do
       unitname:=op(1,op(unitnum[i],OpData));
       # do inlets to the unit
       for el in op(1,op(i,opset)) do
           matname:=el;
           member(matname,op(2,op(unitnum[i],OpData)),'place');
           port2:=op(place,op(6,op(unitnum[i],OpData)));
           if {el} intersect rm <> {} then 
              Streams:=Streams union
                 {[matname,unitname,1,port2,none]}
	   elif vals[el] = 1 then
              Streams:=Streams union
                 {[matname,unitname,5,port2,none]}
           else
              nextunit:=arrivals(el,Gr);
              nextunitname:=op(op(nextunit),LabelsList);
              for k from 1 to nops(OpData)  do
                  if op(1,op(k,OpData)) = nextunitname then
                     if member(matname,op(3,op(k,OpData)),'place')
                        then break  fi;
                  fi;
              od; 
              port1:=op(place,op(7,op(k,OpData)));
              Streams:=Streams union
                 {[nextunitname,unitname,port1,port2,matname]}
           fi;
       od;
       # do outlets of the unit
       for el in op(2,op(i,opset)) do
           matname:=el;
           member(matname,op(3,op(unitnum[i],OpData)),'place');
           port1:=op(place,op(7,op(unitnum[i],OpData)));
           if {el} intersect pr <> {} then 
              Streams:=Streams union
                 {[unitname,matname,port1,2,none]}
           elif vals[el] = 1 then
              Streams:=Streams union
                 {[unitname,matname,port1,3,none]}
           else
              NULL;  # the connection was already made
           fi;
       od;
   od;
   # print the stored stream connectivities,
   for el in Streams do
       fprintf(fileID,`"%s" "%s" "%d" "%d" "%s"\n`,op(el))
   od;
# print the main route through the process followed by sets
# with the other units.
   fprintf(fileID,`Routes:\n`);
   # find all routes from feeds to products,  tag the longest one
   max_els:=0;
   for el in rm do
       Tree:=shortpathtree(Gr,el);
       for el2 in pr  do
           # Note: path returns the vertices in reverse order
           routes[el,el2]:=path([el2,el],Tree);
           # count only the vertices that are part of the diagram
           route_count[el,el2]:=
              nops({op(routes[el,el2])} intersect boxes);
           # find the maximum pathway
           if route_count[el,el2] > max_els then
              route:=el,el2;
              max_els:=route_count[el,el2];
           fi;
       od;
   od;
   # print the maximum route  (it is in reverse order)
   fprintf(fileID,`%d\n`,max_els);
   k:=nops(routes[route]);
   for i from 1 to k  do
       j:=k+1-i;
       if {op(j,routes[route])} intersect boxes <> {} then
          unitprint:=op(j,routes[route]);
          if whattype(unitprint) = integer then  
             fprintf(fileID,`"%s"\n`,op(unitprint,LabelsList))
          else
             fprintf(fileID,`"%s"\n`,unitprint);
          fi;
       fi;
   od;
   # if there are more than 1 box left, then find the longest path
   # through them and output it.
   boxes:=boxes minus {op(routes[route])};
   if nops(boxes) > 1  then
       # get a new graph with only the remaining vertices
       new(Gr2);
       for el in boxes do
           Gr2:=gunion(Gr2,induce(incident(el,Gr),Gr),'SIMPLE');
       od;
       # find the distances for all pairs
       AP:=allpairs(Gr2);
       # find the maximum distance
       route_pair:=0,0;
       max_route:=0;
       table_data:=op(op(AP));
       for i from 1 to nops(table_data) do
           route_value:=op(2,op(i,table_data));
           if route_value > max_route and route_value < infinity then
              max_route:=route_value;
              route_pair:=op(1,op(i,table_data));
           fi;
       od;
       # find the route between the route_pair set
       Tree:=shortpathtree(Gr2,op(1,[route_pair]));
       route:=path([op(2,[route_pair]),op(1,[route_pair])],Tree);
       # print the number of units in the path and their names in reverse
       fprintf(fileID,`%d\n`,nops({op(route)} intersect boxes));
       k:=nops(route);
       if k > 1 then
          for i from 1 to k  do
              el:=op(k+1-i,route);
              if ({el} intersect boxes)  <> {} then 
                 if whattype(el) = integer then 
                   fprintf(fileID,`"%s"\n`,op(el,LabelsList))
                 else
                   fprintf(fileID,`"%s"\n`,el)
                 fi;
              fi;
          od;
       else 
          fprintf(fileID,`0\n`);
          route:=[];
       fi;
   else
       fprintf(fileID,`0\n`)
   fi;
   # if there are more than 1 box left, then find the longest path
   # through them and output it.
   boxes:=boxes minus {op(route)};
   if nops(boxes) > 1  then
       # get a new graph with only the remaining vertices
       new(Gr2);
       for el in boxes do
           Gr2:=gunion(Gr2,induce(incident(el,Gr),Gr),'SIMPLE');
       od;
       # find the distances for all pairs
       AP:=allpairs(Gr2);
       # find the maximum distance
       route_pair:=0,0;
       max_route:=0;
       table_data:=op(op(AP));
       for i from 1 to nops(table_data) do
           route_value:=op(2,op(i,table_data));
           if route_value > max_route and route_value < infinity then
              max_route:=route_value;
              route_pair:=op(1,op(i,table_data));
           fi;
       od;
       # find the route between the route_pair set
       Tree:=shortpathtree(Gr2,op(1,[route_pair]));
       route:=path([op(2,[route_pair]),op(1,[route_pair])],Tree);
       # print the number of units in the path and their names in reverse
       fprintf(fileID,`%d\n`,nops({op(route)} intersect boxes));
       k:=nops(route);
       if k > 1 then
          for i from 1 to k  do
              el:=op(k+1-i,route);
              if ({el} intersect boxes)  <> {} then 
                 if whattype(el) = integer then 
                   fprintf(fileID,`"%s"\n`,op(el,LabelsList))
                 else
                   fprintf(fileID,`"%s"\n`,el)
                 fi;
              fi;
          od;
       else 
          fprintf(fileID,`0\n`);
          route:=[];
       fi;
   else
       fprintf(fileID,`0\n`)
   fi;
   # print the remaining units in the order of 1) raw materials
   # 2) mid materials, 3) unit ops, and 4) products
   boxes:=boxes minus {op(route)};
   fprintf(fileID,`%d\n`,nops(boxes));
   for el in boxes  do
       if ({el} intersect rm) <> {}  then fprintf(fileID,`"%s"\n`,el) fi
   od;
   for el in boxes  do
       if ({el} intersect mid) <> {} then fprintf(fileID,`"%s"\n`,el) fi
   od;
   for el in boxes  do
       if whattype(el) = integer then 
       fprintf(fileID,`"%s"\n`,op(el,LabelsList))   fi
   od;
   for el in boxes  do
       if ({el} intersect pr) <> {}  then fprintf(fileID,`"%s"\n`,el) fi
   od;
   fprintf(fileID,`End:\n`);
   close(fileID);
end:

# calculate the rank value of a decision mapping
rank_value:=proc(mapping)
   local el, i, m, opset, rankings, inset, outset, inset2, outset2,
         unitname, unitrank, PGr, tabledata, sumrank, mat_ins,
         mat_outs, rm, pr;
   global OpData, Material_Rank;
   # rankings is a table for storing the rank of each unit op
   rankings:=table();
   # the inputted decision map is converted to a P graph (m,o)
   PGr:=Pgraph(mapping);
   m:=op(1,PGr);
   opset:=op(2,PGr);
   # check for material inlets and outlets of the structure
   mat_ins:=Psi_minus(opset);
   mat_outs:=Psi_plus(opset);
   # find all raw materials (they are only inlets to unit ops)
   rm:={};
   # find all products (they are only outlets of unit ops)
   pr:={};
   for el in m  do
       if member(el,mat_ins) and member(el,mat_outs) = false then
          rm:=rm union {el}
       fi;
       if member(el,mat_ins) = false and member(el,mat_outs) then
          pr:=pr union {el}
       fi;
   od;
   # cycle through the unit ops within opset to find their location
   # within the OpData array;  extract the name & rank for rankings table
   for el in opset do
       # match the unit op to its row within OpData
       inset:=op(1,el);
       outset:=op(2,el);
       for i from 1 to nops(OpData) do
           inset2:={op(op(2,op(i,OpData)))};
           outset2:={op(op(3,op(i,OpData)))};
           if inset = inset2 and outset = outset2 then break fi;
       od;
       # extract the name and rank
       unitname:=op(1,op(i,OpData));
       unitrank:=op(8,op(i,OpData));
       # store the info within the rankings table;  use the highest
       # ranking that occurs for a given unit name
       if whattype(rankings[unitname]) = indexed  then
           rankings[unitname]:=unitrank
       else
           if rankings[unitname] < unitrank then 
              rankings[unitname]:=unitrank fi;
       fi;
   od;
   # Get the data from inside the rankings table
   tabledata:=op(2,op(1,rankings));
   # Sum the ranks of the unit operations
   sumrank:=convert([seq(rhs(op(i,tabledata)),i=1..nops(tabledata))],`+`);
   # Add the ranks of the raw materials and products within the flowsheet
   m:=rm union pr;
   for el in m  do
       if whattype(Material_Rank[el]) <> indexed then
          sumrank:=sumrank+Material_Rank[el]  fi; 
   od;
   sumrank;
end:

# provide the sort criteria for the flowsheets list; sort by
# increasing rank value.
sort_rank:=proc (a, b)
   #  note:  the rank is stored as the first element of the list
   if op(1,a) < op(1,b) then  true
   else  false
   fi;
end:

# add the rank to each decision mapping and sort it
process_maps:=proc(dm) 
   local j;
   global  flowsheets;
   # create a list containing lists with the flowsheet rank and its
   # decision map.
   flowsheets:=[seq([rank_value(op(j,dm)),op(j,dm)],
                j=1..nops(dm))];
   # sort the list by increasing rank; the sort_rank function provides
   # the boolean function for the sorting task.
   flowsheets:=sort(flowsheets,sort_rank);
   flowsheets;
end:

# use cycles in the flowsheet to determine key material
# balance requirements,
cycles:=proc(first,Gr,rootVertex)
   local el, el2, edgeset, Grund, G2, eset, G3, cyc, G4,
         cycleVertices, chords, numberChords, units,
         cycleFeeds, cycleProducts;
   global LabelsList, Raw_Materials;
   # first, generate an undirected graph from Gr
   edgeset:={};
   for el in edges(Gr) do
       edgeset:=edgeset union {{head(el,Gr),tail(el,Gr)}};
   od;
   Grund:=graph(vertices(Gr),edgeset);
   # find a spanning tree of Grund rooted at a raw material,
   G2:=spantree(Grund,rootVertex);
   # find the chords of the spanning tree,
   chords:=edges(Grund) minus edges(G2);
   # if there are no chords then stop,
   numberChords:=nops(chords);
   if numberChords = 0 then  RETURN(`No loops found.`)  fi;
   # loop through all fundamental cycles,
   for el in chords  do
      # add a chord to the edge set
      eset:={el} union edges(G2);
      # find the fundamental cycle
      cyc:=fundcyc(eset,Grund);
      # create a subgraph from the cycle
      G3:=induce(cyc,Grund);
      # find the vertices in the graph
      cycleVertices:=vertices(G3);
      # delete any vertices that are raw materials or products
      for el2 in cycleVertices do
          if nops(arrivals(el2,Gr)) = 0 then
             cycleVertices:=cycleVertices minus {el2} fi;
          if nops(departures(el2,Gr)) = 0 then
             cycleVertices:=cycleVertices minus {el2} fi; 
      od;
      # create a new graph from the original Gr with
      # the condensation of the cycleVertices,
      G4:=duplicate(Gr);
      shrink(cycleVertices,G4,new_vertex);
      # print results,
      units:={};
      for el2 in cycleVertices do
          if whattype(el2) = integer then
             units:=units union {op(el2,LabelsList)}
          else
             units:=units union {el2}
          fi;
      od;
      cycleFeeds:={};
      for el2 in arrivals(new_vertex,G4) do
          if whattype(el2) = integer  then
             cycleFeeds:=cycleFeeds union {op(el2,LabelsList)}
          else
             cycleFeeds:=cycleFeeds union {el2}
          fi;
      od;
      cycleProducts:={};
      for el2 in departures(new_vertex,G4) do
          if whattype(el2) = integer  then
             cycleProducts:=cycleProducts union {op(el2,LabelsList)}
          else
             cycleProducts:=cycleProducts union {el2}
          fi;
      od;
      lprint(`Cycle found:`);
      lprint(`Fundamental cycle vertices (w/o RM's and Pr's) = `, units);
      lprint(`Feeds from =`,cycleFeeds);
      lprint(`Products to =`,cycleProducts);
   od;
   `# of cycles processed = `.numberChords;
end: