Skip to content

Commit d500a54

Browse files
rileyjmurrayrileyjmurray
committed
initial changes to RawTVDFunction
1 parent 916877c commit d500a54

File tree

2 files changed

+50
-35
lines changed

2 files changed

+50
-35
lines changed

‎pygsti/algorithms/core.py‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1007,11 +1007,13 @@ def _do_runopt(objective, optimizer, printer):
10071007
tm = _time.time()
10081008
nDataParams = objective.num_data_params() # TODO - cache this somehow in term-based calcs...
10091009
profiler.add_time("run_gst_fit: num data params", tm)
1010-
1011-
chi2_k_qty = opt_result.chi2_k_distributed_qty # total chi2 or 2*deltaLogL
10121010
desc = objective.description
1011+
chi2_k_qty = opt_result.chi2_k_distributed_qty # total chi2 or 2*deltaLogL
1012+
if chi2_k_qty > 0:
10131013
# reject GST model if p-value < threshold (~0.05?)
1014-
pvalue = 1.0 - _stats.chi2.cdf(chi2_k_qty, nDataParams - nModelParams)
1014+
pvalue = 1.0 - _stats.chi2.cdf(chi2_k_qty, nDataParams - nModelParams)
1015+
else:
1016+
pvalue = 0.0
10151017
printer.log("%s = %g (%d data params - %d (approx) model params = expected mean of %g; p-value = %g)" %
10161018
(desc, chi2_k_qty, nDataParams, nModelParams, nDataParams - nModelParams, pvalue), 1)
10171019

‎pygsti/objectivefns/objectivefns.py‎

Lines changed: 45 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414
import sys as _sys
1515
import time as _time
1616
import pathlib as _pathlib
17+
import warnings as _warnings
1718

1819
import numpy as _np
1920

@@ -557,6 +558,8 @@ def lsvec(self, probs, counts, total_counts, freqs, intermediates=None):
557558
"""
558559
return _np.sqrt(self.terms(probs, counts, total_counts, freqs, intermediates))
559560

561+
# Infinite loop in evaluation of "dterms" and "dlsvec".
562+
560563
def dterms(self, probs, counts, total_counts, freqs, intermediates=None):
561564
"""
562565
Compute the derivatives of the terms of this objective function.
@@ -591,8 +594,9 @@ def dterms(self, probs, counts, total_counts, freqs, intermediates=None):
591594
"""
592595
if intermediates is None:
593596
intermediates = self._intermediates(probs, counts, total_counts, freqs)
594-
return 2 * self.lsvec(probs, counts, total_counts, freqs, intermediates) \
595-
* self.dlsvec(probs, counts, total_counts, freqs, intermediates)
597+
u = self.lsvec(probs, counts, total_counts, freqs, intermediates)
598+
v = self.dlsvec(probs, counts, total_counts, freqs, intermediates)
599+
return 2 * u * v
596600

597601
def dlsvec(self, probs, counts, total_counts, freqs, intermediates=None):
598602
"""
@@ -2753,35 +2757,37 @@ def zero_freq_hterms(self, total_counts, probs):
27532757
return 2 * _np.ones(len(probs))
27542758

27552759

