text:
# Execution control

export
    CuDim, cudacall

"""
    CuDim3(x)

    CuDim3((x,))
    CuDim3((x, y))
    CuDim3((x, y, x))

A type used to specify dimensions, consisting of 3 integers for respectively the `x`, `y`
and `z` dimension. Unspecified dimensions default to `1`.

Often accepted as argument through the `CuDim` type alias, eg. in the case of
[`cudacall`](@ref) or [`launch`](@ref), allowing to pass dimensions as a plain integer or a
tuple without having to construct an explicit `CuDim3` object.
"""
struct CuDim3
    x::Cuint
    y::Cuint
    z::Cuint
end

CuDim3(dims::Integer)             = CuDim3(dims,    Cuint(1), Cuint(1))
CuDim3(dims::NTuple{1,<:Integer}) = CuDim3(dims[1], Cuint(1), Cuint(1))
CuDim3(dims::NTuple{2,<:Integer}) = CuDim3(dims[1], dims[2],  Cuint(1))
CuDim3(dims::NTuple{3,<:Integer}) = CuDim3(dims[1], dims[2],  dims[3])

# Type alias for conveniently specifying the dimensions
# (e.g. `(len, 2)` instead of `CuDim3((len, 2))`)
const CuDim = Union{Integer,
                    Tuple{Integer},
                    Tuple{Integer, Integer},
                    Tuple{Integer, Integer, Integer}}

# pack arguments in a buffer that CUDA expects
@generated function pack_arguments(f::Function, args...)
    all(isbitstype, args) || throw(ArgumentError("Arguments to kernel should be bitstype."))

    ex = quote
        Base.@_inline_meta
    end

    # If f has N parameters, then kernelParams needs to be an array of N pointers.
    # Each of kernelParams[0] through kernelParams[N-1] must point to a region of memory
    # from which the actual kernel parameter will be copied.

    # put arguments in Ref boxes so that we can get a pointers to them
    arg_refs = Vector{Symbol}(undef, length(args))
    for i in 1:length(args)
        arg_refs[i] = gensym()
        push!(ex.args, :($(arg_refs[i]) = Base.RefValue(args[$i])))
    end

    # generate an array with pointers
    arg_ptrs = [:(Base.unsafe_convert(Ptr{Cvoid}, $(arg_refs[i]))) for i in 1:length(args)]

    append!(ex.args, (quote
        GC.@preserve $(arg_refs...) begin
            kernelParams = [$(arg_ptrs...)]
            f(kernelParams)
        end
    end).args)
    return ex
end

"""
    launch(f::CuFunction; args...; blocks::CuDim=1, threads::CuDim=1,
           cooperative=false, shmem=0, stream=CuDefaultStream())

Low-level call to launch a CUDA function `f` on the GPU, using `blocks` and `threads` as
respectively the grid and block configuration. Dynamic shared memory is allocated according
to `shmem`, and the kernel is launched on stream `stream`.

Arguments to a kernel should either be bitstype, in which case they will be copied to the
internal kernel parameter buffer, or a pointer to device memory.

This is a low-level call, prefer to use [`cudacall`](@ref) instead.
"""
function launch(f::CuFunction, args...; blocks::CuDim=1, threads::CuDim=1,
                cooperative::Bool=false, shmem::Integer=0,
                stream::CuStream=CuDefaultStream())
    blocks = CuDim3(blocks)
    threads = CuDim3(threads)
    (blocks.x>0 && blocks.y>0 && blocks.z>0)    || throw(ArgumentError("Grid dimensions should be non-null"))
    (threads.x>0 && threads.y>0 && threads.z>0) || throw(ArgumentError("Block dimensions should be non-null"))

    pack_arguments(args...) do kernelParams
        if cooperative
            cuLaunchCooperativeKernel(f,
                                      blocks.x, blocks.y, blocks.z,
                                      threads.x, threads.y, threads.z,
                                      shmem, stream, kernelParams)
        else
            cuLaunchKernel(f,
                           blocks.x, blocks.y, blocks.z,
                           threads.x, threads.y, threads.z,
                           shmem, stream, kernelParams, C_NULL)
        end
    end
end

# convert the argument values to match the kernel's signature (specified by the user)
# (this mimics `lower-ccall` in julia-syntax.scm)
@generated function convert_arguments(f::Function, ::Type{tt}, args...) where {tt}
    types = tt.parameters

    ex = quote
        Base.@_inline_meta
    end

    converted_args = Vector{Symbol}(undef, length(args))
    arg_ptrs = Vector{Symbol}(undef, length(args))
    for i in 1:length(args)
        converted_args[i] = gensym()
        arg_ptrs[i] = gensym()
        push!(ex.args, :($(converted_args[i]) = Base.cconvert($(types[i]), args[$i])))
        push!(ex.args, :($(arg_ptrs[i]) = Base.unsafe_convert($(types[i]), $(converted_args[i]))))
    end

    append!(ex.args, (quote
        GC.@preserve $(converted_args...) begin
            f($(arg_ptrs...))
        end
    end).args)

    return ex
end

"""
    cudacall(f::CuFunction, types, values...; blocks::CuDim, threads::CuDim,
             cooperative=false, shmem=0, stream=CuDefaultStream())

`ccall`-like interface for launching a CUDA function `f` on a GPU.

For example:

    vadd = CuFunction(md, "vadd")
    a = rand(Float32, 10)
    b = rand(Float32, 10)
    ad = Mem.upload(a)
    bd = Mem.upload(b)
    c = zeros(Float32, 10)
    cd = Mem.alloc(c)

    cudacall(vadd, (CuPtr{Cfloat},CuPtr{Cfloat},CuPtr{Cfloat}), ad, bd, cd; threads=10)
    Mem.download!(c, cd)

The `blocks` and `threads` arguments control the launch configuration, and should both
consist of either an integer, or a tuple of 1 to 3 integers (omitted dimensions default to
1). The `types` argument can contain both a tuple of types, and a tuple type, the latter
being slightly faster.
"""
cudacall

# FIXME: can we make this infer properly?
cudacall(f::CuFunction, types::Tuple, args...; kwargs...) where {N} =
    cudacall(f, Base.to_tuple_type(types), args...; kwargs...)

function cudacall(f::CuFunction, types::Type, args...; kwargs...)
    convert_arguments(types, args...) do pointers...
        launch(f, pointers...; kwargs...)
    end
end

## attributes

export attributes

struct AttributeDict <: AbstractDict{CUfunction_attribute,Cint}
    f::CuFunction
end

attributes(f::CuFunction) = AttributeDict(f)

@enum_without_prefix CUfunction_attribute CU_

function Base.getindex(dict::AttributeDict, attr::CUfunction_attribute)
    val = Ref{Cint}()
    cuFuncGetAttribute(val, attr, dict.f)
    return val[]
end

Base.setindex!(dict::AttributeDict, val::Integer, attr::CUfunction_attribute) =
    cuFuncSetAttribute(dict.f, attr, val)


text:
using Base.Test
using QuantumOptics

basis = FockBasis(2)

# Test creation
@test basis.Nmin == 0
@test basis.Nmax == 2
@test basis.shape[1] == 3

# Test equality
@test FockBasis(2) == FockBasis(2)
@test FockBasis(2) == FockBasis(0,2)
@test FockBasis(2) != FockBasis(3)
@test FockBasis(1,3) != FockBasis(2,4)

# Test operators
@test number(basis) == SparseOperator(basis, spdiagm(Complex128[0, 1, 2]))
@test destroy(basis) == SparseOperator(basis, sparse(Complex128[0 1 0; 0 0 sqrt(2); 0 0 0]))
@test create(basis) == SparseOperator(basis, sparse(Complex128[0 0 0; 1 0 0; 0 sqrt(2) 0]))
@test number(basis) == dagger(number(basis))
@test create(basis) == dagger(destroy(basis))
@test destroy(basis) == dagger(create(basis))
@test_approx_eq_eps tracedistance(full(create(basis)*destroy(basis)), full(number(basis))) 0. 1e-15

# Test application onto statevectors
@test create(basis)*fockstate(basis, 0) == fockstate(basis, 1)
@test create(basis)*fockstate(basis, 1) == sqrt(2)*fockstate(basis, 2)
@test dagger(fockstate(basis, 0))*destroy(basis) == dagger(fockstate(basis, 1))
@test dagger(fockstate(basis, 1))*destroy(basis) == sqrt(2)*dagger(fockstate(basis, 2))

@test destroy(basis)*fockstate(basis, 1) == fockstate(basis, 0)
@test destroy(basis)*fockstate(basis, 2) == sqrt(2)*fockstate(basis, 1)
@test dagger(fockstate(basis, 1))*create(basis) == dagger(fockstate(basis, 0))
@test dagger(fockstate(basis, 2))*create(basis) == sqrt(2)*dagger(fockstate(basis, 1))

# Test Fock states
b1 = FockBasis(2, 5)
b2 = FockBasis(5)

@test expect(number(b1), fockstate(b1, 3)) == complex(3.)
@test expect(number(b2), fockstate(b2, 3)) == complex(3.)

# Test coherent states
b1 = FockBasis(100)
b2 = FockBasis(2, 5)
alpha = complex(3.)

@test_approx_eq_eps norm(expect(destroy(b1), coherentstate(b1, alpha)) - alpha) 0. 1e-14
@test_approx_eq_eps norm(coherentstate(b1, alpha).data[3:6] - coherentstate(b2, alpha).data) 0. 1e-14

# Test qfunc
b = FockBasis(50)
alpha = complex(1., 2.)
X = [-2:0.1:2;]
Y = [0:0.1:3;]
rho = tensor(coherentstate(b, alpha), dagger(coherentstate(b, alpha)))
Q = qfunc(rho, X, Y)
for (i,x)=enumerate(X), (j,y)=enumerate(Y)
    c = complex(x, y)
    @test_approx_eq_eps Q[i,j] exp(-abs2(c) - abs2(alpha) + 2*real(alpha*conj(c)))/pi 1e-14
end

ext: jl


text:
function contains(string,token)
    return occursin(token,string)
end

function write_combined_program_components_to_disk(file_path::String, set_of_program_components::Set{VLProgramComponent})

    # check - do we have teh file path?
    if (isdir(file_path) == false)
        mkpath(file_path)
    end

    # initialize -
    buffer = ""
    for program_component in set_of_program_components
        component_type = program_component.type
        if (component_type == :buffer)

            # get the data -
            program_buffer = program_component.buffer

            # cache -
            buffer*=program_buffer
            buffer*="\n"
        end
    end

    # all have the same file name, so just grab a component -
    filename = pop!(set_of_program_components).filename

    # build the path -
    path_to_program_file = file_path*"/"*filename

    # write -
    outfile = open(path_to_program_file, "w")
    write(outfile, buffer);
    close(outfile);
end

function write_combined_program_components_to_disk(file_path::String, array_of_program_components::Array{VLProgramComponent,1})

    # check - do we have teh file path?
    if (isdir(file_path) == false)
        mkpath(file_path)
    end

    # initialize -
    buffer = ""
    for program_component in array_of_program_components
        component_type = program_component.type
        if (component_type == :buffer)

            # get the data -
            program_buffer = program_component.buffer

            # cache -
            buffer*=program_buffer
            buffer*="\n"
        end
    end

    # all have the same file name, so just grab a component -
    filename = pop!(array_of_program_components).filename

    # build the path -
    path_to_program_file = file_path*"/"*filename

    # write -
    outfile = open(path_to_program_file, "w")
    write(outfile, buffer);
    close(outfile);
end

function write_program_components_to_disk(file_path::String, set_of_program_components::Set{VLProgramComponent})

  # check - do we have teh file path?
  if (isdir(file_path) == false)
    mkpath(file_path)
  end

  # go through each component, and dump the buffer to disk -
  for program_component in set_of_program_components

    # We switch on type -
    filename = program_component.filename
    component_type = program_component.type
    if (component_type == :buffer)

        # get the data -
        program_buffer = program_component.buffer

        # build the path -
        path_to_program_file = file_path*"/"*filename

        # Write the file -
        outfile = open(path_to_program_file, "w")
        write(outfile,program_buffer);
        close(outfile);

    elseif (component_type == :matrix || component_type == :vector)

        # get the matrix -
        program_matrix = program_component.matrix

        # build the path -
        path_to_program_file = file_path*"/"*filename

        # write the file -
        writedlm(path_to_program_file, program_matrix)

    else
        error("unsupported program component type: $(component_type)")
    end
  end
end

function write_file_to_path(path_to_file::String,buffer::Array{String,1})

    # check, is this a legit dir?
    is_dir_path_ok(path_to_file)

    # write -
    open("$(path_to_file)", "w") do f

        for line_item in buffer
            write(f,"$(line_item)\n")
        end
    end
end

function read_file_from_path(path_to_file::String)::Array{String,1}

    # is this mapping file path legit?
    is_file_path_ok(path_to_file)

    # initialize -
    buffer = String[]

    # Read in the file -
    open("$(path_to_file)", "r") do file
        for line in eachline(file)
            +(buffer,line)
        end
    end

    # return -
    return buffer
