## Phase Retrieval example.
# The Phase Retrieval problem seeks to identify
# missing Phase information in imaging problems.
# Classically, it involves a Fourier transform
# and complex variables, but we can study a simple
# variation on the problem.
#
# Given a full-rank matrix A, we are given the
# result of |Ax|.^2 for some hidden vector x,
# (that is, if y = Ax, then we get z_i = |y_i|^2).
# The goal is to the vector x, or a vector that
# has the same norm.

using Random
using LinearAlgebra
Random.seed!(0)
n = 10000
A = randn(n,4)/sqrt(n)
xtrue = [1,0,-1,0]
y = A*xtrue # we don't learn this
z = y.^2 # we only learn this!
# The goal is to get x from z and A!

10000-element Vector{Float64}:
 1.0858550389561124e-6
 0.0001954686523053184
 2.486424171477468e-10
 0.00019093339447459035
 4.8915919759924644e-5
 0.0005273673719465171
 0.000673464506734742
 0.00026015282118625293
 5.750385715444543e-5
 3.1646434531709566e-5
 3.1304632385335174e-6
 1.2278777259350853e-8
 3.618569775306423e-6
 ⋮
 5.6494539435039535e-6
 6.424156950427331e-5
 8.876882082922743e-9
 1.3865009197848494e-5
 0.0007685433406874259
 3.120626844494615e-5
 4.103836122384246e-5
 0.0001510148418411183
 0.00028624929653821166
 8.035396897070258e-5
 2.3340791123459308e-5
 9.980399569477155e-5

