@@ -33,10 +33,18 @@ def apply_reduction(
3333 net : PetriNet , im : Marking , fm : Marking
3434) -> Tuple [PetriNet , Marking , Marking ]:
3535 """
36- Apply the Murata reduction to an accepting Petri net, removing the structurally redundant places.
36+ Apply the Murata reduction to an accepting Petri net, removing structurally redundant (implicit) places.
3737
38- The implementation follows the Berthelot algorithm as in:
39- https://svn.win.tue.nl/repos/prom/Packages/Murata/Trunk/src/org/processmining/algorithms/BerthelotAlgorithm.java
38+ This implementation follows the (Berthelot) implicit-place check by searching, via ILP,
39+ for a place-flow f = a_p * p - Σ a_q * q (a_p >= 1, a_q >= 0) such that:
40+
41+ (1) f is a P-invariant (flow): a_p*C(p,t) - Σ a_q*C(q,t) = 0 for all transitions t
42+ where C(p,t) = Post(p,t) - Pre(p,t)
43+
44+ (2) a_p*Pre(p,t) - Σ a_q*Pre(q,t) <= k for all transitions t
45+ with k = a_p*M0(p) - Σ a_q*M0(q), and k >= 0
46+
47+ If such a solution exists, the candidate place is implicit and can be removed without changing behavior.
4048
4149 Parameters
4250 ---------------
@@ -50,101 +58,137 @@ def apply_reduction(
5058 Returns
5159 --------------
5260 net
53- Petri net
61+ Reduced Petri net
5462 im
55- Initial marking
63+ Initial marking (unchanged)
5664 fm
57- Final marking
65+ Final marking (unchanged)
5866 """
5967 places = sorted (list (net .places ), key = lambda x : x .name )
60- place_index = {p : i for i , p in enumerate (places )}
61- n_places = len (places )
62-
6368 redundant = set ()
6469
70+ # Choose solver backend
71+ proposed_solver = solver .SCIPY
72+ if importlib .util .find_spec ("pulp" ):
73+ proposed_solver = solver .PULP
74+ else :
75+ if constants .SHOW_INTERNAL_WARNINGS :
76+ warnings .warn (
77+ "solution from scipy may be unstable. Please install PuLP (pip install pulp) for fully reliable results."
78+ )
79+
6580 for place in places :
81+ # Skip already proven redundant places
82+ if place in redundant :
83+ continue
84+
6685 # Skip places in the initial or final markings
6786 if place in im or place in fm :
6887 continue
6988
89+ # Work on the "current" net as if previously found redundant places were removed:
90+ active_places = [p for p in places if p not in redundant ]
91+ active_places = sorted (active_places , key = lambda x : x .name )
92+ place_index = {p : i for i , p in enumerate (active_places )}
93+
94+ n_places = len (active_places )
95+ n_vars = n_places + 1 # last variable is k
96+ k_idx = n_vars - 1
97+
7098 Aeq = []
7199 Aub = []
72100 beq = []
73101 bub = []
74102
75- # first constraint
76- constraint = [0 ] * (n_places + 1 )
77- for p2 in im :
78- if p2 not in redundant :
79- if p2 == place :
80- constraint [place_index [p2 ]] = im [p2 ]
81- else :
82- constraint [place_index [p2 ]] = - im [p2 ]
83- constraint [- 1 ] = - 1
84- Aeq .append (constraint )
103+ # ---------------------------------------------------------------------
104+ # (A) k = a_p*M0(p) - Σ a_q*M0(q)
105+ # => a_p*M0(p) - Σ a_q*M0(q) - k = 0
106+ # ---------------------------------------------------------------------
107+ eq = [0 ] * n_vars
108+ for p2 in active_places :
109+ m0 = im [p2 ] if p2 in im else 0
110+ if m0 == 0 :
111+ continue
112+ if p2 == place :
113+ eq [place_index [p2 ]] += m0
114+ else :
115+ eq [place_index [p2 ]] -= m0
116+ eq [k_idx ] = - 1
117+ Aeq .append (eq )
85118 beq .append (0 )
86119
87- # second constraints
120+ # ---------------------------------------------------------------------
121+ # (B) Flow constraints: for all transitions t,
122+ # a_p*C(p,t) - Σ a_q*C(q,t) = 0
123+ # where C(p,t)=Post(p,t)-Pre(p,t)
124+ # ---------------------------------------------------------------------
88125 for trans in net .transitions :
89- constraint = [0 ] * (n_places + 1 )
126+ pre = {}
127+ post = {}
90128
91129 for arc in trans .in_arcs :
92- p2 = arc .source
93- if p2 not in redundant :
94- if p2 == place :
95- constraint [place_index [p2 ]] = arc .weight
96- else :
97- constraint [place_index [p2 ]] = - arc .weight
98-
99- constraint [- 1 ] = - 1
100- Aub .append (constraint )
101- bub .append (0 )
102-
103- # third constraints (FIXED SIGNS)
104- for trans in net .transitions :
105- constraint = [0 ] * (n_places + 1 )
130+ src = arc .source
131+ if src in place_index :
132+ pre [src ] = pre .get (src , 0 ) + arc .weight
106133
107134 for arc in trans .out_arcs :
108- p2 = arc .target
109- if p2 not in redundant :
110- # FIX: candidate place must have +weight, other places -weight
111- # (the previous implementation inverted these signs)
112- if p2 == place :
113- constraint [place_index [p2 ]] = arc .weight
114- else :
115- constraint [place_index [p2 ]] = - arc .weight
116-
117- Aub .append (constraint )
118- bub .append (0 )
119-
120- # fourth constraint
121- for p2 in places :
122- if p2 not in redundant :
123- constraint = [0 ] * (n_places + 1 )
124- constraint [place_index [p2 ]] = - 1
125- Aub .append (constraint )
135+ tgt = arc .target
136+ if tgt in place_index :
137+ post [tgt ] = post .get (tgt , 0 ) + arc .weight
138+
139+ eq = [0 ] * n_vars
140+ for p2 in active_places :
141+ c_val = post .get (p2 , 0 ) - pre .get (p2 , 0 )
142+ if c_val == 0 :
143+ continue
126144 if p2 == place :
127- bub . append ( - 1 )
145+ eq [ place_index [ p2 ]] += c_val
128146 else :
129- bub .append (0 )
130-
131- # fifth constraint
132- constraint = [0 ] * (n_places + 1 )
133- constraint [- 1 ] = - 1
134- Aub .append (constraint )
135- bub .append (0 )
136-
137- c = [1 ] * (n_places + 1 )
138- integrality = [1 ] * (n_places + 1 )
139-
140- proposed_solver = solver .SCIPY
141- if importlib .util .find_spec ("pulp" ):
142- proposed_solver = solver .PULP
143- else :
144- if constants .SHOW_INTERNAL_WARNINGS :
145- warnings .warn (
146- "solution from scipy may be unstable. Please install PuLP (pip install pulp) for fully reliable results."
147- )
147+ eq [place_index [p2 ]] -= c_val
148+ # k has coefficient 0 here
149+ Aeq .append (eq )
150+ beq .append (0 )
151+
152+ # ---------------------------------------------------------------------
153+ # (C) Implicitness inequalities: for all transitions t,
154+ # a_p*Pre(p,t) - Σ a_q*Pre(q,t) - k <= 0
155+ # ---------------------------------------------------------------------
156+ for trans in net .transitions :
157+ ineq = [0 ] * n_vars
158+ for arc in trans .in_arcs :
159+ p2 = arc .source
160+ if p2 not in place_index :
161+ continue
162+ if p2 == place :
163+ ineq [place_index [p2 ]] += arc .weight
164+ else :
165+ ineq [place_index [p2 ]] -= arc .weight
166+ ineq [k_idx ] = - 1
167+ Aub .append (ineq )
168+ bub .append (0 )
169+
170+ # ---------------------------------------------------------------------
171+ # (D) Variable lower bounds:
172+ # a_p >= 1, a_q >= 0, k >= 0
173+ # Implemented as -a <= b
174+ # ---------------------------------------------------------------------
175+ for p2 in active_places :
176+ ineq = [0 ] * n_vars
177+ ineq [place_index [p2 ]] = - 1
178+ Aub .append (ineq )
179+ if p2 == place :
180+ bub .append (- 1 ) # -a_p <= -1 -> a_p >= 1
181+ else :
182+ bub .append (0 ) # -a_q <= 0 -> a_q >= 0
183+
184+ ineq = [0 ] * n_vars
185+ ineq [k_idx ] = - 1
186+ Aub .append (ineq )
187+ bub .append (0 ) # -k <= 0 -> k >= 0
188+
189+ # Objective (any feasible solution is enough; minimize sum to keep it small)
190+ c = [1 ] * n_vars
191+ integrality = [1 ] * n_vars
148192
149193 xx = solver .apply (
150194 c ,
@@ -161,7 +205,8 @@ def apply_reduction(
161205 ):
162206 redundant .add (place )
163207
164- for place in redundant :
165- net = remove_place (net , place )
208+ # Remove redundant places (deterministic order)
209+ for pl in sorted (list (redundant ), key = lambda x : x .name ):
210+ net = remove_place (net , pl )
166211
167212 return net , im , fm
0 commit comments