end

function move_existing_project_at_path(path_to_existing_project::String)::Bool

    # we are getting called *if* we already know there is a dir conflict -
    # if this is getting called, we have an existing dir where the user wants to write code.
    # we need then create a new dir called *.0, and mv the offending dir to this location?
    # return true if this worked, otherwise false -

    # parent and child dir that we are generating into -
    parent_dir = dirname(path_to_existing_project)
    child_dir = basename(path_to_existing_project)
    destination_path = ""

    # current backup index  -
    current_backup_index = 0

    # do we already have the destination?
    loop_flag = true
    while loop_flag

         # make a destination path -
        destination_path = joinpath(parent_dir,"$(child_dir).$(current_backup_index)")

        # we don't have this dir, we are done -
        if (isdir(destination_path) == false)
            loop_flag = false
        end

        # ok, looks like we already have this dir, update the counter -
        current_backup_index = current_backup_index + 1
    end

    # mv -
    mv(path_to_existing_project, destination_path)

    # check -
    if (isdir(destination_path) == false)
        return false
    end

    return true
end

function extract_section(file_buffer_array::Array{String,1}, start_section_marker::String, end_section_marker::String)::Array{String,1}

    # initialize -
    section_buffer = String[]

    # find the SECTION START AND END -
    section_line_start = 1
    section_line_end = 1
    for (index,line) in enumerate(file_buffer_array)

        if (occursin(start_section_marker,line) == true)
            section_line_start = index
        elseif (occursin(end_section_marker,line) == true)
            section_line_end = index
        end
    end

    for line_index = (section_line_start+1):(section_line_end-1)
        line_item = file_buffer_array[line_index]
        push!(section_buffer,line_item)
    end

    # return -
    return section_buffer
end

function transfer_distribution_file(path_to_distribution_files::String,
                                      input_file_name_with_extension::String,
                                      path_to_output_files::String,
                                      output_file_name_with_extension::String)

    # Load the specific file -
    # create src_buffer -
    src_buffer::Array{String} = String[]

    # check - do we have the file path?
    if (isdir(path_to_output_files) == false)
        mkpath(path_to_output_files)
    end

    # path to distrubtion -
    path_to_src_file = path_to_distribution_files*"/"*input_file_name_with_extension
    open(path_to_src_file,"r") do src_file
        for line in eachline(src_file)

            # need to add a new line for some reason in Julia 0.6
            new_line_with_line_ending = line*"\n"
            push!(src_buffer,new_line_with_line_ending)
        end
    end

    # Write the file to the output -
    path_to_program_file = path_to_output_files*"/"*output_file_name_with_extension
    outfile = open(path_to_program_file, "w")
    write(outfile,src_buffer);
    close(outfile);
end

function transfer_distribution_files(path_to_distribution_files::String,
                                      path_to_output_files::String,
                                      file_extension::String)

    # Search the directory for src files -
    # load the files -
    searchdir(path,key) = filter(x->contains(x,key),readdir(path))

    # build src file list -
    list_of_src_files = searchdir(path_to_distribution_files,file_extension)

    # check - do we have the file path?
    if (isdir(path_to_output_files) == false)
        mkpath(path_to_output_files)
    end

    # go thru the src file list, and copy the files to the output path -
    for src_file in list_of_src_files

        # create src_buffer -
        src_buffer::Array{String,1} = String[]

        # path to distrubtion -
        path_to_src_file = path_to_distribution_files*"/"*src_file
        open(path_to_src_file,"r") do src_file
            for line in eachline(src_file)

                # need to add a new line for some reason in Julia 0.6
                new_line_with_line_ending = line*"\n"
                push!(src_buffer,new_line_with_line_ending)
            end
        end

        # Write the file to the output -
        path_to_program_file = path_to_output_files*"/"*src_file
        open(path_to_program_file, "w") do f
            for line in src_buffer
                write(f,line)
            end
        end
    end
end

function include_function(path_to_src_file::String)::Array{String,1}

    # create src_buffer -
    src_buffer::Array{String,1} = String[]

    # read -
    open(path_to_src_file,"r") do src_file
        for line in eachline(src_file)

            new_line_with_line_ending = line*"\n"
            push!(src_buffer,new_line_with_line_ending)
        end
    end

    # return the raw buffer -
    return src_buffer
end

function include_function(path_to_src_file::String, prefix_pad_string::String)::String

    # create src_buffer -
    src_buffer::Array{String,1} = String[]

    # read -
    open(path_to_src_file,"r") do src_file
        for line in eachline(src_file)

            new_line_with_line_ending = line*"\n"
            push!(src_buffer,new_line_with_line_ending)
        end
    end

    string_value = ""
    for line in src_buffer
        string_value *= prefix_pad_string*line
    end

    return src_buffer
end

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
# source.jl
# this file defines the source term functors for the advection equation,
# similar to the boundary condition functors

@doc """
### AdvectionEquationMod.SRC0

  This is the zero source term.  This is the default of source term
  is specified
"""->
mutable struct SRC0 <: SRCType
end

function (obj::SRC0)(params::ParamTypes, coords::AbstractVector, t)
  return 0
end

@doc """
### AdvectionEquationMod.SRC1

  This source term returns 1 everywhere.
"""->
mutable struct SRC1 <: SRCType
end

function (obj::SRC1)(params::ParamType2, coords::AbstractVector, t)
  return 1
end

@doc """
### AdvectionEquationMod.SRC2

  This source term returns 1 everywhere.
"""->
mutable struct SRC2 <: SRCType
end

function (obj::SRC2)(params::ParamType2, coords::AbstractVector, t)
  return 2
end

@doc """
### AdvectionEquationMod.SRCx

  This source term that returns: f = x
"""->
mutable struct SRCx <: SRCType
end

function (obj::SRCx)(params::ParamType2, coords::AbstractVector, t)
  return coords[1]
end

@doc """
### AdvectionEquationMod.SRC1

  This source term that returns: the derivative of mms1
"""->
mutable struct SRCmms1 <: SRCType
end

function (obj::SRCmms1)(params::ParamType2, coords::AbstractVector, t)
  return params.alpha_x*calc_mms1dx(params, coords, t)
end

@doc """
### AdvectionEquationMod.SRCx4

  This source term that returns: the source term for a manufactured solution
  using a 4th order polynomial
"""->
mutable struct SRCx4 <: SRCType
end

function (obj::SRCx4)(params::ParamType2, coords::AbstractVector, t)
  return params.alpha_x*calc_x4der(params, coords, t)
end

@doc """
### AdvectionEquationMod.SRCp1

  This source term that returns: the source term for a manufactured solution
  using a 1st order polynomial
"""->
mutable struct SRCp1 <: SRCType
end

function (obj::SRCp1)(params::ParamType2, coords::AbstractVector, t)
  return params.alpha_x*calc_p1dx(params, coords, t) + params.alpha_y*calc_p1dy(params, coords, t)

end

function (obj::SRCp1)(params::ParamType3, coords::AbstractVector, t)
  return params.alpha_x*calc_p1dx(params, coords, t) + params.alpha_y*calc_p1dy(params, coords, t) + params.alpha_z*calc_p1dz(params, coords, t)
end

@doc """
### AdvectionEquationMod.SRCp3

  This source term that returns: the source term for a manufactured solution
  using a 2nd order polynomial
"""->
mutable struct SRCp2 <: SRCType
end

function (obj::SRCp2)(params::ParamType2, coords::AbstractVector, t)
  return params.alpha_x*calc_p2dx(params, coords, t) + params.alpha_y*calc_p2dy(params, coords, t)

end

function (obj::SRCp2)(params::ParamType3, coords::AbstractVector, t)
  return params.alpha_x*calc_p2dx(params, coords, t) + params.alpha_y*calc_p2dy(params, coords, t) + params.alpha_z*calc_p2dz(params, coords, t)

end

@doc """
### AdvectionEquationMod.SRCp3

  This source term that returns: the source term for a manufactured solution
  using a 3rd order polynomial
"""->
mutable struct SRCp3 <: SRCType
end

function (obj::SRCp3)(params::ParamType2, coords::AbstractVector, t)
  return params.alpha_x*calc_p3dx(params, coords, t) + params.alpha_y*calc_p3dy(params, coords, t)

end

function (obj::SRCp3)(params::ParamType3, coords::AbstractVector, t)
  return params.alpha_x*calc_p3dx(params, coords, t) + params.alpha_y*calc_p3dy(params, coords, t) + params.alpha_z*calc_p3dz(params, coords, t)

end

@doc """
### AdvectionEquationMod.SRCp4

  This source term that returns: the source term for a manufactured solution
  using a 4th order polynomial
"""->
mutable struct SRCp4 <: SRCType
end

function (obj::SRCp4)(params::ParamType2, coords::AbstractVector, t)
  return params.alpha_x*calc_p4dx(params, coords, t) + params.alpha_y*calc_p4dy(params, coords, t)
end

function (obj::SRCp4)(params::ParamType3, coords::AbstractVector, t)
  return params.alpha_x*calc_p4dx(params, coords, t) + params.alpha_y*calc_p4dy(params, coords, t) + params.alpha_z*calc_p4dz(params, coords, t)
end

@doc """
### AdvectionEquationMod.SRCp5

  This source term that returns: the source term for a manufactured solution
  using a 5th order polynomial
"""->
mutable struct SRCp5 <: SRCType
end

function (obj::SRCp5)(params::ParamType2, coords::AbstractVector, t)
  return params.alpha_x*calc_p5dx(params, coords, t) + params.alpha_y*calc_p5dy(params, coords, t)
end

function (obj::SRCp5)(params::ParamType3, coords::AbstractVector, t)
  return params.alpha_x*calc_p5dx(params, coords, t) + params.alpha_y*calc_p5dy(params, coords, t) + params.alpha_z*calc_p5dz(params, coords, t)
end

@doc """
### AdvectionEquationMod.SRCexp_xplusy

  This is a source term that returns a source term for e^(x+y)
"""->
mutable struct SRCexp_xplusy <: SRCType
end

function (obj::SRCexp_xplusy)(params::ParamType2, coords::AbstractVector, t)
  u = calc_exp_xplusy(params, coords, t)
  return params.alpha_x*u + params.alpha_y*u
end

@doc """
### AdvectionEquationMod.SRCx5plusy5

  This is a source term that returns a source term for e^(x+y)
"""->
mutable struct SRCx5plusy5 <: SRCType
end

function (obj::SRCx5plusy5)(params::ParamType2, coords::AbstractVector, t)
  # u = calc_x5plusy5(params, coords, t)
  x = coords[1]
  y = coords[2]

  return params.alpha_x*5*(x^4) + params.alpha_y*5*(y^4)
end

@doc """
### AdvectionEquationMod.SRCexp5xplus4yplus2

Calculates the source term for q = exp(5*x + 4*y +2)

"""->
mutable struct SRCexp5xplus4yplus2 <: SRCType
end

function (obj::SRCexp5xplus4yplus2)(params::ParamType2, coords::AbstractVector, t)

  u = calc_exp5xplus4yplus2(params, coords, t)
  return params.alpha_x*5*u + params.alpha_y*4*u
end

@doc """
### AdvectionEquationMod.SRCexp5xplusy

Calculates the source term for q = exp(5*x + y)
"""->
mutable struct SRCexp5xplusy <: SRCType
end

function (obj::SRCexp5xplusy)(params::ParamType2, coords::AbstractVector, t)
  u = calc_exp5xplusy(params, coords, t)
  return params.alpha_x*5*u + params.alpha_y*u
end

@doc """
### AdvectionEquationMod.SRCexp3xplusy

Calculates the source term for q = exp(3*x + y)
"""->
mutable struct SRCexp3xplusy <: SRCType
end

function (obj::SRCexp3xplusy)(params::ParamType2, coords::AbstractVector, t)
  u = calc_exp3xplusy(params, coords, t)
  return params.alpha_x*3*u + params.alpha_y*u
end

@doc """
### AdvectionEquationMod.SRCexp2xplus2y

Calculates the source term for q = exp(2*x + 2*y)
"""->
mutable struct SRCexp2xplus2y <: SRCType
end

function (obj::SRCexp2xplus2y)(params::ParamType2, coords::AbstractVector, t)
  u = calc_exp2xplus2y(params, coords, t)
  return params.alpha_x*2*u + params.alpha_y*2*u
end

@doc """
### AdvectionEquationMod.SRCexp_xy

Calculates the source term for q = exp(x*y)
"""->

mutable struct SRCexp_xy <: SRCType
end

function (obj::SRCexp_xy)(params::ParamType2, coords::AbstractVector, t)
  x = coords[1]
  y = coords[2]
  u = calc_exp_xy(params, coords, t)
  return params.alpha_x*y*u + params.alpha_y*x*u
end

@doc """
### AdvectionEquationMod.SRCxplusy

calculates the source term for q = x + y

"""->
mutable struct SRCxplusy <: SRCType
end

