/
model.py
2094 lines (1763 loc) · 95 KB
/
model.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
""" Defines the Model class and supporting functionality."""
from __future__ import division, print_function, absolute_import, unicode_literals
#***************************************************************************************************
# Copyright 2015, 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
# Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights
# in this software.
# Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
# in compliance with the License. You may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory.
#***************************************************************************************************
import numpy as _np
import scipy as _scipy
import itertools as _itertools
import collections as _collections
import warnings as _warnings
import time as _time
import uuid as _uuid
import bisect as _bisect
import copy as _copy
from ..tools import matrixtools as _mt
from ..tools import optools as _gt
from ..tools import slicetools as _slct
from ..tools import likelihoodfns as _lf
from ..tools import jamiolkowski as _jt
from ..tools import compattools as _compat
from ..tools import basistools as _bt
from ..tools import listtools as _lt
from ..tools import symplectic as _symp
from . import modelmember as _gm
from . import circuit as _cir
from . import operation as _op
from . import spamvec as _sv
from . import povm as _povm
from . import instrument as _instrument
from . import labeldicts as _ld
from . import gaugegroup as _gg
from . import matrixforwardsim as _matrixfwdsim
from . import mapforwardsim as _mapfwdsim
from . import termforwardsim as _termfwdsim
from . import explicitcalc as _explicitcalc
from ..baseobjs import VerbosityPrinter as _VerbosityPrinter
from ..baseobjs import Basis as _Basis
from ..baseobjs import BuiltinBasis as _BuiltinBasis
from ..baseobjs import Label as _Label
class Model(object):
"""
A predictive model for a Quantum Information Processor (QIP).
The main function of a `Model` object is to compute the outcome
probabilities of :class:`Circuit` objects based on the action of the
model's ideal operations plus (potentially) noise which makes the
outcome probabilities deviate from the perfect ones.
"""
def __init__(self, state_space_labels):
"""
Creates a new Model. Rarely used except from derived classes
`__init__` functions.
Parameters
----------
state_space_labels : StateSpaceLabels or list or tuple
The decomposition (with labels) of (pure) state-space this model
acts upon. Regardless of whether the model contains operators or
superoperators, this argument describes the Hilbert space dimension
and imposed structure. If a list or tuple is given, it must be
of a from that can be passed to `StateSpaceLabels.__init__`.
"""
if isinstance(state_space_labels, _ld.StateSpaceLabels):
self._state_space_labels = state_space_labels
else:
self._state_space_labels = _ld.StateSpaceLabels(state_space_labels)
self._hyperparams = {}
self._paramvec = _np.zeros(0, 'd')
self._paramlbls = None # a placeholder for FUTURE functionality
self.uuid = _uuid.uuid4() # a Model's uuid is like a persistent id(), useful for hashing
@property
def state_space_labels(self):
""" State space labels """
return self._state_space_labels
@property
def hyperparams(self):
""" Dictionary of hyperparameters associated with this model """
return self._hyperparams # Note: no need to set this param - just set/update values
def num_params(self):
"""
Return the number of free parameters when vectorizing
this model.
Returns
-------
int
the number of model parameters.
"""
return len(self._paramvec)
def to_vector(self):
"""
Returns the model vectorized according to the optional parameters.
Returns
-------
numpy array
The vectorized model parameters.
"""
return self._paramvec
def from_vector(self, v, reset_basis=False):
"""
The inverse of to_vector. Loads values of gates and rho and E vecs from
from the vector `v`. Note that `v` does not specify the number of
gates, etc., and their labels: this information must be contained in
this `Model` prior to calling `from_vector`. In practice, this just
means you should call the `from_vector` method using the same `Model`
that was used to generate the vector `v` in the first place.
"""
assert(len(v) == self.num_params())
self._paramvec = v.copy()
def probs(self, circuit, clipTo=None):
"""
Construct a dictionary containing the probabilities of every spam label
given a operation sequence.
Parameters
----------
circuit : Circuit or tuple of operation labels
The sequence of operation labels specifying the operation sequence.
clipTo : 2-tuple, optional
(min,max) to clip probabilities to if not None.
Returns
-------
probs : dictionary
A dictionary such that
probs[SL] = pr(SL,circuit,clipTo)
for each spam label (string) SL.
"""
raise NotImplementedError("Derived classes should implement this!")
def dprobs(self, circuit, returnPr=False, clipTo=None):
"""
Construct a dictionary containing the probability derivatives of every
spam label for a given operation sequence.
Parameters
----------
circuit : Circuit or tuple of operation labels
The sequence of operation labels specifying the operation sequence.
returnPr : bool, optional
when set to True, additionally return the probabilities.
clipTo : 2-tuple, optional
(min,max) to clip returned probability to if not None.
Only relevant when returnPr == True.
Returns
-------
dprobs : dictionary
A dictionary such that
dprobs[SL] = dpr(SL,circuit,gates,G0,SPAM,SP0,returnPr,clipTo)
for each spam label (string) SL.
"""
#Finite difference default?
raise NotImplementedError("Derived classes should implement this!")
def hprobs(self, circuit, returnPr=False, returnDeriv=False, clipTo=None):
"""
Construct a dictionary containing the probability derivatives of every
spam label for a given operation sequence.
Parameters
----------
circuit : Circuit or tuple of operation labels
The sequence of operation labels specifying the operation sequence.
returnPr : bool, optional
when set to True, additionally return the probabilities.
returnDeriv : bool, optional
when set to True, additionally return the derivatives of the
probabilities.
clipTo : 2-tuple, optional
(min,max) to clip returned probability to if not None.
Only relevant when returnPr == True.
Returns
-------
hprobs : dictionary
A dictionary such that
hprobs[SL] = hpr(SL,circuit,gates,G0,SPAM,SP0,returnPr,returnDeriv,clipTo)
for each spam label (string) SL.
"""
raise NotImplementedError("Derived classes should implement this!")
def bulk_evaltree_from_resources(self, circuit_list, comm=None, memLimit=None,
distributeMethod="default", subcalls=[],
dataset=None, verbosity=0):
raise NotImplementedError("Derived classes should implement this!")
#return circuit_list # MORE?
def bulk_evaltree(self, circuit_list, minSubtrees=None, maxTreeSize=None,
numSubtreeComms=1, dataset=None, verbosity=0):
raise NotImplementedError("Derived classes should implement this!")
#return circuit_list # MORE?
#def uses_evaltrees(self):
# """
# Whether or not this model uses evaluation trees to compute many
# (bulk) probabilities and their derivatives.
#
# Returns
# -------
# bool
# """
# return False
def bulk_probs(self, circuit_list, clipTo=None, check=False,
comm=None, memLimit=None, dataset=None, smartc=None):
raise NotImplementedError("Derived classes should implement this!")
def bulk_dprobs(self, circuit_list, returnPr=False, clipTo=None,
check=False, comm=None, wrtBlockSize=None, dataset=None):
raise NotImplementedError("Derived classes should implement this!")
def bulk_hprobs(self, circuit_list, returnPr=False, returnDeriv=False,
clipTo=None, check=False, comm=None,
wrtBlockSize1=None, wrtBlockSize2=None, dataset=None):
raise NotImplementedError("Derived classes should implement this!")
def bulk_fill_probs(self, mxToFill, evalTree, clipTo=None, check=False, comm=None):
raise NotImplementedError("Derived classes should implement this!")
def bulk_fill_dprobs(self, mxToFill, evalTree, prMxToFill=None, clipTo=None,
check=False, comm=None, wrtBlockSize=None,
profiler=None, gatherMemLimit=None):
raise NotImplementedError("Derived classes should implement this!")
def bulk_fill_hprobs(self, mxToFill, evalTree=None,
prMxToFill=None, derivMxToFill=None,
clipTo=None, check=False, comm=None,
wrtBlockSize1=None, wrtBlockSize2=None,
gatherMemLimit=None):
raise NotImplementedError("Derived classes should implement this!")
def bulk_hprobs_by_block(self, evalTree, wrtSlicesList,
bReturnDProbs12=False, comm=None):
raise NotImplementedError("Derived classes should implement this!")
def _init_copy(self, copyInto):
"""
Copies any "tricky" member of this model into `copyInto`, before
deep copying everything else within a .copy() operation.
"""
copyInto.uuid = _uuid.uuid4() # new uuid for a copy (don't duplicate!)
def copy(self):
"""
Copy this model.
Returns
-------
Model
a (deep) copy of this model.
"""
#Avoid having to reconstruct everything via __init__;
# essentially deepcopy this object, but give the
# class opportunity to initialize tricky members instead
# of letting deepcopy do it.
newModel = type(self).__new__(self.__class__) # empty object
#first call _init_copy to initialize any tricky members
# (like those that contain references to self or other members)
self._init_copy(newModel)
for attr, val in self.__dict__.items():
if not hasattr(newModel, attr):
assert(attr != "uuid"), "Should not be copying UUID!"
setattr(newModel, attr, _copy.deepcopy(val))
return newModel
def __str__(self):
pass
def __hash__(self):
if self.uuid is not None:
return hash(self.uuid)
else:
raise TypeError('Use digest hash')
class OpModel(Model):
"""
A Model containing operators (i.e. "members") which are independently
(sort of) parameterized and can be thought to have dense representations
(even if they're not actually stored that way). This gives rise to the
model having `basis` and `evotype` members.
Secondly, attached to an `OpModel` is the idea of "circuit simplification"
whereby the operators (preps, operations, povms, instruments) within
a circuit get simplified to things corresponding to a single outcome
probability, i.e. pseudo-circuits containing just preps, operations,
and POMV effects.
Thirdly, an `OpModel` is assumed to use a *layer-by-layer* evolution, and,
because of circuit simplification process, the calculaton of circuit
outcome probabilities has been pushed to a :class:`ForwardSimulator`
object which just deals with the forward simulation of simplified circuits.
Furthermore, instead of relying on a static set of operations a forward
simulator queries a :class:`LayerLizard` for layer operations, making it
possible to build up layer operations in an on-demand fashion from pieces
within the model.
"""
#Whether to perform extra parameter-vector integrity checks
_pcheck = False
def __init__(self, state_space_labels, basis, evotype, simplifier_helper, sim_type="auto"):
"""
Creates a new OpModel. Rarely used except from derived classes
`__init__` functions.
Parameters
----------
state_space_labels : StateSpaceLabels or list or tuple
The decomposition (with labels) of (pure) state-space this model
acts upon. Regardless of whether the model contains operators or
superoperators, this argument describes the Hilbert space dimension
and imposed structure. If a list or tuple is given, it must be
of a from that can be passed to `StateSpaceLabels.__init__`.
basis : Basis
The basis used for the state space by dense operator representations.
evotype : {"densitymx", "statevec", "stabilizer", "svterm", "cterm"}
The evolution type of this model, describing how states are
represented, allowing compatibility checks with (super)operator
objects.
simplifier_helper : SimplifierHelper
Provides a minimal interface for compiling circuits for forward
simulation.
sim_type : {"auto", "matrix", "map", "termorder:X"}
The type of forward simulator this model should use. `"auto"`
tries to determine the best type automatically.
"""
self._evotype = evotype
self.set_state_space(state_space_labels, basis)
#sets self._state_space_labels, self._basis, self._dim
self.set_simtype(sim_type)
#sets self._calcClass, self._sim_type, self._sim_args
self._shlp = simplifier_helper
self._need_to_rebuild = True # whether we call _rebuild_paramvec() in to_vector() or num_params()
self.dirty = False # indicates when objects and _paramvec may be out of sync
super(OpModel, self).__init__(self.state_space_labels)
##########################################
## Get/Set methods
##########################################
@property
def simtype(self):
""" Forward simulation type """
return self._sim_type
@property
def evotype(self):
""" Evolution type """
return self._evotype
@property
def basis(self):
""" The basis used to represent dense (super)operators of this model """
return self._basis
@basis.setter
def basis(self, basis):
if isinstance(basis, _Basis):
assert(basis.dim == self.state_space_labels.dim), \
"Cannot set basis w/dim=%d when sslbls dim=%d!" % (basis.dim, self.state_space_labels.dim)
self._basis = basis
else: # create a basis with the proper structure & dimension
self._basis = _Basis.cast(basis, self.state_space_labels)
def set_simtype(self, sim_type, calc_cache=None, max_cache_size=None):
"""
Reset the forward simulation type of this model.
Parameters
----------
sim_type : {"auto", "matrix", "map", "termorder:X"}
The type of forward simulator this model should use. `"auto"`
tries to determine the best type automatically.
calc_cache : dict or None
A cache of pre-computed values used in Taylor-term-based forward
simulation.
Returns
-------
None
"""
#Calculator selection based on simulation type
if sim_type == "auto":
default_param = self.operations.default_param # assume the same for other dicts
if _gt.is_valid_lindblad_paramtype(default_param) and \
_gt.split_lindblad_paramtype(default_param)[1] in ("svterm", "cterm"):
sim_type = "termorder:1"
else:
d = self._dim if (self._dim is not None) else 0
sim_type = "matrix" if d <= 16 else "map"
simtype_and_args = sim_type.split(":")
sim_type = simtype_and_args[0]
if sim_type == "matrix": c = _matrixfwdsim.MatrixForwardSimulator
elif sim_type == "map": c = _mapfwdsim.MapForwardSimulator
elif sim_type in ("termorder", "termgap", "termdirect"): c = _termfwdsim.TermForwardSimulator
else: raise ValueError("Invalid `sim_type` (%s)" % sim_type)
self._calcClass = c
self._sim_type = sim_type
self._sim_args = list(simtype_and_args[1:])
if sim_type.startswith("term"):
#cache = calc_cache if (calc_cache is not None) else {} # make a temp cache if none is given
cache = calc_cache # allow None cache to indicate *direct* computation of terms (no polys)
self._sim_args.append(cache) # add calculation cache as another argument
elif sim_type == "map":
self._sim_args.append(max_cache_size) # add cache size as another argument
#TODO REMOVE
#def reset_basis(self):
# """
# "Forgets" the current basis, so that
# self.basis becomes a dummy Basis w/name "unknown".
# """
# self._basis = _BuiltinBasis('unknown', 0)
def set_state_space(self, lbls, basis="pp"):
"""
Sets labels for the components of the Hilbert space upon which
the gates of this Model act.
Parameters
----------
lbls : list or tuple or StateSpaceLabels object
A list of state-space labels (can be strings or integers), e.g.
`['Q0','Q1']` or a :class:`StateSpaceLabels` object.
basis : Basis or str
A :class:`Basis` object or a basis name (like `"pp"`), specifying
the basis used to interpret the operators in this Model. If a
`Basis` object, then its dimensions must match those of `lbls`.
Returns
-------
None
"""
if isinstance(lbls, _ld.StateSpaceLabels):
self._state_space_labels = lbls
else:
self._state_space_labels = _ld.StateSpaceLabels(lbls, evotype=self._evotype)
self.basis = basis # invokes basis setter to set self._basis
#Operator dimension of this Model
self._dim = self.state_space_labels.dim
#e.g. 4 for 1Q (densitymx) or 2 for 1Q (statevec)
@property
def dim(self):
"""
The dimension of the model, which equals d when the gate
matrices have shape d x d and spam vectors have shape d x 1.
Returns
-------
int
model dimension
"""
return self._dim
def get_dimension(self):
"""
Get the dimension of the model, which equals d when the gate
matrices have shape d x d and spam vectors have shape d x 1.
Equivalent to model.dim.
Returns
-------
int
model dimension
"""
return self._dim
####################################################
## Parameter vector maintenance
####################################################
def num_params(self):
"""
Return the number of free parameters when vectorizing
this model.
Returns
-------
int
the number of model parameters.
"""
self._clean_paramvec()
return len(self._paramvec)
def _iter_parameterized_objs(self):
raise NotImplementedError("Derived Model classes should implement _iter_parameterized_objs")
#return # default is to have no parameterized objects
def _check_paramvec(self, debug=False):
if debug: print("---- Model._check_paramvec ----")
TOL = 1e-8
for lbl, obj in self._iter_parameterized_objs():
if debug: print(lbl, ":", obj.num_params(), obj.gpindices)
w = obj.to_vector()
msg = "None" if (obj.parent is None) else id(obj.parent)
assert(obj.parent is self), "%s's parent is not set correctly (%s)!" % (lbl, msg)
if obj.gpindices is not None and len(w) > 0:
if _np.linalg.norm(self._paramvec[obj.gpindices] - w) > TOL:
if debug: print(lbl, ".to_vector() = ", w, " but Model's paramvec = ",
self._paramvec[obj.gpindices])
raise ValueError("%s is out of sync with paramvec!!!" % lbl)
if not self.dirty and obj.dirty:
raise ValueError("%s is dirty but Model.dirty=False!!" % lbl)
def _clean_paramvec(self):
""" Updates _paramvec corresponding to any "dirty" elements, which may
have been modified without out knowing, leaving _paramvec out of
sync with the element's internal data. It *may* be necessary
to resolve conflicts where multiple dirty elements want different
values for a single parameter. This method is used as a safety net
that tries to insure _paramvec & Model elements are consistent
before their use."""
#print("Cleaning Paramvec (dirty=%s, rebuild=%s)" % (self.dirty, self._need_to_rebuild))
if self._need_to_rebuild:
self._rebuild_paramvec()
self._need_to_rebuild = False
if self.dirty: # if any member object is dirty (ModelMember.dirty setter should set this value)
TOL = 1e-8
#Note: lbl args used *just* for potential debugging - could strip out once
# we're confident this code always works.
def clean_single_obj(obj, lbl): # sync an object's to_vector result w/_paramvec
if obj.dirty:
w = obj.to_vector()
chk_norm = _np.linalg.norm(self._paramvec[obj.gpindices] - w)
#print(lbl, " is dirty! vec = ", w, " chk_norm = ",chk_norm)
if (not _np.isfinite(chk_norm)) or chk_norm > TOL:
self._paramvec[obj.gpindices] = w
obj.dirty = False
def clean_obj(obj, lbl): # recursive so works with objects that have sub-members
for i, subm in enumerate(obj.submembers()):
clean_obj(subm, _Label(lbl.name + ":%d" % i, lbl.sslbls))
clean_single_obj(obj, lbl)
def reset_dirty(obj): # recursive so works with objects that have sub-members
for i, subm in enumerate(obj.submembers()): reset_dirty(subm)
obj.dirty = False
for lbl, obj in self._iter_parameterized_objs():
clean_obj(obj, lbl)
#re-update everything to ensure consistency ~ self.from_vector(self._paramvec)
#print("DEBUG: non-trivially CLEANED paramvec due to dirty elements")
for _, obj in self._iter_parameterized_objs():
obj.from_vector(self._paramvec[obj.gpindices])
reset_dirty(obj) # like "obj.dirty = False" but recursive
#object is known to be consistent with _paramvec
if OpModel._pcheck: self._check_paramvec()
def _mark_for_rebuild(self, modified_obj=None):
#re-initialze any members that also depend on the updated parameters
self._need_to_rebuild = True
for _, o in self._iter_parameterized_objs():
if o._obj_refcount(modified_obj) > 0:
o.clear_gpindices() # ~ o.gpindices = None but works w/submembers
# (so params for this obj will be rebuilt)
self.dirty = True
#since it's likely we'll set at least one of our object's .dirty flags
# to True (and said object may have parent=None and so won't
# auto-propagate up to set this model's dirty flag (self.dirty)
def _print_gpindices(self):
print("PRINTING MODEL GPINDICES!!!")
for lbl, obj in self._iter_parameterized_objs():
print("LABEL ", lbl)
obj._print_gpindices()
def _rebuild_paramvec(self):
""" Resizes self._paramvec and updates gpindices & parent members as needed,
and will initialize new elements of _paramvec, but does NOT change
existing elements of _paramvec (use _update_paramvec for this)"""
v = self._paramvec; Np = len(self._paramvec) # NOT self.num_params() since the latter calls us!
off = 0; shift = 0
#print("DEBUG: rebuilding...")
#Step 1: remove any unused indices from paramvec and shift accordingly
used_gpindices = set()
for _, obj in self._iter_parameterized_objs():
if obj.gpindices is not None:
assert(obj.parent is self), "Member's parent is not set correctly (%s)!" % str(obj.parent)
used_gpindices.update(obj.gpindices_as_array())
else:
assert(obj.parent is self or obj.parent is None)
#Note: ok for objects to have parent == None if their gpindices is also None
indices_to_remove = sorted(set(range(Np)) - used_gpindices)
if len(indices_to_remove) > 0:
#print("DEBUG: Removing %d params:" % len(indices_to_remove), indices_to_remove)
v = _np.delete(v, indices_to_remove)
def get_shift(j): return _bisect.bisect_left(indices_to_remove, j)
memo = set() # keep track of which object's gpindices have been set
for _, obj in self._iter_parameterized_objs():
if obj.gpindices is not None:
if id(obj) in memo: continue # already processed
if isinstance(obj.gpindices, slice):
new_inds = _slct.shift(obj.gpindices,
-get_shift(obj.gpindices.start))
else:
new_inds = []
for i in obj.gpindices:
new_inds.append(i - get_shift(i))
new_inds = _np.array(new_inds, _np.int64)
obj.set_gpindices(new_inds, self, memo)
# Step 2: add parameters that don't exist yet
# Note that iteration order (that of _iter_parameterized_objs) determines
# parameter index ordering, so "normally" an object that occurs before
# another in the iteration order will have gpindices which are lower - and
# when new indices are allocated we try to maintain this normal order by
# inserting them at an appropriate place in the parameter vector.
# off : holds the current point where new params should be inserted
# shift : holds the amount existing parameters that are > offset (not in `memo`) should be shifted
# Note: Adding more explicit "> offset" logic may obviate the need for the memo arg?
memo = set() # keep track of which object's gpindices have been set
for lbl, obj in self._iter_parameterized_objs():
if shift > 0 and obj.gpindices is not None:
if isinstance(obj.gpindices, slice):
obj.set_gpindices(_slct.shift(obj.gpindices, shift), self, memo)
else:
obj.set_gpindices(obj.gpindices + shift, self, memo) # works for integer arrays
if obj.gpindices is None or obj.parent is not self:
#Assume all parameters of obj are new independent parameters
num_new_params = obj.allocate_gpindices(off, self, memo)
objvec = obj.to_vector() # may include more than "new" indices
if num_new_params > 0:
new_local_inds = _gm._decompose_gpindices(obj.gpindices, slice(off, off + num_new_params))
assert(len(objvec[new_local_inds]) == num_new_params)
v = _np.insert(v, off, objvec[new_local_inds])
# print("objvec len = ",len(objvec), "num_new_params=",num_new_params,
# " gpinds=",obj.gpindices) #," loc=",new_local_inds)
#obj.set_gpindices( slice(off, off+obj.num_params()), self )
#shift += obj.num_params()
#off += obj.num_params()
shift += num_new_params
off += num_new_params
#print("DEBUG: %s: alloc'd & inserted %d new params. indices = " \
# % (str(lbl),obj.num_params()), obj.gpindices, " off=",off)
else:
inds = obj.gpindices_as_array()
M = max(inds) if len(inds) > 0 else -1; L = len(v)
#print("DEBUG: %s: existing indices = " % (str(lbl)), obj.gpindices, " M=",M," L=",L)
if M >= L:
#Some indices specified by obj are absent, and must be created.
w = obj.to_vector()
v = _np.concatenate((v, _np.empty(M + 1 - L, 'd')), axis=0) # [v.resize(M+1) doesn't work]
shift += M + 1 - L
for ii, i in enumerate(inds):
if i >= L: v[i] = w[ii]
#print("DEBUG: --> added %d new params" % (M+1-L))
if M >= 0: # M == -1 signifies this object has no parameters, so we'll just leave `off` alone
off = M + 1
self._paramvec = v
#print("DEBUG: Done rebuild: %d params" % len(v))
def _init_virtual_obj(self, obj):
"""
Initializes a "virtual object" - an object (e.g. LinearOperator) that *could* be a
member of the Model but won't be, as it's just built for temporary
use (e.g. the parallel action of several "base" gates). As such
we need to fully initialize its parent and gpindices members so it
knows it belongs to this Model BUT it's not allowed to add any new
parameters (they'd just be temporary). It's also assumed that virtual
objects don't need to be to/from-vectored as there are already enough
real (non-virtual) gates/spamvecs/etc. to accomplish this.
"""
if obj.gpindices is not None:
assert(obj.parent is self), "Virtual obj has incorrect parent already set!"
return # if parent is already set we assume obj has already been init
#Assume all parameters of obj are new independent parameters
num_new_params = obj.allocate_gpindices(self.num_params(), self)
assert(num_new_params == 0), "Virtual object is requesting %d new params!" % num_new_params
def _obj_refcount(self, obj):
""" Number of references to `obj` contained within this Model """
cnt = 0
for _, o in self._iter_parameterized_objs():
cnt += o._obj_refcount(obj)
return cnt
def to_vector(self):
"""
Returns the model vectorized according to the optional parameters.
Returns
-------
numpy array
The vectorized model parameters.
"""
self._clean_paramvec() # will rebuild if needed
return self._paramvec
def from_vector(self, v):
"""
The inverse of to_vector. Loads values of gates and rho and E vecs from
from the vector `v`. Note that `v` does not specify the number of
gates, etc., and their labels: this information must be contained in
this `Model` prior to calling `from_vector`. In practice, this just
means you should call the `from_vector` method using the same `Model`
that was used to generate the vector `v` in the first place.
"""
assert(len(v) == self.num_params())
self._paramvec = v.copy()
for _, obj in self._iter_parameterized_objs():
obj.from_vector(v[obj.gpindices])
obj.dirty = False # object is known to be consistent with _paramvec
#if reset_basis:
# self.reset_basis()
# assume the vector we're loading isn't producing gates & vectors in
# a known basis.
if OpModel._pcheck: self._check_paramvec()
######################################
## Compilation
######################################
def _layer_lizard(self):
""" Return a layer lizard for this model """
raise NotImplementedError("Derived Model classes should implement this!")
def _fwdsim(self):
""" Create & return a forward-simulator ("calculator") for this model """
self._clean_paramvec()
layer_lizard = self._layer_lizard()
kwargs = {}
if self._sim_type == "termorder":
assert(len(self._sim_args) == 1 + 1), "termorder must have <order> arg, e.g. 'termorder:1'"
kwargs['mode'] = "taylor-order"
kwargs['max_order'] = int(self._sim_args[0])
kwargs['cache'] = self._sim_args[-1] # always the last argument
if self._sim_type in ("termgap", "termdirect"):
assert(len(self._sim_args) >= 3 + 1), \
"%s must have <max-order>, <gap> and <min> args, e.g. '%s:3:0.1:0.01'" \
% (self._sim_type, self._sim_type)
kwargs['mode'] = "pruned" if (self._sim_type == "termgap") else "direct"
kwargs['max_order'] = int(self._sim_args[0])
kwargs['pathmag_gap'] = float(self._sim_args[1])
kwargs['min_term_mag'] = float(self._sim_args[2])
kwargs['opt_mode'] = bool(self._sim_args[3]) if len(self._sim_args) > 3 else False
# indicates fwdsim is being used within an optimization loop (only recomp paths on deriv evals)
kwargs['cache'] = self._sim_args[-1] # always the last argument
if self._sim_type == "map":
kwargs['max_cache_size'] = self._sim_args[0] if len(self._sim_args) > 0 else None # backward compat
assert(self._calcClass is not None), "Model does not have a calculator setup yet!"
return self._calcClass(self._dim, layer_lizard, self._paramvec, **kwargs) # fwdsim class
def split_circuit(self, circuit, erroron=('prep', 'povm')):
"""
Splits a operation sequence into prepLabel + opsOnlyString + povmLabel
components. If `circuit` does not contain a prep label or a
povm label a default label is returned if one exists.
Parameters
----------
circuit : Circuit
A operation sequence, possibly beginning with a state preparation
label and ending with a povm label.
erroron : tuple of {'prep','povm'}
A ValueError is raised if a preparation or povm label cannot be
resolved when 'prep' or 'povm' is included in 'erroron'. Otherwise
`None` is returned in place of unresolvable labels. An exception
is when this model has no preps or povms, in which case `None`
is always returned and errors are never raised, since in this
case one usually doesn't expect to use the Model to compute
probabilities (e.g. in germ selection).
Returns
-------
prepLabel : str or None
opsOnlyString : Circuit
povmLabel : str or None
"""
if len(circuit) > 0 and self._shlp.is_prep_lbl(circuit[0]):
prep_lbl = circuit[0]
circuit = circuit[1:]
elif self._shlp.get_default_prep_lbl() is not None:
prep_lbl = self._shlp.get_default_prep_lbl()
else:
if 'prep' in erroron and self._shlp.has_preps():
raise ValueError("Cannot resolve state prep in %s" % circuit)
else: prep_lbl = None
if len(circuit) > 0 and self._shlp.is_povm_lbl(circuit[-1]):
povm_lbl = circuit[-1]
circuit = circuit[:-1]
elif self._shlp.get_default_povm_lbl() is not None:
povm_lbl = self._shlp.get_default_povm_lbl()
else:
if 'povm' in erroron and self._shlp.has_povms():
raise ValueError("Cannot resolve POVM in %s" % circuit)
else: povm_lbl = None
return prep_lbl, circuit, povm_lbl
def simplify_circuits(self, circuits, dataset=None):
"""
Simplifies a list of :class:`Circuit`s.
Circuits must be "simplified" before probabilities can be computed for
them. Each string corresponds to some number of "outcomes", indexed by an
"outcome label" that is a tuple of POVM-effect or instrument-element
labels like "0". Compiling creates maps between operation sequences and their
outcomes and the structures used in probability computation (see return
values below).
Parameters
----------
circuits : list of Circuits
The list to simplify.
dataset : DataSet, optional
If not None, restrict what is simplified to only those
probabilities corresponding to non-zero counts (observed
outcomes) in this data set.
Returns
-------
raw_spamTuples_dict : collections.OrderedDict
A dictionary whose keys are raw operation sequences (containing just
"simplified" gates, i.e. not instruments), and whose values are
lists of (preplbl, effectlbl) tuples. The effectlbl names a
"simplified" effect vector; preplbl is just a prep label. Each tuple
corresponds to a single "final element" of the computation, e.g. a
probability. The ordering is important - and is why this needs to be
an ordered dictionary - when the lists of tuples are concatenated (by
key) the resulting tuple orderings corresponds to the final-element
axis of an output array that is being filled (computed).
elIndices : collections.OrderedDict
A dictionary whose keys are integer indices into `circuits` and
whose values are slices and/or integer-arrays into the space/axis of
final elements. Thus, to get the final elements corresponding to
`circuits[i]`, use `filledArray[ elIndices[i] ]`.
outcomes : collections.OrderedDict
A dictionary whose keys are integer indices into `circuits` and
whose values are lists of outcome labels (an outcome label is a tuple
of POVM-effect and/or instrument-element labels). Thus, to obtain
what outcomes the i-th operation sequences's final elements
(`filledArray[ elIndices[i] ]`) correspond to, use `outcomes[i]`.
nTotElements : int
The total number of "final elements" - this is how big of an array
is need to hold all of the probabilities `circuits` generates.
"""
# model.simplify -> odict[raw_gstr] = spamTuples, elementIndices, nElements
# dataset.simplify -> outcomeLabels[i] = list_of_ds_outcomes, elementIndices, nElements
# simplify all gsplaq strs -> elementIndices[(i,j)],
circuits = [(opstr if isinstance(opstr, _cir.Circuit) else _cir.Circuit(opstr))
for opstr in circuits] # cast to Circuits
#Indexed by raw operation sequence
raw_spamTuples_dict = _collections.OrderedDict() # final
raw_opOutcomes_dict = _collections.OrderedDict()
raw_offsets = _collections.OrderedDict()
#Indexed by parent index (an integer)
elIndicesByParent = _collections.OrderedDict() # final
outcomesByParent = _collections.OrderedDict() # final
elIndsToOutcomesByParent = _collections.OrderedDict()
def resolveSPAM(circuit):
""" Determines spam tuples that correspond to circuit
and strips any spam-related pieces off """
prep_lbl, circuit, povm_lbl = \
self.split_circuit(circuit)
if prep_lbl is None or povm_lbl is None:
spamtups = [None] # put a single "dummy" spam-tuple placeholder
# so that there's a single "element" for each simplified string,
# which means that the usual "lookup" or "elIndices" will map
# original circuit-list indices to simplified-string, i.e.,
# evalTree index, which is useful when computing products
# (often the case when a Model has no preps or povms,
# e.g. in germ selection)
else:
if dataset is not None:
#Then we don't need to consider *all* possible spam tuples -
# just the ones that are observed, i.e. that correspond to
# a final element in the "full" (tuple) outcome labels that
# were observed.
observed_povm_outcomes = sorted(set(
[full_out_tup[-1] for full_out_tup in dataset[circuit].outcomes]))
spamtups = [(prep_lbl, povm_lbl + "_" + oout)
for oout in observed_povm_outcomes]
# elbl = oout[-1] -- the last element corresponds
# to the POVM (earlier ones = instruments)
else:
spamtups = [(prep_lbl, povm_lbl + "_" + elbl)
for elbl in self._shlp.get_effect_labels_for_povm(povm_lbl)]
return circuit, spamtups
def process(s, spamtuples, observed_outcomes, elIndsToOutcomes,
op_outcomes=(), start=0):
"""
Implements recursive processing of a string. Separately
implements two different behaviors:
"add" : add entries to raw_spamTuples_dict and raw_opOutcomes_dict
"index" : adds entries to elIndicesByParent and outcomesByParent
assuming that raw_spamTuples_dict and raw_opOutcomes_dict
are already build (and won't be modified anymore).
"""
sub = s if start == 0 else s[start:]
for i, op_label in enumerate(sub, start=start):
# OLD: now allow "gate-level" labels which can contain
# multiple (parallel) instrument labels
#if op_label in self.instruments:
# #we've found an instrument - recurse!
# for inst_el_lbl in self.instruments[op_label]:
# simplified_el_lbl = op_label + "_" + inst_el_lbl
# process(s[0:i] + _cir.Circuit((simplified_el_lbl,)) + s[i+1:],
# spamtuples, elIndsToOutcomes, op_outcomes + (inst_el_lbl,), i+1)
# break
if any([self._shlp.is_instrument_lbl(sub_gl) for sub_gl in op_label.components]):
# we've found an instrument - recurse!
sublabel_tups_to_iter = [] # one per label component (may be only 1)
for sub_gl in op_label.components:
if self._shlp.is_instrument_lbl(sub_gl):
sublabel_tups_to_iter.append(
[(sub_gl, inst_el_lbl)
for inst_el_lbl in self._shlp.get_member_labels_for_instrument(sub_gl)])
else:
sublabel_tups_to_iter.append([(sub_gl, None)]) # just a single element
for sublabel_tups in _itertools.product(*sublabel_tups_to_iter):
sublabels = [] # the sub-labels of the overall operation label to add
outcomes = [] # the outcome tuple associated with this overall label
for sub_gl, inst_el_lbl in sublabel_tups:
if inst_el_lbl is not None:
sublabels.append(sub_gl + "_" + inst_el_lbl)
outcomes.append(inst_el_lbl)
else:
sublabels.append(sub_gl)
simplified_el_lbl = _Label(sublabels)
simplified_el_outcomes = tuple(outcomes)