Script to reproduce https://github.com/JuliaNLSolvers/Optim.jl/issues/953#issuecomment-951543625

simple_state_to_state.jl

# Put this script in an empty folder, an inside the folder, run
#
#     julia --project=.
#     julia> include("simple_state_to_state.jl")
#
using Pkg
Pkg.add(url="https://github.com/JuliaQuantumControl/QuantumPropagators.jl.git", rev="da1dde1b1b3c7caa94be4aff5401669b2c7c8f2a")
Pkg.add(url="https://github.com/JuliaQuantumControl/QuantumControlBase.jl.git", rev="12ce4061fe8b67fe9e419e016b2631fbfbc1ba19")
Pkg.add(url="https://github.com/JuliaQuantumControl/GRAPE.jl.git", rev="d92c6e162a094c4590cbd668328ec349bfff0a0d")
Pkg.add(url="https://github.com/JuliaQuantumControl/QuantumControl.jl.git", rev="5b59c7b5c96280576d77d1d38697affa921c973d")
Pkg.add(url="https://github.com/JuliaQuantumControl/GRAPELinesearchAnalysis.jl.git", rev="9ce2830e503af1c38034e08f3b1204addebbcc77")
Pkg.add("Optim")
Pkg.add("LBFGSB")
Pkg.add("LineSearches")
Pkg.add("UnicodePlots")
Pkg.add("PyPlot")
Pkg.instantiate()

using Printf
using QuantumControl
using LinearAlgebra
using Optim
using LBFGSB
using GRAPE
using QuantumControlBase: chain_infohooks
using GRAPELinesearchAnalysis
using LineSearches
using UnicodePlots
using PyPlot: matplotlib
matplotlib.use("Agg")

using Test

ϵ(t) = 0.2 * QuantumControl.shapes.flattop(t, T = 5, t_rise = 0.3, func = :blackman);

"""Two-level-system Hamiltonian."""
function hamiltonian(Ω = 1.0, ϵ = ϵ)
    σ̂_z = ComplexF64[1 0; 0 -1]
    σ̂_x = ComplexF64[0 1; 1 0]
    Ĥ₀ = -0.5 * Ω * σ̂_z
    Ĥ₁ = σ̂_x
    return (Ĥ₀, (Ĥ₁, ϵ))
end;

H = hamiltonian();
@test length(H) == 2

tlist = collect(range(0, 5, length = 500));

function ket(label)
    result = Dict("0" => Vector{ComplexF64}([1, 0]), "1" => Vector{ComplexF64}([0, 1]))
    return result[string(label)]
end;

objectives = [Objective(initial_state = ket(0), generator = H, target_state = ket(1))]

problem = ControlProblem(
    objectives = objectives,
    tlist = tlist,
    pulse_options=Dict(),
    iter_stop = 500,
    J_T = QuantumControl.functionals.J_T_sm,
    gradient=QuantumControl.functionals.grad_J_T_sm!,
    check_convergence = res -> begin
        ((res.J_T < 1e-3) && (res.converged = true) && (res.message = "J_T < 10⁻³"))
    end,
);

guess_dynamics = propagate(
    objectives[1],
    problem.tlist;
    storage = true,
    observables = (Ψ -> abs.(Ψ) .^ 2,),
)


opt_result = optimize_grape(
        problem,
        info_hook=chain_infohooks(
            GRAPELinesearchAnalysis.plot_linesearch(@__DIR__),
            GRAPE.print_table,
        ),
        #=optimizer=LBFGSB.L_BFGS_B(length(tlist)-1, 10),=#
        optimizer=Optim.LBFGS(),
);

display(opt_result)
display(opt_result.optim_res)
@test opt_result.J_T < 1e-3

println(lineplot(tlist, opt_result.optimized_controls[1]))

opt_dynamics = propagate(
    objectives[1],
    problem.tlist;
    controls_map = IdDict(ϵ => opt_result.optimized_controls[1]),
    storage = true,
    observables = (Ψ -> abs.(Ψ) .^ 2,),
)

@test opt_dynamics[2,end] > 0.99