function (obj::SRCxplusy)(params::ParamType2, coords::AbstractVector, t)

  return params.alpha_x + params.alpha_y
end

"""
  Source term for unsteady mms
"""
mutable struct SRCunsteadymms <: SRCType
end

function (obj::SRCunsteadymms)(params::ParamType2, coords::AbstractVector, t)

  alpha_x = params.alpha_x
  alpha_y = params.alpha_y

  dudt = calc_unsteadymmsdt(params, coords, t)
  dudx = calc_unsteadymmsdx(params, coords, t)
  dudy = calc_unsteadymmsdy(params, coords, t)
  return  dudt + alpha_x*dudx + alpha_y*dudy
end

function (obj::SRCunsteadymms)(params::ParamType3, coords::AbstractVector, t)

  alpha_x = params.alpha_x
  alpha_y = params.alpha_y
  alpha_z = params.alpha_z

  dudt = calc_unsteadymmsdt(params, coords, t)
  dudx = calc_unsteadymmsdx(params, coords, t)
  dudy = calc_unsteadymmsdy(params, coords, t)
  dudz = calc_unsteadymmsdz(params, coords, t)
  return  dudt + alpha_x*dudx + alpha_y*dudy + alpha_z*dudz
end

"""
  Source term for unsteady poly
"""
mutable struct SRCunsteadypoly <: SRCType
end

function (obj::SRCunsteadypoly)(params::ParamType2, coords::AbstractVector, t)

  alpha_x = params.alpha_x
  alpha_y = params.alpha_y

  dudt = calc_unsteadypolydt(params, coords, t)
  dudx = calc_unsteadypolydx(params, coords, t)
  dudy = calc_unsteadypolydy(params, coords, t)

  return dudt + alpha_x*dudx + alpha_y*dudy
end

@doc """
### AdvectionEquationMod.SRCDict

  It stores all the possible source term dictionary options. Whenever a
  new source is created, it should get added to SRCDict.

  All functors must have the signature:

  src_func(params, coords::ParamType, t)

  where coords is the vector of length 2 containing the x and y coordinates
  of the node, params.alpha_x and params.alpha_y are the advection velocities in the x an y
  directions, and t is the current time
"""->
global const SRCDict = Dict{String, SRCType}(
"SRC0" => SRC0(),
"SRC1" => SRC1(),
"SRC2" => SRC2(),
"SRCx" => SRCx(),
"SRCmms1" => SRCmms1(),
"SRCx4" => SRCx4(),
"SRCp0" => SRC0(),  # this *should* avoid calculating teh source term?
"SRCp1" => SRCp1(),
"SRCp2" => SRCp2(),
"SRCp3" => SRCp3(),
"SRCp4" => SRCp4(),
"SRCp5" => SRCp5(),
"SRCexp_xplusy" => SRCexp_xplusy(),
"SRCx5plusy5" => SRCx5plusy5(),
"SRCexp5xplus4yplus2" => SRCexp5xplus4yplus2(),
"SRCexp5xplusy" => SRCexp5xplusy(),
"SRCexp3xplusy" => SRCexp3xplusy(),
"SRCexp2xplus2y" => SRCexp2xplus2y(),
"SRCexp_xy" => SRCexp_xy(),
"SRCxplusy" => SRCxplusy(),
"SRCunsteadymms" => SRCunsteadymms(),
"SRCunsteadypoly" => SRCunsteadypoly(),
)

@doc """

  This function gets the functor specified by opts["SRCname"] and stores
  it to the equation object.  Currently one 1 source functor is allowed.

"""->
function getSRCFunctors(mesh::AbstractMesh, sbp::AbstractOperator,
                        eqn::AdvectionData, opts)

  # currently we only allow 1 source functor
  eqn.src_func = SRCDict[opts["SRCname"]]
#  println("using source term functor ", eqn.src_func)
  return nothing
end


text:
##################################################
# Compute the NR local diffusion coefficients in (E,L)
##################################################

include("../sources/julia/Main.jl")

r = 1.0 # Radius
E = -0.5 # Energy
L = 0.1 # Angular momentum
m = 0.1 # Mass of field star

println("--------------------")

println("Plummer potential")

println("q = ", qCalc)
println("m = ", m)
println("r = ", r)
println("E = ", E)
println("L = ", L)

println("--------------------")

println("Lc(E) = ", Lc(_tE(E)))
println("Coulomb log = ", logCoulomb)

println("--------------------")

@time dE, dE2, dL, dEdL, dL2 = localOrbitChange(r,E,L,m,100,m)

println("D_E  =  ", dE)
println("D_EE =  ", dE2)
println("D_L  =  ", dL)
println("D_LL =  ", dEdL)
println("D_EL =  ", dL2)

println("--------------------")

max_stars_repo_name: KerwannTEP/CAT

max_issues_repo_name: KerwannTEP/CAT

max_forks_repo_name: KerwannTEP/CAT


text:
@testset "Algebra and Elements" begin
    A = [:a, :b, :c]
    b = StarAlgebras.Basis{UInt8}(words(A, radius = 2))
    l = length(b)

    RG = StarAlgebra(one(first(b)), b, (4, 4))

    a = rand(l)

    @test AlgebraElement(a, RG) isa AlgebraElement
    @test all(RG(g) isa AlgebraElement{typeof(RG)} for g in b)

    @test typeof(zero(RG)) == typeof(RG(0))
    @test typeof(one(RG)) == typeof(RG(1))

    @test_throws AssertionError AlgebraElement([1, 2, 3], RG)
    @test AlgebraElement([1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], RG) isa AlgebraElement

    x = AlgebraElement([1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], RG)
    @test AlgebraElement{Float64}(x) isa AlgebraElement
    y = AlgebraElement{Float64}(x)
    @test parent(x) === parent(y)
    @test x == y

    p = Word(A, [2, 3])
    a = RG(p)
    @test coeffs(a)[b[p]] == 1
    @test coeffs(a) isa SparseVector
    @test all(coeffs(a)[i] == 0 for i in 1:length(b) if i ≠ b[p])
    @test a(p) == 1
    @test all(a(g) == 0 for g in b if g != p)

    @test sprint(show, zero(RG)) == "0·(id)"
    @test sprint(show, one(RG)) == "1·(id)"
    @test isone(one(a))
    @test iszero(zero(a))
    @test sprint(show, a) == "1·b·c"
    @test sprint(show, -a) == "-1·b·c"

    @test hash(a) == hash(one(RG) + RG(p) - one(RG))

    z = zeros(l)
    z[b[p]] = 1
    @test AlgebraElement(z, RG) == a

    @test StarAlgebras.supp(a) == [p]
    @test StarAlgebras.supp_ind(a) == [b[p]]

    s = one(first(b))
    @test a(s) == 0

    a[s] = 2

    @test coeffs(a)[b[s]] == 2
    @test a[b[s]] == 2
    @test a(s) == 2

    dense_a = AlgebraElement(Vector(coeffs(a)), RG)
    @test a == dense_a
    @test hash(a) == hash(dense_a)

    @test StarAlgebras.supp_ind(a) == [b[s], b[p]] == StarAlgebras.supp_ind(dense_a)
    @test supp(a) == [s, p] == StarAlgebras.supp(dense_a)

    aa = a - RG(p)
    dense_aa = dense_a - RG(p)
    @test StarAlgebras.supp_ind(aa) == [b[s]] == StarAlgebras.supp_ind(dense_aa)
    @test supp(aa) == [s] == StarAlgebras.supp(dense_aa)

    @test sprint(show, a) == "2·(id) +1·b·c"
    @test sprint(show, -a) == "-2·(id) -1·b·c"
    z[b[s]] = 2
    @test AlgebraElement(z, RG) == a
    @test sprint(show, AlgebraElement(z, RG)) == "2.0·(id) +1.0·b·c"

    @test LinearAlgebra.norm(a, 1) == 3
end


text:
@testset "EvalDfs" begin
    @testset "triggering" begin
        d = Dag()

        rn1 = node!(d, Int64)
        rn2 = node!(d, Int64)

        an = add!(d, [
            ("left", rn1, latest(0), true),
            ("right", rn2, latest(0), false),
        ])

        res = Dict()
        onfire!(d, rn1.nid, v -> push!(get!(res, "rn1", []), v))
        onfire!(d, rn2.nid, v -> push!(get!(res, "rn2", []), v))
        onfire!(d, an.nid, v -> push!(get!(res, "an", []), v))

        eve = EvalDfs()
        fire(eve, d, rn1.nid, 42)
        @test res["rn1"] == [42]
        # TODO check throw @test res["rn1"] == [42]
        @test res["an"] == [42]

        fire(eve, d, rn2.nid, 24)  # does not trigger `an` display
        @test res["rn1"] == [42]
        @test res["rn2"] == [24]
        @test res["an"] == [42]

        fire(eve, d, rn1.nid, 642)
        @test res["rn1"] == [42, 642]
        @test res["rn2"] == [24]
        @test res["an"] == [42, 666]
    end

    @testset "not eliding" begin
        d = Dag()

        rn1 = node!(d, Tuple{Int64, Int64})
        rn11 = first!(d, rn1)
        rn12 = last!(d, rn1)

        an = add!(d, [
            ("left", rn11, latest(0), true),
            ("right", rn12, latest(0), true),
        ])

        res = Dict()
        onfire!(d, an.nid, v -> push!(get!(res, "an", []), v))

        eve = EvalDfs()
        fire(eve, d, rn1.nid, (1, 2))
        @test res["an"] == [1, 3]
    end
end


text:
# FWI on the 2D Overthrust model with stochastic gradient descent
# Author: pwitte.slim@gmail.com
# Date: December 2018
#

using Statistics, Random, LinearAlgebra
using JUDI, SlimOptim, HDF5, SegyIO, PyPlot

# Load starting model
if ~isfile("$(JUDI.JUDI_DATA)/overthrust_model_2D.h5")
    ftp_data("ftp://slim.gatech.edu/data/SoftwareRelease/WaveformInversion.jl/2DFWI/overthrust_model_2D.h5")
end
n, d, o, m0 = read(h5open("$(JUDI.JUDI_DATA)/overthrust_model_2D.h5", "r"), "n", "d", "o", "m0")
model0 = Model((n[1], n[2]), (d[1], d[2]), (o[1], o[2]), m0)

# Bound constraints
v0 = sqrt.(1f0 ./ model0.m)
vmin = ones(Float32, model0.n) .* 1.3f0
vmax = ones(Float32, model0.n) .* 6.5f0
vmin[:, 1:21] .= v0[:, 1:21]   # keep water column fixed
vmax[:, 1:21] .= v0[:, 1:21]

# Slowness squared [s^2/km^2]
mmin = vec((1f0 ./ vmax).^2)
mmax = vec((1f0 ./ vmin).^2)

# Load data and create data vector
if ~isfile("$(JUDI.JUDI_DATA)/overthrust_2D.segy")
    ftp_data("ftp://slim.gatech.edu/data/SoftwareRelease/WaveformInversion.jl/2DFWI/overthrust_2D.segy")
end
block = segy_read("$(JUDI.JUDI_DATA)/overthrust_2D.segy")
d_obs = judiVector(block)

# Set up wavelet and source vector
src_geometry = Geometry(block; key = "source", segy_depth_key = "SourceDepth")
wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], 0.008f0) # 8 Hz wavelet
q = judiVector(src_geometry, wavelet)

########################################### FWI ####################################################

# Optimization parameters
maxiter = 20
batchsize = 20
fhistory_SGD = zeros(Float32, maxiter)
fTerm = 1e3
gTerm = 1e1

# Projection operator for bound constraints
proj(x) = reshape(median([vec(mmin) vec(x) vec(mmax)], dims=2), model0.n)

# Define misfit function
function fwi_misfit(model::Model, q::judiVector, d::judiVector; misfit = "L2", compute_gradient = true)

    # Set up operators
    M = judiModeling(model, q.geometry, d.geometry)
    J = judiJacobian(M, q)

    # Data residual, function value and gradient
    if misfit == "L2"
        r = M*q - d
        f = .5f0*norm(r)^2
        compute_gradient == true && (g = adjoint(J)*r)    # gradient not necessary for line search
    elseif misifit == "huber"
        r = M*q - d
        f = eps^2*sqrt(1f0 + dot(r, r)/eps^2) - eps^2
        compute_gradient == true && (g = adjoint(J)*r/sqrt(1f0 + dot(r, r)/eps^2))
    else
        throw("Wrong misfit method specified")
    end
    if compute_gradient == true
        return f, g
    else
        return f
    end
end

# Main loop
for j = 1: maxiter

    # select current subset of shots
    i = randperm(d_obs.nsrc)[1:batchsize]
    f, g = fwi_misfit(model0, q[i], d_obs[i])
    println("FWI iteration no: ", j, "; function value: ", f)
    fhistory_SGD[j] = f

    # linesearch
    step = backtracking_linesearch(model0, q[i], d_obs[i], f, g,proj, fwi_misfit; alpha=1f0)

    # Update model and bound projection
    model0.m = proj(model0.m + reshape(step, model0.n))

    if f <= fTerm || norm(g) <= gTerm
        break
    end
