-
Notifications
You must be signed in to change notification settings - Fork 0
/
correlation.jl
1987 lines (1576 loc) · 58.1 KB
/
correlation.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
### A Pluto.jl notebook ###
# v0.16.1
using Markdown
using InteractiveUtils
# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
quote
local el = $(esc(element))
global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
el
end
end
# ╔═╡ e501c230-d8b5-11eb-37a3-e7c9ecf7bd80
begin
using Plots
using LinearAlgebra
using Kronecker
using PlutoUI
using SparseArrays
#https://juliapackages.com/p/kronecker
end
# ╔═╡ 95a5da6c-08b1-45ab-ac60-957f6671f1f3
TableOfContents()
# ╔═╡ eb3906c7-d1e5-479b-93c4-5b59707c55f1
md""" # Test Kronecker package"""
# ╔═╡ 736b4743-bcb3-42c7-8b4c-067f6a3f18b3
begin
n = 2
m = 3
A = randn(n, n);
B = rand(m, m);
#v = rand(n*m);
K = A ⊗ B
v = rand(size(K,1),1)
K * v # equivalent with vec(B * reshape(v, 50, 100) * A')
# trace(K)
end
# ╔═╡ db75c98c-afe2-478b-9a93-9f87765e43cc
vec(K)
# ╔═╡ 27796834-1dfc-48b6-8005-c4dee17f3490
collect(vec(K))
# ╔═╡ 5c77981b-893a-4d77-b3db-3377afac0b63
begin
function density_matrix(n, neg = 0, real_valued = false)
#trace 1 remains preserved!
dia = rand(n)
if neg > 0
dia[rand(1:n, neg)] .*= -1
end
dia = dia/sum(dia)
if real_valued
temp = rand(n,n)
else
temp = rand(n,n) + rand(n,n) * im
end
ew, ev = eigen(Hermitian(temp'*temp))
d_mat = Hermitian(ev* Diagonal(dia) * ev')
return d_mat
end
function unitary(n)
A = (rand(n,n) .-0.5) + im*(rand(n,n) .-0.5)
return eigvecs(A'+A)
end
function unity(i,n)
E = spzeros(Bool,n,n)
E[i] = true
return E
end
function partial_trace_matrix(dim_1, dim_2, traced_system)
# traces out the system indicated by traced_system
# for C = kron(A,B), tracing over the first system (traced_system = 1);
# yields: reshape(p1*C[:], 4,4)/tr(A) == B
#
# in colex notation:
# traced_system == 1 represents a trace over index 2 and 4 and
# traced_system == 2 represents a trace over index 1 and 3 and
if traced_system == 1
T = [kron(sparse(I, dim_1,dim_1), unity(i,dim_2))[:] for i = 1:dim_2^2]
col = reduce(vcat,[findnz(T[i])[1] for i = 1:dim_2^2])
row = kron(1:dim_2^2, ones(Int64,dim_1))
return sparse(row, col, trues(dim_2^2*dim_1))
else
T = [kron(unity(i, dim_1), sparse(I, dim_2, dim_2))[:] for i = 1:dim_1^2]
col = reduce(vcat, [findnz(T[i])[1] for i = 1:dim_1^2])
row = kron(1:dim_1^2, ones(Int64,dim_2))
return sparse(row, col, trues(dim_1^2*dim_2))
#return Tr_B = [reshape(kron(reshape(unity(i,dim_1), dim_1, dim_1), I(dim_2)), 1, dim_2^2 * dim_1^2)[1,j] for i = 1:dim_1^2, j =1:dim_1^2*dim_2^2]
end
end
function choi(A, dim_1)
mm, nn = size(A)
n = mm ÷ dim_1
return reshape(permutedims(reshape(A,n,dim_1,n,dim_1),(1,3,2,4)), mm,mm)
end
function choi_rect(A,dim_1, dim_2, dim_3, dim_4)
return reshape(permutedims(reshape(A, dim_1, dim_2, dim_3, dim_4), (1,3,2,4)), dim_1*dim_3, dim_2*dim_4)
end
function choi_index(i, n)
return (i.-1) .% n[1] .+ 1 .+
n[1] .* ((i.-1) .% prod(n[1:3]) .÷ prod(n[1:2]) .+
n[2] .* ((i.-1) .% prod(n[1:2]) .÷ n[1] .+
n[3] .* ((i.-1) .÷ prod(n[1:3]))))
end
function ptr(ρ, dim_1, traced_system = 1)
# get partial trace superoperator
N,M = size(ρ)
dim_2 = N÷dim_1
PT = partial_trace_matrix(dim_1, dim_2, traced_system)
new_dim = traced_system == 1 ? dim_2 : dim_1
# apply superoperator to vectorized density matrix and reshape to new dimensions
σ = reshape(PT * ρ[:], new_dim, new_dim)
return σ
end
function factorize(A, dim_1, dim_2)
σ = ptr(A, dim_1, 2)
ρ = ptr(A, dim_1, 1)
A_tensor_prod = kron(σ, ρ)
rho_c = A - A_tensor_prod
return A_tensor_prod, rho_c
end
function normed_vector(k)
v = (rand(k) .- 0.5) + im*(rand(k) .- 0.5)
return v/norm(v)
end
function prob_dist(k)
p = rand(k)
return p/sum(p)
end
md"""
# Functions
"""
end
# ╔═╡ ca92b1da-20cd-40b8-81ea-72c49f8898cd
function PPT(rho,dim_1)
dim_2 = size(rho,1)÷dim_1
if dim_1 *dim_2 > 6
error("positive partial trace does not work for too large dimensions")
end
return eigvals(reshape(permutedims(reshape(rho,dim_1,dim_2,dim_1,dim_2), [1,4,3,2]),dim_1*dim_2, dim_1*dim_2))
end
# ╔═╡ 3a9e432b-7c7b-463b-bf13-e30c017b5880
function sep_state(k=1,n=2)
VVV = sum([kron(normed_vector(n), normed_vector(n)) for _=1:k] .* (prob_dist(k)))
return VVV / norm(VVV)
end
# ╔═╡ 627cb3ec-704b-4d9a-8a93-253b4c68971d
begin
#test for k-block positivity
# strategy: generate random matrices in 2x2 and test whether Choi matrix maps to positive operators, if true for all, then positive.
# alternatively: check k-positivity: k as optional parameter. Create states with Schmidt rank k
#Schmidt rank:
function pos_test(rho, k = 1, n = 2)
runs = 10000
test = zeros(runs)
vector = zeros(ComplexF64, n^2)
for i = 1:runs
#sum([density_matrix(2) for _=1:k] * p_rand/sum(p_rand))
VVN = sep_state(k,n)
test[i] = real(VVN'* rho *VVN)
if real(test[i]) < 0
vector = VVN
end
end
return sum(test .> 0) / runs
end
end
# ╔═╡ 5450c94b-f78b-4bf2-b025-f035a28850cd
md"""
# Create a nice initial correlation term
"""
# ╔═╡ 2e11a379-2639-4cb9-a892-aa610680e671
md"
### Detecting entanglement with positive partial transpose
"
# ╔═╡ b54087af-d08c-402c-a4ff-52899b4787d1
begin
# get r_c an arbitrary correlation term
rr = density_matrix(4)
r_t, r_c = factorize(rr,2,2)
#partial trace of correlation term is zero
#partial transpose is equal to swapping two indices in the four index picture
r_t2, r_c2 = factorize(density_matrix(4),2,2)
aaa, bbb = factorize(1/2*r_t + 1/2*r_t2,2,2)
#convex combination of product states turn to be not product states
if any(PPT(rr,2) .< 0)
md"entanglement detected"
else
md"separable state"
end
end
# ╔═╡ 29fbbc0a-3c97-4b15-9f20-3f66a562e2b7
md"
### get fully entangled state
"
# ╔═╡ cbb8399c-44c0-4def-8458-b6f34c94be66
begin
dim_n = 2
#note: check what for makes with the scope if combined with a markup element!
let
global ent = zeros(dim_n^2, dim_n^2)
for k = 1:dim_n^2
ent += kron(unity(k,dim_n), unity(k,dim_n))
end
ent = ent/tr(ent)
end
#apply a random weight:
#sum(reshape(reduce(hcat, [kron(unity(k,dim_n), unity(k,dim_n)) for k = 1:4]),4,4,4) .* reshape([0.25, 0.3, 0.3, 0.2], 1,1,4), dims=3)
md"get full entagled state"
end
# ╔═╡ 508b3c1e-1835-4575-a59e-bfc4f2504860
begin
R = ent#0.5*kron(density_matrix(2), density_matrix(2)) + 0.5*kron(density_matrix(2), density_matrix(2))
eigen(choi(R, 2))
C = choi(R,2)
end
# ╔═╡ aca27f55-197e-4a27-98d5-be27981c3077
md"
### Create correlated initial state
"
# ╔═╡ 1cd9f410-b75a-4567-aeac-0809110be2d7
md" c23: $(@bind c23 Scrubbable(-1:0.1:1))"
# ╔═╡ 00647334-6d64-4cc3-95ab-f51196300893
md"""
# Time evolution and dynamical map"""
# ╔═╡ cdde600f-d922-4376-81e8-949339378024
begin
diagonal_index = 1:(n+1):n^2
# generate a unit basis of density operators
rho_basis = zeros(ComplexF64,n^2,n^2)
for k = 1:n
for j = 1:n
compl_factor = j>=k ? 1 : im
temp = unity((k-1)*n + j, n)* compl_factor + unity(diagonal_index[k],n) + unity(diagonal_index[j],n)
rho_basis[:,(k-1)*n+j] =((temp + temp')./(2*tr(temp)))[:]
end
end
end
# ╔═╡ 6724ccf7-dbd8-44cb-a96d-a7074b914a81
begin
#
md"""
## quantum--classical correlated initial state of the form
$\rho(t_0) = \sum_j p_j \Pi_j \otimes \rho_{B,j}$
"""
end
# ╔═╡ 4a1ccfb2-bae6-464d-9790-0d09a4dcfbab
md""" $(@bind new_density Button("random density operator")) """
# ╔═╡ 7656ee8c-0b8c-402b-84b0-5517b39aac3a
begin
new_density
rho_a = density_matrix(2)# [0.5 0.1im; -0.1im 0.5]
rho_b = density_matrix(2)
#r1, r2 = factorize(density_matrix(4),2,2)
#r2 does not produce any entanglement!!!
r2 = [0 0.1 -0.1 -0.2; 0.1 0 -0.2 0.1; -0.1 -0.2 0 0.1; -0.2 0.1 0.1 0]
#r2 = [0.25 0.0 0.0 0.5; 0.0 -0.25 0.0 0.0; 0.0 0.0 -0.25 0.0; 0.5 0.0 0.0 0.25]
# for 0.7*r2 always non complete positive
#rho = density_matrix(4) #kron(rho_a, rho_b) + r2*0 #
#rho = 1/4*(kron(I(2), I(2)) + kron( τ₁ * σ₁ + τ₂ * σ₂ + τ₃ * σ₃, I(2)) - c23*kron(σ₂, σ₃))
# Create quantum-classical correlated state
Dim = 2
RHO_b = [density_matrix(Dim) for _=1:Dim]
RHO_s = density_matrix(Dim)
evals, Pro = eigen(RHO_s)
prob = prob_dist(Dim)
rho = sum([kron(Pro[:,i] * Pro[:,i]', RHO_b[i]) for i = 1:Dim] .*prob)
end
# ╔═╡ 07199f9a-1180-4901-848e-3aca182cdf52
r2
# ╔═╡ 228654fc-8422-4f64-964c-c8b69c7fd175
eigvals(rho)
# ╔═╡ cc77f9ed-73ef-4129-9e6c-1a58dff2141d
if any(PPT(rho,2) .< 0)
md"entanglement detected 🧬"
else
md"separable state 📎📎"
end
# ╔═╡ 4e33f911-b190-4771-a1f4-2b15547fc84c
ptr(rho,Dim,2)
# ╔═╡ 56e4959e-586a-42b5-be90-8da1d719d6e5
RHO_b
# ╔═╡ 1289264a-dfaf-4889-a81e-2991295fde84
Pro[:,2]
# ╔═╡ 0c3439c3-948e-4cbd-965c-6b91bd10e11a
kron((Pro[1,:]*Pro[1,:]'),I(3))
# ╔═╡ 36c5e484-275b-4325-8a45-89e23d9f48f8
Pro[:,1] * Pro[:,1]'
# ╔═╡ 519bfd98-a91b-496b-a484-8fa861aef467
eigvals(rho)
# ╔═╡ a421f989-e88d-45ab-a9b8-e79fa321e2d7
if any(PPT(rho,Dim) .< 0)
md"entanglement detected 🧬"
else
md"separable state 📎📎"
end
# ╔═╡ b57f0f09-831b-4c99-ba6a-7b49af111c9b
tr(sqrt(choi(rho,Dim)' * choi(rho,Dim))) # cross norm criteria or realignment criteria (with trace norm)
# ╔═╡ c1906f46-ffd6-4c26-909b-d952066813d9
md"""
``\boldsymbol{K}_i = \sqrt{\varepsilon_{k'}} \langle b_k | U(t,t_0) | b_{k'} \rangle ``
"""
# ╔═╡ bc6dea59-2db2-4d9a-aca9-3ae77159a070
md"""
Time: $(@bind time Slider(0:0.1:40, show_value=true, default = 0.5))
"""
# ╔═╡ d8b26c46-2c3b-4e8c-bfad-f81a0406ac6e
time
# ╔═╡ ec05b0c0-92d4-44ef-80be-ef3d3e53e688
time_vector = 0:0.1:50
# ╔═╡ deed3f43-ddbf-47fc-abb7-3fdfe8016ba9
begin
md"""
Energy of first system site $(@bind eps_1 Slider(-4:0.01:4, show_value=true, default = 1))
Second system: $(@bind eps_2 Slider(-4:0.01:4, show_value=true, default = 2))
Hopping: $(@bind t Slider(-4:0.01:4, show_value=true, default = 0.5))
"""
end
# ╔═╡ 6a3574a9-9797-49b3-891f-d04b4919614f
begin
eps_3 = 3
eps = [0, eps_1, eps_2, eps_3]
if Dim == 2
H_1 = [0 0 0 0; 0 eps[2] 0 0; 0 0 eps[3] 0; 0 0 0 eps[2] + eps[3]]
H_2 = [0 0 0 0; 0 0 t 0; 0 t 0 0; 0 0 0 0]
H_full= H_1 + H_2
U = exp(-im * H_full * time)
else
U = unitary(Dim^2)
end
md" time slider: $(time)"
end
# ╔═╡ 90a2efc4-1706-4882-a1d4-c41f42f9c230
begin
bath = 2 # means colex 2 and 4
# bath == 2 -> choi is defined differently (add transpose)
# choose rather colexicographical order also for tensor product!
# A ⊗ B -> kron(B,A)
# TrB -> ptr(. ,1)
# factorize??? -> swap arguments in kron?? should work perfectly as it is!
#
system = 1 # means coles 1 and 3
#decompose initial:
sigma_rho_0, rho_c_0 = factorize(rho,Dim,bath)
# time evolution
aim_density = ptr(U * rho * U',Dim,bath)
sigma_c = ptr(U * rho_c_0 * U',Dim,bath)
#not necessary anymore
#=
Lambda = repeat(sigma_c[:],1,n^2) * inv(rho_basis)
u_svd, z_svd, v_svd = svd(Lambda)
#lv_inv = lv^-1
temp_vec = u_svd[:,real(z_svd) .== maximum(real.(z_svd))]
temp_vec_2 = v_svd[:,real(z_svd) .== maximum(real.(z_svd))]
Proj = temp_vec * temp_vec_2' * maximum(z_svd)
=#
Choi_c = kron(I(Dim), sigma_c)#choi(Proj,2)
# compare UDM
sigma = ptr(U * sigma_rho_0 * U',Dim,bath)
sigma_full = ptr(U * rho * U',Dim,bath)
# create Krausoperators
rho_b_0 = ptr(sigma_rho_0,Dim,system)
sigma_0 = ptr(sigma_rho_0,Dim,bath)
bw, bv = eigen(rho_b_0)
Kr = [zeros(ComplexF64,Dim,Dim) for _=1:Dim^2]
Up = [zeros(ComplexF64,Dim,Dim) for _=1:Dim,_=1:Dim]
for k = 1:Dim
for j = 1:Dim
if bath == 1
Up[k,j] = U[(1:Dim).+(k-1)*Dim, (1:Dim).+(j-1)*Dim] #U[k:2:end,j:2:end] #
else
Up[k,j] = U[k:Dim:end,j:Dim:end]
end
end
end
for k = 1:Dim, j=1:Dim, g=1:Dim, h=1:Dim
Kr[(k-1)*Dim + j] += bv[g,j]' * Up[g,h] * bv[h,k] * sqrt(bw[k])
end
#create Choi of UDM
Choi_t = sum([Kr[k][:] * Kr[k][:]' for k = 1:Dim^2])
#alternative
if bath == 2
Choi_t_2 = transpose(choi(U,Dim)) * kron(rho_b_0, I(Dim)) * transpose(choi(U,Dim))'
else
Choi_t_2 = choi(U,Dim) * kron(rho_b_0, I(Dim)) * choi(U,Dim)'
end
sigma_dyn_t = reshape(choi(Choi_t + Choi_c,Dim) * sigma_0[:],Dim,Dim)
end
# ╔═╡ 09519cbd-1821-4453-ac28-db81f0adc701
norm(sigma_full - sigma_dyn_t) #choi superoperator from inhomogeneous and homogeneous part vs full time evolution
# ╔═╡ 1f7e4796-c90f-4801-b98c-12518199e1c3
norm(Choi_t_2 - Choi_t)
# ╔═╡ 586e00d9-6144-4fec-b0e6-77465483a78a
eigvals(Choi_t + Choi_c)
# ╔═╡ 256c69ba-d43e-4571-8269-54885157d67d
round.(rho_c_0,digits=5)
# ╔═╡ 5c15961d-339a-4511-a441-f11cb9bc7399
begin
choi_evals = real(eigvals(Choi_t + Choi_c))
if any(choi_evals .< 0)
#eigvals(sum([Kr[k][:] * Kr[k][:]' for k = 1:4]) + Choi_c)
md"not completely positive 😈 $(string(round.(choi_evals,digits=4))) "
else
md" completely positivie 😇 $(string(round.(choi_evals,digits=4)))"
end
end
# ╔═╡ 41d63a63-5dff-4b99-b343-a23124922e07
round.(choi_evals,digits=5)
# ╔═╡ e5b9cb77-da33-45e5-827e-49f768cd1f07
pos_test(Choi_t + Choi_c,1,Dim)
# ╔═╡ be0fb311-f271-40ca-856d-184d98c8c0e7
real(eigvals(reshape(permutedims(reshape(Choi_t + Choi_c,Dim,Dim,Dim,Dim), [3,2,1,4]), Dim^2,Dim^2)))
# ╔═╡ aa527940-4029-48ce-abed-fc2695bd3900
real.(Choi_t)
# ╔═╡ 80a80124-14c2-4cdb-bbd3-16951a06e2f4
sigma_dyn_t
# ╔═╡ 0b4788d7-0c22-4821-b355-19b31aadfdf1
sigma_full
# ╔═╡ 8bcc43e5-00dd-400e-ad15-db0358c9dce2
sigma
# ╔═╡ 331505b8-6195-4ba3-a8cc-abe217c98867
sigma_0
# ╔═╡ 71efd276-b9ed-4d9c-a902-1e028ace8114
begin
rho_try = zeros(ComplexF64,2,2)
for k = 1:4
rho_try += Kr[k] * sigma_0 * Kr[k]'
end
test = zeros(ComplexF64,2,2)
for k = 1:4
test += Kr[k]' * Kr[k]
end
test
end
# ╔═╡ cc741913-140c-4efb-ab3a-27520fb35563
rho_try
# ╔═╡ 2fe62982-40af-4684-a519-1a520892a41c
aim_density
# ╔═╡ 574eda40-33f0-4a03-8d23-103db03614d3
begin
Pro
# try to find the composition of the Choi operator for the quantum--classical correlated initial state
RHO_b
KK = [[zeros(ComplexF64,Dim,Dim) for _=1:Dim^2] for _=1:length(RHO_b)]
for h = 1:length(RHO_b)
local bw, bv = eigen(RHO_b[h])
local Kr = [zeros(ComplexF64,Dim,Dim) for _=1:Dim^2]
local Up = [zeros(ComplexF64,Dim,Dim) for _=1:Dim,_=1:Dim]
for k = 1:Dim
for j = 1:Dim
if bath == 1
Up[k,j] = U[(1:Dim).+(k-1)*Dim, (1:Dim).+(j-1)*Dim] #U[k:2:end,j:2:end] #
else
Up[k,j] = U[k:Dim:end,j:Dim:end]
end
end
end
for k = 1:Dim, j=1:Dim, g=1:Dim, h=1:Dim
Kr[(k-1)*Dim + j] += bv[g,j]' * Up[g,h] * bv[h,k] * sqrt(bw[k])
end
KK[h] = Kr
end
test_1 = sum([sum([(KK[h][i]*Pro[:,h]*Pro[:,h]')[:] * (KK[g][i]*Pro[:,g]*Pro[:,g]')[:]' for h=1:3, g=1:3]) for i=1:9])
KKK = [sum([KK[h][i] * (Pro[:,h]*Pro[:,h]') for h=1:3]) for i=1:9]
test_2 = sum([KKK[i][:] * KKK[i][:]' for i = 1:9])
end
# ╔═╡ 40777c5b-f031-404a-8c68-4dc7199245eb
norm(reshape(choi(test_2,Dim) * sigma_0[:],Dim,Dim) - sigma_dyn_t)
# ╔═╡ 87abb02f-d6cf-469d-af1f-e0433b19e572
test_2
# ╔═╡ 0a661954-8cd4-4dbd-a2b6-d3a759dc15a1
test_1
# ╔═╡ 9cca546c-2773-47f5-a45e-831ffa0d886c
sum([ kron( Pro[:,i] * Pro[:,i]', I(Dim)) * choi(U,Dim) * kron(RHO_b[i], RHO_b[j]) * choi(U,Dim)' * kron(I(Dim), Pro[:,j] * Pro[:,j]' ) for i = 1:Dim, j=1:Dim])
# ╔═╡ ab40cfbc-51e4-4d5a-bf97-db4ef3e5c002
begin
# relation of partial trace for initial product states - works! :D
abs.((choi((U),2) * kron(rho_b_0,I(2)) * choi(U,2)') - (Choi_t))
end
# ╔═╡ b105e3f9-7387-4999-a3af-cbdda12eb9ed
rho_t = [exp(-im * H_full * i) * sigma_rho_0 * exp(im * H_full * i) for i = time_vector]
# ╔═╡ d8c75acc-e6c7-443d-9ad8-e804aec54349
rho_r = reshape(reduce(hcat, rho_t),4,4,length(time_vector));
# ╔═╡ f449c516-435e-428c-afad-1e4315cf7e99
time_slider = @bind tt Slider(1:length(time_vector), show_value = true, default = 5)
# ╔═╡ fc9b2490-233d-44b9-8ddf-6368c4f74bae
time_slider
# ╔═╡ a141a91d-5e12-4633-8318-b7f91dc19c46
begin
sys_1 = reshape(reduce(hcat,[ptr(rho_r[:,:,i],2,2) for i = 1:length(time_vector)]),2,2,length(time_vector))
plot(time_vector, real(sys_1[1,1,:]))
plot!(time_vector, real(sys_1[2,2,:]))
plot!(time_vector, real(sys_1[1,2,:]))
plot!(time_vector, imag(sys_1[1,2,:]))
end
# ╔═╡ c94aefb0-42e0-46eb-9274-055378adeae2
sys_1[:,:,1]
# ╔═╡ 6ed3ba14-19bf-4e35-933e-5740bbaa227b
sys_1[:,:,tt]
# ╔═╡ 46e923aa-60a2-4458-b26f-47edaa4b53e0
begin
factors = [factorize(rho_r[:,:,i],2,2) for i = 1:length(time_vector)]
c = [sum(abs.(factors[i][2]).^2) for i = 1:length(time_vector)]
end
# ╔═╡ 6ed477ad-a5f8-4682-b171-323c096c9748
begin
plot(time_vector, real.(rho_r[2,2,:]))
plot!(time_vector, real.(rho_r[3,3,:]))
plot!(time_vector, real([tr(ptr(rho_r[:,:,i],2,2)^2) for i = 1:length(time_vector)]), ylim = [0,0.7], label = "purity of sys 1")
plot!(time_vector, c)
end
# ╔═╡ 484a0ec1-a352-473e-9826-e7d38d1de570
sys_1[:]
# ╔═╡ 970f944e-e57a-43cd-9e3e-a129d4ffb53d
size(time_vector)
# ╔═╡ 15493335-e1b7-489e-95c4-5f74a2e23534
size(rho_r[3,3,:])
# ╔═╡ 6df2ffbf-f921-477f-b6bf-b8d2070f9cc2
time_vector
# ╔═╡ 7b5f91e5-f5c4-4971-b62f-7ecf5c5393c6
md" # Partial trace"
# ╔═╡ 733ad294-2c04-4bbb-bf5f-778a5bed34ca
begin
dim_1 = 3
a = rand(dim_1,dim_1)
b = rand(2,2)
G = kron(a,b)
end
# ╔═╡ b4f086ea-5f62-42e9-a7d0-02dff480fdcf
begin
function partial_trace(K, dim_1, which_space)
dim_2 = size(K,1) ÷ dim_1
# check QuantumOptics.jl
#https://qojulia.org
if which_space == 1
L = kron(I(dim_1), trues(dim_2, dim_2))
P = dropdims(sum(reshape(K[L],dim_2, dim_2, dim_1),dims=3),dims=3)
else
L = kron(trues(dim_1,dim_1), I(dim_2))
P = dropdims(sum(reshape(K[L], dim_1, dim_2, dim_1),dims=2),dims=2)
end
return P
end
end
# ╔═╡ b50bde64-a115-4abe-a952-88983c9fa511
@time partial_trace(G,2,2)
# ╔═╡ 992f454a-4e30-4080-86c0-82aeec7d9597
partial_trace(G,3,1)/tr(a) - b
# ╔═╡ a2d4e774-6c12-4c2c-a213-44cedc6c1667
begin
md"""
# Issue of how to sum over a 3rd dimension
There are three versions:
1. dropdims(sum(A,dims=3),dims=3)
2. sum(eachslice(A,dims=3))
3. Using a function (here called f)
"""
end
# ╔═╡ 7ce9332f-3639-465f-bea1-417308914c47
function f(x::AbstractArray{<:Number, 3})
nrows, ncols, n = size(x)
result = zeros(eltype(x), nrows, ncols)
for k in 1:n
for j in 1:ncols
for i in 1:nrows
result[i,j] += x[i,j,k]
end
end
end
result
end
# ╔═╡ b77e9684-fb8a-4d71-9820-f11654508850
begin
x = rand(3,3,100000)
@time dropdims(sum(x, dims=3), dims=3)
@time sum(eachslice(x,dims=3))
@time f(x)
md""" Check result in the console"""
end
# ╔═╡ 4ed1a61c-a916-46f2-abb2-f95837db3149
dim_2 = size(K,2)÷dim_1
# ╔═╡ bf12a345-7f36-41cf-8567-84819a99e17a
md"""
# Parametrizing a qubit
"""
# ╔═╡ 7fa6e288-20bd-468f-8fb4-bb50cb4ee548
p = 1
# ╔═╡ 8a6325b6-8c55-4418-8163-72210a9ef3b4
md"""
# Tensor product of superoperators
``\mathcal{T} \otimes \mathcal{S}``
"""
# ╔═╡ 4ae7f11e-728c-47b6-a810-39ce5bccd1ea
begin
Tsup = rand(4,4) + 1im*rand(4,4)
Ssup = rand(4,4) + 1im * rand(4,4)
X_T = [rand(2,2) for _=1:5]
X_S = [rand(2,2) for _=1:5]
X = sum([kron(X_S[i], X_T[i]) for i = 1:5])
#elementwise
final_sol = sum([kron(reshape(Ssup * X_S[i][:],2,2),reshape(Tsup * X_T[i][:],2,2)) for i = 1:5])
# final version with
CC = sparse(1:16, choi_index(1:16,[2,2,2,2]),ones(16))
# full operator
T_total = kron(Ssup, Tsup)
T_total = T_total[choi_index(1:16,[2,2,2,2]),choi_index(1:16,[2,2,2,2])]
final_vergleich = reshape(T_total*X[:], 4,4)
final_2 = reshape(CC*kron(Ssup, Tsup) * CC * X[:],4,4)
end
# ╔═╡ 7b80e86a-3009-42d3-9ea9-8e9e245fecd5
final_sol
# ╔═╡ 103372c5-49aa-44d7-ae79-5af08892c5b1
CC * kron(I(4), Tsup) * CC
# ╔═╡ 770d537d-dc21-45fa-8ab5-55554201e479
Tsup
# ╔═╡ 146252a2-834d-4ead-9358-f3e520dc54a1
reshape( kron(I(4), Tsup) * I(4)[:],4,4)
# ╔═╡ ff6fdcad-1196-468a-bfe3-52e406c5decb
norm(final_sol - final_vergleich)
# ╔═╡ a4fd9c40-4e99-4f5c-b9b4-6078611407dc
md""" # Examine the tensor decomposition of a Hamiltonian"""
# ╔═╡ 89bce87e-b4db-4726-bb9e-e67dab205089
begin
H = rand(6,6) .+5
H = H+H'
end
# ╔═╡ 403d03b1-e6d4-423a-94b4-05dd6fa37ebc
tr(H)
# ╔═╡ 733acbdd-c57f-4e05-a522-be97f7ef2ae6
@bind trb Slider(-10:0.1:20, show_value = true)
# ╔═╡ 6834d85e-271f-4994-932c-919ca5cd36f8
tr(H)/4
# ╔═╡ e18bdab5-458d-4aaa-bbf4-06b07d2071f7
begin
d_1 = 3
d_2 = size(K,1) ÷ d_1
end
# ╔═╡ dff95946-ce40-469b-83a4-0c3459840d98
#@bind trs Slider(-10:0.1:10, show_value = true)
trs = 1/d_2 * (tr(H) - d_1 * trb)
# ╔═╡ 30d18139-5c12-4ac1-b196-63bcd90ef5d7
H_S = (partial_trace(H,d_1,2) - I(d_1)*trb)./d_2
# ╔═╡ cda2c9cd-935e-4ca0-bd4a-a99caf68b5a7
tr(H_S) - trs
# ╔═╡ da36ec94-fa8a-4927-8b37-756c4a0bdba8
H_B = (partial_trace(H,d_1,1) - I(d_2) * trs)./d_1
# ╔═╡ a9f16b3c-bbc5-4f0d-84cf-a9ee59c986b4
tr(H_B) - trb
# ╔═╡ bd40b6dd-3441-4f76-9983-73b6fea43a62
H_I = H- kron(H_S, I(d_2)) - kron(I(d_1), H_B)
# ╔═╡ 250fcce2-7a9d-417a-8aaa-f407450ffb09
tr(H_I)
# ╔═╡ 7d157865-23f7-46a9-b04f-7e61ce3085e5
rank(H_I)
# ╔═╡ c31f3bc1-f78d-49e7-99a3-e21fcf2abb2a
H_I * H - H * H_I
# ╔═╡ c330dc0f-b941-4ece-ba1c-0a8301aa7d33
partial_trace(H_I, d_1,1)
# ╔═╡ e48108d1-1507-4ed5-bb35-6d5b18f83c76
begin
trb_ = -41:0.1:40
trs_ = -40:0.11:40
z = zeros(length(trb_), length(trs_))
for i = 1:length(trb_)
for j = 1:length(trs_)
trs = trs_[j]
trb = trb_[i]
H_S = (partial_trace(H,d_1,2) - I(d_1)*trb)./d_2
H_B = (partial_trace(H,d_1,1) - I(d_2) * trs)./d_1
z[i,j] = abs(tr(H_B) - trb) + abs(tr(H_S) - trs)
end
end
end
# ╔═╡ 901ac58c-24e1-4242-9a80-c1925cb36184
partial_trace(H,dim_1,1)
# ╔═╡ ec72a849-17df-40c7-aafc-390083fdfd6e
begin
contour(trs_, trb_, z)
plot!(trs_, trs_ .* -d_2 ./d_1 .+ tr(H)/d_1)
end
# ╔═╡ 62fbffea-699c-44e8-a1fe-6b957c641a61
# ╔═╡ 75ff52ac-8832-4a47-aa56-81d9083348bb
minimum(z[:])
# ╔═╡ 19d8a1f4-8bf2-4f18-a8cd-1f1e4d12d143
tr(H)
# ╔═╡ fa193035-ff5c-4334-b65f-695445db4160
partial_trace(H,2,1)
# ╔═╡ 676bb8c8-5e13-4c49-bed4-62fd92e44779
partial_trace(H,2,2)
# ╔═╡ 7b1ad86e-1ca0-4666-8807-32f60b1e6191
begin
scrub_tau_1 = @bind τ₁ Scrubbable(-1:0.1:1)
scrub_tau_2 = @bind τ₂ Scrubbable(-1:0.1:1)
scrub_tau_3 = @bind τ₃ Scrubbable(-1:0.1:1)
end
# ╔═╡ 9ce64a5b-05fe-4bed-b2c4-044053255a77
begin
md"""
``\tau_1``: $(scrub_tau_1)
``\tau_2``: $(scrub_tau_2)
``\tau_3``: $(scrub_tau_3)
"""
end
# ╔═╡ 359dc820-3ee9-4633-99e4-4e4eb76b6548
begin
md"""
``\tau_1``: $(scrub_tau_1)
``\tau_2``: $(scrub_tau_2)
``\tau_3``: $(scrub_tau_3)
To gain a valid density matrix the parameters ``τᵢ`` have to lie within a sqhere of radius 0.5.
Those parameter sets on that sphere yield pure states.
"""
end
# ╔═╡ 6b666e20-e9f2-4c3c-adb4-9a556139dc20
begin
σ₁ = [0 1; 1 0]
σ₂ = [0 -im; im 0]
σ₃ = [1 0; 0 -1]
qubit = Diagonal([0.5, 0.5]) + τ₁ * σ₁ + τ₂ * σ₂ + τ₃ * σ₃
end
# ╔═╡ 69ef6ce0-cd38-41a8-9463-1326c42d29e4
begin
pauli = collect([I(2)[:] σ₁[:] σ₂[:] σ₃[:]])
choi(repeat(σ₃[:], 1,4) * rho_basis^-1,2)
end
# ╔═╡ f2e6d34f-13c4-4e6e-b1c7-09d63d660c36
kron(I(2),σ₂)
# ╔═╡ 92755d15-ea82-4313-8e67-945d777080c0
σ₁
# ╔═╡ a09f4eea-50ff-4307-96a9-811ec1148949
begin
σₚ = (I(2) + σ₃)/2
σₘ = (I(2) - σ₃)/2
#note, that D is the transposition operator (non CP)
D = -p*(im* σ₂[:]) * (im * σ₂[:])'/2 +(1+p)/2*( σₚ[:] * σₚ[:]' + σₘ[:] * σₘ[:]' + σ₁[:] * σ₁[:]'/2 )
eigvals(D)
#p3 = sqrt(0.5^2 - τ₁^2 - τ₂^2 )
D2 = kron(I(2), τ₁ * σ₁ + τ₂ * σ₂ + τ₃ * σ₃)
# get affine part!!
end
# ╔═╡ a9cb62af-091f-472f-892e-c0aca842e5b1
begin
#check whether \sum K' * K = 1
(1+p)/2 * (σ₁' * σ₁/2 + σₘ' * σₘ + σₚ' * σₚ) - p*(im * σ₂)' * (im* σ₂)/2
end
# ╔═╡ 2b41dc34-d3d9-4736-8511-17049f87e604
sum(D, dims = 1)
# ╔═╡ 7dff4447-10fa-46a5-a1ac-0a64962449c3
tr(reshape(choi(D2,2) * density_matrix(2)[:],2,2))
# ╔═╡ e080aaf6-1453-4142-91cc-7ee26dd68f13