1+ import numpy as np
2+ from qiskit import QuantumCircuit , QuantumRegister , ClassicalRegister , transpile
3+ from qiskit_aer import AerSimulator
4+ import matplotlib .pyplot as plt
5+ from itertools import product
6+ from collections import defaultdict
7+ from sklearn .linear_model import Ridge
8+ from reservoirpy .datasets import narma ,logistic_map
9+
10+
11+
12+
13+ def fit_model (model , res_states , series , WARMUP , timeplex = 1 ):
14+ warmup = int (len (series ) * WARMUP )
15+
16+ X = res_states [warmup :- 1 ]
17+ y = series [warmup + 1 :]
18+
19+ model .fit (X , y )
20+
21+ return model , X , y
22+
23+
24+ def run_prediction (model , res_states , timeplex = 1 ):
25+
26+ X = np .copy (res_states )
27+
28+ X = X [- 1 ,:]
29+ X = X .reshape ((1 , - 1 ))
30+ return model .predict (X )
31+
32+
33+
34+ def compute_z_expectations (memory , num_measurement , num_steps ):
35+ timestep_counts = [defaultdict (int ) for _ in range (num_steps )]
36+
37+ for shot in memory :
38+ steps = shot .strip ().split () # e.g., ['01', '11', '00', ...] — one per timestep
39+ for step_idx , bits in enumerate (steps ):
40+ # Reverse bitstring if using little-endian (rightmost qubit is qubit 0)
41+ bitstring = bits [::- 1 ][:num_measurement ]
42+ timestep_counts [step_idx ][bitstring ] += 1
43+
44+ # Convert defaultdicts to regular dicts (optional)
45+ return [dict (d ) for d in timestep_counts ]
46+
47+ def pauli_z_expectation (counts , pauli_str ):
48+ total_counts = sum (counts .values ())
49+ expectation = 0.0
50+
51+ for bitstring , count in counts .items ():
52+ # Compute product of eigenvalues for this bitstring
53+ val = 1
54+ # bitstring assumed little-endian (bit 0 is rightmost), adjust if needed
55+ for i , p in enumerate (pauli_str ):
56+ if p == 'Z' :
57+ # bit at position i in little endian (rightmost = qubit 0)
58+ bit = bitstring [::- 1 ][i ]
59+ val *= 1 if bit == '0' else - 1
60+ elif p == 'I' :
61+ val *= 1
62+ else :
63+ raise ValueError ("Pauli string must be 'I' or 'Z' only" )
64+ expectation += val * count
65+
66+ return expectation / total_counts
67+
68+ def compute_expectations_for_all_timesteps (timestep_counts , pauli_strings ):
69+ all_expectations = []
70+
71+ for counts in timestep_counts :
72+ expectations = []
73+ for p_str in pauli_strings :
74+ exp_val = pauli_z_expectation (counts , p_str )
75+ expectations .append (exp_val )
76+ all_expectations .append (expectations )
77+
78+ return all_expectations [::- 1 ]
79+
80+ def generate_ZI_pauli_strings (n ):
81+ pauli_strings = ['' .join (p ) for p in product ('ZI' , repeat = n )]
82+ return pauli_strings
83+
84+
85+
86+
87+ def encode_circuit (t ,num_qubits ):
88+ qr = QuantumRegister (num_qubits )
89+ qc = QuantumCircuit (qr )
90+
91+ for k in range (num_qubits ):
92+ beta = 2 ** (- (k )/ num_qubits )
93+ #beta=1
94+ #qc.rx(t * np.pi * beta, num_qubits-k-1)
95+ qc .rx (t * np .pi * beta , k )
96+ return qc
97+
98+ def quantum_part_circuit (timeseries , num_qubits , num_measurement , reservoir_circuit ,type = 'quantum_part' ):
99+ if type == 'quantum_part' :
100+ qr = QuantumRegister (num_qubits , "q" )
101+ crs = [ClassicalRegister (num_measurement , f"c{ i } ) for i in range (len (timeseries ))]
102+ qc = QuantumCircuit (qr , * crs )
103+
104+ for step_idx , t in enumerate (timeseries ):
105+ #reservoir = get_Ising_circuit(num_qubits, isingparams)
106+ qc .compose (encode_circuit (t , num_qubits ), qubits = qr , inplace = True )
107+ qc .compose (reservoir_circuit , qubits = qr , inplace = True )
108+
109+ # Mid measurement: store each step's result in a different ClassicalRegister
110+ for i in range (num_measurement ):
111+ qbit = num_qubits - num_measurement + i
112+ qc .measure (qr [qbit ], crs [step_idx ][i ])
113+ qc .reset (qr [qbit ]) # Reset qubit to |0> after measurement
114+
115+ return qc
116+ else :
117+ #res = Stabilizer(num_qubits, num_measurement, backend = AerSimulator(noise_model=None),\
118+ #degree=num_measurement, num_reservoirs=1, isingparams=ising_params,decode=True)
119+ qc = reservoir_circuit .circuit (timeseries )
120+ return qc
121+
122+ def execute_reservoir (timeseries ,num_measurement , num_shots , circuit ):
123+ backend = AerSimulator ()
124+
125+
126+ result = backend .run (circuit , shots = num_shots , memory = True ).result ()
127+
128+ # Access full memory results (per shot)
129+ memory = result .get_memory ()
130+ #print(memory)
131+ z_expectations = compute_z_expectations (
132+ memory = memory ,
133+ num_measurement = num_measurement ,
134+ num_steps = len (timeseries )
135+ )
136+ return z_expectations
137+
138+
139+ def sine (n ):
140+ ts = []
141+ for i in range (1 ,n ):
142+ ts .append ((np .sin (7.35 * i * np .pi / n )+ 1 )/ 2 )
143+ #ts.append(0.5)
144+ return np .array (ts )
145+
146+ def henon1d (n , a = 1.4 , b = 0.3 ):
147+ ts = [0 ,0 ]
148+ for i in range (2 ,n + 2 ):
149+ ts .append (1 - a * ts [i - 1 ]** 2 + b * ts [i - 2 ])
150+ return np .array (ts [2 :])
151+
152+ def narma_task (n ,order ):
153+ rng = np .random .default_rng (seed = 2341 )
154+ u = rng .uniform (0 , 0.5 , size = (n + order , 1 ))
155+ y = narma (n , order = order , u = u )
156+ print (y )
157+ return np .array (y .flatten ()[order :])
158+
159+ def logistic (n ):
160+ y = logistic_map (n_timesteps = n )
161+ print (y )
162+ return np .array (y .flatten ())
163+
164+
165+
166+ # Prediction function
167+ def predict_one_step_ahead (model , timeseries_aux , num_qubits , num_measurement , reservoirs ,backend , method , ising_params = None , shots = 10000 ):
168+ ts_list = timeseries_aux .flatten ().tolist ()
169+ k = 1
170+ obs = []
171+ num_tot_obs = len (generate_ZI_pauli_strings (num_measurement ))
172+ for i ,reservoir_circuit in enumerate (reservoirs ):
173+ circuit = reservoir_circuit .circuit (ts_list )
174+ circuit = transpile (circuit , backend )
175+ counts = execute_reservoir (ts_list ,num_measurement ,shots , circuit )
176+ obs .append (compute_expectations_for_all_timesteps (counts , generate_ZI_pauli_strings (num_measurement )[1 :num_tot_obs - 1 ]))
177+
178+
179+ Xi = []
180+ for j in range (k ):
181+ for r in range (0 , len (reservoirs )):
182+ Xi .extend (obs [r ][len (obs [0 ])- j - 1 ]) # take last k observations
183+
184+ Xi = np .array (Xi [::- 1 ]).reshape (1 ,- 1 )
185+ y_pred = model .predict (Xi )
186+ if k == 1 :
187+ return y_pred [- 1 ]
188+ else :
189+ return y_pred [- 1 ][- 1 ]
0 commit comments