end

# Save results
h5open("result_2D_overthrust_sgd.h5", "w") do file
    write(file, "x", sqrt.(1f0 ./ reshape(model0.m, model0.n)), "fhistory", fhistory_SGD)
end

# Plot convergence and final result
figure(); plot(1:maxiter, fhistory_SGD/norm(fhistory_SGD, Inf));
xlabel("Iteration no."); ylabel("Normalized residual"); title("Convergence of FWI w/ SGD")
figure(); imshow(sqrt.(1f0./adjoint(model0.m)), extent=(0, 20.0, 5.15, 0)); title("FWI with SGD");
xlabel("Lateral position [km]"); ylabel("Depth [km]")


text:
@with_kw mutable struct Strategy{S}
    sim::S
    最大持仓::Int = 0
    持仓天数::Int = 1
    最多交易次数::Int = 1000
    禁止平今::Bool = false
    夜盘最早开仓时间::Float32 = -Inf
    夜盘最晚开仓时间::Float32 = Inf
    夜盘最晚平仓时间::Float32 = Inf
    早盘最早开仓时间::Float32 = -Inf
    早盘最晚开仓时间::Float32 = Inf
    早盘最晚平仓时间::Float32 = Inf
end

function simulate end

function transition end

function threshgrid end

select(::Any) = false


text:
function main()
    println("Before")
    exit(0)
    println("After")
end

main()


text:
using sem, Feather, ModelingToolkit, Statistics, LinearAlgebra,
    Optim, SparseArrays, Test, Zygote, LineSearches, ForwardDiff

## Observed Data
three_path_dat = Feather.read("test/comparisons/three_path_dat.feather")
three_path_dat_loadeq = Feather.read("test/comparisons/three_path_loadeq_dat.feather")
three_path_par = Feather.read("test/comparisons/three_path_par.feather")
three_path_par_alleq = Feather.read("test/comparisons/three_path_2_par.feather")
three_path_par_loadeq = Feather.read("test/comparisons/three_path_loadeq_par.feather")
three_path_start = Feather.read("test/comparisons/three_path_start.feather")
testdat = Feather.read("test/comparisons/testdat.feather")

## Model definition
@variables x[1:31]

S =[x[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
    0     x[2]  0     0     0     0     0     0     0     0     0     0     0     0
    0     0     x[3]  0     0     0     0     0     0     0     0     0     0     0
    0     0     0     x[4]  0     0     0     x[15] 0     0     0     0     0     0
    0     0     0     0     x[5]  0     x[16] 0     x[17] 0     0     0     0     0
    0     0     0     0     0     x[6]  0     0     0     x[18] 0     0     0     0
    0     0     0     0     x[16] 0     x[7]  0     0     0     x[19] 0     0     0
    0     0     0     x[15] 0     0     0     x[8]  0     0     0     0     0     0
    0     0     0     0     x[17] 0     0     0     x[9]  0     x[20] 0     0     0
    0     0     0     0     0     x[18] 0     0     0     x[10] 0     0     0     0
    0     0     0     0     0     0     x[19] 0     x[20] 0     x[11] 0     0     0
    0     0     0     0     0     0     0     0     0     0     0     x[12] 0     0
    0     0     0     0     0     0     0     0     0     0     0     0     x[13] 0
    0     0     0     0     0     0     0     0     0     0     0     0     0     x[14]]

F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0]

A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[21] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[22] 0     0
    0  0  0  0  0  0  0  0  0  0  0     0     1     0
    0  0  0  0  0  0  0  0  0  0  0     0     x[23] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[24] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[25] 0
    0  0  0  0  0  0  0  0  0  0  0     0     0     1
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[26]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[27]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[28]
    0  0  0  0  0  0  0  0  0  0  0     0     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[29] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[30] x[31] 0]

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

start_val = vcat(
    vec(var(Matrix(three_path_dat), dims = 1))./2,
    fill(0.05, 3),
    fill(0.0, 6),
    fill(1.0, 8),
    fill(0, 3)
    )

semobserved = SemObsCommon(data = Matrix(three_path_dat))

diff_fin = SemFiniteDiff(
    BFGS(
        ;alphaguess = LineSearches.InitialStatic(;scaled = true),
        linesearch = LineSearches.Static()),
    Optim.Options(;g_tol = 0.001)
    )

# diff_fin = SemFiniteDiff(
#     LBFGS(;alphaguess = LineSearches.InitialStatic(;scaled = true),
#         linesearch = LineSearches.Static()),
#     Optim.Options(;g_tol = 0.001)
#     )

loss = Loss([SemML(semobserved, [0.0], similar(start_val))])

imply = ImplySymbolic(A, S, F, x, start_val)

rowind = [1:25, collect(26:50), collect(51:75)]

rowind_loadeq = [1:40, 41:80, 81:150]

differ_group = [
    fill(1, 31),
    vcat(fill(1, 14), fill(2, 17)),
    vcat(fill(1, 7), fill(2, 7), fill(3, 10), fill(2, 7))
    ]

differ_group_2 = [
    fill(1, 31),
    fill(1, 31),
    fill(1, 31)
    ]

differ_group_loadeq = [
    fill(1, 31),
    vcat(fill(2, 20), fill(1, 11)),
    vcat(fill(3, 20), fill(1, 11))
    ]

obs_list = get_observed(rowind, three_path_dat, SemObsCommon)
obs_list_loadeq = get_observed(rowind_loadeq, three_path_dat_loadeq, SemObsCommon)

sem_mg, start_val_mg = MGSem(three_path_dat, start_val, differ_group_2,
    rowind, obs_list, imply, loss, diff_fin)

sem_mg_loadeq, start_val_loadeq = MGSem(three_path_dat_loadeq, start_val,
    differ_group_loadeq, rowind_loadeq, obs_list_loadeq, imply,
    loss, diff_fin)

solution = sem_fit(sem_mg, start_val_mg)

solution_loadeq = sem_fit(sem_mg_loadeq, start_val_loadeq)

##
three_path_par_alleq = three_path_par_alleq[1:34, :]

par_order_alleq = [collect(21:34); collect(15:20); 2; 3; 5; 6; 7; collect(9:14)]

three_path_par_alleq.est[par_order_alleq]

solution.minimizer

all(
    abs.(solution.minimizer .- three_path_par_alleq.est[par_order_alleq]
        ) .< 0.05*abs.(three_path_par_alleq.est[par_order_alleq]))

##
par_old = copy(par_loadeq)

par_loadeq == par_old

par_loadeq =
    solution_loadeq.minimizer[vcat(sem_mg_loadeq.par_subsets...)]

est_loadeq = three_path_par_loadeq.est[vcat(par_order_alleq, 34 .+ par_order_alleq,
    68 .+ par_order_alleq)]

all(
    abs.(par_loadeq .- est_loadeq
        ) .< 0.05*abs.(est_loadeq))


text:
""" Find the Best Matching Unit (BMU) given an array dist of distances."""
function bestMatchingUnit(dist::Array{Float64,1})
    return argmin(abs.(dist))
end

function neighbourhood(bmu::Int,s::Int,smax::Int,nodes::Array{Float64,1},w::Float64,t::Float64,r::Float64)

    decay = αexp(s,t)
    shrink = αexp(s,r)

    # Gaussian neighbourhood function
    d = decay*exp.(-(((nodes .- nodes[bmu])/(w*shrink)).^2)./2)
    return d
end


text:
using MacroTools: @capture, @q, shortdef
using ZygoteRules: named, typeless, isvararg
using Base: tail

drop(x, n) = n == 0 ? x : :(tail($(drop(x, n-1))))
drop(n) = x -> drop(x, n)

# TODO: move to ZygoteRules
function tangent end

function gradm(ex)
  @capture(shortdef(ex), (name_(args__) = body_) |
                         (name_(args__) where {Ts__} = body_)) || error("Need a function definition")
  kw = length(args) > 1 && isexpr(args[1], :parameters) ? esc(popfirst!(args)) : nothing
  isclosure = isexpr(name, :(::)) && length(name.args) > 1
  f, T = isexpr(name, :(::)) ?
    (length(name.args) == 1 ? (esc(gensym()), esc(name.args[1])) : esc.(name.args)) :
    (esc(gensym()), :(Core.Typeof($(esc(name)))))
  kT = :(Core.kwftype($T))
  Ts == nothing && (Ts = [])
  args = named.(args)
  argnames = Any[typeless(arg) for arg in args]
  !isempty(args) && isvararg(args[end]) && (argnames[end] = :($(argnames[end])...,))
  args = esc.(args)
  argnames = esc.(argnames)
  Ts = esc.(Ts)
  fargs = kw == nothing ? [:($f::$T), args...] : [kw, :($f::$T), args...]
  dropg  = isclosure ? identity : drop(1)
  dropkw = isclosure ?  drop(2) : drop(3)
  adj = @q @inline Zygote.Forward.tangent($(fargs...)) where $(Ts...) = $(esc(body))
  quote
    $adj
    @inline function Zygote.Forward._pushforward(partials, $f::$T, $(args...)) where $(Ts...)
      y, forw = tangent($f, $(argnames...))
      return y, forw($(dropg(:partials))...)
    end
    @inline function Zygote.Forward._pushforward(dargs, ::$kT, kw, $f::$T, $(args...)) where $(Ts...)
      y, forw = tangent($f, $(argnames...))
      return y, forw($(dropkw(:partials))...)
    end
    nothing
  end
end

macro tangent(ex)
  gradm(ex)
end

pushforward(f, x...) = (ẋ...) -> _pushforward((zerolike(f), ẋ...), f, x...)[2]


text:
module DistributedDictVectors

using Rimu
using MPI
using Folds
using KrylovKit
using LinearAlgebra

include("storage.jl")
include("tvec.jl")
include("workingmemory.jl")
include("eigsolve.jl")

export Storage, TVec, WorkingMemory

end # module


text:
module PhotonTest

@info "Running Unit tests"

using Photon
using Test

@testset "Photon" begin

@info "Testing Basics"
include("./basic.jl")
include("./dense.jl")
include("./conv.jl")
include("./recurrent.jl")

@info "Testing Large models"
include("./modelzoo.jl")

@info "Testing training of models"
include("./training.jl")
include("./optimise.jl")
include("./loss.jl")

@info "Testing metrics and meters"
include("./metrics.jl")

@info "Testing complex training scenarios"
include("./complex.jl")

@info "Testing dataset and dataloader"
include("./data.jl")

end

end


text:
using DiffEqBayes, OrdinaryDiffEq, ParameterizedFunctions, RecursiveArrayTools
using Stan, Mamba, Test

ProjDir = dirname(@__FILE__)

#=
Experimental link to DiffEqBayes examples

Possible issues to investigate:

1. Sampling results vary too much. chains often get stuck
2. @ode_def_nohes needs to be kept outside local scope block
3. Minimize string replacements in Stan language model to prevent recompilations

=#

f1 = @ode_def_nohes LotkaVolterraTest1 begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a=>1.5 b=1.0 c=3.0 d=1.0

cd(ProjDir) do

  isdir("tmp") && rm("tmp", recursive=true)
  isfile("LotkaVolterra1-1.pdf") &&
    rm("LotkaVolterra1-1.pdf")

    println("One parameter case")
    u0 = [1.0,1.0]
    tspan = (0.0,10.0)
    prob1 = ODEProblem(f1,u0,tspan)
    sol = solve(prob1,Tsit5())
    t = collect(range(1, stop=10, length=10))
    randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
    data = convert(Array,randomized)
    priors = [Truncated(Normal(1.5,0.1),0,2)]

    bayesian_result = stan_inference(prob1,t,data,priors;num_samples=300,
      num_warmup=500,likelihood=Normal,vars =("StanODEData",InverseGamma(3,2)))

    global sim1 = bayesian_result.chain_results[1:end,["theta.1", "sigma1.1", "sigma1.2"], :]

  ## Plotting
  p = plot(sim1, [:trace, :mean, :density, :autocor], legend=true);
  draw(p, ncol=4, filename="LotkaVolterra1", fmt=:pdf)

  theta1 = bayesian_result.chain_results[:,["theta.1"],:]
  @test mean(theta1.value) ≈ 1.5 atol=3e-1
end