2756-
# The log(Likelihood) within the Poisson picture is: # noqa
2757-
# # noqa
2758-
# L = prod_{i,sl} lambda_{i,sl}^N_{i,sl} e^{-lambda_{i,sl}} / N_{i,sl}! # noqa
2759-
# # noqa
2760-
# Where lamba_{i,sl} := p_{i,sl}*N[i] is a rate, i indexes the operation sequence, # noqa
2761-
# and sl indexes the spam label. N[i] is the total counts for the i-th circuit, and # noqa
2762-
# so sum_{sl} N_{i,sl} == N[i]. We can ignore the p-independent N_j! and take the log: # noqa
2763-
# # noqa
2764-
# log L = sum_{i,sl} N_{i,sl} log(N[i]*p_{i,sl}) - N[i]*p_{i,sl} # noqa
2765-
# = sum_{i,sl} N_{i,sl} log(p_{i,sl}) - N[i]*p_{i,sl} (where we ignore the p-independent log(N[i]) terms) # noqa
2766-
# # noqa
2767-
# The objective function computes the negative log(Likelihood) as a vector of leastsq # noqa
2768-
# terms, where each term == sqrt( N_{i,sl} * -log(p_{i,sl}) + N[i] * p_{i,sl} ) # noqa
2769-
# # noqa
2770-
# See LikelihoodFunctions.py for details on patching # noqa
2771-
# The log(Likelihood) within the standard picture is:
2772-
#
2773-
# L = prod_{i,sl} p_{i,sl}^N_{i,sl}
2774-
#
2775-
# Where i indexes the operation sequence, and sl indexes the spam label.
2776-
# N[i] is the total counts for the i-th circuit, and
2777-
# so sum_{sl} N_{i,sl} == N[i]. We take the log:
2778-
#
2779-
# log L = sum_{i,sl} N_{i,sl} log(p_{i,sl})
2780-
#
2781-
# The objective function computes the negative log(Likelihood) as a vector of leastsq
2782-
# terms, where each term == sqrt( N_{i,sl} * -log(p_{i,sl}) )
2783-
#
2784-
# See LikelihoodFunction.py for details on patching
2760+
"""
2761+
The log(Likelihood) within the Poisson picture is: # noqa
2762+
# noqa
2763+
L = prod_{i,sl} lambda_{i,sl}^N_{i,sl} e^{-lambda_{i,sl}} / N_{i,sl}! # noqa
2764+
# noqa
2765+
Where lamba_{i,sl} := p_{i,sl}*N[i] is a rate, i indexes the operation sequence, # noqa
2766+
and sl indexes the spam label. N[i] is the total counts for the i-th circuit, and # noqa
2767+
so sum_{sl} N_{i,sl} == N[i]. We can ignore the p-independent N_j! and take the log: # noqa
2768+
# noqa
2769+
log L = sum_{i,sl} N_{i,sl} log(N[i]*p_{i,sl}) - N[i]*p_{i,sl} # noqa
2770+
= sum_{i,sl} N_{i,sl} log(p_{i,sl}) - N[i]*p_{i,sl} (where we ignore the p-independent log(N[i]) terms) # noqa
2771+
# noqa
2772+
The objective function computes the negative log(Likelihood) as a vector of leastsq # noqa
2773+
terms, where each term == sqrt( N_{i,sl} * -log(p_{i,sl}) + N[i] * p_{i,sl} ) # noqa
2774+
# noqa
2775+
See LikelihoodFunctions.py for details on patching # noqa
2776+
The log(Likelihood) within the standard picture is:
2777+
2778+
L = prod_{i,sl} p_{i,sl}^N_{i,sl}
2779+
2780+
Where i indexes the operation sequence, and sl indexes the spam label.
2781+
N[i] is the total counts for the i-th circuit, and
2782+
so sum_{sl} N_{i,sl} == N[i]. We take the log:
2783+
2784+
log L = sum_{i,sl} N_{i,sl} log(p_{i,sl})
2785+
2786+
The objective function computes the negative log(Likelihood) as a vector of leastsq
2787+
terms, where each term == sqrt( N_{i,sl} * -log(p_{i,sl}) )
2788+
2789+
See LikelihoodFunction.py for details on patching
2790+
"""
27852791
class RawPoissonPicDeltaLogLFunction(RawObjectiveFunction):
27862792
"""
27872793
The function `N*f*log(f/p) - N*(f-p)`.
@@ -4018,6 +4024,9 @@ def __init__(self, regularization=None,
40184024
resource_alloc=None, name='tvd', description="Total Variational Distance (TVD)", verbosity=0):
40194025
super().__init__(regularization, resource_alloc, name, description, verbosity)
40204026

4027+
def chi2k_distributed_qty(self, objective_function_value):
4028+
return -1
4029+
40214030
def terms(self, probs, counts, total_counts, freqs, intermediates=None):
40224031
"""
40234032
Compute the terms of the objective function.
@@ -4083,7 +4092,11 @@ def dterms(self, probs, counts, total_counts, freqs, intermediates=None):
40834092
numpy.ndarray
40844093
A 1D array of length equal to that of each array argument.
40854094
"""
4086-
raise NotImplementedError("Derivatives not implemented for TVD yet!")
4095+
_warnings.warn('This derivative is discontinuous and does not return a full subgradient.')
4096+
t = self.terms(probs, counts, total_counts, freqs, intermediates)
4097+
d = 0.5*_np.ones_like(t)
4098+
d[t < 0] *= -1
4099+
return d
40874100

40884101
def hterms(self, probs, counts, total_counts, freqs, intermediates=None):
40894102
"""

0 commit comments

Comments
 (0)