1616
1717
1818def F (curve , w_t ):
19- """ The F(γ) operator, defined as F(γ) = W(γ)/L(γ) .
19+ """ The F(γ) operator, minimization target in the insertion step .
2020
2121 Parameters
2222 ----------
23- curve : DGCG.classes.curve class
24- w_t : DGCG.classes.dual_variable
23+ curve : :py:class:`src.classes.curve`
24+ Curve γ where the F operator is evaluated.
25+ w_t : :py:class:`src.classes.dual_variable`
26+ Dual variable that defines the F operator.
2527
2628 Returns
2729 -------
28- double number.
30+ float
31+
2932
3033 Notes
3134 -----
32- When solving the insertion step, this is the main energy to minimize.
35+ The F operator is defined via the dual variable as
36+
37+ .. math::
38+ F(\\ gamma) = -\\ frac{a_{\\ gamma}}{T+1} \\ sum_{t=0}^T w_t(\\ gamma(t))
39+
3340 """
3441 assert isinstance (curve , classes .curve ) and \
3542 isinstance (w_t , classes .dual_variable )
@@ -41,19 +48,23 @@ def grad_F(curve, w_t):
4148
4249 Parameters
4350 ----------
44- curve : DGCG.classes.curve class
45- w_t : DGCG.classes.dual_variable
51+ curve : :py:class:`src.classes.curve`
52+ Curve γ where the F operator is evaluated.
53+ w_t : :py:class:`src.classes.dual_variable`
54+ Dual variable that defines the F operator.
4655
4756 Returns
4857 -------
49- double number.
58+ :py:class:`src.classes.curve`
5059
5160 Notes
5261 -----
53- We use the gradient to minimize F(γ).
62+ The F operator is defined on the Hilbert space of curves, therefore the
63+ gradient should be a curve.
5464 """
5565 assert isinstance (curve , classes .curve ) and \
5666 isinstance (w_t , classes .dual_variable )
67+ # F = W/L
5768 L_gamma = curve .energy ()
5869 W_gamma = - curve .integrate_against (w_t )
5970 # ∇L(γ) computation
@@ -68,15 +79,13 @@ def grad_F(curve, w_t):
6879 w_t_curve = lambda t : w_t .grad_eval (t , curve .eval_discrete (t )).reshape (2 )
6980 T = config .T
7081 grad_W_gamma = - np .array ([t_weigh [t ]* w_t_curve (t ) for t in range (T )])
71- # grad_W_gamma = -np.array([config.time_weights[t] *
72- # w_t.grad_eval(t, curve.eval_discrete(t)).reshape(2)
73- # for t in range(config.T)])
7482 # (L(γ)∇W(γ)-W(γ)∇L(γ))/L(γ)²
7583 pos_gradient = (L_gamma * grad_W_gamma - W_gamma * grad_L_gamma )/ L_gamma ** 2
7684 gradient_curve = classes .curve (pos_gradient )
7785 return gradient_curve
7886
79- def after_optimization_sparsifier (current_measure , energy_curves = None ):
87+
88+ def after_optimization_sparsifier (current_measure ):
8089 """ Trims a sparse measure by merging atoms that are too close.
8190
8291 Given a measure composed of atoms, it will look for the atoms that are
@@ -85,10 +94,8 @@ def after_optimization_sparsifier(current_measure, energy_curves=None):
8594
8695 Parameters
8796 ----------
88- current_measure : DGCG.classes.measure class
89- energy_curves : numpy.ndarray, optional
90- vector indicating the energy of the curves of the measure. To
91- accelerate the comparisons.
97+ current_measure : :py:class:`src.classes.measure`
98+ Target measure to trim.
9299
93100 Returns
94101 -------
@@ -97,9 +104,13 @@ def after_optimization_sparsifier(current_measure, energy_curves=None):
97104 Notes
98105 -----
99106 This method is required because the quadratic optimization step is realized
100- by an interior point method. Therefore, it is likely to find minimums in
101- between two identical items instead of selecting one and discarding the
102- other.
107+ by an interior point method. Therefore, in the case that there are repeated
108+ (or very close to repeated) atoms in the current measure, the quadratic
109+ optimization step can give positive weights to both of them.
110+
111+ This is not desirable, since besides incrementing the computing power for
112+ the sliding step, we would prefer each atom numerically represented only
113+ once.
103114 """
104115 output_measure = copy .deepcopy (current_measure )
105116 id1 = 0
@@ -140,6 +151,8 @@ def after_optimization_sparsifier(current_measure, energy_curves=None):
140151 return output_measure
141152
142153def solve_quadratic_program (current_measure ):
154+ """
155+ """
143156 assert isinstance (current_measure , classes .measure )
144157 # Build the quadratic system of step 5 and then use some generic python
145158 # solver to get a solution.
@@ -178,32 +191,6 @@ def solve_quadratic_program(current_measure):
178191 logger .status ([1 , 2 , 2 ], coefficients )
179192 return curves_list , coefficients
180193
181- def to_positive_semidefinite (Q ):
182- """ Takes a symmetric matrix and returns a positive semidefinite projection
183-
184- Parameters
185- ----------
186- Q : numpy.ndarray
187- symmetric matrix
188-
189- Returns
190- -------
191- numpy.ndarray, symmetric positive semidefinite matrix.
192- """
193- min_eigval = 0
194- eigval , eigvec = np .linalg .eigh (Q )
195- if min (eigval ) < min_eigval :
196- # truncate
197- print ("Negative eigenvalues: " ,
198- [eig for eig in eigval if eig < min_eigval ])
199- eigval = np .maximum (eigval , min_eigval )
200- # Recompute Q = VΣV^(-1)
201- Q2 = np .linalg .solve (eigvec .T , np .diag (eigval )@eigvec .T ).T
202- print ("PSD projection relative norm difference:" ,
203- np .linalg .norm (Q - Q2 )/ np .linalg .norm (Q ))
204- return Q2
205- return Q
206-
207194def weight_optimization_step (current_measure ):
208195 config .logger .status ([1 , 2 , 1 ])
209196 # optimizes the coefficients for the current_measure
0 commit comments