text:
# ---
# title: 438. Find All Anagrams in a String
# id: problem438
# author: Indigo
# date: 2021-06-21
# difficulty: Medium
# categories: Hash Table
# link: <https://leetcode.com/problems/find-all-anagrams-in-a-string/description/>
# hidden: true
# ---
#
# Given a string **s** and a **non-empty** string **p** , find all the start
# indices of **p** 's anagrams in **s**.
#
# Strings consists of lowercase English letters only and the length of both
# strings **s** and **p** will not be larger than 20,100.
#
# The order of output does not matter.
#
# **Example 1:**
#
#
#
#     Input:
#     s: "cbaebabacd" p: "abc"
#
#     Output:
#     [0, 6]
#
#     Explanation:
#     The substring with start index = 0 is "cba", which is an anagram of "abc".
#     The substring with start index = 6 is "bac", which is an anagram of "abc".
#
#
# **Example 2:**
#
#
#
#     Input:
#     s: "abab" p: "ab"
#
#     Output:
#     [0, 1, 2]
#
#     Explanation:
#     The substring with start index = 0 is "ab", which is an anagram of "ab".
#     The substring with start index = 1 is "ba", which is an anagram of "ab".
#     The substring with start index = 2 is "ab", which is an anagram of "ab".
#
#
#
## @lc code=start
using LeetCode

function find_anagrams(s::String, p::String)
    len = length(p)
    cntp = fill(0, 26)
    cnts = cntp[:]
    res = Int[]
    for c in p
        cntp[c - 'a' + 1] += 1
    end
    for i in 1:len-1
        cnts[s[i] - 'a' + 1] += 1
    end
    for i in len:length(s)
        cnts[s[i] - 'a' + 1] += 1
        cnts == cntp && push!(res, i - len + 1)
        cnts[s[i - len + 1] - 'a' + 1] -= 1
    end
    return res
end
## @lc code=end


text:
# SLAMinDB service functions
using Dates

function getcredentials(;nparticles=true, drawdepth=true, multisession=false, drawedges=true)
  cloudGraph, addrdict = standardcloudgraphsetup(nparticles=nparticles, drawdepth=drawdepth, drawedges=drawedges, multisession=multisession)
  return addrdict
end

function slamindbsavejld(fgl::G, session::AbstractString, itercount::Int) where G <: AbstractDFG
  dt = Base.Dates.now()
  filenamejld = "$(session)_$(Dates.format(dt, "dduyy-HH:MM:SS"))_slamindb_$(itercount).jld"
  println("------Save fg to file: $(filenamejld)------")
  savejld(fgl,file=filenamejld)
  nothing
end

"""
    $(SIGNATURES)

Runs SlamInDb for given number of iterations against a specific session.
"""
function runSlamInDbOnSession(
            caesarConfig::CaesarConfig,
            cloudGraph::CloudGraph,
            userId::String,
            robotId::String,
            sessionId::String,
            iterations::Int64,
            isRecursiveSolver::Bool,
            solverStatus::SolverStatus,
            iterationCompleteCallback)::Nothing
    #
    N = caesarConfig.numParticles

    # TODO: Constants to refactor
    drawbayestree = false

    # Update the status parameters
    solverStatus.isAttached = true
    solverStatus.attachedSessionTimestamp = string(unix2datetime(time()))
    solverStatus.userId = userId
    solverStatus.robotId = robotId
    solverStatus.sessionId = sessionId
    solverStatus.result = "IN PROGRESS"

    itercount = 0
    while ((iterations > 0 || iterations == -1) && solverStatus.isAttached)
        iterationStats = IterationStatistics()
        try
          iterations = iterations == -1 ? iterations : iterations-1 # stop at 0 or continue indefinitely if -1

          startns = time_ns()
          solverStatus.iteration = itercount

          println("===================CONVERT===================")
          solverStatus.currentStep = "Prep_Convert"
          fg = initfg(sessionname=sessionId, robotname=robotId, username=userId, cloudgraph=cloudGraph)
          updatenewverts!(fg, N=N)
          println()

          println("=============ENSURE INITIALIZED==============")
          solverStatus.currentStep = "Prep_EnsureInitialized"
          ensureAllInitialized!(fg)
          println()

          println("================MULTI-SESSION================")
          solverStatus.currentStep = "Prep_MultiSession"
          rmInstMultisessionPriors!(cloudGraph, session=sessionId, robot=robotId, user=userId, multisessions=caesarConfig.multiSession)
          println()

          println("====================SOLVE====================")
          solverStatus.currentStep = "Init_Solve"
          fg = initfg(sessionname=sessionId, robotname=robotId, username=userId, cloudgraph=cloudGraph)

          setBackendWorkingSet!(cloudGraph.neo4j.connection, sessionId, robotId, userId)

          println("Get local copy of graph")

          solverStatus.currentStep = "Init_LocalGraphCopy"
          if fullLocalGraphCopy!(fg)
            iterationStats.numNodes = length(fg.IDs) + length(fg.fIDs) #Variables and factors
            # (savejlds && itercount == 0) ? slamindbsavejld(fg, sessionId, itercount) : nothing
            itercount += 1

            println("-------------Ensure Initialization-----------")
            solverStatus.currentStep = "Solve_EnsureInitialized"
            ensureAllInitialized!(fg)

            println("------------Bayes (Junction) Tree------------")
            solverStatus.currentStep = "Solve_BuildBayesTree"
            tree = wipeBuildNewTree!(fg,drawpdf=drawbayestree)

            solverStatus.currentStep = "Solve_InferOverTree"
            if !isRecursiveSolver
              inferOverTree!(fg, tree, N=N)
            else
              inferOverTreeR!(fg, tree, N=N)
            end

          else
            sleep(0.2)
          end

          # Notify iteration update.
          solverStatus.lastIterationDurationSeconds = (time_ns() - startns) / 1e9
          solverStatus.currentStep = "Idle"
          iterationStats.result = "GOOD"
          solverStatus.result = "GOOD"
      catch ex
          io = IOBuffer()
          showerror(io, ex, catch_backtrace())
          err = String(take!(io))
          msg = "ERROR\r\n$err"
          iterationStats.result = msg
          solverStatus.result = msg
      finally
          iterationStats.endTimestamp = Dates.now()
          iterationCompleteCallback(iterationStats)
      end
    end

    solverStatus.isAttached = false
    solverStatus.detachedSessionTimestamp = string(unix2datetime(time()))
    solverStatus.userId = ""
    solverStatus.robotId = ""
    solverStatus.sessionId = ""
    return nothing
end

"""
    $(SIGNATURES)

Low-level call to iterate the SlamInDb solver for given number of iterations against a specific session and keyword parameters.
"""
function runDbSolver(cloudGraph::CloudGraphs.CloudGraph,
            robotId::A,
            sessionId::A,
            userId::A;
            N::Int=100,
            loopctrl::Vector{Bool}=Bool[true],
            iterations::Int=-1,
            multisession::Vector{A}=Sting[""],
            savejlds::Bool=false,
            recursivesolver::Bool=false,
            drawbayestree::Bool=false  ) where {A <: AbstractString}

  itercount = 0
  while loopctrl[1] && (iterations > 0 || iterations == -1) # loopctrl for future use
    iterations = iterations == -1 ? iterations : iterations-1 # stop at 0 or continue indefinitely if -1
    println("===================CONVERT===================")
    fgl = initfg(sessionname=sessionId, robotname=robotId, username=userId, cloudgraph=cloudGraph)
    updatenewverts!(fgl, N=N)
    println()

    println("=============ENSURE INITIALIZED==============")
    ensureAllInitialized!(fgl)
    println()

    println("================MULTI-SESSION================")
    rmInstMultisessionPriors!(cloudGraph, session=sessionId, robot=robotId, user=userId, multisessions=multisession)
    println()

    println("====================SOLVE====================")
    fgl = initfg(sessionname=sessionId, robotname=robotId, username=userId, cloudgraph=cloudGraph)

    setBackendWorkingSet!(cloudGraph.neo4j.connection, sessionId, robotId, userId)

    println("get local copy of graph")

    if fullLocalGraphCopy!(fgl)
      (savejlds && itercount == 0) ? slamindbsavejld(fgl, sessionName, itercount) : nothing
      itercount += 1

      println("-------------Ensure Initialization-----------")
      ensureAllInitialized!(fgl)

      println("------------Bayes (Junction) Tree------------")
      tree = wipeBuildNewTree!(fgl,drawpdf=drawbayestree)
      if !recursivesolver
        inferOverTree!(fgl, tree, N=N)
      else
        inferOverTreeR!(fgl, tree, N=N)
      end
      savejlds ? slamindbsavejld(fgl, sessionName, itercount) : nothing
    else
      sleep(0.2)
    end
  end
end

"""
Manually call the SLAMinDB solver to perform inference on a specified session given the keyword parameter settings.
"""
function slamindb(;addrdict=nothing,
            N::Int=-1,
            loopctrl::Vector{Bool}=Bool[true],
            iterations::Int=-1,
            multisession::Bool=false,
            savejlds::Bool=false,
            recursivesolver::Bool=false,
            drawbayestree::Bool=false  )

  nparticles = false
  if N > 0
    addrdict["num particles"] = string(N)
  else
    nparticles = true
  end

  cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict,
          nparticles=nparticles, multisession=multisession)

  N = parse(Int, addrdict["num particles"])
  session = addrdict["session"]
  robot = addrdict["robot"]
  user = addrdict["user"]

  if !haskey(addrdict, "multisession")
    addrdict["multisession"]=String[]
  end

  runDbSolver(cloudGraph,
              addrdict["robotId"],
              session,
              N=N,
              loopctrl=loopctrl,
              iterations=iterations,
              multisession=addrdict["multisession"],
              savejlds=savejlds,
              recursivesolver=recursivesolver,
              drawbayestree=drawbayestree)
  nothing
end

function convertdb(;addrdict=nothing,
      N::Int=-1  )
  #
  nparticles = false
  if N > 0
    addrdict["num particles"] = string(N)
  else
    nparticles = true
  end
  cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict, nparticles=nparticles)
  N = parse(Int, addrdict["num particles"])
  fg = initfg(sessionname=addrdict["session"], cloudgraph=cloudGraph)
  updatenewverts!(fg, N=N)
end

function resetconvertdb(;addrdict=nothing)
  cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict)
  println("Clearing slamindb data, leaving front-end data, session: $(addrdict["session"])")
  resetentireremotesession(cloudGraph.neo4j.connection, addrdict["session"])
end

#


text:
using MeinBeispiel
using Test

@test meineFunktion(2, 1) == 3

@testset "MeinBeispiel.jl" begin
    @test meineFunktion(2, 1) == 3
    @test meineFunktion(4, 1) == 5
     # Write your tests here.
end


text:
# Implementation of data structures that automatically
#  perform copy-on-write after task copying.

# If current_task is an existing key in s, then return s[current_task].
# Otherwise, return s[current_task] = s[last_task].
immutable TArray{T,N} <: DenseArray{T,N}
  ref :: Symbol  # object_id
  TArray() = new(gensym())
end

call{T}(::Type{TArray{T,1}}, d::Integer)  = TArray(T,  d)
call{T}(::Type{TArray{T}}, d::Integer...) = TArray(T, convert(Tuple{Vararg{Int}}, d))
call{T,N}(::Type{TArray{T,N}}, d::Integer...) = length(d)==N ? TArray(T, convert(Tuple{Vararg{Int}}, d)) : error("malformed dims")
call{T,N}(::Type{TArray{T,N}}, dim::NTuple{N,Int}) = TArray(T, dim)

function TArray(T::Type, dim)
  res = TArray{T,length(dim)}();
  t = current_task()
  d = Array{T}(dim)
  task_local_storage(res.ref, (t,d))
  res
end

# pass through getindex and setindex!
# duplicate TArray if task id does not match current_task
function Base.getindex(S::TArray, i::Real)
  t, d = task_local_storage(S.ref)
  newd = d
#   ct = current_task()
#   if t != ct
#     # println("[getindex]: $(S.ref ) copying data")
#     newd = deepcopy(d)
#     task_local_storage(S.ref, (ct, newd))
#   end
  getindex(newd, i)
end

function Base.setindex!(S::TArray, x, i::Real)
  t, d = task_local_storage(S.ref)
  ct = current_task()
  newd = d
  if t != ct
    # println("[setindex!]: $(S.ref) copying data")
    newd = deepcopy(d)
    task_local_storage(S.ref, (ct, newd))
  end
  setindex!(newd, x, i)
end

function Base.push!(S::TArray, x)
  t, d = task_local_storage(S.ref)
  ct = current_task()
  newd = d
  if t != ct
    newd = deepcopy(d)
    task_local_storage(S.ref, (ct, newd))
  end
  push!(newd, x)
end

function Base.pop!(S::TArray)
  t, d = task_local_storage(S.ref)
  ct = current_task()
  newd = d
  if t != ct
    newd = deepcopy(d)
    task_local_storage(S.ref, (ct, newd))
  end
  pop!(d)
end

function Base.convert(::Type{TArray}, x::Array)
  res = TArray{typeof(x[1]),ndims(x)}();
  t = current_task()
  task_local_storage(res.ref, (t,x))
  res
end

Base.show(io::IO, S::TArray) = Base.show(io::IO, task_local_storage(S.ref)[2])
Base.display(io::IO, S::TArray) = Base.display(io::IO, task_local_storage(S.ref)[2])
Base.size(S::TArray) = Base.size(task_local_storage(S.ref)[2])
Base.ndims(S::TArray) = Base.ndims(task_local_storage(S.ref)[2])

