From 368675b8937f63feebd9ff04a51a72db7de19a24 Mon Sep 17 00:00:00 2001
From: Katharine Hyatt <67932820+kshyatt-aws@users.noreply.github.com>
Date: Thu, 25 Apr 2024 11:24:47 -0400
Subject: [PATCH] feature: Local simulator support (#72)

feature: Initial implementation of pure Julia local statevector simulator
---
 Project.toml                                  |  18 +-
 PyBraket/CondaPkg.toml                        |   4 +-
 PyBraket/Project.toml                         |   8 +-
 PyBraket/src/py_type_conversion.jl            |   2 +-
 PyBraket/src/pycircuit.jl                     |   2 +-
 PyBraket/test/ahs.jl                          |   4 +-
 PyBraket/test/circuits.jl                     |  30 +--
 PyBraket/test/gates.jl                        |  14 +-
 .../integ_tests/local_braket_simulator.jl     |   8 +-
 PyBraket/test/noise.jl                        |  16 +-
 docs/make.jl                                  |   2 +-
 docs/src/local_simulator.md                   |  13 ++
 src/Braket.jl                                 |  24 +-
 src/circuit.jl                                |  31 +--
 src/gate_applicators.jl                       |  38 ++--
 src/gates.jl                                  |  23 +-
 src/local_simulator.jl                        | 197 +++++++++++++++++
 src/moments.jl                                |  48 ++--
 src/observables.jl                            | 207 ++++++++++--------
 src/operators.jl                              |  37 +++-
 src/qubit_set.jl                              |   3 +-
 src/raw_schema.jl                             |   6 +-
 src/raw_task_result_types.jl                  |   2 +-
 src/results.jl                                |  18 +-
 src/schemas.jl                                |  40 ++--
 src/task.jl                                   |  20 +-
 src/utils.jl                                  |   8 +
 test/observables.jl                           |  23 +-
 28 files changed, 593 insertions(+), 253 deletions(-)
 create mode 100644 docs/src/local_simulator.md
 create mode 100644 src/local_simulator.jl

diff --git a/Project.toml b/Project.toml
index 2e5e80c9..821525ac 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Braket"
 uuid = "19504a0f-b47d-4348-9127-acc6cc69ef67"
 authors = ["Katharine Hyatt <hyatkath@amazon.com>"]
-version = "0.8.2"
+version = "0.8.3"
 
 [deps]
 AWS = "fbe9abb3-538b-5e4e-ba9e-bc94f4f92ebc"
@@ -26,6 +26,7 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NamedTupleTools = "d9ec5142-1e00-5aa0-9d6a-321866360f50"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StructTypes = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
 Tar = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
@@ -35,19 +36,19 @@ UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 AWS = "=1.90.4"
 AWSS3 = "=0.11.2"
 Aqua = "=0.8"
-Base64 = "1.6"
 AxisArrays = "=0.4.7"
-CSV = "=0.10.13"
+Base64 = "1.6"
+CSV = "=0.10.14"
 CodeTracking = "=1.3.5"
 Compat = "=4.14.0"
-DataStructures = "=0.18.18"
+DataStructures = "=0.18.20"
 Dates = "1.6"
 DecFP = "=1.3.2"
-Distributions = "=0.25.76"
 Distributed = "1.6"
+Distributions = "=0.25.76"
 Downloads = "1"
-Graphs = "=1.9.0"
-HTTP = "=1.10.3"
+Graphs = "=1.10.0"
+HTTP = "=1.10.6"
 InteractiveUtils = "1.6"
 JLD2 = "=0.4.46"
 JSON3 = "=1.14.0"
@@ -59,8 +60,9 @@ NamedTupleTools = "=0.14.3"
 OrderedCollections = "=1.6.3"
 Pkg = "1.6"
 Random = "1.6"
-Statistics = "1.6"
 SparseArrays = "1.6"
+StaticArrays = "=1.9.3"
+Statistics = "1.6"
 StructTypes = "=1.10.0"
 Tar = "1.9.3"
 Test = "1.6"
diff --git a/PyBraket/CondaPkg.toml b/PyBraket/CondaPkg.toml
index 4477a755..87c0f149 100644
--- a/PyBraket/CondaPkg.toml
+++ b/PyBraket/CondaPkg.toml
@@ -1,8 +1,8 @@
 [deps]
-python = ">=3.7,<=3.9"
+python = ">=3"
 pydantic = ""
 scipy = ""
 numpy = ""
 
 [pip.deps]
-amazon-braket-sdk = ">=1.35.0"
+amazon-braket-sdk = ">=1.70.0"
diff --git a/PyBraket/Project.toml b/PyBraket/Project.toml
index aa6bdd36..cdf76eba 100644
--- a/PyBraket/Project.toml
+++ b/PyBraket/Project.toml
@@ -1,7 +1,7 @@
 name = "PyBraket"
 uuid = "e85266a6-1825-490b-a80e-9b9469c53660"
 authors = ["Katharine Hyatt <hyatkath@amazon.com>"]
-version = "0.8.2"
+version = "0.8.3"
 
 [deps]
 Braket = "19504a0f-b47d-4348-9127-acc6cc69ef67"
@@ -14,11 +14,11 @@ StructTypes = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
 
 [compat]
 Aqua = "=0.8"
-Braket = "=0.8.2"
+Braket = "=0.8.3"
 CondaPkg = "=0.2.22"
-DataStructures = "=0.18.18"
+DataStructures = "=0.18.20"
 LinearAlgebra = "1.6"
-PythonCall = "=0.9.15"
+PythonCall = "=0.9.19"
 Statistics = "1"
 StructTypes = "=1.10.0"
 Test = "1.6"
diff --git a/PyBraket/src/py_type_conversion.jl b/PyBraket/src/py_type_conversion.jl
index 7d10fa19..4a12360d 100644
--- a/PyBraket/src/py_type_conversion.jl
+++ b/PyBraket/src/py_type_conversion.jl
@@ -54,7 +54,7 @@ function jl_convert_attr(n, t, attr)
             return pyconvert(t, attr)
         end
     else
-        PythonCall.pyisnone(attr) && return nothing
+        PythonCall.pyis(attr, PythonCall.pybuiltins.None) && return nothing
         return union_convert(t, attr)
     end
 end
diff --git a/PyBraket/src/pycircuit.jl b/PyBraket/src/pycircuit.jl
index ed349cc3..6724c6e1 100644
--- a/PyBraket/src/pycircuit.jl
+++ b/PyBraket/src/pycircuit.jl
@@ -6,7 +6,7 @@ mutable struct PyCircuit
 end
 (c::PyCircuit)(arg::Number; kwargs...)   = PyCircuit(Py(c)(arg; kwargs...))
 (c::PyCircuit)(; kwargs...)              = PyCircuit(Py(c)(; kwargs...))
-Base.:(==)(c1::PyCircuit, c2::PyCircuit) = PythonCall.pyisTrue(Py(c1) == Py(c2))
+Base.:(==)(c1::PyCircuit, c2::PyCircuit) = pyconvert(Bool, Py(c1) == Py(c2))
 Py(c::PyCircuit) = getfield(c, :o)
 Base.getproperty(c::PyCircuit, s::Symbol)             = pyconvert(Any, getproperty(Py(c), s))
 Base.getproperty(c::PyCircuit, s::AbstractString)     = pyconvert(Any, getproperty(Py(c), s))
diff --git a/PyBraket/test/ahs.jl b/PyBraket/test/ahs.jl
index 00ee1075..da68e540 100644
--- a/PyBraket/test/ahs.jl
+++ b/PyBraket/test/ahs.jl
@@ -43,7 +43,7 @@ using Braket: AtomArrangement, AtomArrangementItem, TimeSeries, DrivingField, Aw
     drive                   = DrivingField(Ω, ϕ, Δ)
     ahs_program             = AnalogHamiltonianSimulation(register, drive)
 
-    ahs_local    = LocalSimulator("braket_ahs")
+    ahs_local    = PyBraket.LocalSimulator("braket_ahs")
     local_result = result(run(ahs_local, ahs_program, shots=1_000))
     @test length(local_result.measurements) == 1_000
 
@@ -82,7 +82,7 @@ end
     drive                   = DrivingField(Ω, ϕ, Δ)
     ahs_program             = AnalogHamiltonianSimulation(register, drive)
 
-    ahs_local    = LocalSimulator("braket_ahs")
+    ahs_local    = PyBraket.LocalSimulator("braket_ahs")
     local_result = result(run(ahs_local, ahs_program, shots=1_000))
     
     g_count = 0
diff --git a/PyBraket/test/circuits.jl b/PyBraket/test/circuits.jl
index 3e8f4fb1..c2e8a2e2 100644
--- a/PyBraket/test/circuits.jl
+++ b/PyBraket/test/circuits.jl
@@ -1,11 +1,11 @@
 using Braket, PyBraket, Test, LinearAlgebra, PythonCall
-using PythonCall: pyconvert, Py, pyisTrue, pyisinstance
+using PythonCall: pyconvert, Py, pyisinstance
 @testset "PyBraket circuits" begin
     @testset for ir_type in (:JAQCD, :OpenQASM)
         Braket.IRType[] = ir_type 
         @testset "Expectation" begin
             c = Circuit([(H, 0), (CNot, 0, 1), (Expectation, Braket.Observables.Z(), 0)])
-            dev = LocalSimulator()
+            dev = PyBraket.LocalSimulator()
             braket_task = run(dev, c, shots=0)
             res = result(braket_task)
             @test res.values[1] ≈ 0.0 atol=1e-12
@@ -13,7 +13,7 @@ using PythonCall: pyconvert, Py, pyisTrue, pyisinstance
 
         @testset "Variance" begin
             c = Circuit() |> (ci->H(ci, 0)) |> (ci->CNot(ci, 0, 1)) |> (ci->Variance(ci, Braket.Observables.Z(), 0))
-            dev = LocalSimulator()
+            dev = PyBraket.LocalSimulator()
             braket_task = run(dev, c, shots=0)
             res = result(braket_task)
             @test res.values[1] ≈ 1.0 atol=1e-12
@@ -21,7 +21,7 @@ using PythonCall: pyconvert, Py, pyisTrue, pyisinstance
 
         @testset "Probability" begin
             c = Circuit() |> (ci->H(ci, 0)) |> (ci->CNot(ci, 0, 1)) |> (ci->Probability(ci))
-            dev = LocalSimulator()
+            dev = PyBraket.LocalSimulator()
             braket_task = run(dev, c, shots=0)
             res = result(braket_task)
             @test res.values[1] ≈ [0.5, 0.0, 0.0, 0.5] atol=1e-12
@@ -29,7 +29,7 @@ using PythonCall: pyconvert, Py, pyisTrue, pyisinstance
 
         @testset "DensityMatrix" begin
             c = Circuit() |> (ci->H(ci, 0)) |> (ci->CNot(ci, 0, 1)) |> (ci->DensityMatrix(ci))
-            dev = LocalSimulator()
+            dev = PyBraket.LocalSimulator()
             braket_task = run(dev, c, shots=0)
             res = result(braket_task)
             ρ = zeros(ComplexF64, 4, 4)
@@ -44,7 +44,7 @@ using PythonCall: pyconvert, Py, pyisTrue, pyisinstance
 
         @testset "Results type" begin
             c = Circuit() |> (ci->H(ci, 0)) |> (ci->CNot(ci, 0, 1)) |> (ci->Expectation(ci, Braket.Observables.Z(), 0))
-            dev = LocalSimulator()
+            dev = PyBraket.LocalSimulator()
             braket_task = run(dev, c, shots=0)
             res = result(braket_task)
             @test res.values[1] ≈ 0.0 atol=1e-12
@@ -53,7 +53,7 @@ using PythonCall: pyconvert, Py, pyisTrue, pyisinstance
 
         @testset "Observable on all qubits" begin
             c = Circuit() |> (ci->H(ci, 0)) |> (ci->CNot(ci, 0, 1)) |> (ci->Expectation(ci, Braket.Observables.Z()))
-            dev = LocalSimulator()
+            dev = PyBraket.LocalSimulator()
             res = result(run(dev, c, shots=0))
             @test pyconvert(Vector{Float64}, res.values[1]) ≈ [0.0, 0.0] atol=1e-12
         end
@@ -66,7 +66,7 @@ using PythonCall: pyconvert, Py, pyisTrue, pyisinstance
             Unitary(c, [0, 1, 2], ccz_mat)
             H(c, [0, 1, 2])
             StateVector(c)
-            dev = LocalSimulator()
+            dev = PyBraket.LocalSimulator()
             braket_task = run(dev, c, shots=0)
             res = result(braket_task)
             @test res.values[1] ≈ ComplexF64[0.75, 0.25, 0.25, -0.25, 0.25, -0.25, -0.25, 0.25] atol=1e-12
@@ -93,7 +93,7 @@ using PythonCall: pyconvert, Py, pyisTrue, pyisinstance
                     obs = rt(o, [0])
                     py_obs = Py(obs)
                     ir_obs = ir(obs)
-                    @test pyisTrue(py_obs.to_ir() == Py(ir_obs))
+                    @test pyconvert(Bool, py_obs.to_ir() == Py(ir_obs))
                     @test pyconvert(ir_rt, Py(ir_obs)) == ir_obs
                 end
                 @testset for (rt, ir_rt) in ((Braket.Probability, Braket.IR.Probability),
@@ -101,27 +101,27 @@ using PythonCall: pyconvert, Py, pyisTrue, pyisinstance
                     obs = rt([0])
                     py_obs = Py(obs)
                     ir_obs = ir(obs)
-                    @test pyisTrue(py_obs.to_ir() == Py(ir_obs))
+                    @test pyconvert(Bool, py_obs.to_ir() == Py(ir_obs))
                     @test pyconvert(ir_rt, Py(ir_obs)) == ir_obs
 
                     obs = rt()
                     py_obs = Py(obs)
                     ir_obs = ir(obs)
-                    @test pyisTrue(py_obs.to_ir() == Py(ir_obs))
+                    @test pyconvert(Bool, py_obs.to_ir() == Py(ir_obs))
                     @test pyconvert(ir_rt, Py(ir_obs)) == ir_obs
                 end
                 @testset "(rt, ir_rt) = (Amplitude, IR.Amplitude)" begin
                     obs = Braket.Amplitude(["0000"])
                     py_obs = Py(obs)
                     ir_obs = ir(obs)
-                    @test pyisTrue(py_obs.to_ir() == Py(ir_obs))
+                    @test pyconvert(Bool, py_obs.to_ir() == Py(ir_obs))
                     @test pyconvert(Braket.IR.Amplitude, Py(ir_obs)) == ir_obs
                 end
                 @testset "(rt, ir_rt) = (StateVector, IR.StateVector)" begin
                     obs = Braket.StateVector()
                     py_obs = Py(obs)
                     ir_obs = ir(obs)
-                    @test pyisTrue(py_obs.to_ir() == Py(ir_obs))
+                    @test pyconvert(Bool, py_obs.to_ir() == Py(ir_obs))
                     @test pyconvert(Braket.IR.StateVector, Py(ir_obs)) == ir_obs
                 end
                 @testset "HermitianObservable and TensorProduct" begin
@@ -136,7 +136,7 @@ using PythonCall: pyconvert, Py, pyisTrue, pyisinstance
                     o = 3.0 * Braket.Observables.Z()
                     py_obs = Py(o)
                     @test pyisinstance(py_obs, PyBraket.braketobs.Z)
-                    @test pyisTrue(py_obs.coefficient == 3.0)
+                    @test pyconvert(Bool, py_obs.coefficient == 3.0)
                 end
                 @testset "Sum" begin
                     m = [1. -im; im -1.]
@@ -168,7 +168,7 @@ using PythonCall: pyconvert, Py, pyisTrue, pyisinstance
             py_c2 = PyCircuit(non_para_circ)
             @test py_c2 == py_c1(1.0)
             @testset "running with inputs" begin 
-                dev = LocalSimulator()
+                dev = PyBraket.LocalSimulator()
                 oq3_circ = ir(circ, Val(:OpenQASM))
                 braket_task = run(dev, oq3_circ, shots=0, inputs=Dict(string(α)=>1.0, string(θ)=>2.0))
                 res = result(braket_task)
diff --git a/PyBraket/test/gates.jl b/PyBraket/test/gates.jl
index 392e5644..686ecf70 100644
--- a/PyBraket/test/gates.jl
+++ b/PyBraket/test/gates.jl
@@ -1,5 +1,5 @@
 using Test, PyBraket, Braket, Braket.IR, PythonCall
-using PythonCall: Py, pyconvert, pyisTrue
+using PythonCall: Py, pyconvert
 
 using Braket: I
 
@@ -20,7 +20,7 @@ using Braket: I
         g    = gate()
         ir_g = Braket.ir(g, 0)
         py_g = Py(g)
-        @test pyisTrue(py_g.to_ir([0]) == Py(ir_g))
+        @test pyconvert(Bool, py_g.to_ir([0]) == Py(ir_g))
         @test pyconvert(ir_gate, Py(ir_g)) == ir_g
     end
     @testset for (gate, ir_gate) in ((PhaseShift, IR.PhaseShift),
@@ -32,7 +32,7 @@ using Braket: I
         g     = gate(angle)
         ir_g  = Braket.ir(g, 0)
         py_g  = Py(g)
-        @test pyisTrue(py_g.to_ir([0]) == Py(ir_g))
+        @test pyconvert(Bool, py_g.to_ir([0]) == Py(ir_g))
         @test pyconvert(ir_gate, Py(ir_g)) == ir_g
     end
     @testset for (gate, ir_gate) in ((CNot, IR.CNot),
@@ -46,7 +46,7 @@ using Braket: I
         g    = gate()
         ir_g = Braket.ir(g, [0, 1])
         py_g = Py(g)
-        @test pyisTrue(py_g.to_ir([0, 1]) == Py(ir_g))
+        @test pyconvert(Bool, py_g.to_ir([0, 1]) == Py(ir_g))
         @test pyconvert(ir_gate, Py(ir_g)) == ir_g
     end
     @testset for (gate, ir_gate) in ((CPhaseShift, IR.CPhaseShift),
@@ -62,14 +62,14 @@ using Braket: I
         g     = gate(angle)
         ir_g  = Braket.ir(g, [0, 1])
         py_g  = Py(g)
-        @test pyisTrue(py_g.to_ir([0, 1]) == Py(ir_g))
+        @test pyconvert(Bool, py_g.to_ir([0, 1]) == Py(ir_g))
         @test pyconvert(ir_gate, Py(ir_g)) == ir_g
     end
     @testset for (gate, ir_gate) in ((CCNot, IR.CCNot), (CSwap, IR.CSwap))
         g    = gate()
         ir_g = Braket.ir(g, [0, 1, 2])
         py_g = Py(g)
-        @test pyisTrue(py_g.to_ir([0, 1, 2]) == Py(ir_g))
+        @test pyconvert(Bool, py_g.to_ir([0, 1, 2]) == Py(ir_g))
         @test pyconvert(ir_gate, Py(ir_g)) == ir_g
     end
     @testset "(gate, ir_gate) = (Unitary, IR.Unitary)" begin
@@ -77,7 +77,7 @@ using Braket: I
         n    = Unitary(mat)
         ir_n = Braket.ir(n, 0)
         py_n = Py(n)
-        @test pyisTrue(py_n.to_ir([0]) == Py(ir_n))
+        @test pyconvert(Bool, py_n.to_ir([0]) == Py(ir_n))
         @test pyconvert(IR.Unitary, Py(ir_n)) == ir_n
     end
     Braket.IRType[] = :OpenQASM
diff --git a/PyBraket/test/integ_tests/local_braket_simulator.jl b/PyBraket/test/integ_tests/local_braket_simulator.jl
index 440ff08e..93998fcd 100644
--- a/PyBraket/test/integ_tests/local_braket_simulator.jl
+++ b/PyBraket/test/integ_tests/local_braket_simulator.jl
@@ -1,8 +1,8 @@
 using Braket, Braket.Observables, LinearAlgebra, PyBraket, Statistics, Test
 using Braket: I
 
-PURE_DEVICE = LocalSimulator()
-NOISE_DEVICE = LocalSimulator("braket_dm")
+PURE_DEVICE = PyBraket.LocalSimulator()
+NOISE_DEVICE = PyBraket.LocalSimulator("braket_dm")
 SHOTS = 8000
 
 get_tol(shots::Int) = return (shots > 0 ? Dict("atol"=> 0.1, "rtol"=>0.15) : Dict("atol"=>0.01, "rtol"=>0))
@@ -29,7 +29,7 @@ end
     @testset for (backend, device_name) in [("default", "StateVectorSimulator"),
                                             ("braket_sv", "StateVectorSimulator"),
                                             ("braket_dm", "DensityMatrixSimulator")]
-        local_simulator_device = LocalSimulator(backend)
+        local_simulator_device = PyBraket.LocalSimulator(backend)
         @test name(local_simulator_device) == device_name
     end
     @testset for DEVICE in (PURE_DEVICE, NOISE_DEVICE)
@@ -357,4 +357,4 @@ end
             end
         end
     end
-end
\ No newline at end of file
+end
diff --git a/PyBraket/test/noise.jl b/PyBraket/test/noise.jl
index e4423483..6a9eb1d2 100644
--- a/PyBraket/test/noise.jl
+++ b/PyBraket/test/noise.jl
@@ -1,12 +1,12 @@
 using Test, PyBraket, Braket, Braket.IR, PythonCall
-using PythonCall: Py, pyconvert, pyisTrue
+using PythonCall: Py, pyconvert
 
 @testset "Noise" begin
     circ  = CNot(H(Circuit(), 0), 0, 1)
     noise = BitFlip(0.1)
     circ  = Braket.apply_gate_noise!(circ, noise)
 
-    device = LocalSimulator("braket_dm")
+    device = PyBraket.LocalSimulator("braket_dm")
     # run the circuit on the local simulator
     task = run(device, circ, shots=1000)
 
@@ -37,14 +37,14 @@ using PythonCall: Py, pyconvert, pyisTrue
             n = noise(0.1)
             ir_n = Braket.ir(n, 0)
             py_n = Py(n)
-            @test pyisTrue(py_n.to_ir([0]) == Py(ir_n))
+            @test pyconvert(Bool, py_n.to_ir([0]) == Py(ir_n))
             @test pyconvert(ir_noise, Py(ir_n)) == ir_n
         end
         @testset "(noise, ir_noise) = (PauliChannel, IR.PauliChannel)" begin
             n = PauliChannel(0.1, 0.2, 0.1)
             ir_n = Braket.ir(n, 0)
             py_n = Py(n)
-            @test pyisTrue(py_n.to_ir([0]) == Py(ir_n))
+            @test pyconvert(Bool, py_n.to_ir([0]) == Py(ir_n))
             @test pyconvert(IR.PauliChannel, Py(ir_n)) == ir_n
         end
         @testset "(noise, ir_noise) = (MultiQubitPauliChannel{1}, IR.MultiQubitPauliChannel)" begin
@@ -56,7 +56,7 @@ using PythonCall: Py, pyconvert, pyisTrue
             n = TwoQubitPauliChannel(Dict("XX"=>0.1, "YY"=>0.2))
             ir_n = Braket.ir(n, [0, 1])
             py_n = Py(n)
-            @test pyisTrue(py_n.to_ir([0, 1]) == Py(ir_n))
+            @test pyconvert(Bool, py_n.to_ir([0, 1]) == Py(ir_n))
             @test pyconvert(IR.MultiQubitPauliChannel, Py(ir_n)) == ir_n
         end
         @testset for (noise, ir_noise) in ((TwoQubitDephasing, IR.TwoQubitDephasing),
@@ -64,14 +64,14 @@ using PythonCall: Py, pyconvert, pyisTrue
             n = noise(0.4)
             ir_n = Braket.ir(n, [0, 1])
             py_n = Py(n)
-            @test pyisTrue(py_n.to_ir([0, 1]) == Py(ir_n))
+            @test pyconvert(Bool, py_n.to_ir([0, 1]) == Py(ir_n))
             @test pyconvert(ir_noise, Py(ir_n)) == ir_n
         end
         @testset "(noise, ir_noise) = (GeneralizedAmplitudeDamping, IR.GeneralizedAmplitudeDamping)" begin
             n = GeneralizedAmplitudeDamping(0.1, 0.2)
             ir_n = Braket.ir(n, 0)
             py_n = Py(n)
-            @test pyisTrue(py_n.to_ir([0]) == Py(ir_n))
+            @test pyconvert(Bool, py_n.to_ir([0]) == Py(ir_n))
             @test pyconvert(IR.GeneralizedAmplitudeDamping, Py(ir_n)) == ir_n
         end
         @testset "(noise, ir_noise) = (Kraus, IR.Kraus)" begin
@@ -79,7 +79,7 @@ using PythonCall: Py, pyconvert, pyisTrue
             n = Kraus([mat])
             ir_n = Braket.ir(n, 0)
             py_n = Py(n)
-            @test pyisTrue(py_n.to_ir([0]) == Py(ir_n))
+            @test pyconvert(Bool, py_n.to_ir([0]) == Py(ir_n))
             @test pyconvert(IR.Kraus, Py(ir_n)) == ir_n
         end 
         Braket.IRType[] = :OpenQASM
diff --git a/docs/make.jl b/docs/make.jl
index d046272b..e1ebc472 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -4,4 +4,4 @@ using Documenter, Braket
 
 makedocs(sitename="Braket.jl")
 
-deploydocs(repo="github.com/awslabs/braket.jl.git",)
+deploydocs(repo="github.com/amazon-braket/braket.jl.git",)
diff --git a/docs/src/local_simulator.md b/docs/src/local_simulator.md
new file mode 100644
index 00000000..701686d2
--- /dev/null
+++ b/docs/src/local_simulator.md
@@ -0,0 +1,13 @@
+# Local Simulators
+
+Local simulators allow you to *classically* simulate your [`Circuit`](@ref) or `OpenQasmProgram`
+on local hardware, rather than sending it to the cloud for on-demand execution. Local simulators
+can run task *batches* in parallel using multithreading.
+
+```@docs
+Braket._simulator_devices
+Braket.LocalQuantumTask
+Braket.LocalQuantumTaskBatch
+LocalSimulator
+simulate
+```
diff --git a/src/Braket.jl b/src/Braket.jl
index a56e810a..42bb0c7c 100644
--- a/src/Braket.jl
+++ b/src/Braket.jl
@@ -1,7 +1,7 @@
 module Braket
 
 export Circuit, QubitSet, Qubit, Device, AwsDevice, AwsQuantumTask, AwsQuantumTaskBatch
-export metadata, status, Observable, Result, FreeParameter, Job, AwsQuantumJob, LocalQuantumJob
+export metadata, status, Observable, Result, FreeParameter, Job, AwsQuantumJob, LocalQuantumJob, LocalSimulator
 export Tracker, simulator_tasks_cost, qpu_tasks_cost
 export arn, cancel, state, result, results, name, download_result, id, ir, isavailable, search_devices, get_devices
 export provider_name, properties, type
@@ -9,6 +9,7 @@ export apply_gate_noise!, apply
 export logs, log_metric, metrics, @hybrid_job
 export depth, qubit_count, qubits, ir, IRType, OpenQASMSerializationProperties
 export OpenQasmProgram
+export simulate
 export QueueDepthInfo, QueueType, Normal, Priority, queue_depth, queue_position
 
 export AdjointGradient, Expectation, Sample, Variance, Amplitude, Probability, StateVector, DensityMatrix, Result
@@ -30,10 +31,10 @@ using Downloads
 using DecFP
 using Graphs
 using HTTP
+using StaticArrays
 using JSON3, StructTypes
 using LinearAlgebra
 using DataStructures
-
 using NamedTupleTools
 using OrderedCollections
 using Tar
@@ -68,6 +69,18 @@ const IRType = Ref{Symbol}()
 include("tracker.jl")
 const Prices = Ref{Pricing}()
 const GlobalTrackerContext = Ref{TrackerContext}()
+abstract type AbstractBraketSimulator end
+
+"""
+    _simulator_devices
+
+A `Ref{Dict}` which records the `id`s of registered
+[`LocalSimulator`](@ref) backends so that they can be
+retrieved by this `id`. A new simulator backend
+should register itself in `_simulator_devices`
+in its module's `__init__` function.
+"""
+const _simulator_devices = Ref{Dict}()
 
 function __init__()
     downloader = Downloads.Downloader()
@@ -77,8 +90,11 @@ function __init__()
     IRType[] = :OpenQASM
     Prices[] = Pricing([])
     GlobalTrackerContext[] = TrackerContext()
+    _simulator_devices[] = Dict()
 end
 
+qubit_count(o) = throw(MethodError(qubit_count, (o,)))
+
 """
     Device
 
@@ -87,7 +103,7 @@ Abstract type representing a generic device which tasks (local or managed) may b
 abstract type Device end
 
 
-qubit_count(t) = throw(MethodError(qubit_count, t))
+#qubit_count(t) = throw(MethodError(qubit_count, t))
 chars(t) = (string(t),) 
 
 include("qubit_set.jl")
@@ -130,6 +146,7 @@ include("noise_model.jl")
 include("ahs.jl")
 include("queue_information.jl")
 include("device.jl")
+include("local_simulator.jl")
 include("gate_applicators.jl")
 include("noise_applicators.jl")
 include("jobs.jl")
@@ -138,5 +155,4 @@ include("aws_jobs.jl")
 include("local_jobs.jl")
 include("task.jl")
 include("task_batch.jl")
-
 end # module
diff --git a/src/circuit.jl b/src/circuit.jl
index 61042703..0860b68c 100644
--- a/src/circuit.jl
+++ b/src/circuit.jl
@@ -14,9 +14,9 @@ See:
 """
 mutable struct Circuit
     moments::Moments
-    instructions::Vector{Instruction}
+    instructions::Vector{Instruction{<:Operator}}
     result_types::Vector{Result}
-    basis_rotation_instructions::Vector{Instruction}
+    basis_rotation_instructions::Vector{Instruction{<:Operator}}
     qubit_observable_mapping::Dict{Int, Observables.Observable}
     qubit_observable_target_mapping::Dict{Int, Tuple}
     qubit_observable_set::Set{Int}
@@ -40,8 +40,8 @@ a vector of [`Result`](@ref)s, and a vector of basis rotation instructions.
 function Circuit(m::Moments, ixs::Vector, rts::Vector{Result}, bri::Vector)
     c = Circuit()
     c.moments = deepcopy(m)
-    c.instructions = deepcopy(ixs)
-    c.basis_rotation_instructions = deepcopy(bri)
+    c.instructions = Instruction.(deepcopy(ixs))
+    c.basis_rotation_instructions = Instruction.(deepcopy(bri))
     foreach(rt->add_result_type!(c, rt), rts)
     return c
 end
@@ -174,9 +174,9 @@ end
 
 (c::Circuit)(::Type{T}, args...) where {T<:Gate}  = apply_gate!(T, c, args...)
 (c::Circuit)(::Type{T}, args...) where {T<:Noise} = apply_noise!(T, c, args...)
-(c::Circuit)(::Type{T}) where {T<:CompilerDirective} = add_instruction!(c, Instruction(T()))
-(c::Circuit)(g::QuantumOperator, args...) = add_instruction!(c, Instruction(g, args...))
-(c::Circuit)(g::CompilerDirective) = add_instruction!(c, Instruction(g))
+(c::Circuit)(::Type{T}) where {T<:CompilerDirective} = add_instruction!(c, Instruction{T}(T()))
+(c::Circuit)(g::QO, args...) where {QO<:QuantumOperator} = add_instruction!(c, Instruction(g, args...))
+(c::Circuit)(g::CD) where {CD<:CompilerDirective} = add_instruction!(c, Instruction{CD}(g))
 (c::Circuit)(v::AbstractVector) = foreach(vi->c(vi...), v)
 (c::Circuit)(rt::Result, args...) = add_result_type!(c, rt, args...)
 (c::Circuit)(::Type{T}, args...) where {T<:Result} = T(c, args...)
@@ -264,7 +264,12 @@ julia> qubits(c)
 QubitSet(0, 1)
 ```
 """
-qubits(c::Circuit)      = (qs = union!(copy(c.moments._qubits), c.qubit_observable_set); QubitSet(qs))
+qubits(c::Circuit) = (qs = union!(copy(c.moments._qubits), c.qubit_observable_set); QubitSet(qs))
+function qubits(p::Program)
+    inst_qubits = mapreduce(ix->ix.target, union, p.instructions, init=Set{Int}())
+    res_qubits  = mapreduce(ix->(hasproperty(ix, :targets) && !isnothing(ix.targets)) ? reduce(vcat, ix.targets) : Set{Int}(), union, p.results, init=Set{Int}())
+    return union(inst_qubits, res_qubits)
+end
 """
     qubit_count(c::Circuit) -> Int
 
@@ -283,7 +288,7 @@ julia> qubit_count(c)
 ```
 """
 qubit_count(c::Circuit) = length(qubits(c))
-qubit_count(p::Program) = length(mapreduce(ix->ix.target, union, p.instructions))
+qubit_count(p::Program) = length(qubits(p))
 
 (rt::Result)(c::Circuit) = add_result_type!(c, rt)
 
@@ -388,7 +393,7 @@ basis_rotation_gates(o::Observables.TensorProduct) = tuple(reduce(vcat, basis_ro
 basis_rotation_gates(o::Observables.HermitianObservable) = (Unitary(Matrix(adjoint(eigvecs(o.matrix)))),)
 
 
-_observable_to_instruction(observable::Observables.Observable, target_list)::Vector{Instruction} = [Instruction(gate, target_list) for gate in basis_rotation_gates(observable)]
+_observable_to_instruction(observable::Observables.Observable, target_list)::Vector{Instruction{<:Operator}} = [Instruction(gate, target_list) for gate in basis_rotation_gates(observable)]
 
 """
     basis_rotation_instructions!(c::Circuit)
@@ -500,7 +505,7 @@ end
 add_result_type!(c::Circuit, rt::Result, target) = add_result_type!(c, remap(rt, target))
 add_result_type!(c::Circuit, rt::Result, target_mapping::Dict{<:Integer, <:Integer}) = add_result_type!(c, remap(rt, target_mapping))
 
-function add_instruction!(c::Circuit, ix::Instruction)
+function add_instruction!(c::Circuit, ix::Instruction{O}) where {O<:Operator}
     to_add = [ix]
     if ix.operator isa QuantumOperator && Parametrizable(ix.operator) == Parametrized()
         for param in parameters(ix.operator)
@@ -512,7 +517,7 @@ function add_instruction!(c::Circuit, ix::Instruction)
     return c
 end
 
-function add_instruction!(c::Circuit, ix::Instruction, target)
+function add_instruction!(c::Circuit, ix::Instruction{O}, target) where {O<:Operator}
     to_add = Instruction[]
     if qubit_count(ix.operator) == 1
         to_add = [remap(ix, q) for q in target]
@@ -523,7 +528,7 @@ function add_instruction!(c::Circuit, ix::Instruction, target)
     return c
 end
 
-function add_instruction!(c::Circuit, ix::Instruction, target_mapping::Dict{<:Integer, <:Integer})
+function add_instruction!(c::Circuit, ix::Instruction{O}, target_mapping::Dict{<:Integer, <:Integer}) where {O<:Operator}
     to_add = [remap(ix, target_mapping)]
     foreach(ix->add_instruction!(c, ix), to_add)
     return c
diff --git a/src/gate_applicators.jl b/src/gate_applicators.jl
index d14af113..0a3359d7 100644
--- a/src/gate_applicators.jl
+++ b/src/gate_applicators.jl
@@ -2,27 +2,27 @@
 (::Type{G})(c::Circuit, args...) where {G<:Gate} = apply_gate!(Val(n_angles(G)), IR.Control(ir_typ(G)), IR.Target(ir_typ(G)), G, c, args...)
 apply_gate!(::Type{G}, c::Circuit, args...) where {G<:Gate} = apply_gate!(Val(n_angles(G)), IR.Control(ir_typ(G)), IR.Target(ir_typ(G)), G, c, args...)
 
-apply_gate!(::Val{N}, ::IR.NoControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, arg::IntOrQubit, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction(G(angle), arg))
-apply_gate!(::Val{N}, ::IR.NoControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, v::VecOrQubitSet, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = (foreach(i->add_instruction!(c, Instruction(G(angle), i)), v); return c)
-apply_gate!(::Val{N}, ::IR.NoControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, v::NTuple{Nq, Ti}, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N, Nq, Ti} = (foreach(i->add_instruction!(c, Instruction(G(angle), i)), v); return c)
-apply_gate!(::Val{N}, ::IR.NoControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, args...) where {G<:Gate, N} = (foreach(i->add_instruction!(c, Instruction(G(args[end-(N-1):end]), i)), args[1:end-N]); return c)
-apply_gate!(::Val{0}, ::IR.NoControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, qs::Vararg{IntOrQubit}) where {G<:Gate} = (foreach(i->add_instruction!(c, Instruction(G(), i)), qs); return c)
+apply_gate!(::Val{N}, ::IR.NoControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, arg::IntOrQubit, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction{G}(G(angle), arg))
+apply_gate!(::Val{N}, ::IR.NoControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, v::VecOrQubitSet, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = (foreach(i->add_instruction!(c, Instruction{G}(G(angle), i)), v); return c)
+apply_gate!(::Val{N}, ::IR.NoControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, v::NTuple{Nq, Ti}, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N, Nq, Ti} = (foreach(i->add_instruction!(c, Instruction{G}(G(angle), i)), v); return c)
+apply_gate!(::Val{N}, ::IR.NoControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, args...) where {G<:Gate, N} = (foreach(i->add_instruction!(c, Instruction{G}(G(args[end-(N-1):end]), i)), args[1:end-N]); return c)
+apply_gate!(::Val{0}, ::IR.NoControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, qs::Vararg{IntOrQubit}) where {G<:Gate} = (foreach(i->add_instruction!(c, Instruction{G}(G(), i)), qs); return c)
 
-apply_gate!(::Val{N}, ::IR.NoControl, ::IR.DoubleTarget, ::Type{G}, c::Circuit, t1::IntOrQubit, t2::IntOrQubit, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction(G(angle), [t1, t2]))
-apply_gate!(::Val{N}, ::IR.NoControl, ::IR.DoubleTarget, ::Type{G}, c::Circuit, args::VecOrQubitSet, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction(G(angle), args[1:2]))
-apply_gate!(::Val{N}, ::IR.NoControl, ::IR.DoubleTarget, ::Type{G}, c::Circuit, args::NTuple{2, Ti}, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N, Ti} = add_instruction!(c, Instruction(G(angle), args[1:2]))
+apply_gate!(::Val{N}, ::IR.NoControl, ::IR.DoubleTarget, ::Type{G}, c::Circuit, t1::IntOrQubit, t2::IntOrQubit, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction{G}(G(angle), [t1, t2]))
+apply_gate!(::Val{N}, ::IR.NoControl, ::IR.DoubleTarget, ::Type{G}, c::Circuit, args::VecOrQubitSet, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction{G}(G(angle), args[1:2]))
+apply_gate!(::Val{N}, ::IR.NoControl, ::IR.DoubleTarget, ::Type{G}, c::Circuit, args::NTuple{2, Ti}, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N, Ti} = add_instruction!(c, Instruction{G}(G(angle), args[1:2]))
 
-apply_gate!(::Val{N}, ::IR.SingleControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, ci::IntOrQubit, ti::IntOrQubit, angles::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction(G(angles), [ci, ti]))
-apply_gate!(::Val{N}, ::IR.SingleControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, args::VecOrQubitSet, angles::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction(G(angles), args[1:2]))
-apply_gate!(::Val{N}, ::IR.SingleControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, args::NTuple{2, Ti}, angles::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N, Ti} = add_instruction!(c, Instruction(G(angles), [args...]))
+apply_gate!(::Val{N}, ::IR.SingleControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, ci::IntOrQubit, ti::IntOrQubit, angles::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction{G}(G(angles), [ci, ti]))
+apply_gate!(::Val{N}, ::IR.SingleControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, args::VecOrQubitSet, angles::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction{G}(G(angles), args[1:2]))
+apply_gate!(::Val{N}, ::IR.SingleControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, args::NTuple{2, Ti}, angles::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N, Ti} = add_instruction!(c, Instruction{G}(G(angles), [args...]))
 
-apply_gate!(::Val{N}, ::IR.SingleControl, ::IR.DoubleTarget, ::Type{G}, c::Circuit, ci::IntOrQubit, t1::IntOrQubit, t2::IntOrQubit, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction(G(angle), [ci, t1, t2]))
-apply_gate!(::Val{N}, ::IR.SingleControl, ::IR.DoubleTarget, ::Type{G}, c::Circuit, args::VecOrQubitSet, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction(G(angle), args[1:3]))
-apply_gate!(::Val{N}, ::IR.SingleControl, ::IR.DoubleTarget, ::Type{G}, c::Circuit, args::NTuple{3, Ti}, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N, Ti} = add_instruction!(c, Instruction(G(angle), [args...]))
+apply_gate!(::Val{N}, ::IR.SingleControl, ::IR.DoubleTarget, ::Type{G}, c::Circuit, ci::IntOrQubit, t1::IntOrQubit, t2::IntOrQubit, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction{G}(G(angle), [ci, t1, t2]))
+apply_gate!(::Val{N}, ::IR.SingleControl, ::IR.DoubleTarget, ::Type{G}, c::Circuit, args::VecOrQubitSet, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction{G}(G(angle), args[1:3]))
+apply_gate!(::Val{N}, ::IR.SingleControl, ::IR.DoubleTarget, ::Type{G}, c::Circuit, args::NTuple{3, Ti}, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N, Ti} = add_instruction!(c, Instruction{G}(G(angle), [args...]))
 
-apply_gate!(::Val{N}, ::IR.DoubleControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, c1::IntOrQubit, c2::IntOrQubit, ti::IntOrQubit, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction(G(angle), [c1, c2, ti]))
-apply_gate!(::Val{N}, ::IR.DoubleControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, args::VecOrQubitSet, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction(G(angle), args[1:3]))
-apply_gate!(::Val{N}, ::IR.DoubleControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, args::NTuple{3, Ti}, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N, Ti} = add_instruction!(c, Instruction(G(angle), [args...]))
+apply_gate!(::Val{N}, ::IR.DoubleControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, c1::IntOrQubit, c2::IntOrQubit, ti::IntOrQubit, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction{G}(G(angle), [c1, c2, ti]))
+apply_gate!(::Val{N}, ::IR.DoubleControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, args::VecOrQubitSet, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N} = add_instruction!(c, Instruction{G}(G(angle), args[1:3]))
+apply_gate!(::Val{N}, ::IR.DoubleControl, ::IR.SingleTarget, ::Type{G}, c::Circuit, args::NTuple{3, Ti}, angle::Vararg{<:Union{Float64, FreeParameter}}) where {G<:Gate, N, Ti} = add_instruction!(c, Instruction{G}(G(angle), [args...]))
 
-apply_gate!(::Val{0}, ::IR.NoControl, ::IR.MultiTarget, ::Type{Unitary}, c::Circuit, v::VecOrQubitSet, m::Matrix{ComplexF64}) = add_instruction!(c, Instruction(Unitary(m), v))
-apply_gate!(::Val{0}, ::IR.NoControl, ::IR.MultiTarget, ::Type{Unitary}, c::Circuit, args...) = add_instruction!(c, Instruction(Unitary(args[end]), [args[1:end-1]...]))
+apply_gate!(::Val{0}, ::IR.NoControl, ::IR.MultiTarget, ::Type{Unitary}, c::Circuit, v::VecOrQubitSet, m::Matrix{ComplexF64}) = add_instruction!(c, Instruction{Unitary}(Unitary(m), v))
+apply_gate!(::Val{0}, ::IR.NoControl, ::IR.MultiTarget, ::Type{Unitary}, c::Circuit, args...) = add_instruction!(c, Instruction{Unitary}(Unitary(args[end]), [args[1:end-1]...]))
diff --git a/src/gates.jl b/src/gates.jl
index 8ea1b7fe..cefd8d6f 100644
--- a/src/gates.jl
+++ b/src/gates.jl
@@ -45,14 +45,14 @@ for gate_def in (
         """
         struct $G <: AngledGate{$n_angle}
             angle::NTuple{$n_angle, Union{Float64, FreeParameter}}
-            $G(angle::NTuple{$n_angle, Union{Float64, FreeParameter}}) = new(angle)
+            $G(angle::T) where {T<:NTuple{$n_angle, Union{Float64, FreeParameter}}} = new(angle)
         end
         $G(angles::Vararg{Union{Float64, FreeParameter}}) = $G(tuple(angles...))
+        $G(angles::Vararg{Number}) = $G((Float64(a) for a in angles)...)
         chars(::Type{$G}) = $c
         qubit_count(::Type{$G}) = $qc
     end
 end
-(::Type{G})(angle::Union{Float64, FreeParameter}) where {G<:AngledGate} = G((angles,))
 
 for gate_def in (
 	(:H, :1, ("H",)),
@@ -91,16 +91,17 @@ end
 label(::Type{G}) where {G<:Gate} = lowercase(string(G))
 (::Type{G})(x::Tuple{}) where {G<:Gate} = G()
 (::Type{G})(x::Tuple{}) where {G<:AngledGate} = throw(ArgumentError("angled gate must be constructed with at least one angle."))
-(::Type{G})(x::AbstractVector) where {G<:AngledGate} = G(x...) 
-qubit_count(g::G) where {G<:Gate} = qubit_count(G)
-angles(g::G) where {G<:Gate} = ()
+(::Type{G})(x::AbstractVector) where {G<:AngledGate} = G(x...)
+(::Type{G})(x::T) where {G<:AngledGate{1}, T<:Union{Float64, FreeParameter}} = G((x,))
+qubit_count(g::G) where {G<:Gate}  = qubit_count(G)
+angles(g::G) where {G<:Gate}       = ()
 angles(g::AngledGate{N}) where {N} = g.angle
-chars(g::G) where {G<:Gate} = map(char->replace(string(char), "ang"=>join(string.(angles(g)), ", ")), chars(G))
-ir_typ(::Type{G}) where {G<:Gate} = getproperty(IR, Symbol(G))
-ir_typ(g::G) where {G<:Gate} = ir_typ(G)
-label(g::G) where {G<:Gate} = label(G)
+chars(g::G) where {G<:Gate}        = map(char->replace(string(char), "ang"=>join(string.(angles(g)), ", ")), chars(G))
+ir_typ(::Type{G}) where {G<:Gate}  = getproperty(IR, Symbol(G))
+ir_typ(g::G) where {G<:Gate}       = ir_typ(G)
+label(g::G) where {G<:Gate}        = label(G)
 ir_str(g::G) where {G<:AngledGate} = label(g) * "(" * join(string.(angles(g)), ", ") * ")"
-ir_str(g::G) where {G<:Gate} = label(g)
+ir_str(g::G) where {G<:Gate}       = label(g)
 targets_and_controls(g::G, target::QubitSet) where {G<:Gate} = IR._generate_control_and_target(IR.ControlAndTarget(ir_typ(g))..., target)
 function ir(g::G, target::QubitSet, ::Val{:JAQCD}; kwargs...) where {G<:Gate}
     t_c = targets_and_controls(g, target)
@@ -125,7 +126,7 @@ struct Unitary <: Gate
     matrix::Matrix{ComplexF64}
     Unitary(matrix::Matrix{<:Number}) = new(ComplexF64.(matrix))
 end
-Unitary(mat::Vector{Vector{Vector{Float64}}}) = Unitary(complex_matrix_from_ir(mat))
+Unitary(mat::Vector{Vector{Vector{T}}}) where {T} = Unitary(complex_matrix_from_ir(mat))
 Base.:(==)(u1::Unitary, u2::Unitary) = u1.matrix == u2.matrix
 qubit_count(g::Unitary) = convert(Int, log2(size(g.matrix, 1)))
 chars(g::Unitary)       = ntuple(i->"U", qubit_count(g))
diff --git a/src/local_simulator.jl b/src/local_simulator.jl
new file mode 100644
index 00000000..6982f484
--- /dev/null
+++ b/src/local_simulator.jl
@@ -0,0 +1,197 @@
+"""
+    LocalQuantumTask(id::String, result::GateModelQuantumTaskResult)
+
+A quantum task which has been run *locally* using a [`LocalSimulator`](@ref).
+The `state` of a `LocalQuantumTask` is always `"COMPLETED"` as the task object
+is only created once the loca simulation has finished.
+"""
+struct LocalQuantumTask
+    id::String
+    result::GateModelQuantumTaskResult
+end
+state(b::LocalQuantumTask)  = "COMPLETED"
+id(b::LocalQuantumTask)     = b.id
+result(b::LocalQuantumTask) = b.result
+
+"""
+    LocalQuantumTaskBatch(ids::Vector{String}, results::Vector{GateModelQuantumTaskResult})
+
+A *batch* of quantum tasks which has been run *locally* using a
+[`LocalSimulator`](@ref). 
+"""
+struct LocalQuantumTaskBatch
+    ids::Vector{String}
+    results::Vector{GateModelQuantumTaskResult}
+end
+results(lqtb::LocalQuantumTaskBatch) = lqtb.results
+
+"""
+    LocalSimulator(backend::Union{String, AbstractBraketSimulator})
+
+A quantum simulator which runs *locally* rather than sending the circuit(s)
+to the cloud to be executed on-demand. A `LocalSimulator` must be created with
+a backend -- either a handle, in the form of a `String`, which uniquely identifies
+a simulator backend registered in [`_simulator_devices`](@ref Braket._simulator_devices), or an already instantiated
+simulator object.
+
+`LocalSimulator`s should implement their own method for [`simulate`](@ref) if needed. They
+can process single tasks or task batches.
+"""
+struct LocalSimulator <: Device
+    backend::String
+    _delegate::AbstractBraketSimulator
+    function LocalSimulator(backend::Union{String, AbstractBraketSimulator})
+        backend_name = device_id(backend)
+        haskey(_simulator_devices[], backend_name) && return new(backend_name, _simulator_devices[][backend_name](0,0))
+        !isdefined(Main, Symbol(backend_name)) && throw(ArgumentError("no local simulator with name $backend_name is loaded!"))
+        return new(backend_name, Symbol(backend_name)(0, 0))
+    end
+end
+name(s::LocalSimulator) = name(s._delegate)
+device_id(s::String) = s
+
+"""
+    simulate(d::LocalSimulator, task_spec::Union{Circuit, AbstractProgram}, args...; shots::Int=0, inputs::Dict{String, Float64} = Dict{String, Float64}(), kwargs...)
+
+Simulate the execution of `task_spec` using the backend of [`LocalSimulator](@ref) `d`.
+`args` are additional arguments to be provided to the backend. `inputs` is used to set
+the value of any [`FreeParameter`](@ref) in `task_spec` and will override the existing
+`inputs` field of an `OpenQasmProgram`. Other `kwargs` will be passed to the backend
+simulator. Returns a [`LocalQuantumTask`](@ref Braket.LocalQuantumTask).
+"""
+function simulate(d::LocalSimulator, task_spec::Union{Circuit, AbstractProgram}, args...; shots::Int=0, inputs::Dict{String, Float64} = Dict{String, Float64}(), kwargs...)
+    sim = d._delegate
+    @debug "Single task. Starting run..."
+    stats = @timed _run_internal(sim, task_spec, args...; inputs=inputs, shots=shots, kwargs...)
+    @debug "Single task. Time to run internally: $(stats.time). GC time: $(stats.gctime)."
+    local_result = stats.value
+    return LocalQuantumTask(local_result.task_metadata.id, local_result)
+end
+
+"""
+    simulate(d::LocalSimulator, task_specs::Vector{T}, args...; kwargs...) where {T}
+
+Simulate the execution of a *batch* of tasks specified by `task_specs`
+using the backend of [`LocalSimulator](@ref) `d`.
+`args` are additional arguments to be provided to the backend.
+
+`kwargs` used by the `LocalSimulator` are:
+  - `shots::Int` - the number of shots to run for all tasks in `task_specs`. Default is `0`.
+  - `max_parallel::Int` - the maximum number of simulations to execute simultaneously. Default is `32`.
+  - `inputs::Union{Vector{Dict{String, Float64}}, Dict{String, Float64}}` - used to set
+    the value of any [`FreeParameter`](@ref) in each task specification. It must either be a
+    `Dict` or a single-element `Vector` (in which case the same parameter values are used for
+    all elements of `task_specs` *or* of the same length as `task_specs` (in which case the `i`-th 
+    specification is paired with the `i`-th input dictionary). Default is an empty dictionary. 
+
+Other `kwargs` are passed through to the backend simulator. Returns a [`LocalQuantumTaskBatch`](@ref Braket.LocalQuantumTaskBatch).
+
+!!! note
+    Because Julia uses dynamic threading and `Task`s can migrate between threads, each simulation is a `Task`
+    which itself can spawn many more `Task`s, as the internal implementation of [`LocalSimulator`](@ref)'s backend
+    may use threading where appropriate. On systems with many CPU cores, spawning too many `Task`s may overwhelm
+    the Julia scheduler and degrade performance. "Too many" depends on the particulars of your hardware, so on
+    many-core systems you may need to tune this value for best performance.
+"""
+function simulate(d::LocalSimulator, task_specs::Vector{T}, args...; shots::Int=0, max_parallel::Int=-1, inputs::Union{Vector{Dict{String, Float64}}, Dict{String, Float64}} = Dict{String, Float64}(), kwargs...) where {T}
+    is_single_task  = length(task_specs) == 1
+    is_single_input = inputs isa Dict{String, Float64} || length(inputs) == 1
+    if is_single_input
+        if is_single_task 
+            return d(task_specs[1], args...; shots=shots, inputs=inputs, kwargs...)
+        elseif inputs isa Dict{String, Float64}
+            inputs = [deepcopy(inputs) for ix in 1:length(task_specs)]
+        else
+            inputs = [deepcopy(inputs[1]) for ix in 1:length(task_specs)]
+        end
+    end
+    !is_single_task && !is_single_input && length(task_specs) != length(inputs) && throw(ArgumentError("number of inputs ($(length(inputs))) and task specifications ($(length(task_specs))) must be equal."))
+    # let each thread pick up an idle simulator
+    #sims     = Channel(Inf)
+    #foreach(i -> put!(sims, copy(d._delegate)), 1:Threads.nthreads())
+
+    tasks_and_inputs = zip(1:length(task_specs), task_specs, inputs)
+    # is this actually faster?
+    todo_tasks_ch = Channel(length(tasks_and_inputs))
+    for (ix, spec, input) in tasks_and_inputs
+        if spec isa Circuit
+            param_names = Set(string(param.name) for param in spec.parameters)
+            unbounded_params = setdiff(param_names, collect(keys(input)))
+            !isempty(unbounded_params) && throw(ErrorException("cannot execute circuit with unbound parameters $unbounded_params."))
+        end
+        put!(todo_tasks_ch, (ix, spec, input))
+    end
+    @debug "batch size is $(length(task_specs)). Starting run..."
+    n_task_threads = 32
+    stats = @timed begin
+	    done_tasks_ch = Channel(length(tasks_and_inputs)) do ch
+		function task_processor(input_tup)
+		    ix, spec, input = input_tup
+		    sim = copy(d._delegate)
+		    res = _run_internal(sim, spec, args...; inputs=input, shots=shots, kwargs...)
+		    put!(ch, (ix, res))
+		end
+		Threads.foreach(task_processor, todo_tasks_ch, ntasks=n_task_threads)
+	    end
+	    r_ix = 1
+	    results = Vector{GateModelQuantumTaskResult}(undef, length(task_specs))
+	    while r_ix <= length(task_specs)
+            ix, res = take!(done_tasks_ch)
+            results[ix] = res
+            r_ix += 1
+	    end
+    end
+    @debug "batch size is $(length(task_specs)). Time to run internally: $(stats.time). GC time: $(stats.gctime)."
+    return LocalQuantumTaskBatch([local_result.task_metadata.id for local_result in results], results)
+end
+(d::LocalSimulator)(args...; kwargs...) = simulate(d, args...; kwargs...)
+
+function _run_internal(simulator, circuit::Circuit, args...; shots::Int=0, inputs::Dict{String, Float64}=Dict{String, Float64}(), kwargs...)
+    if haskey(properties(simulator).action, "braket.ir.openqasm.program")
+        validate_circuit_and_shots(circuit, shots)
+        program      = ir(circuit, Val(:OpenQASM))
+        full_inputs  = isnothing(program.inputs) ? inputs : merge(program.inputs, inputs)
+        full_program = OpenQasmProgram(program.braketSchemaHeader, program.source, full_inputs) 
+        r            = simulate(simulator, full_program, args...; shots=shots, kwargs...)
+        return format_result(r) 
+    elseif haskey(properties(simulator).action, "braket.ir.jaqcd.program")
+        validate_circuit_and_shots(circuit, shots) 
+        program = ir(circuit, Val(:JAQCD))
+        qubits  = qubit_count(circuit)
+        r       = simulate(simulator, program, qubits, args...; shots=shots, inputs=inputs, kwargs...)
+        return format_result(r)
+    else
+        throw(ErrorException("$(typeof(simulator)) does not support qubit gate-based programs."))
+    end
+end
+function _run_internal(simulator, program::OpenQasmProgram, args...; shots::Int=0, inputs::Dict{String, Float64}=Dict{String, Float64}(), kwargs...)
+    if haskey(properties(simulator).action, "braket.ir.openqasm.program")
+        stats = @timed begin
+            simulate(simulator, program; shots=shots, inputs=inputs, kwargs...)
+        end
+        @debug "Time to invoke simulator: $(stats.time)"
+        r     = stats.value
+        stats = @timed format_result(r)
+        @debug "Time to format results: $(stats.time)"
+        return stats.value
+    else
+        throw(ErrorException("$(typeof(simulator)) does not support qubit gate-based programs."))
+    end
+end
+function _run_internal(simulator, program::Program, args...; shots::Int=0, inputs::Dict{String, Float64}=Dict{String, Float64}(), kwargs...)
+    if haskey(properties(simulator).action, "braket.ir.jaqcd.program")
+        stats   = @timed qubit_count(program)
+        @debug "Time to get qubit count: $(stats.time)"
+        qubits  = stats.value
+        stats = @timed begin
+            simulate(simulator, program, qubits, args...; shots=shots, inputs=inputs, kwargs...)
+        end
+        @debug "Time to invoke simulator: $(stats.time)"
+        r     = stats.value
+        stats = @timed format_result(r)
+        @debug "Time to format results: $(stats.time)"
+        return stats.value
+    else
+        throw(ErrorException("$(typeof(simulator)) does not support qubit gate-based programs."))
+    end
+end
diff --git a/src/moments.jl b/src/moments.jl
index 9ed812d2..20306a21 100644
--- a/src/moments.jl
+++ b/src/moments.jl
@@ -21,7 +21,7 @@ mutable struct Moments <: AbstractDict{MomentKey,Instruction}
 end
 Base.get(m::Moments, mk::MomentKey, default) = get(m._moments, mk, default)
 
-function Moments(ixs::Vector)
+function Moments(ixs::Vector{Instruction})
     m = Moments()
     foreach(ix->push!(m, ix), ixs)
     return m
@@ -36,32 +36,34 @@ function Base.show(io::IO, m::Moments)
     println(io, "Max times: ", m._max_times)
 end
 
-function Base.push!(m::Moments, ix::Instruction, mt::MomentType_enum=mtGateNoise, noise_index::Int=0)
-    op = ix.operator
-    if op isa CompilerDirective
-        qubits = QubitSet(keys(m._max_times))
-        time   = _update_qubit_time!(m, qubits)
-        key    = MomentKey(time, QubitSet(), mtCompilerDirective, 0)
-        m._moments[key] = ix
-        m._depth = max(m._depth, time+1)
-        m._time_all_qubits = time
-    elseif op isa Noise
-        add_noise!(m, ix, "noise", noise_index)
-    else # it's a Gate
-        qubits = QubitSet(ix.target)
-        time   = _update_qubit_time!(m, qubits)
-        key    = MomentKey(time, qubits, mtGate, noise_index)
-        m._moments[key] = ix
-        union!(m._qubits, qubits)
-        m._depth = max(m._depth, time+1)
-    end
+function Base.push!(m::Moments, ix::Instruction{O}, mt::MomentType_enum=mtGateNoise, noise_index::Int=0) where {O<:CompilerDirective}
+    qubits = QubitSet(keys(m._max_times))
+    time   = _update_qubit_time!(m, qubits)
+    key    = MomentKey(time, QubitSet(), mtCompilerDirective, 0)
+    m._moments[key] = ix
+    m._depth = max(m._depth, time+1)
+    m._time_all_qubits = time
+end
+
+function Base.push!(m::Moments, ix::Instruction{O}, mt::MomentType_enum=mtGateNoise, noise_index::Int=0) where {O<:Noise}
+    add_noise!(m, ix, "noise", noise_index)
+    return m
+end
+
+function Base.push!(m::Moments, ix::Instruction{O}, mt::MomentType_enum=mtGateNoise, noise_index::Int=0) where {O}
+    qubits = QubitSet(ix.target)
+    time   = _update_qubit_time!(m, qubits)
+    key    = MomentKey(time, qubits, mtGate, noise_index)
+    m._moments[key] = ix
+    union!(m._qubits, qubits)
+    m._depth = max(m._depth, time+1)
     return m
 end
-Base.push!(m::Moments, ixs::Vector{Instruction}, mt::MomentType_enum=mtGateNoise, noise_index::Int=0) = (foreach(ix->push!(m, ix, mt, noise_index), ixs); return m)
+Base.push!(m::Moments, ixs::Vector{Instruction{O}}, mt::MomentType_enum=mtGateNoise, noise_index::Int=0) where {O} = (foreach(ix->push!(m, ix, mt, noise_index), ixs); return m)
 
-function add_noise!(m::Moments, ix::Instruction, input_type::String="noise", noise_index::Int=0)
+function add_noise!(m::Moments, ix::Instruction{O}, input_type::String="noise", noise_index::Int=0) where {O}
     qubits = QubitSet(ix.target)
-    time = max(0, prod(get(m._max_times, q, -1) for q in qubits))
+    time   = max(0, prod(get(m._max_times, q, -1) for q in qubits, init=1))
     if MomentType_toenum[input_type] == mtInitializationNoise
         time = 0
     end
diff --git a/src/observables.jl b/src/observables.jl
index b07d1e28..8703be15 100644
--- a/src/observables.jl
+++ b/src/observables.jl
@@ -2,7 +2,7 @@ module Observables
 
 using JSON3, StructTypes, LinearAlgebra
 
-import ..Braket: ir, qubit_count, pauli_eigenvalues, complex_matrix_from_ir, complex_matrix_to_ir, Operator, SerializationProperties, format_qubits, format, format_matrix, OpenQASMSerializationProperties, IRType, QubitSet, Qubit, IntOrQubit, IRObservable, chars
+import ..Braket: ir, qubit_count, complex_matrix_from_ir, complex_matrix_to_ir, Operator, SerializationProperties, format_qubits, format, format_matrix, OpenQASMSerializationProperties, IRType, QubitSet, Qubit, IntOrQubit, IRObservable, chars, PauliEigenvalues
 export Observable, TensorProduct, HermitianObservable, Sum
 
 """
@@ -14,11 +14,14 @@ have `eigvals` defined.
 See also: [`H`](@ref), [`I`](@ref), [`X`](@ref), [`Y`](@ref), [`Z`](@ref), [`TensorProduct`](@ref), [`HermitianObservable`](@ref).
 """
 abstract type Observable <: Operator end
+abstract type NonCompositeObservable <: Observable end
+abstract type StandardObservable <: NonCompositeObservable end
+
 LinearAlgebra.ishermitian(o::Observable) = true
 
 coef(o::O) where {O<:Observable} = o.coefficient
 
-for (typ, label) in ((:H, "h"), (:X, "x"), (:Y, "y"), (:Z, "z"), (:I, "i"))
+for (typ, label, super) in ((:H, "h", :StandardObservable), (:X, "x", :StandardObservable), (:Y, "y", :StandardObservable), (:Z, "z", :StandardObservable), (:I, "i", :NonCompositeObservable))
     @eval begin
         @doc """
             $($typ) <: Observable
@@ -27,7 +30,7 @@ for (typ, label) in ((:H, "h"), (:X, "x"), (:Y, "y"), (:Z, "z"), (:I, "i"))
         Struct representing a `$($typ)` observable in a measurement. The observable
         may be scaled by `coeff`.
         """
-        struct $typ <: Observable
+        struct $typ <: $super
             coefficient::Float64
             $typ(coef::Float64=1.0) = new(coef)
         end
@@ -49,9 +52,53 @@ for (typ, label) in ((:H, "h"), (:X, "x"), (:Y, "y"), (:Z, "z"), (:I, "i"))
         chars(o::$typ) = (uppercase($label),)
     end
 end
-const StandardObservable = Union{H, X, Y, Z}
-LinearAlgebra.eigvals(o::I) = o.coefficient .* [1.0, 1.0]
-LinearAlgebra.eigvals(o::StandardObservable) = o.coefficient .* [1.0, -1.0]
+LinearAlgebra.eigvals(o::I) = [o.coefficient, o.coefficient]
+LinearAlgebra.eigvals(o::StandardObservable) = [o.coefficient, -o.coefficient]
+
+"""
+    HermitianObservable <: Observable
+    HermitianObservable(matrix::Matrix) -> HermitianObservable
+
+Struct representing an observable of an arbitrary complex Hermitian matrix.
+
+# Examples
+```jldoctest
+julia> ho = Braket.Observables.HermitianObservable([0 1; 1 0])
+Braket.Observables.HermitianObservable(Complex{Int64}[0 + 0im 1 + 0im; 1 + 0im 0 + 0im])
+```
+"""
+struct HermitianObservable <: NonCompositeObservable
+    matrix::Matrix{<:Complex}
+    coefficient::Float64
+    function HermitianObservable(mat::Matrix{<:Number})
+        ishermitian(mat) || throw(ArgumentError("input matrix to HermitianObservable must be Hermitian."))
+        new(complex(mat), 1.0)
+    end
+end
+HermitianObservable(v::Vector{Vector{Vector{T}}}) where {T<:Number} = HermitianObservable(complex_matrix_from_ir(v))
+Base.copy(o::HermitianObservable) = HermitianObservable(copy(o.matrix))
+StructTypes.StructType(::Type{HermitianObservable}) = StructTypes.CustomStruct()
+StructTypes.lower(x::HermitianObservable) = Union{String, Vector{Vector{Vector{Float64}}}}[complex_matrix_to_ir(ComplexF64.(x.matrix))]
+Base.:(==)(h1::HermitianObservable, h2::HermitianObservable) = (size(h1.matrix) == size(h2.matrix) && h1.matrix ≈ h2.matrix)
+qubit_count(o::HermitianObservable) = convert(Int, log2(size(o.matrix, 1)))
+LinearAlgebra.eigvals(o::HermitianObservable) = eigvals(Hermitian(o.matrix))
+unscaled(o::HermitianObservable) = o
+function ir(ho::HermitianObservable, target::QubitSet, ::Val{:OpenQASM}; serialization_properties=OpenQASMSerializationProperties())
+    function c_str(x::Number)
+        iszero(real(x)) && iszero(imag(x)) && return "0im"
+        iszero(real(x)) && return string(imag(x))*"im"
+        iszero(imag(x)) && return string(real(x))*"+0im"
+        imag_sign = imag(x) < 0 ? "" : "+"
+        return string(real(x)) * imag_sign * string(imag(x)) * "im"
+    end
+    m = "[" * join(["["*join(c_str.(ho.matrix[i, :]), ", ")*"]" for i in 1:size(ho.matrix, 1)], ", ") * "]"
+    t = isempty(target) ? "all" : format_qubits(target, serialization_properties)
+    return "hermitian($m) $t"
+end
+Base.:(*)(o::HermitianObservable, n::Real) = HermitianObservable(Float64(n) .* o.matrix)
+Base.show(io::IO, ho::HermitianObservable) = print(io, "HermitianObservable($(size(ho.matrix)))")
+ir(ho::HermitianObservable, target::Nothing, ::Val{:OpenQASM}; kwargs...) = ir(ho, QubitSet(), Val(:OpenQASM); kwargs...)
+chars(o::HermitianObservable) = ("Hermitian",)
 
 """
     TensorProduct <: Observable
@@ -71,31 +118,49 @@ julia> Braket.Observables.TensorProduct([ho, Braket.Observables.Z()])
 Braket.Observables.TensorProduct(Braket.Observables.Observable[Braket.Observables.HermitianObservable(Complex{Int64}[0 + 0im 1 + 0im; 1 + 0im 0 + 0im]), Braket.Observables.Z()])
 ```
 """
-struct TensorProduct <: Observable
-    factors::Vector{<:Observable}
+struct TensorProduct{O} <: Observable where {O<:Observable}
+    factors::Vector{O}
     coefficient::Float64
-    function TensorProduct(v::Vector{<:Observable}, coefficient::Float64=1.0)
-        flattened_v = Observable[]
-        for o in v
-            if o isa TensorProduct
-                foreach(o_->push!(flattened_v, o_), o.factors)
-            elseif o isa Sum
-                throw(ArgumentError("Sum observable not allowed in TensorProduct."))
-            else
-                push!(flattened_v, o)
+    function TensorProduct{O}(v::Vector{O}, coefficient::Float64=1.0) where {O<:NonCompositeObservable}
+        coeff = coefficient
+        flattened_v = Vector{O}(undef, length(v))
+        for (oi, o) in enumerate(v)
+            flattened_v[oi] = unscaled(o)
+            coeff *= coef(o)
+        end
+        return new(flattened_v, coeff)
+    end
+    function TensorProduct{O}(v::Vector{O}, coefficient::Float64=1.0) where {O<:Observable}
+        any(v_->v_ isa Sum, v) && throw(ArgumentError("Sum observable not allowed in TensorProduct."))
+        coeff = prod(coef, v, init=coefficient)
+        if VERSION >= v"1.9.0"
+            flattened_v = Iterators.flatmap(v) do o 
+                if o isa TensorProduct
+                    return (unscaled(o) for o in o.factors)
+                else
+                    return (unscaled(o),)
+                end
+            end
+        else
+            mapped_v = map(v) do o 
+                if o isa TensorProduct
+                    return (unscaled(o) for o in o.factors)
+                else
+                    return (unscaled(o),)
+                end
             end
+            flattened_v = Iterators.flatten(mapped_v)
         end
-        coeff = mapreduce(coef, *, flattened_v, init=1.0) * coefficient
-        unscaled_factors = Observable[unscaled(f) for f in flattened_v]
-        return new(unscaled_factors, coeff)
+        flat_v = collect(flattened_v)
+        return new(flat_v, coeff)
     end
 end
-Base.:(*)(o::TensorProduct, n::Real) = TensorProduct(deepcopy(o.factors), Float64(n*o.coefficient))
+TensorProduct(o::Vector{T}, coefficient::Float64=1.0) where {T<:Union{String, Observable}} = TensorProduct{T}(o, coefficient)
 function TensorProduct(o::Vector{String}, coefficient::Float64=1.0)
-    sts = StructTypes.subtypes(Observable)
-    ops = Observable[sts[Symbol(lowercase(oi))]() for oi in o]
-    return TensorProduct(ops, coefficient)
+    return TensorProduct([StructTypes.constructfrom(Observable, s) for s in o], coefficient)
 end
+
+Base.:(*)(o::TensorProduct, n::Real) = TensorProduct(deepcopy(o.factors), Float64(n*o.coefficient))
 function ir(tp::TensorProduct, target::QubitSet, ::Val{:OpenQASM}; serialization_properties::SerializationProperties=OpenQASMSerializationProperties())
     factors = []
     use_qubits = collect(target)
@@ -114,15 +179,13 @@ end
 ir(tp::TensorProduct; kwargs...) = ir(tp, Val(IRType[]); kwargs...)
 unscaled(o::TensorProduct) = TensorProduct(o.factors, 1.0)
 qubit_count(o::TensorProduct) = sum(qubit_count.(o.factors))
-StructTypes.StructType(::Type{TensorProduct}) = StructTypes.CustomStruct()
-StructTypes.lower(x::TensorProduct) = Union{String, Vector{Vector{Vector{Float64}}}}[convert(Union{String, Vector{Vector{Vector{Float64}}}}, o) for o in mapreduce(StructTypes.lower, vcat, x.factors)]
+StructTypes.StructType(::Type{TensorProduct{O}}) where {O<:Observable} = StructTypes.CustomStruct()
+StructTypes.lower(x::TensorProduct{O}) where {O<:Observable} = Union{String, Vector{Vector{Vector{Float64}}}}[convert(Union{String, Vector{Vector{Vector{Float64}}}}, o) for o in mapreduce(StructTypes.lower, vcat, x.factors)]
 Base.:(==)(t1::TensorProduct, t2::TensorProduct) = t1.factors == t2.factors && t1.coefficient ≈ t2.coefficient
 Base.copy(t::TensorProduct) = TensorProduct(deepcopy(t.factors), t.coefficient)
-function LinearAlgebra.eigvals(o::TensorProduct)
-    all(o_ isa StandardObservable for o_ in o.factors) && return o.coefficient .* pauli_eigenvalues(length(o.factors))
-    evs = mapfoldl(eigvals, kron, o.factors, init=[1.0])
-    return o.coefficient .* evs
-end
+_evs(os::Vector{O}, coeff::Float64=1.0) where {O<:StandardObservable} = PauliEigenvalues(Val(length(os)), coeff)
+_evs(os::Vector{O}, coeff::Float64=1.0) where {O<:Observable} = mapfoldl(eigvals, kron, os, init=[coeff])::Vector{Float64}
+LinearAlgebra.eigvals(o::TensorProduct{O}) where {O} = return _evs(o.factors, o.coefficient)
 function Base.show(io::IO, o::TensorProduct)
     coef_str = isone(o.coefficient) ? "" : string(o.coefficient) * " * "
     print(io, coef_str)
@@ -135,49 +198,6 @@ function Base.show(io::IO, o::TensorProduct)
 end
 chars(o::TensorProduct) = (sprint(show, o),)
 
-"""
-    HermitianObservable <: Observable
-    HermitianObservable(matrix::Matrix) -> HermitianObservable
-
-Struct representing an observable of an arbitrary complex Hermitian matrix.
-
-# Examples
-```jldoctest
-julia> ho = Braket.Observables.HermitianObservable([0 1; 1 0])
-Braket.Observables.HermitianObservable(Complex{Int64}[0 + 0im 1 + 0im; 1 + 0im 0 + 0im])
-```
-"""
-struct HermitianObservable <: Observable
-    matrix::Matrix{<:Complex}
-    coefficient::Float64
-    function HermitianObservable(mat::Matrix{<:Number})
-        ishermitian(mat) || throw(ArgumentError("input matrix to HermitianObservable must be Hermitian."))
-        new(complex(mat), 1.0)
-    end
-end
-HermitianObservable(v::Vector{Vector{Vector{T}}}) where {T<:Number} = HermitianObservable(complex_matrix_from_ir(v))
-Base.copy(o::HermitianObservable) = HermitianObservable(copy(o.matrix))
-StructTypes.StructType(::Type{HermitianObservable}) = StructTypes.CustomStruct()
-StructTypes.lower(x::HermitianObservable) = Union{String, Vector{Vector{Vector{Float64}}}}[complex_matrix_to_ir(ComplexF64.(x.matrix))]
-Base.:(==)(h1::HermitianObservable, h2::HermitianObservable) = (size(h1.matrix) == size(h2.matrix) && h1.matrix ≈ h2.matrix)
-qubit_count(o::HermitianObservable) = convert(Int, log2(size(o.matrix, 1)))
-LinearAlgebra.eigvals(o::HermitianObservable) = eigvals(Hermitian(o.matrix))
-unscaled(o::HermitianObservable) = o
-function ir(ho::HermitianObservable, target::QubitSet, ::Val{:OpenQASM}; serialization_properties=OpenQASMSerializationProperties())
-    function c_str(x::Number)
-        iszero(real(x)) && iszero(imag(x)) && return "0im"
-        iszero(real(x)) && return string(imag(x))*"im"
-        iszero(imag(x)) && return string(real(x))*"+0im"
-        imag_sign = imag(x) < 0 ? "" : "+"
-        return string(real(x)) * imag_sign * string(imag(x)) * "im"
-    end
-    m = "[" * join(["["*join(c_str.(ho.matrix[i, :]), ", ")*"]" for i in 1:size(ho.matrix, 1)], ", ") * "]"
-    t = isempty(target) ? "all" : format_qubits(target, serialization_properties)
-    return "hermitian($m) $t"
-end
-Base.:(*)(o::HermitianObservable, n::Real) = HermitianObservable(Float64(n) .* o.matrix)
-ir(ho::HermitianObservable, target::Nothing, ::Val{:OpenQASM}; kwargs...) = ir(ho, QubitSet(), Val(:OpenQASM); kwargs...)
-chars(o::HermitianObservable) = ("Hermitian",)
 
 """
     Sum <: Observable
@@ -198,16 +218,18 @@ Braket.Observables.Sum()
 struct Sum <: Observable
     summands::Vector{Observable}
     coefficient::Float64
-    function Sum(observables)
-        flattened_observables = Observable[]
-        for obs in observables
-            if obs isa Sum
-                append!(flattened_observables, obs.summands)
-            else
-                push!(flattened_observables, obs)
+    function Sum(v)
+        if VERSION >= v"1.9.0"
+            flattened_v = Iterators.flatmap(v) do obs
+                obs isa Sum ? obs.summands : (obs,)
             end
+        else
+            mapped_v = map(v) do obs
+                obs isa Sum ? obs.summands : (obs,)
+            end
+            flattened_v = Iterators.flatten(mapped_v)
         end
-        new(flattened_observables, 1.0)
+        new(collect(flattened_v), 1.0)
     end
 end
 Base.length(s::Sum) = length(s.summands)
@@ -236,6 +258,7 @@ function Base.show(io::IO, s::Sum)
     print(io, ")")
     return
 end
+StructTypes.lower(s::Sum) = [StructTypes.lower(summand) for summand in s.summands]
 ir(s::Sum, target::Vector{QubitSet}, ::Val{:JAQCD}; kwargs...) = throw(ErrorException("Sum observables are not supported in JAQCD."))
 ir(s::Sum, target::Vector{<:IntOrQubit}, ::Val{:JAQCD}; kwargs...) = throw(ErrorException("Sum observables are not supported in JAQCD."))
 function ir(s::Sum, target::Vector{QubitSet}, ::Val{:OpenQASM}; kwargs...)
@@ -254,13 +277,23 @@ chars(s::Sum) = ("Sum",)
 StructTypes.StructType(::Type{Observable}) = StructTypes.AbstractType()
 StructTypes.subtypes(::Type{Observable}) = (i=I, h=H, x=X, y=Y, z=Z)
 StructTypes.lowertype(::Type{O}) where {O<:Observable} = IRObservable
-function StructTypes.constructfrom(::Type{Observable}, obj::IRObservable)
-    obj isa String && return StructTypes.subtypes(Observable)[Symbol(obj)]()
-    length(obj) == 1 && obj[1] isa Vector{Vector{Vector{Float64}}} && return HermitianObservable(obj[1])
-    return TensorProduct([StructTypes.constructfrom(Observable, (o isa String ? o : convert(Vector{Union{String, Vector{Vector{Vector{Float64}}}}}, [o]))) for o in obj])
+function StructTypes.constructfrom(::Type{Observable}, obj::String)
+    (obj == "i" || obj == "I") && return I()
+    (obj == "x" || obj == "X") && return X()
+    (obj == "y" || obj == "Y") && return Y()
+    (obj == "z" || obj == "Z") && return Z()
+    (obj == "h" || obj == "H") && return H()
+    throw(ArgumentError("Observable of type \"$obj\" provided, only \"i\", \"x\", \"y\", \"z\", and \"h\" are valid."))
+end
+StructTypes.constructfrom(::Type{Observable}, obj::Vector{String}) = length(obj) == 1 ? StructTypes.constructfrom(Observable, first(obj)) : TensorProduct(obj)
+StructTypes.constructfrom(::Type{Observable}, obj::Vector{Vector{Vector{Float64}}}) = HermitianObservable(obj)
+StructTypes.constructfrom(::Type{Observable}, obj::Vector{Vector{Vector{Vector{Float64}}}}) = length(obj) == 1 ? HermitianObservable(obj[1]) : TensorProduct(obj)
+function StructTypes.constructfrom(::Type{Observable}, obj::Vector{T}) where {T}
+    length(obj) == 1 && return StructTypes.constructfrom(Observable, obj[1])
+    return TensorProduct([StructTypes.constructfrom(Observable, o) for o in obj])
 end
 
-Base.:(*)(o1::Observable, o2::Observable) = TensorProduct([o1, o2])
+Base.:(*)(o1::O1, o2::O2) where {O1<:Observable, O2<:Observable} = TensorProduct([o1, o2], 1.0)
 Base.:(*)(n::Real, o::Observable) = o*n 
 Base.:(+)(o1::Observable, o2::Observable) = Sum([o1, o2])
 Base.:(-)(o1::Observable, o2::Observable) = Sum([o1, -1.0 * o2])
diff --git a/src/operators.jl b/src/operators.jl
index 4225927b..05d59b1f 100644
--- a/src/operators.jl
+++ b/src/operators.jl
@@ -20,7 +20,38 @@ ir(o::Operator, t::QubitSet; kwargs...)  = ir(o, t, Val(IRType[]); kwargs...)
 ir(o::Operator, t::IntOrQubit; kwargs...)          = ir(o, QubitSet(t), Val(IRType[]); kwargs...)
 ir(o::Operator, t::IntOrQubit, args...; kwargs...) = ir(o, QubitSet(t), args...; kwargs...)
 
-function pauli_eigenvalues(n::Int)
-    n == 1 && return [1.0, -1.0]
-    return reduce(vcat, [pauli_eigenvalues(n - 1), -1 .* pauli_eigenvalues(n - 1)])
+struct PauliEigenvalues{N}
+    coeff::Float64
+    PauliEigenvalues{N}(coeff::Float64=1.0) where {N} = new(coeff)
 end
+PauliEigenvalues(::Val{N}, coeff::Float64=1.0) where {N} = PauliEigenvalues{N}(coeff)
+Base.length(p::PauliEigenvalues{N}) where {N} = 2^N
+function Base.iterate(p::PauliEigenvalues{N}, ix::Int=1) where {N}
+    return ix <= length(p) ? (p[ix], ix+1) : nothing
+end
+
+Base.getindex(p::PauliEigenvalues{1}, i::Int)::Float64 = getindex((p.coeff, -p.coeff), i)
+function Base.getindex(p::PauliEigenvalues{N}, i::Int)::Float64 where N
+    i_block = div(i-1, 2)
+    split = div(2^(N-1)-1, 2)
+    if N < 5
+        total_evs = 2^N
+        is_front = !isodd(mod(i-1, 2))
+        ev = is_front ? p.coeff : -p.coeff
+        mi = mod(i_block, 2)
+        di = div(i_block, 2)
+        if i_block <= split
+            return isodd(mi) ⊻ isodd(di) ? -ev : ev
+        else
+            mi = mod(i_block - split - 1, 2)
+            di = div(i_block - split - 1, 2)
+            return isodd(mi) ⊻ isodd(di) ? ev : -ev
+        end
+    else
+        new_i = i > 2^(N-1) ? i - 2^(N-1) : i
+	one_down_pe = PauliEigenvalues(Val(N-1))
+	one_down = one_down_pe[new_i]
+        return i_block <= split ? one_down : -one_down
+    end
+end
+Base.getindex(p::PauliEigenvalues{N}, ix::Vector{Int}) where {N} = [p[i] for i in ix]
diff --git a/src/qubit_set.jl b/src/qubit_set.jl
index 8b34de91..82641fcd 100644
--- a/src/qubit_set.jl
+++ b/src/qubit_set.jl
@@ -106,7 +106,8 @@ function Base.show(io::IO, qs::QubitSet)
     print(io, join([sprint(show, q) for q in qs], ", "))
     print(io, ")")
 end
-
+Base.convert(::Type{QubitSet}, q::Integer) = QubitSet(q)
+Base.convert(::Type{QubitSet}, v::Vector{<:Integer}) = QubitSet(v)
 Base.sort(qs::QubitSet; kwargs...) = QubitSet(sort(collect(qs); kwargs...))
 
 const VecOrQubitSet = Union{Vector, QubitSet}
diff --git a/src/raw_schema.jl b/src/raw_schema.jl
index 1f5e8b85..872ed386 100644
--- a/src/raw_schema.jl
+++ b/src/raw_schema.jl
@@ -167,6 +167,8 @@ const IRObservable = Union{Vector{Union{String, Vector{Vector{Vector{Float64}}}}
 Base.convert(::Type{IRObservable}, v::Vector{String}) = convert(Vector{Union{String, Vector{Vector{Vector{Float64}}}}}, v)
 Base.convert(::Type{IRObservable}, s::String)         = s
 Base.convert(::Type{IRObservable}, v::Vector{Vector{Vector{Vector{Float64}}}}) = convert(Vector{Union{String, Vector{Vector{Vector{Float64}}}}}, v)
+Base.:(==)(s::String, iro::IRObservable) = all(iro .== s)
+Base.:(==)(iro::IRObservable, s::String) = all(iro .== s)
 
 abstract type CompilerDirective <: AbstractIR end
 StructTypes.StructType(::Type{CompilerDirective}) = StructTypes.AbstractType()
@@ -240,9 +242,9 @@ StructTypes.defaults(::Type{T}) = Dict{Symbol, Any}(:type => "t")
 
 struct Program <: AbstractProgram
     braketSchemaHeader::braketSchemaHeader
-    instructions::Vector{Any}
+    instructions::Vector{<:Any}
     results::Union{Nothing, Vector{AbstractProgramResult}}
-    basis_rotation_instructions::Union{Nothing, Vector{Any}}
+    basis_rotation_instructions::Union{Nothing, Vector{<:Any}}
 end
 StructTypes.StructType(::Type{Program}) = StructTypes.UnorderedStruct()
 StructTypes.defaults(::Type{Program}) = Dict{Symbol, Any}(:braketSchemaHeader => braketSchemaHeader("braket.ir.jaqcd.program", "1"))
diff --git a/src/raw_task_result_types.jl b/src/raw_task_result_types.jl
index efa1544d..528ad6e1 100644
--- a/src/raw_task_result_types.jl
+++ b/src/raw_task_result_types.jl
@@ -39,7 +39,7 @@ struct GateModelQuantumTaskResult <: AbstractQuantumTaskResult
     additional_metadata::AdditionalMetadata
     result_types::Union{Nothing, Vector{ResultTypeValue}}
     values::Union{Nothing, Vector{Any}}
-    measurements::Union{Nothing, Array}
+    measurements::Union{Nothing, Matrix{Int}, Vector{Vector{Int}}}
     measured_qubits::Union{Nothing, Vector{Int}}
     measurement_counts::Union{Nothing, Accumulator}
     measurement_probabilities::Union{Nothing, Dict{String, Float64}}
diff --git a/src/results.jl b/src/results.jl
index 3ccad38a..54c0e144 100644
--- a/src/results.jl
+++ b/src/results.jl
@@ -41,9 +41,15 @@ for (typ, ir_typ, label) in ((:Expectation, :(Braket.IR.Expectation), "expectati
           - An `Int` or `Qubit`
           - Absent, in which case the observable `o` will be applied to all qubits provided it is a single qubit observable.
         """ $typ(o, targets) = $typ(o, QubitSet(targets))
-        $typ(o::String, targets::QubitSet) = $typ(Observables.TensorProduct([o for ii in 1:length(targets)]), targets)
+        function $typ(o::String, targets::QubitSet)
+            return if !isempty(targets)
+                $typ(Observables.TensorProduct([o for ii in 1:length(targets)]), targets)
+            else
+                $typ(StructTypes.constructfrom(Observable, o), targets)
+            end
+        end
         $typ(o::Vector{String}, targets::QubitSet) = $typ(Observables.TensorProduct(o), targets)
-        $typ(o) = $typ(o, QubitSet())
+        $typ(o::Union{Observables.Observable, String, Vector{String}}) = $typ(o, QubitSet())
         Base.:(==)(e1::$typ, e2::$typ) = (e1.observable == e2.observable && e1.targets == e2.targets)
         StructTypes.StructType(::Type{$typ}) = StructTypes.CustomStruct()
         StructTypes.lower(x::$typ) = $ir_typ(StructTypes.lower(x.observable), (isempty(x.targets) ? nothing : Int.(x.targets)), $label)
@@ -52,7 +58,7 @@ for (typ, ir_typ, label) in ((:Expectation, :(Braket.IR.Expectation), "expectati
             obs_ir = ir(r.observable, r.targets, Val(:OpenQASM); serialization_properties=serialization_properties)
             return "#pragma braket result " * $label * " $obs_ir"
         end
-        $typ(x::union_obs_typ) = $typ(StructTypes.constructfrom(Observables.Observable, x.observable), x.targets)
+        $typ(x::union_obs_typ) = $typ(StructTypes.constructfrom(Observables.Observable, x.observable), QubitSet(x.targets))
         chars(x::$typ) = (uppercasefirst($label) * "(" * chars(x.observable)[1] * ")",) 
     end
 end
@@ -187,7 +193,7 @@ function ir(ag::AdjointGradient, ::Val{:OpenQASM}; serialization_properties::Ser
 end
 StructTypes.StructType(::Type{AdjointGradient}) = StructTypes.CustomStruct()
 function StructTypes.lower(x::AdjointGradient)
-    lowered_obs = ir(x.observable, Val(:JAQCD))
+    lowered_obs     = StructTypes.lower(x.observable)
     lowered_targets = (isempty(x.targets) ? nothing : convert(Vector{Vector{Int}}, x.targets))
     Braket.IR.AdjointGradient(x.parameters, lowered_obs, lowered_targets, "adjoint_gradient")
 end
@@ -196,7 +202,7 @@ function AdjointGradient(x::@NamedTuple{parameters::Union{Nothing, Vector{String
     params = isnothing(x.parameters) || isempty(x.parameters) ? ["all"] : x.parameters
     obs = StructTypes.constructfrom(Observables.Observable, x.observable)
     targets = isnothing(x.targets) || isempty(x.targets) ? QubitSet[] : [QubitSet(t) for t in x.targets]
-    return AdjointGradient(obs, targets, parameters)
+    return AdjointGradient(obs, targets, params)
 end
 chars(x::AdjointGradient) = ("AdjointGradient(H)",) 
 
@@ -259,7 +265,7 @@ StateVector(x::@NamedTuple{type::String}) = StateVector()
 ir(r::Result, ::Val{:JAQCD}; kwargs...) = StructTypes.lower(r)
 ir(r::Result; kwargs...) = ir(r, Val(IRType[]); kwargs...)
 StructTypes.StructType(::Type{Result}) = StructTypes.AbstractType()
-StructTypes.subtypes(::Type{Result}) = (amplitude=Amplitude, expectation=Expectation, probability=Probability, statevector=StateVector, variance=Variance, sample=Sample, densitymatrix=DensityMatrix)
+StructTypes.subtypes(::Type{Result}) = (amplitude=Amplitude, expectation=Expectation, probability=Probability, statevector=StateVector, variance=Variance, sample=Sample, densitymatrix=DensityMatrix, adjoint_gradient=AdjointGradient)
 
 function StructTypes.constructfrom(::Type{R}, obj) where {R<:Result}
     return R((getproperty(obj, fn) for fn in fieldnames(R))...)
diff --git a/src/schemas.jl b/src/schemas.jl
index 949b2e5b..25b1b582 100644
--- a/src/schemas.jl
+++ b/src/schemas.jl
@@ -20,30 +20,35 @@ julia> Instruction(StartVerbatimBox(), QubitSet())
 Braket.Instruction(StartVerbatimBox(), QubitSet())
 ```
 """
-struct Instruction
-    operator::Operator
+struct Instruction{O<:Operator}
+    operator::O
     target::QubitSet
-    Instruction(o::Operator, target) = new(o, QubitSet(target))
 end
-Instruction(cd::CompilerDirective) = Instruction(cd, Int[])
-operator(ix::Instruction) = ix.operator
+Instruction(o::O, target...) where {O<:Operator} = Instruction{O}(o, QubitSet(target...))
+Instruction(o::O, target) where {O<:Operator} = Instruction{O}(o, QubitSet(target...))
+Instruction{CD}(cd::CD) where {CD<:CompilerDirective} = Instruction{CD}(cd, QubitSet(Int[]))
+Instruction(cd::CD) where {CD<:CompilerDirective} = Instruction{CD}(cd, Int[])
+operator(ix::Instruction{O}) where {O<:Operator} = ix.operator
+StructTypes.StructType(::Type{Instruction{O}}) where {O} = StructTypes.CustomStruct()
 StructTypes.StructType(::Type{Instruction}) = StructTypes.CustomStruct()
-StructTypes.lower(x::Instruction) = isempty(x.target) ? ir(x.operator, Val(:JAQCD)) : ir(x.operator, x.target, Val(:JAQCD))
-ir(x::Instruction, ::Val{:OpenQASM}; kwargs...) = isempty(x.target) ? ir(x.operator, Val(:OpenQASM); kwargs...) : ir(x.operator, x.target, Val(:OpenQASM); kwargs...)
+StructTypes.lower(x::Instruction{O}) where {O<:Operator} = isempty(x.target) ? ir(x.operator, Val(:JAQCD)) : ir(x.operator, x.target, Val(:JAQCD))
+ir(x::Instruction{O}, ::Val{:OpenQASM}; kwargs...) where {O<:Operator} = isempty(x.target) ? ir(x.operator, Val(:OpenQASM); kwargs...) : ir(x.operator, x.target, Val(:OpenQASM); kwargs...)
 
 conc_types = filter(Base.isconcretetype, vcat(subtypes(AbstractIR), subtypes(CompilerDirective)))
 nt_dict = merge([Dict(zip(fieldnames(t), (Union{Nothing, x} for x in fieldtypes(t)))) for t in conc_types]...)
 ks = tuple(:angles, keys(nt_dict)...)
 vs = Tuple{Union{Nothing, Vector{Float64}}, values(nt_dict)...}
 inst_typ = NamedTuple{ks, vs}
+StructTypes.lowertype(::Type{Instruction{O}}) where {O} = inst_typ
 StructTypes.lowertype(::Type{Instruction}) = inst_typ
+Instruction(x::Instruction{O}) where {O<:Braket.Operator} = x
 function Instruction(x)
     sts    = merge(StructTypes.subtypes(Gate), StructTypes.subtypes(Noise), StructTypes.subtypes(CompilerDirective))
     o_type = sts[Symbol(x[:type])]
     (o_type <: CompilerDirective) && return Instruction(o_type(), Int[])
     if o_type <: AngledGate{1}
-        o_fns  = (:angles, :angle) 
-        args   = (x[:angle],)
+        o_fns  = (:angles, :angle)
+        args   = (Float64(x[:angle]),)
         op     = o_type(args)
     elseif o_type <: AngledGate{3}
         o_fns  = (:angles, :angle1, :angle2, :angle3) 
@@ -63,14 +68,21 @@ function Instruction(x)
     target = reduce(vcat, raw_target)
     return Instruction(op, target)
 end
+Instruction(x::Dict{String, <:Any}) = Instruction(Dict(Symbol(k)=>v for (k,v) in x))
 StructTypes.constructfrom(::Type{Instruction}, x::Union{Dict{Symbol, Any}, NamedTuple}) = Instruction(x)
-Base.:(==)(ix1::Instruction, ix2::Instruction) = (ix1.operator == ix2.operator && ix1.target == ix2.target)
+Base.:(==)(ix1::Instruction{O}, ix2::Instruction{O}) where {O<:Operator} = (ix1.operator == ix2.operator && ix1.target == ix2.target)
+qubits(ix::Instruction{O}) where {O<:Operator} = ix.target
+qubit_count(ix::Instruction{O}) where {O<:Operator} = length(ix.target)
+qubits(ixs::Vector{Instruction{O}}) where {O<:Operator}  = mapreduce(ix->ix.target, union!, ixs, init=Set{Int}())
+qubits(ixs::Vector{Instruction}) = mapreduce(ix->ix.target, union!, ixs, init=Set{Int}())
+qubit_count(ixs::Vector{Instruction{O}}) where {O<:Operator} = length(qubits(ixs)) 
+qubit_count(ixs::Vector{Instruction}) = length(qubits(ixs)) 
 
-bind_value!(ix::Instruction, param_values::Dict{Symbol, Number}) = Instruction(bind_value!(ix.operator, param_values), ix.target)
+bind_value!(ix::Instruction{O}, param_values::Dict{Symbol, Number}) where {O<:Operator} = Instruction{O}(bind_value!(ix.operator, param_values), ix.target)
 
-remap(ix::Instruction, mapping::Dict{<:Integer, <:Integer}) = Instruction(copy(ix.operator), [mapping[q] for q in ix.target])
-remap(ix::Instruction, target::VecOrQubitSet) = Instruction(copy(ix.operator), target[1:length(ix.target)])
-remap(ix::Instruction, target::IntOrQubit)    = Instruction(copy(ix.operator), target)
+remap(ix::Instruction{O}, mapping::Dict{<:Integer, <:Integer}) where {O<:Operator} = Instruction{O}(copy(ix.operator), [mapping[q] for q in ix.target])
+remap(ix::Instruction{O}, target::VecOrQubitSet) where {O<:Operator} = Instruction{O}(copy(ix.operator), target[1:length(ix.target)])
+remap(ix::Instruction{O}, target::IntOrQubit) where {O<:Operator}    = Instruction{O}(copy(ix.operator), target)
 
 function StructTypes.constructfrom(::Type{Program}, obj)
     new_obj = copy(obj)
diff --git a/src/task.jl b/src/task.jl
index a500947f..16bb21af 100644
--- a/src/task.jl
+++ b/src/task.jl
@@ -375,7 +375,7 @@ function measurement_probabilities(measurement_counts::Accumulator)
 end
 
 function _measurements(probs::Dict{String, Float64}, shots::Int)
-    measurements = []
+    measurements = Vector{Int}[]
     for (bitstring, prob) in probs
         int_list = [tryparse(Int, string(b)) for b in bitstring]
         measurement = [int_list for ii in 1:convert(Int, round(prob*shots))]
@@ -443,13 +443,13 @@ end
 
 function calculate_result_types(ir_obj, measurements, measured_qubits::Vector{Int})
     result_types = ResultTypeValue[]
-    (!haskey(ir_obj, "results") || isnothing(ir_obj["results"])) && return result_types
-    for result_type in ir_obj["results"]
-        ir_obs = get(result_type, "observable", nothing)
+    (!haskey(ir_obj, :results) || isnothing(ir_obj[:results])) && return result_types
+    for result_type in ir_obj[:results]
+        ir_obs = get(result_type, :observable, nothing)
         rt     = JSON3.read(JSON3.write(result_type), AbstractProgramResult)
         obs    = isnothing(ir_obs) ? nothing : StructTypes.constructfrom(Observables.Observable, rt.observable)
         targs  = rt.targets
-        typ    = result_type["type"]
+        typ    = result_type[:type]
         if typ == "probability"
             val = _calculate_for_targets(measurements, measured_qubits, targs, Val(Symbol(typ)))
             push!(result_types, ResultTypeValue(rt, val))
@@ -486,11 +486,15 @@ function computational_basis_sampling(::Type{GateModelQuantumTaskResult}, r::Gat
     end
     measured_qubits = r.measuredQubits
     if isnothing(r.resultTypes) || isempty(r.resultTypes)
-        result_types = calculate_result_types(JSON3.read(JSON3.write(addl_mtd.action)), measurements, measured_qubits)
+        # do this in case we have custom gates in the IR
+        json_          = hasproperty(addl_mtd.action, :results) ? Dict("results"=>[rt for rt in addl_mtd.action.results]) : addl_mtd.action
+        rebuilt_action = JSON3.read(JSON3.write(json_))
+        result_types   = calculate_result_types(rebuilt_action, measurements, measured_qubits)
     else
         if !isempty(r.resultTypes) && !isnothing(r.resultTypes[1])
-            json_ = Dict("results"=>[rt.type for rt in r.resultTypes])
-            result_types = calculate_result_types(JSON3.read(JSON3.write(json_)), measurements, measured_qubits)
+            json_          = Dict("results"=>[rt.type for rt in r.resultTypes])
+            rebuilt_action = JSON3.read(JSON3.write(json_))
+            result_types   = calculate_result_types(rebuilt_action, measurements, measured_qubits)
         else
             result_types = r.resultTypes
         end
diff --git a/src/utils.jl b/src/utils.jl
index 2d092742..be944239 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -99,3 +99,11 @@ function complex_matrix_from_ir(mat::Vector{Vector{Vector{T}}}) where {T<:Number
     end
     return m
 end
+
+function complex_matrix_from_ir(mat::Vector{Vector{Vector{Any}}})
+    m = zeros(ComplexF64, length(mat), length(mat))
+    for ii in 1:length(mat), jj in 1:length(mat)
+        m[ii,jj] = complex(mat[ii][jj][1], mat[ii][jj][2])
+    end
+    return m
+end
diff --git a/test/observables.jl b/test/observables.jl
index 58f80231..5dfd3a0a 100644
--- a/test/observables.jl
+++ b/test/observables.jl
@@ -1,21 +1,26 @@
 using Braket, Braket.Observables, Test, JSON3, StructTypes, LinearAlgebra
-using Braket: VIRTUAL, PHYSICAL, OpenQASMSerializationProperties, pauli_eigenvalues, IRObservable
+using Braket: VIRTUAL, PHYSICAL, OpenQASMSerializationProperties, PauliEigenvalues, IRObservable
 using LinearAlgebra: eigvals
 
 @testset "Observables" begin
     Braket.IRType[] = :JAQCD
     @testset "pauli eigenvalues" begin
         z = [1.0 0.0; 0.0 -1.0]
-        @test pauli_eigenvalues(1) == diag(z)
-        @test pauli_eigenvalues(2) == diag(kron(z,z))
-        @test pauli_eigenvalues(3) == diag(kron(z,z,z,))
+        for n in 2:6
+            pe  = PauliEigenvalues(Val(n))
+            mat = kron(ntuple(i->diag(z), n)...)
+            for ix in 1:2^n
+                @test pe[ix] == mat[ix]
+            end
+        end
     end
     @testset "Hermitian" begin
         m = [1. -im; im -1.]
         o = Observables.HermitianObservable(m)
         @test qubit_count(o) == 1
         rt = Expectation(o, [0])
-        @test JSON3.read(JSON3.write(rt), Braket.Result) == rt
+        read_rt = JSON3.read(JSON3.write(rt), Braket.Result)
+        @test read_rt == rt
         @test ir(o) isa IRObservable
         mult_h = 2.0 * o
         @test mult_h.matrix == 2.0 * m
@@ -100,10 +105,12 @@ using LinearAlgebra: eigvals
             (Observables.TensorProduct(["x", "y", "i"]),[1, 1, -1, -1, -1, -1, 1, 1]),
             (Observables.TensorProduct([Observables.X(), ho, Observables.Y()]),[-1, 1, -1, 1, 1, -1, 1, -1, 1, -1, 1, -1, -1, 1, -1, 1])
         )
-            @test eigvals(tp) == convert(Vector{Float64}, evs)
+            for ix in 1:length(evs)
+                @test eigvals(tp)[ix] == convert(Float64, evs[ix])
+            end
         end
         @testset for typ in (Observables.H(), Observables.X(), Observables.Y(), Observables.Z())
-            @test eigvals(typ) == [1.0, -1.0]
+            @test eigvals(typ)[[1, 2]] == [1.0, -1.0]
         end
         @testset for (mat, evs) in [([1.0 0.0; 0.0 1.0], [1, 1]), ([0 -im; im 0], [-1.0, 1.0]), ([1 1-im; 1+im -1], [-sqrt(3), sqrt(3)])]
             @test eigvals(Observables.HermitianObservable(mat)) ≈ evs
@@ -111,7 +118,7 @@ using LinearAlgebra: eigvals
     end
     @testset "OpenQASM" begin
         @testset for ir_bolus in [
-            (Observables.I(), OpenQASMSerializationProperties(qubit_reference_type=VIRTUAL), [3], "i(q[3])"),
+            ( Observables.I(), OpenQASMSerializationProperties(qubit_reference_type=VIRTUAL), [3], "i(q[3])"),
             ( Observables.I(), OpenQASMSerializationProperties(qubit_reference_type=PHYSICAL), [3], "i(\$3)"),
             ( Observables.I(), OpenQASMSerializationProperties(qubit_reference_type=VIRTUAL), nothing, "i all"),
             ( Observables.X(), OpenQASMSerializationProperties(qubit_reference_type=VIRTUAL), [3], "x(q[3])"),