## There is a very good idea that results from using
# the largest eigenvector (this is the eigenvector
# associated with the largest eigenvalue)
# of the matrix A'*Diag(z)*A
#
x = vec(eigen(Symmetric(A'*Diagonal(z)*A)).vectors[:,end])

4-element Vector{Float64}:
  0.7126030254459261
 -0.015344188020373961
 -0.7013937814041298
 -0.002871833373255558

##
ntrial = 50
ns = 10.0.^(range(1;length=ntrial,stop=5)) # log-spaced values of n
c = zeros(0)
for n in ns
  n = ceil(Int,n)
  A = randn(n,4)/sqrt(n)
  y = A*xtrue
  z = y.^2
  x = vec(eigen(Symmetric(A'*Diagonal(z)*A)).vectors[:,end])
  @show x
  push!(c, abs(x'*xtrue)/norm(xtrue))
end

x = [0.3159872772456981, 0.012922281759764953, -0.9289432207356257, 0.19248259116682395]
x = [0.5117248044362032, -0.2272992173450204, -0.7958251032156578, 0.23051072732273908]
x = [0.3784461512878178, 0.004827518437893419, -0.7389520970944593, 0.5574091888737308]
x = [0.23263033491196716, 0.4172267272908101, -0.87486882412086, 0.08005951470206277]
x = [0.13941366305366332, 0.08501315056041105, -0.9591049822288574, 0.23120170382057115]
x = [0.22817583455795515, 0.25281291171369624, -0.8104675485057382, -0.47661700873313634]
x = [0.41082208893629235, -0.48601269760660726, -0.7700261090276349, -0.04557038976057226]
x = [0.5605409022429553, -0.5276500014110821, -0.22091648142406184, 0.598811557302254]
x = [0.4365991817689517, 0.17748842337966808, -0.8762002833945071, -0.10075751795407971]
x = [-0.7668250827931897, -0.3340409572224251, 0.38290303656358726, 0.3921494560592027]
x = [0.6372119665203615, 0.03813972072266936, -0.7697344039388089, -0.003901130530648852]
x = [0.7951789685089613, 0.02851723100381065, -0.6012142474381326, 0.07361116935831785]
x = [0.8571910675154513, -0.2238834139038147, -0.458103539142117, -0.07239363353196553]
x = [0.8339507016057736, 0.15030117441889812, -0.5174389630700911, 0.11913313458684506]
x = [0.6615884228842419, 0.09001164463201594, -0.7405436478496649, -0.0761168060660281]
x = [-0.8301804545124758, -0.020035530365501708, 0.5570319298252601, 0.010696710876665237]
x = [0.7442609984346825, -0.06697387099436508, -0.6644888496073845, -0.006680984994347441]
x = [0.7231042510481134, 0.21358369623528678, -0.6341516096818898, -0.17133004043800504]
x = [0.8142423814229958, 0.03980878091496154, -0.5684663208293663, -0.11077295400682265]
x = [0.6279159070648598, 0.005701817846919699, -0.771769345940274, -0.100305431533274]
x = [0.6825743506114914, 0.04763632918022653, -0.724334154308068, -0.08463491556296253]
x = [-0.6956771446214938, 0.05283423031556589, 0.7162086094790162, 0.016943502187585363]
x = [0.6749324970044808, -0.13714831372907055, -0.7222784336457123, -0.06301054529890313]
x = [0.6684949760748954, 0.02059377933040779, -0.7359007344106078, 0.10554843584555565]
x = [-0.7119038345569589, -0.0059324357770442065, 0.7020255030171322, -0.017830582222464017]
x = [0.7379101822618372, -0.05341186979156973, -0.6713562472813283, -0.04377812599920836]
x = [0.6640993218530125, -0.003509462987244416, -0.7466868095688022, -0.037664078375465364]
x = [0.7348279813708875, 0.003881170328870847, -0.6712619265071795, 0.09705771650484543]
x = [-0.7015993671368232, 0.0009432214205285505, 0.7125615986430611, -0.0036614895611536548]
x = [0.6856610415135276, 0.021365882048781537, -0.7249157624855854, -0.0625265746287453]
x = [-0.7216860881160887, -0.0485057236450677, 0.6905016518896594, -0.004884028154807285]
x = [0.7124117922643148, -0.07444499918652134, -0.6977958941741482, -0.0028758324997803374]
x = [-0.7297581304748607, 0.01383604873348146, 0.6835566073845574, 0.0034639951205372448]
x = [0.687521283083389, -0.026675326143255657, -0.7251651775516125, 0.027173103412381217]
x = [0.6878992749656291, 0.016186388370919992, -0.7256091310004937, 0.00489666641461741]
x = [0.7086000199063986, -0.009039008872181953, -0.7055498348280207, 0.0019335670180241386]
x = [-0.7278224132765602, -0.006898296185256392, 0.685705662945635, 0.005889995442467302]
x = [0.712401704372281, 0.010419916087145029, -0.7016593673659977, 0.00702631793717526]
x = [0.6996348687237911, -0.028266466366450138, -0.7136233114833859, -0.021303207563314067]
x = [0.7015870564517942, 0.03940599697929986, -0.7114374461428856, -0.00891795083317859]
x = [0.7047566077155747, -0.013733174397024338, -0.7093070720090352, 0.0036057454410978384]
x = [-0.7222906461923017, -0.014495275793312997, 0.6914188300363551, 0.0051098800319966265]
x = [0.7159453160207139, 0.005502409007516917, -0.6981246508337463, -0.003741638862021182]
x = [0.7041376751355655, -0.009863817160119845, -0.70987326114716, -0.013145062723629182]
x = [0.7064910888425463, 0.009932442684318477, -0.7076400021902404, -0.0041611619282304795]
x = [-0.7108317417085707, 0.013118778424154445, 0.7032328084330526, -0.0031224630628151084]
x = [0.723390129064396, -0.00589501675375631, -0.6903672686719593, 0.008062524175952821]
x = [0.7050470693513484, -9.267621352227612e-5, -0.709123363482435, -0.0072578766559227015]
x = [0.7032408677211578, 0.008393794314768388, -0.7108936934466996, -0.003461618336057757]
x = [0.7025147009029871, -0.010583696607921356, -0.7115864452969182, 0.002410653644520807]

##
using Plots
pyplot()
myc = c
scatter(ns, myc; xscale=:log10)
#gui()

ArgumentError: Package PyPlot not found in current path.
- Run `import Pkg; Pkg.add("PyPlot")` to install the PyPlot package.

Stacktrace:
  [1] macro expansion
    @ ./loading.jl:2296 [inlined]
  [2] macro expansion
    @ ./lock.jl:273 [inlined]
  [3] __require(into::Module, mod::Symbol)
    @ Base ./loading.jl:2271
  [4] #invoke_in_world#3
    @ ./essentials.jl:1089 [inlined]
  [5] invoke_in_world
    @ ./essentials.jl:1086 [inlined]
  [6] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:2260
  [7] top-level scope
    @ ~/.julia/packages/Plots/FFuQi/src/backends.jl:883
  [8] eval
    @ ./boot.jl:430 [inlined]
  [9] _initialize_backend(pkg::Plots.PyPlotBackend)
    @ Plots ~/.julia/packages/Plots/FFuQi/src/backends.jl:882
 [10] backend(pkg::Plots.PyPlotBackend)
    @ Plots ~/.julia/packages/Plots/FFuQi/src/backends.jl:245
 [11] pyplot()
    @ Plots ~/.julia/packages/Plots/FFuQi/src/backends.jl:86
 [12] top-level scope
    @ In[4]:3