Base.get(t::Task, S) = S
Base.get(t::Task, S::TArray) = (t.storage[S.ref][2])

## convenience constructors ##

"""
     tzeros(dims, ...)
Construct a distributed array of zeros.
Trailing arguments are the same as those accepted by `TArray`.
"""
function tzeros(T::Type, dim)
  res = TArray{T,length(dim)}();
  t = current_task()
  d = zeros(T,dim)
  task_local_storage(res.ref, (t,d))
  res
end
tzeros{T}(::Type{T}, d1::Integer, drest::Integer...) = tzeros(T, convert(Dims, tuple(d1, drest...)))
tzeros(d1::Integer, drest::Integer...) = tzeros(Float64, convert(Dims, tuple(d1, drest...)))
tzeros(d::Dims) = tzeros(Float64, d)


text:
# Copyright 2021 The FirstOrderLp Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
Parameters of the Malitsky and Pock lineseach algorithm
(https://arxiv.org/pdf/1608.08883.pdf).
 """
struct MalitskyPockStepsizeParameters

  """
  Factor by which the step size is multiplied for in the inner loop.
  Valid values: interval (0, 1). Corresponds to mu in the paper.
  """
  downscaling_factor::Float64

  # Y. Malitsky notes that while the theory requires the value to be strictly
  # less than 1, a value of 1 should work fine in practice.
  """
  Breaking factor that defines the stopping criteria of the linesearch.
  Valid values: interval (0, 1]. Corresponds to delta in the paper.
  """
  breaking_factor::Float64

  """
  Interpolation coefficient to pick next step size. The next step size can be
  picked within an interval [a, b] (See Step 2 of Algorithm 1). The solver uses
  a + interpolation_coefficient * (b - a). Valid values: interval [0, 1].
  """
  interpolation_coefficient::Float64
end

"""
Parameters used for the adaptive stepsize policy.

At each inner iteration we update the step size as follows
Our step sizes are a factor
 1 - (iteration + 1)^(-reduction_exponent)
smaller than they could be as a margin to reduce rejected steps.
From the first term when we have to reject a step, the step_size
decreases by a factor of at least 1 - (iteration + 1)^(-reduction_exponent).
From the second term we increase the step_size by a factor of at most
1 + (iteration + 1)^(-growth_exponent)
Therefore if more than order
(iteration + 1)^(reduction_exponent - growth_exponent)
fraction of the iterations have a rejected step we overall decrease the
step_size. When the step_size is below the inverse of the max singular
value we stop having rejected steps.
"""
struct AdaptiveStepsizeParams
  reduction_exponent::Float64
  growth_exponent::Float64
end

"""
Empty placeholder for the parameters of a constant step size policy.
"""
struct ConstantStepsizeParams end

"""
A PdhgParameters struct specifies the parameters for solving the saddle
point formulation of an problem using primal-dual hybrid gradient.
Quadratic Programming Problem (see quadratic_programming.jl):
minimize 1/2 * x' * objective_matrix * x + objective_vector' * x
         + objective_constant
s.t. constraint_matrix[1:num_equalities, :] * x =
     right_hand_side[1:num_equalities]
     constraint_matrix[(num_equalities + 1):end, :] * x >=
     right_hand_side[(num_equalities + 1):end, :]
     variable_lower_bound <= x <= variable_upper_bound
We use notation from Chambolle and Pock, "On the ergodic convergence rates of a
first-order primal-dual algorithm"
(http://www.optimization-online.org/DB_FILE/2014/09/4532.pdf).
That paper doesn't explicitly use the terminology "primal-dual hybrid gradient"
but their Theorem 1 is analyzing PDHG. In this file "Theorem 1" without further
reference refers to that paper.
Our problem is equivalent to the saddle point problem:
    min_x max_y L(x, y)
where
    L(x, y) = y' K x + f(x) + g(x) - h*(y)
    K = -constraint_matrix
    f(x) = objective_constant + objective_vector' x + 1/2*x' objective_matrix x
    g(x) = 0 if variable_lower_bound <= x <= variable_upper_bound
           otherwise infinity
    h*(y) = -right_hand_side' y if y[(num_equalities + 1):end] >= 0
                                otherwise infinity
Note that the places where g(x) and h*(y) are infinite effectively limits the
domain of the min and max. Therefore there's no infinity in the code.

The code uses slightly different notation from the Chambolle and Pock paper.
The primal and dual step sizes are parameterized as:
  tau = primal_step_size = step_size / primal_weight
  sigma = dual_step_size = step_size * primal_weight.
The primal_weight factor is named as such because this parameterization is
equivalent to defining the Bregman divergences as:
D_x(x, x bar) = 0.5 * primal_weight ||x - x bar||_2^2, and
D_y(y, y bar) = 0.5 / primal_weight ||y - y bar||_2^2.

The parameter primal_weight is adjusted smoothly at each restart; to balance the
primal and dual distances traveled since the last restart; see
compute_new_primal_weight().

The adaptive rule adjusts step_size to be as large as possible without
violating the condition assumed in Theorem 1. Adjusting the step size
unfortunately seems to invalidate that Theorem (unlike the case of mirror prox)
but this step size adjustment heuristic seems to work fine in practice.
See comments in the code for details.

TODO: compare the above step size scheme with the scheme by Goldstein
et al (https://arxiv.org/pdf/1305.0546.pdf).

TODO: explore PDHG variants with tuning parameters, e.g. the
overrelaxed and intertial variants in Chambolle and Pock and the algorithm in
"An Algorithmic Framework of Generalized Primal-Dual Hybrid Gradient Methods for
Saddle Point Problems" by Bingsheng He, Feng Ma, Xiaoming Yuan
(http://www.optimization-online.org/DB_FILE/2016/02/5315.pdf).
"""
struct PdhgParameters
  """
  Number of L_infinity Ruiz rescaling iterations to apply to the constraint
  matrix. Zero disables this rescaling pass.
  """
  l_inf_ruiz_iterations::Int

  """
  If true, applies L2 norm rescaling after the Ruiz rescaling.
  """
  l2_norm_rescaling::Bool

  """
  If not `nothing`, runs Pock-Chambolle rescaling with the given alpha exponent
  parameter.
  """
  pock_chambolle_alpha::Union{Float64,Nothing}

  """
  Used to bias the initial value of the primal/dual balancing parameter
  primal_weight. Must be positive. See also
  scale_invariant_initial_primal_weight.
  """
  primal_importance::Float64

  """
  If true, computes the initial primal weight with a scale-invariant formula
  biased by primal_importance; see select_initial_primal_weight() for more
  details. If false, primal_importance itself is used as the initial primal
  weight.
  """
  scale_invariant_initial_primal_weight::Bool

  """
  If >= 4 a line of debugging info is printed during some iterations. If >= 2
  some info is printed about the final solution.
  """
  verbosity::Int64

  """
  Whether to record an IterationStats object.
  """
  record_iteration_stats::Bool

  """
  Check for termination with this frequency (in iterations).
  """
  termination_evaluation_frequency::Int32

  """
  The termination criteria for the algorithm.
  """
  termination_criteria::TerminationCriteria

  """
  Parameters that control when the algorithm restarts and whether it resets
  to the average or the current iterate. Also, controls the primal weight
  updates.
  """
  restart_params::RestartParameters

  """
  Parameters of the step size policy. There are three step size policies
  implemented: Adaptive, Malitsky and Pock, and constant step size.
  """
  step_size_policy_params::Union{
    MalitskyPockStepsizeParameters,
    AdaptiveStepsizeParams,
    ConstantStepsizeParams,
  }
end

"""
A PdhgSolverState struct specifies the state of the solver.  It is used to
pass information among the main solver function and other helper functions.
"""
mutable struct PdhgSolverState
  current_primal_solution::Vector{Float64}

  current_dual_solution::Vector{Float64}

  """
  Current primal delta. That is current_primal_solution - previous_primal_solution.
  """
  delta_primal::Vector{Float64}

  """
  Current dual delta. That is current_dual_solution - previous_dual_solution.
  """
  delta_dual::Vector{Float64}

  """
  A cache of constraint_matrix' * current_dual_solution.
  """
  current_dual_product::Vector{Float64}

  solution_weighted_avg::SolutionWeightedAverage

  step_size::Float64

  primal_weight::Float64

  """
  True only if the solver was unable to take a step in the previous
  iterations because of numerical issues, and must terminate on the next step.
  """
  numerical_error::Bool

  """
  Number of KKT passes so far.
  """
  cumulative_kkt_passes::Float64

  """
  Total number of iterations. This includes inner iterations.
  """
  total_number_iterations::Int64

  """
  Latest required_ratio. This field is only used with the adaptive step size.
  The proof of Theorem 1 requires 1 >= required_ratio.
  """
  required_ratio::Union{Float64,Nothing}

  """
  Ratio between the last two step sizes: step_size(n)/step_size(n-1).
  It is only saved while using Malitsky and Pock linesearch.
  """
  ratio_step_sizes::Union{Float64,Nothing}
end

"""
Defines the primal norm and dual norm using the norms of matrices, step_size
and primal_weight. This is used only for interacting with general utilities
in saddle_point.jl. PDHG utilities these norms implicitly.
"""
function define_norms(
  primal_size::Int64,
  dual_size::Int64,
  step_size::Float64,
  primal_weight::Float64,
)
  # TODO: Should these norms include the step size?
  primal_norm_params = 1 / step_size * primal_weight * ones(primal_size)
  dual_norm_params = 1 / step_size / primal_weight * ones(dual_size)

  return primal_norm_params, dual_norm_params
end

"""
Logging while the algorithm is running.
"""
function pdhg_specific_log(
  problem::QuadraticProgrammingProblem,
  iteration::Int64,
  current_primal_solution::Vector{Float64},
  current_dual_solution::Vector{Float64},
  step_size::Float64,
  required_ratio::Union{Float64,Nothing},
  primal_weight::Float64,
)
  Printf.@printf(
    "   %5d norms=(%9g, %9g) inv_step_size=%9g ",
    iteration,
    norm(current_primal_solution),
    norm(current_dual_solution),
    1 / step_size,
  )
  if !isnothing(required_ratio)
    Printf.@printf(
      "   primal_weight=%18g dual_obj=%18g  inverse_ss=%18g\n",
      primal_weight,
      corrected_dual_obj(
        problem,
        current_primal_solution,
        current_dual_solution,
      ),
      required_ratio
    )
  else
    Printf.@printf(
      "   primal_weight=%18g dual_obj=%18g\n",
      primal_weight,
      corrected_dual_obj(
        problem,
        current_primal_solution,
        current_dual_solution,
      )
    )
  end
end

"""
Logging for when the algorithm terminates.
"""
function pdhg_final_log(
  problem::QuadraticProgrammingProblem,
  avg_primal_solution::Vector{Float64},
  avg_dual_solution::Vector{Float64},
  verbosity::Int64,
  iteration::Int64,
  termination_reason::TerminationReason,
  last_iteration_stats::IterationStats,
)

  if verbosity >= 2
    infeas = max_primal_violation(problem, avg_primal_solution)
    primal_obj_val = primal_obj(problem, avg_primal_solution)
    dual_stats =
      compute_dual_stats(problem, avg_primal_solution, avg_dual_solution)
    println("Avg solution:")
    Printf.@printf(
      "  pr_infeas=%12g pr_obj=%15.10g dual_infeas=%12g dual_obj=%15.10g\n",
      infeas,
      primal_obj_val,
      norm(dual_stats.dual_residual, Inf),
      dual_stats.dual_objective
    )
    Printf.@printf(
      "  primal norms: L1=%15.10g, L2=%15.10g, Linf=%15.10g\n",
      norm(avg_primal_solution, 1),
      norm(avg_primal_solution),
      norm(avg_primal_solution, Inf)
    )
    Printf.@printf(
      "  dual norms:   L1=%15.10g, L2=%15.10g, Linf=%15.10g\n",
      norm(avg_dual_solution, 1),
      norm(avg_dual_solution),
      norm(avg_dual_solution, Inf)
    )
  end

  generic_final_log(
    problem,
    avg_primal_solution,
    avg_dual_solution,
    last_iteration_stats,
    verbosity,
    iteration,
    termination_reason,
  )
end

"""
Estimate the probability that the power method, after k iterations, has relative
error > epsilon.  This is based on Theorem 4.1(a) (on page 13) from
"Estimating the Largest Eigenvalue by the Power and Lanczos Algorithms with a
Random Start"
https://pdfs.semanticscholar.org/2b2e/a941e55e5fa2ee9d8f4ff393c14482051143.pdf
"""
function power_method_failure_probability(
  dimension::Int64,
  epsilon::Float64,
  k::Int64,
)
  if k < 2 || epsilon <= 0.0
    # The theorem requires epsilon > 0 and k >= 2.
    return 1.0
  end
  return min(0.824, 0.354 / (epsilon * (k - 1))) *
         sqrt(dimension) *
         (1.0 - epsilon)^(k - 1 / 2)
end

"""
Estimate the maximum singular value using power method
https://en.wikipedia.org/wiki/Power_iteration, returning a result with
desired_relative_error with probability at least 1 - probability_of_failure.

Note that this will take approximately log(n / delta^2)/(2 * epsilon) iterations
as per the discussion at the bottom of page 15 of

"Estimating the Largest Eigenvalue by the Power and Lanczos Algorithms with a
Random Start"
https://pdfs.semanticscholar.org/2b2e/a941e55e5fa2ee9d8f4ff393c14482051143.pdf

For lighter reading on this topic see
https://courses.cs.washington.edu/courses/cse521/16sp/521-lecture-13.pdf
which does not include the failure probability.

# Output
A tuple containing:
- estimate of the maximum singular value
- the number of power iterations required to compute it
"""
function estimate_maximum_singular_value(
  matrix::SparseMatrixCSC{Float64,Int64};
  probability_of_failure = 0.01::Float64,
  desired_relative_error = 0.1::Float64,
  seed::Int64 = 1,
)
  # Epsilon is the relative error on the eigenvalue of matrix' * matrix.
  epsilon = 1.0 - (1.0 - desired_relative_error)^2
  # Use the power method on matrix' * matrix
  x = randn(Random.MersenneTwister(seed), size(matrix, 2))

  number_of_power_iterations = 0
  while power_method_failure_probability(
    size(matrix, 2),
    epsilon,
    number_of_power_iterations,
  ) > probability_of_failure
    x = x / norm(x, 2)
    x = matrix' * (matrix * x)
    number_of_power_iterations += 1
  end

  # The singular value is the square root of the maximum eigenvalue of
  # matrix' * matrix
  return sqrt(dot(x, matrix' * (matrix * x)) / norm(x, 2)^2),
  number_of_power_iterations
end

function compute_next_primal_solution(
  problem::QuadraticProgrammingProblem,
  current_primal_solution::Vector{Float64},
  current_dual_product::Vector{Float64},
  step_size::Float64,
  primal_weight::Float64,
)
  # The next lines compute the primal portion of the PDHG algorithm:
  # argmin_x [gradient(f)(current_primal_solution)'x + g(x)
  #          + current_dual_solution' K x
  #          + (1 / step_size) * D_x(x, current_primal_solution)]
  # See Sections 2-3 of Chambolle and Pock and the comment above
  # PdhgParameters.
  # This minimization is easy to do in closed form since it can be separated
  # into independent problems for each of the primal variables. The
  # projection onto the primal feasibility set comes from the closed form
  # for the above minimization and the cases where g(x) is infinite - there
  # isn't officially any projection step in the algorithm.
  primal_gradient = compute_primal_gradient_from_dual_product(
    problem,
    current_primal_solution,
    current_dual_product,
  )

  next_primal =
    current_primal_solution .- (step_size / primal_weight) .* primal_gradient
  project_primal!(next_primal, problem)
  return next_primal
end

function compute_next_dual_solution(
  problem::QuadraticProgrammingProblem,
  current_primal_solution::Vector{Float64},
  next_primal::Vector{Float64},
  current_dual_solution::Vector{Float64},
  step_size::Float64,
  primal_weight::Float64;
  extrapolation_coefficient::Float64 = 1.0,
)
  # The next two lines compute the dual portion:
  # argmin_y [H*(y) - y' K (next_primal + extrapolation_coefficient*(next_primal - current_primal_solution)
  #           + 0.5*norm_Y(y-current_dual_solution)^2]
  dual_gradient = compute_dual_gradient(
    problem,
    next_primal .+
    extrapolation_coefficient .* (next_primal - current_primal_solution),
  )
  next_dual =
    current_dual_solution .+ (primal_weight * step_size) .* dual_gradient
  project_dual!(next_dual, problem)
  next_dual_product = problem.constraint_matrix' * next_dual
  return next_dual, next_dual_product
end

"""
Updates the solution fields of the solver state with the arguments given.
The function modifies the first argument: solver_state.
"""
function update_solution_in_solver_state(
  solver_state::PdhgSolverState,
  next_primal::Vector{Float64},
  next_dual::Vector{Float64},
  next_dual_product::Vector{Float64},
)
  solver_state.delta_primal = next_primal - solver_state.current_primal_solution
  solver_state.delta_dual = next_dual - solver_state.current_dual_solution
  solver_state.current_primal_solution = next_primal
  solver_state.current_dual_solution = next_dual
  solver_state.current_dual_product = next_dual_product

  weight = solver_state.step_size
  add_to_solution_weighted_average(
    solver_state.solution_weighted_avg,
    solver_state.current_primal_solution,
    solver_state.current_dual_solution,
    weight,
  )
end

"""
Computes the interaction and movement of the new iterates.
The movement is used to check if there is a numerical error (movement == 0.0)
and based on the theory (Theorem 1) the algorithm only moves if
interaction / movement <= step_size.
"""
function compute_interaction_and_movement(
  solver_state::PdhgSolverState,
  problem::QuadraticProgrammingProblem,
  next_primal::Vector{Float64},
  next_dual::Vector{Float64},
  next_dual_product::Vector{Float64},
)
  delta_primal = next_primal .- solver_state.current_primal_solution
  delta_dual = next_dual .- solver_state.current_dual_solution
  if iszero(problem.objective_matrix)
    primal_objective_interaction = 0.0
  else
    primal_objective_interaction =
      0.5 * (delta_primal' * problem.objective_matrix * delta_primal)
  end
  primal_dual_interaction =
    delta_primal' * (next_dual_product .- solver_state.current_dual_product)
  interaction = abs(primal_dual_interaction) + abs(primal_objective_interaction)
  movement =
    0.5 * solver_state.primal_weight * norm(delta_primal)^2 +
    (0.5 / solver_state.primal_weight) * norm(delta_dual)^2
  return interaction, movement
end

"""
Takes a step using Malitsky and Pock linesearch.
It modifies the third arguement: solver_state.
"""
function take_step(
  step_params::MalitskyPockStepsizeParameters,
  problem::QuadraticProgrammingProblem,
  solver_state::PdhgSolverState,
)
  if !is_linear_programming_problem(problem)
    error(
      "Malitsky and Pock linesearch is only supported for linear" *
      " programming problems.",
    )
  end

  step_size = solver_state.step_size
  ratio_step_sizes = solver_state.ratio_step_sizes
  done = false
  iter = 0

  next_primal = compute_next_primal_solution(
    problem,
    solver_state.current_primal_solution,
    solver_state.current_dual_product,
    step_size,
    solver_state.primal_weight,
  )
  solver_state.cumulative_kkt_passes += 0.5
  step_size =
    step_size +
    step_params.interpolation_coefficient *
    (sqrt(1 + ratio_step_sizes) - 1) *
    step_size

  max_iter = 60
  while !done && iter < max_iter
    iter += 1
    solver_state.total_number_iterations += 1
    ratio_step_sizes = step_size / solver_state.step_size

    # TODO: Get rid of the extra multiply by the constraint matrix (see Remark 1 in
    # https://arxiv.org/pdf/1608.08883.pdf)
    next_dual, next_dual_product = compute_next_dual_solution(
      problem,
      solver_state.current_primal_solution,
      next_primal,
      solver_state.current_dual_solution,
      step_size,
      solver_state.primal_weight;
      extrapolation_coefficient = ratio_step_sizes,
    )
    delta_dual = next_dual .- solver_state.current_dual_solution
    delta_dual_product = next_dual_product .- solver_state.current_dual_product
    # This is the ideal count. The current version of the code incurs in 1.0 kkt
    # pass. See TODO before the next_dual update above.
    solver_state.cumulative_kkt_passes += 0.5

    # The primal weight does not play a role in this condition. As noted in the
    # paper (See second paragraph of Section 2 in https://arxiv.org/pdf/1608.08883.pdf)
    # the coefficient on left-hand-side is equal to
    # sqrt(<primal_step_size> * <dual_step_size>) = step_size.
    # where the equality follows since the primal_weight in the primal and dual step
    # sizes cancel out.
    if step_size * norm(delta_dual_product) <=
       step_params.breaking_factor * norm(delta_dual)
      # Malitsky and Pock guarantee uses a nonsymmetric weighted average, the
      # primal variable average involves the initial point, while the dual
      # doesn't. See Theorem 2 in https://arxiv.org/pdf/1608.08883.pdf for
      # details.
      if solver_state.solution_weighted_avg.sum_primal_solutions_count == 0
        add_to_primal_solution_weighted_average(
          solver_state.solution_weighted_avg,
          solver_state.current_primal_solution,
          step_size * ratio_step_sizes,
        )
      end

      update_solution_in_solver_state(
        solver_state,
        next_primal,
        next_dual,
        next_dual_product,
      )
      done = true
    else
      step_size *= step_params.downscaling_factor
    end
  end
  if iter == max_iter && !done
    solver_state.numerical_error = true
    return
  end
  solver_state.step_size = step_size
  solver_state.ratio_step_sizes = ratio_step_sizes

end

"""
Takes a step using the adaptive step size.
It modifies the third argument: solver_state.
"""
function take_step(
  step_params::AdaptiveStepsizeParams,
  problem::QuadraticProgrammingProblem,
  solver_state::PdhgSolverState,
)
  step_size = solver_state.step_size
  done = false
  iter = 0

  while !done
    iter += 1
    solver_state.total_number_iterations += 1

    next_primal = compute_next_primal_solution(
      problem,
      solver_state.current_primal_solution,
      solver_state.current_dual_product,
      step_size,
      solver_state.primal_weight,
    )

    next_dual, next_dual_product = compute_next_dual_solution(
      problem,
      solver_state.current_primal_solution,
      next_primal,
      solver_state.current_dual_solution,
      step_size,
      solver_state.primal_weight,
    )
    interaction, movement = compute_interaction_and_movement(
      solver_state,
      problem,
      next_primal,
      next_dual,
      next_dual_product,
    )
    solver_state.cumulative_kkt_passes += 1

    if movement == 0.0
      # The algorithm will terminate at the beginning of the next iteration
      solver_state.numerical_error = true
      break
    end
    # The proof of Theorem 1 requires movement / step_size >= interaction.
    if interaction > 0
      step_size_limit = movement / interaction
    else
      step_size_limit = Inf
    end

    if step_size <= step_size_limit
      update_solution_in_solver_state(
        solver_state,
        next_primal,
        next_dual,
        next_dual_product,
      )
      done = true
    end

    first_term =
      (
        1 -
        (
          solver_state.total_number_iterations + 1
        )^(-step_params.reduction_exponent)
      ) * step_size_limit
    second_term =
      (
        1 +
        (
          solver_state.total_number_iterations + 1
        )^(-step_params.growth_exponent)
      ) * step_size
    step_size = min(first_term, second_term)
  end
  solver_state.step_size = step_size
end

"""
Takes a step with constant step size.
It modifies the third argument: solver_state.
"""
function take_step(
  step_params::ConstantStepsizeParams,
  problem::QuadraticProgrammingProblem,
  solver_state::PdhgSolverState,
)
  next_primal = compute_next_primal_solution(
    problem,
    solver_state.current_primal_solution,
    solver_state.current_dual_product,
    solver_state.step_size,
    solver_state.primal_weight,
  )

  next_dual, next_dual_product = compute_next_dual_solution(
    problem,
    solver_state.current_primal_solution,
    next_primal,
    solver_state.current_dual_solution,
    solver_state.step_size,
    solver_state.primal_weight,
  )

  solver_state.cumulative_kkt_passes += 1

  update_solution_in_solver_state(
    solver_state,
    next_primal,
    next_dual,
    next_dual_product,
  )
end

"""
`optimize(params::PdhgParameters,
          original_problem::QuadraticProgrammingProblem)`

Solves a quadratic program using primal-dual hybrid gradient.

# Arguments
- `params::PdhgParameters`: parameters.
- `original_problem::QuadraticProgrammingProblem`: the QP to solve.

# Returns
A SaddlePointOutput struct containing the solution found.
"""
function optimize(
  params::PdhgParameters,
  original_problem::QuadraticProgrammingProblem,
)
  validate(original_problem)
  qp_cache = cached_quadratic_program_info(original_problem)
  scaled_problem = rescale_problem(
    params.l_inf_ruiz_iterations,
    params.l2_norm_rescaling,
    params.pock_chambolle_alpha,
    params.verbosity,
    original_problem,
  )
  problem = scaled_problem.scaled_qp

  primal_size = length(problem.variable_lower_bound)
  dual_size = length(problem.right_hand_side)
  if params.primal_importance <= 0 || !isfinite(params.primal_importance)
    error("primal_importance must be positive and finite")
  end

  # TODO: Correctly account for the number of kkt passes in
  # initialization
  solver_state = PdhgSolverState(
    zeros(primal_size),  # current_primal_solution
    zeros(dual_size),    # current_dual_solution
    zeros(primal_size),  # delta_primal
    zeros(dual_size),    # delta_dual
    zeros(primal_size),  # current_dual_product
    initialize_solution_weighted_average(primal_size, dual_size),
    0.0,                 # step_size
    1.0,                 # primal_weight
    false,               # numerical_error
    0.0,                 # cumulative_kkt_passes
    0,                   # total_number_iterations
    nothing,             # required_ratio
    nothing,             # ratio_step_sizes
  )

  if params.step_size_policy_params isa AdaptiveStepsizeParams
    solver_state.cumulative_kkt_passes += 0.5
    solver_state.step_size = 1.0 / norm(problem.constraint_matrix, Inf)
  elseif params.step_size_policy_params isa MalitskyPockStepsizeParameters
    solver_state.cumulative_kkt_passes += 0.5
    solver_state.step_size = 1.0 / norm(problem.constraint_matrix, Inf)
    solver_state.ratio_step_sizes = 1.0
  else
    desired_relative_error = 0.2
    maximum_singular_value, number_of_power_iterations =
      estimate_maximum_singular_value(
        problem.constraint_matrix,
        probability_of_failure = 0.001,
        desired_relative_error = desired_relative_error,
      )
    solver_state.step_size =
      (1 - desired_relative_error) / maximum_singular_value
    solver_state.cumulative_kkt_passes += number_of_power_iterations
  end

  # Idealized number of KKT passes each time the termination criteria and
  # restart scheme is run. One of these comes from evaluating the gradient at
  # the average solution and evaluating the gradient at the current solution.
  # In practice this number is four.
  KKT_PASSES_PER_TERMINATION_EVALUATION = 2.0

  if params.scale_invariant_initial_primal_weight
    solver_state.primal_weight = select_initial_primal_weight(
      problem,
      ones(primal_size),
      ones(dual_size),
      params.primal_importance,
      params.verbosity,
    )
  else
    solver_state.primal_weight = params.primal_importance
  end

  primal_weight_update_smoothing =
    params.restart_params.primal_weight_update_smoothing

  iteration_stats = IterationStats[]
  start_time = time()
  # Basic algorithm refers to the primal and dual steps, and excludes restart
  # schemes and termination evaluation.
  time_spent_doing_basic_algorithm = 0.0

  # This variable is used in the adaptive restart scheme.
  last_restart_info = create_last_restart_info(
    problem,
    solver_state.current_primal_solution,
    solver_state.current_dual_solution,
  )

  # For termination criteria:
  termination_criteria = params.termination_criteria
  iteration_limit = termination_criteria.iteration_limit
  termination_evaluation_frequency = params.termination_evaluation_frequency

  # This flag represents whether a numerical error occurred during the algorithm
  # if it is set to true it will trigger the algorithm to terminate.
  solver_state.numerical_error = false
  display_iteration_stats_heading(params.verbosity)

  iteration = 0
  while true
    iteration += 1

    # Evaluate the iteration stats at frequency
    # termination_evaluation_frequency, when the iteration_limit is reached,
    # or if a numerical error occurs at the previous iteration.
    if mod(iteration - 1, termination_evaluation_frequency) == 0 ||
       iteration == iteration_limit + 1 ||
       iteration <= 10 ||
       solver_state.numerical_error
      # TODO: Experiment with evaluating every power of two iterations.
      # This ensures that we do sufficient primal weight updates in the initial
      # stages of the algorithm.
      solver_state.cumulative_kkt_passes +=
        KKT_PASSES_PER_TERMINATION_EVALUATION
      # Compute the average solution since the last restart point.
      if solver_state.numerical_error ||
         solver_state.solution_weighted_avg.sum_primal_solutions_count == 0 ||
         solver_state.solution_weighted_avg.sum_dual_solutions_count == 0
        avg_primal_solution = solver_state.current_primal_solution
        avg_dual_solution = solver_state.current_dual_solution
      else
        avg_primal_solution, avg_dual_solution =
          compute_average(solver_state.solution_weighted_avg)
      end

      current_iteration_stats = evaluate_unscaled_iteration_stats(
        scaled_problem,
        qp_cache,
        params.termination_criteria,
        params.record_iteration_stats,
        avg_primal_solution,
        avg_dual_solution,
        iteration,
        time() - start_time,
        solver_state.cumulative_kkt_passes,
        termination_criteria.eps_optimal_absolute,
        termination_criteria.eps_optimal_relative,
        solver_state.step_size,
        solver_state.primal_weight,
        POINT_TYPE_AVERAGE_ITERATE,
      )
      method_specific_stats = current_iteration_stats.method_specific_stats
      method_specific_stats["time_spent_doing_basic_algorithm"] =
        time_spent_doing_basic_algorithm

      primal_norm_params, dual_norm_params = define_norms(
        primal_size,
        dual_size,
        solver_state.step_size,
        solver_state.primal_weight,
      )
      update_objective_bound_estimates(
        current_iteration_stats.method_specific_stats,
        problem,
        avg_primal_solution,
        avg_dual_solution,
        primal_norm_params,
        dual_norm_params,
      )
      if params.record_iteration_stats
        push!(iteration_stats, current_iteration_stats)
      end

      # Check the termination criteria.
      termination_reason = check_termination_criteria(
        termination_criteria,
        qp_cache,
        current_iteration_stats,
      )
      if solver_state.numerical_error && termination_reason == false
        termination_reason = TERMINATION_REASON_NUMERICAL_ERROR
      end

      # Print table.
      if print_to_screen_this_iteration(
        termination_reason,
        iteration,
        params.verbosity,
        termination_evaluation_frequency,
      )
        display_iteration_stats(current_iteration_stats, params.verbosity)
      end

      if termination_reason != false
        # ** Terminate the algorithm **
        # This is the only place the algorithm can terminate. Please keep it
        # this way.
        pdhg_final_log(
          problem,
          avg_primal_solution,
          avg_dual_solution,
          params.verbosity,
          iteration,
          termination_reason,
          current_iteration_stats,
        )
        return unscaled_saddle_point_output(
          scaled_problem,
          avg_primal_solution,
          avg_dual_solution,
          termination_reason,
          iteration - 1,
          iteration_stats,
        )
      end

      current_iteration_stats.restart_used = run_restart_scheme(
        problem,
        solver_state.solution_weighted_avg,
        solver_state.current_primal_solution,
        solver_state.current_dual_solution,
        last_restart_info,
        iteration - 1,
        primal_norm_params,
        dual_norm_params,
        solver_state.primal_weight,
        params.verbosity,
        params.restart_params,
      )

      if current_iteration_stats.restart_used != RESTART_CHOICE_NO_RESTART
        solver_state.primal_weight = compute_new_primal_weight(
          last_restart_info,
          solver_state.primal_weight,
          primal_weight_update_smoothing,
          params.verbosity,
        )
        solver_state.ratio_step_sizes = 1.0
      end
      if current_iteration_stats.restart_used ==
         RESTART_CHOICE_RESTART_TO_AVERAGE
        solver_state.current_dual_product =
          problem.constraint_matrix' * solver_state.current_dual_solution
      end
    end

    time_spent_doing_basic_algorithm_checkpoint = time()

    if params.verbosity >= 6 && print_to_screen_this_iteration(
      false, # termination_reason
      iteration,
      params.verbosity,
      termination_evaluation_frequency,
    )
      pdhg_specific_log(
        problem,
        iteration,
        solver_state.current_primal_solution,
        solver_state.current_dual_solution,
        solver_state.step_size,
        solver_state.required_ratio,
        solver_state.primal_weight,
      )
    end

    take_step(params.step_size_policy_params, problem, solver_state)

    time_spent_doing_basic_algorithm +=
      time() - time_spent_doing_basic_algorithm_checkpoint
  end
end


text:
"""
    TypeF(var_1::TypeB)

Lorem ipsum dolor sit amet.
"""
struct TypeF

  var_1::TypeB

end


text:
using CUDex

SIZE = 100000

x32 = DexArray(rand(Float32,SIZE))
x64 = DexArray(rand(Float64,SIZE))
s32 = rand(Float32)
s64 = rand(Float64)

# sqrt
@time sqrt(x32)
@time sqrt(x64)

# exp
@time exp(x32)
@time exp(x64)

# log
@time log(x32)
@time log(x64)

# tanh
@time tanh(x32)
@time tanh(x64)

# .*
@time x32.*s32
@time x64.*s64


text:
# TODO: implement RSA assymetric cryptography. (A sepperate package) https://en.wikipedia.org/wiki/RSA_(cryptosystem).
# Also other systems like ElGammal, XTR could be implemented

struct RSASignature{H,P} <: AbstractSignature
    hash::H # hash could have its own type
    secret
    pubkey::P # in the same way as public key would have its own type
end

verify(s::RSASignature) = decrypt(s.secret,s.pubkey)==s.hash.hash

### For the system hash function which noone is supposed to use
verify(s::RSASignature{UInt64,P}) where P <: Any = decrypt(s.secret,s.pubkey)==s.hash

function RSASignature(data,hash,s::Signer{K,T}) # where K<:AbstractEncryptionKey ...
    h = hash(data)
    return RSASignature(h,encrypt(h,s.privkey),s.pubkey)
end

# ### An exception when one uses Paillier for the job. But with high probably this is shit.
# import Paillier

# encrypt(data::Integer,priv::Paillier.PublicKey) = Paillier.encrypt(priv,data)
# #encrypt(data::Hash,priv::Paillier.PublicKey) = encrypt(data.hash,priv) # data.hash must be integer

# decrypt(data,pub::Paillier.PrivateKey) = Paillier.decrypt(pub,data)


text:
module FastRecurFlux

using LinearAlgebra
using Flux
using LoopVectorization
using PaddedMatrices
using Requires

using Flux: LSTMCell, GRUCell, gate

SubMatrix{T} = SubArray{T, 2, <: AbstractArray{T}, Tuple{Vararg{AbstractUnitRange, N} where N}, true}
MatTypesC{T} = Union{Matrix{T}, SubMatrix{T}}
MatTypesR{T} = Union{Adjoint{T, <: MatTypesC{T}}, Transpose{T, <: MatTypesC{T}}}
MatTypes{T} = Union{MatTypesC{T}, MatTypesR{T}}
SubVector{T} = SubArray{T, 1, <: AbstractArray{T}, Tuple{Vararg{AbstractUnitRange, N} where N}, true}
VecTypes{T} = Union{Vector{T}, SubVector{T}}

function mul(A::AbstractArray{T}, B::AbstractArray{T}) where T
    M, N, K = size(A, 1), size(A, 2), size(B, 2)
    C = Matrix{T}(undef, M, K)
    if M * N * K < 256^3
        PaddedMatrices.jmul!(C, A, B)
    else
        LinearAlgebra.mul!(C, A, B)
    end
    return C
end

@static if Sys.ARCH == :x86_64
    Base.:*(A::TA, B::TB) where {T <: AbstractFloat, TA <: Matrix{T}, TB <: Matrix{T}} = mul(A, B)
    Base.:*(A::TA, B::TB) where {T <: AbstractFloat, TA <: MatTypes{T}, TB <: MatTypes{T}} = mul(A, B)
    Base.:*(A::TA, B::TB) where {T <: AbstractFloat, TA <: MatTypes{T}, TB <: VecTypes{T}} = vec(mul(A, reshape(B, :, 1)))
end

function (m::LSTMCell)((h, c), x)
    b, o = m.b, size(h, 1)
    g = m.Wi * x .+ m.Wh * h .+ b
    input = σ.(gate(g, o, 1))
    forget = σ.(gate(g, o, 2))
    cell = tanh.(gate(g, o, 3))
    output = σ.(gate(g, o, 4))
    c = forget .* c .+ input .* cell
    h′ = output .* (c_ = tanh.(c))
    return (h′, c), h′
end

function (m::GRUCell)(h, x)
    b, o = m.b, size(h, 1)
    gx, gh = m.Wi * x, m.Wh * h
    r = σ.((r_ = gate(gx, o, 1) .+ gate(gh, o, 1) .+ gate(b, o, 1);))
    z = σ.((z_ = gate(gx, o, 2) .+ gate(gh, o, 2) .+ gate(b, o, 2);))
    h̃ = tanh.((h̃_ = gate(gx, o, 3) .+ r .* gate(gh, o, 3) .+ gate(b, o, 3);))
    h′ = (1 .- z) .* h̃ .+ z .* h
    return h′, h′
end

function Flux.NNlib.σ(x::LoopVectorization.SLEEFPirates.FloatType)
    t = exp(-abs(x))
    LoopVectorization.ifelse(x ≥ 0, inv(one(t) + t), t / (one(t) + t))
end

function __init__()
    for f in (:σ, :tanh)
        @eval Base.broadcasted(::typeof($f), x::AbstractArray{T, N}) where {T, N} = vmap($f, x)
    end
    @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" include("tracker.jl")
    @require ReverseDiff="37e2e3b7-166d-5795-8a7a-e32c996b4267" include("reversediff.jl")
end

end