text:
import SparseArrays.spzeros
import SparseArrays.rowvals
import SparseArrays.nzrange
import SparseArrays.nonzeros
# using LinearAlgebraicRepresentation
# Lar = LinearAlgebraicRepresentation

# @everywhere begin

    """
    Filter sparse matrix. Set new values in matrix.
    """
    function sparse_filter!(sparse, what_to_find, true_value, false_value)
        A = sparse
        vals = nonzeros(A)
        for j=1:length(vals)
        # rows = rowvals(A)
        # # m, n = size(A)
        # n = size(A, 2)
        # for i = 1:n
        #    for j in nzrange(A, i)
        #       # row = rows[j]
              # val = vals[j]
              if vals[j] == what_to_find
                  newval = true_value
              else
                  newval = false_value
              end
              vals[j] = newval
        #    end
        end
        return A
    end

    function sparse_filter_old!(sparse, what_to_find, true_value, false_value)
        A = sparse
        rows = rowvals(A)
        vals = nonzeros(A)
        m, n = size(A)
        for i = 1:n
           for j in nzrange(A, i)
              row = rows[j]
              val = vals[j]
              if val == what_to_find
                  newval = true_value
              else
                  newval = false_value
              end
              vals[j] = newval
           end
        end
        return A
    end

	function unsafe_resize(sp::SparseMatrixCSC,m,n)
		newcolptr = sp.colptr
		resize!(newcolptr,n+1)
		for i=sp.n+2:n+1
	    	newcolptr[i] = sp.colptr[sp.n+1]
	 	end
	  	return SparseMatrixCSC(m,n,newcolptr,sp.rowval,sp.nzval)
	end

	"""
	From indexed base matrix to sparse array with 0 and 1
	Original implementation by Alberto Paoluzzi
	"""
	function characteristicMatrix_push( FV::LarSurf.Lar.Cells)::LarSurf.Lar.ChainOp
		I,J,V = Int64[],Int64[],Int8[]
		for f=1:length(FV)
			for k in FV[f]
			push!(I,f)
			push!(J,k)
			push!(V,1)
			end
		end
		M_2 = sparse(I,J,V)
		return M_2
	end

	"""
	From indexed base matrix to sparse array with 0 and 1
	"""
	function characteristicMatrix_set( FV::LarSurf.Lar.Cells)::LarSurf.Lar.ChainOp
	    sz = length(FV) * length(FV[1])
	    I = Array{Int64,1}(undef, sz)
	    J = Array{Int64,1}(undef, sz)
	    V = Array{Int8,1}(undef, sz)
	    i = 1
		for f=1:length(FV)
			for k in FV[f]
	            I[i] = f
	            J[i] = k
	            V[i] = 1
	            i = i + 1
			end
		end
		M_2 = sparse(I,J,V)
		return M_2
	end

	"""
	From indexed base matrix to sparse array with 0 and 1
	"""
	function characteristicMatrix_push( FV::Array{Int64,2})
	    I,J,V = Int64[],Int64[],Int8[]
	    for f=1:size(FV, 1)
	        for k in FV[f, :]
	        push!(I,f)
	        push!(J,k)
	        push!(V,1)
	        end
	    end
	    M_2 = sparse(I,J,V)
	    return M_2
	end

	function characteristicMatrix_set( FV::Array{Int64,2})
	    sz = length(FV)
	    I = Array{Int64,1}(undef, sz)
	    J = Array{Int64,1}(undef, sz)
	    V = Array{Int8,1}(undef, sz)
	    i = 1
	    for f=1:size(FV, 1)
	        for k in FV[f, :]
	            I[i] = f
	            J[i] = k
	            V[i] = 1
	            i = i + 1

	        end
	    end
	    M_2 = sparse(I,J,V)
	    return M_2
	end

	# function characteristicMatrix_parallel( FV::Array{Int64,2})
	# #     sz = size(FV, 1)
	#     sz = length(FV)
	#     I = SharedArray{Int64,1}(sz)
	#     J = SharedArray{Int64,1}(sz)
	#     V = SharedArray{Int8,1}(sz)
	#     ii = SharedArray{Int64,1}(1)
	# #     I,J,V = Int64[],Int64[],Int8[]
	#     println("sz $sz")
	#     println(size(FV))
	#     println(length(FV))
	#     ii[1] = 1
	#     @sync @distributed for f=1:size(FV, 1)
	# #     for f=1:size(FV, 1)
	#         for k in FV[f, :]
	#             I[ii[1]] = f
	#             J[ii[1]] = k
	#             V[ii[1]] = 1
	#             ii[1] = ii[1] + 1
	#
	# #         push!(I,f)
	# #         push!(J,k)
	# #         push!(V,1)
	#         end
	#     end
	#     println("i=${ii[1]}")
	#     M_2 = sparse(I,J,V)
	#     return M_2
	# end
	characteristicMatrix(FV::Array{Int64,2}) = characteristicMatrix_set(FV::Array{Int64,2})
	characteristicMatrix(FV::LarSurf.Lar.Cells) = characteristicMatrix_set(FV::LarSurf.Lar.Cells)

	function characteristicMatrix(FV::Array{Int64,2}, nvertices)
    	mat = characteristicMatrix(FV)
    	return unsafe_resize(mat, size(mat, 1), nvertices)
	end

    """
    Expand indexes into sparse horizontally.
    Convert CV with vertex indexes (Int64) to sparse CV with just zeros and ones.
	Slow implementation.

    It should be the same as lar.CharacteristicMatrix()
    """
    function _CVi_to_CV01(CVi, nvertices)
        nc = size(CVi)[1]
        nvertices_per_C = size(CVi)[2]
        CV01 = spzeros(Int8, nc, nvertices)

        for k in 1:nc
            for i in 1:nvertices_per_C
                # print("vert k: ", k,", i: ", i)
                vertind = CVi[k, i]
                # print(" vert ind: ", vertind, "\n")
                CV01[k, vertind] = 1
            end
        end
        return CV01
    end

    # CV01 = CVi_to_CV01(CVi, nvertices)

    # VF to sparse
    # print("VF to sparse")

    """
    Expand indexes into sparse verticaly.
    Convert CV with vertex indexes (Int64) to sparse CV with just zeros and ones.
	Slow implementation.

    """
    function _VFi_to_VF01(VFi, nvertices)
        nfaces = size(VFi)[2]
        VF01 = spzeros(Int8, nvertices, nfaces)

        for k in 1:size(VFi)[1]
            for i in 1:(size(VFi)[2])
                # print("vert k: ", k,", i: ", i)
                vertind = VFi[k, i]
                # print(" vert ind: ", vertind, "\n")
                VF01[vertind, i] = 1
            end
        end
        return VF01
    end
    # ind_to_sparse_v = VFi_to_VF01

    # VF01 = ind_to_sparse_v(VFi, nvertices)
    """
    Based on matrix with indexes expand the index into sparse array containing one
    and zero. Deprecated. Use LarSurf.Lar.characteristicMatrix().
    """
    function characteristicMatrix(CVi, nvertices, axis=2)
        if axis == 1
            CV01 = _VFi_to_VF01(CVi, nvertices)
        elseif axis == 2
			CV01 = characteristicMatrix(CVi, nvertices)
            # CV01 = _CVi_to_CV01(CVi, nvertices)
        else
            error("Axis should be 1 or 2")
        end
        return CV01
    end

    function characteristicMatrix_for_loop(CVi, nvertices, axis=2)
        if axis == 1
            CV01 = _VFi_to_VF01(CVi, nvertices)
        elseif axis == 2
            CV01 = _CVi_to_CV01(CVi, nvertices)
        else
            error("Axis should be 1 or 2")
        end
        return CV01
    end

    ind_to_sparse = characteristicMatrix

	"""
	Check two arrays. Return true if their elements are same.
	"""
	function match_arr(arr1, arr2)
    for i =1:length(arr1)
        if arr1[i] != arr2[i]
            println("Element i=$i is not in match. $(arr1[i]) != $(arr2[i])")
            return false
        end
    end
    return true
end
    """
    Input is list of lists
    Output is 2D Array
    """
    function ll2array(CVill::Array)
        nc = size(CVill)[1]
        CVi = Array{typeof(CVill[1][1])}(undef, size(CVill)[1], size(CVill[1])[1])
        for k in 1:nc
            CVi[k, :] = CVill[k]
        end
        return CVi
    end

	"""
	Convert 2D array to list of lists.
	"""
	function array2ll(CVarr::Union{Array{Int64,2}, Array{Int8,2}})
    na, nb = size(CVarr)
	arr = Array{Array{eltype(CVarr),1},1}(undef, na)
    for i=1:na
        arr[i] = CVarr[i, :]
    end
    return arr
	end

    """
    Create VV matrix representing identity matrix in list of lists format with indices.

    julia> createVVll(5)

    [[1], [2], [3], [4], [5]]
    """
    function createVVll(number_of_vertices)

        VVm32 = [[i] for i=1:size_of_vertices]
        return VVm32
    end

    """
    Convert size in tuple to size in array
    """
    function size_as_array(block_size)
        if typeof(block_size) == Tuple{Int64,Int64,Int64}
            block_size = [block_size[1], block_size[2], block_size[3]]
        end
        return block_size
    end

    function hard_max!(stop, data_size)
        for i=1:length(stop)
            if data_size[i] < stop[i]
                stop[i] = data_size[i]
            end
        end
        return stop
    end

    # function check_LARmodel(larmodel::Lar.LARmodel)
    function check_LARmodel(larmodel::Lar.LARmodel)
        V, topology = larmodel
        n_nodes = size(V)[2]

        for i=1:length(topology)
            mx = maximum(maximum(topology[i]))
            if n_nodes < mx
                error("Maximal node ID is $n_nodes but in $i-th topology is maximum $mx")
                return false
            end
        end
        return true

    end

# end


text:
"""
main.jl est le fichier principal.
"""

println("Precompiling packages...")

using JuMP, GLPK, LinearAlgebra, PyPlot

@enum Sense forth back
@enum State departure arrival

include("structs.jl")
include("dataManager.jl")
include("JuMPModels.jl")
include("tools.jl")

const COST_OPEN_STATION = 1900000
const COST_BRIDGE = 28600000
const COST_NEW_LINE_PER_KM = 5900000
const COST_VEHICLES = 250000
const COST_MAINTENANCE_VEHICLES_PER_YEAR = 35000
const MAX_VEHICLES = 2
const MAX_CAPAVEHICLE = 20
const SPEED_VEHICLE = 20

const STOP_TIME_STATION = 3

function formul(lat1::Float64, long1::Float64, lat2::Float64, long2::Float64)
    r = 6371008
    distAng = acos( sin(lat1 * pi/180) * sin(lat2 * pi/180) + cos(lat1 * pi/180) * cos(lat2 * pi/180) * cos((long1 - long2)*pi/180) )

    return r * distAng
end

function calculDistLatLong(stations::Vector{Station})

    nbStations = length(stations)
    distStations = Array{Float64, 2}(undef, nbStations, nbStations)

    for iter = 1:(nbStations-1)
		dist = formul(stations[iter].latitude, stations[iter].longitude, stations[iter+1].latitude, stations[iter+1].longitude)
        distStations[iter, iter+1] = dist
		distStations[iter+1, iter] = dist
    end

	for iterRow = 1:nbStations
		for iterCol = (iterRow+2):nbStations
			distStations[iterRow, iterCol] = 0.
			distStations[iterCol, iterRow] = 0.
			for iter = 1:(iterCol-iterRow)
				distStations[iterRow, iterCol] += distStations[iter, iter+1]
				distStations[iterCol, iterRow] += distStations[iter, iter+1]
			end
		end
	end
    return distStations
end

function mainPart1()
    stations = getStations("stations.dat")
    distStations = calculDistLatLong(stations)
    townsites = getTownsites("townsites.dat",stations)
    h = getDensities(townsites)
    Z = Vector{Int}(undef,length(stations))
    couts = zeros(Float64,length(stations))
    for p in 1:length(stations)
        x,z = MCP_model_Lucas(stations,townsites,p,false)

        Z[p] = z

        couts[p] = 1.9 * p
        if x[end].ind >= 10
        	couts[p] += 5.9 * distStations[10, x[end].ind]/1000
        end

    end
    plot(Z, couts, marker=".",label="Coût d'installation en fonction du nombre d'habitants qui ont au moins une station accessible")
    grid()
    title("Représentation pour chaque nombre de station ouverte, le cout en fonction de la densité de population touchée")
    xlabel("Densité de population")
    ylabel("Cout d'installation en Millions d'€")
    legend()
end

function jules()
    stations = getStations("stations.dat")
    println(stations)
    townsites = getTownsites("townsites.dat",stations)
    println("townsites = ",townsites)

    h = getDensities(townsites)
    println("h = ",h)
    x,z = MCP_model_Lucas(stations,townsites,10,true)
    println("z = ",z)
    println("x = ",x)
end

function lucas()
    stations = getStations("stations.dat")
    townsites = getTownsites("townsites.dat", stations)

    return stations, townsites
end

function distDidact()
    distStations = Array{Float64, 2}(undef, 3, 3)

    distStations[1,2] = 750
    distStations[2,3] = 1000
    distStations[2,1] = 750
    distStations[3,2] = 1000
    return distStations
end

function runPESP(;stations::Vector{Station}, distStations::Array{Float64, 2} = calculDistLatLong(stations), T::Int = 3600, nbShuttles::Int = 2)

    nbStations = length(stations)

    T, E, A_run, A_dwell, A_thr, A_head, A_reg, L, U = parserPESP(T, nbStations, nbShuttles, stations, distStations)
    m = PESP_model(T, E, A_run, A_dwell, A_thr, A_head, A_reg, L, U)

    if JuMP.termination_status(m) != MOI.OPTIMAL
        println("La periode est trop courte pour ce nombre de Navette, veuillez diminuer le nombre de Navette ou bien augmenter la periode.")
        return m
    else

		grid()
		title("Représentation graphique de la planification des horaires pour $nbShuttles navette(s)")
		xlabel("Temps en secondes")
		ylabel("Indice de la station actuelle")

		for iterStation in 1:nbStations
			station = stations[nbStations - iterStation + 1]
			plot([], [], "k", label = "$(nbStations - iterStation + 1) : $(station.name)")
		end

        nbNodes = length(E)
        nbNodesPerShuttle = Int(nbNodes/nbShuttles)

        zValues = JuMP.value.(JuMP.all_variables(m)[(nbNodes+1):end]) #Acces a z[i, j] : zValues[i + (j-1)*nbNodes] et Acces z[i-1, i] = zValues[ (i - 1) * (nbNodes + 1)]

        iterNode = 0

        for shuttle = 1:nbShuttles

            values = JuMP.value.(JuMP.all_variables(m)[1 + (shuttle - 1)*(4 * (nbStations - 1)):(Int(length(E)/nbShuttles) + (shuttle - 1)*(4 * (nbStations - 1)))])

            Y = Vector{Int}(undef, Int(length(E)/nbShuttles)) # Idée : Juste les vider genre empty!(X)
            X = Vector{Int}(undef, Int(length(E)/nbShuttles))

            nbAjoutT = 0

            for iter = 1:nbNodesPerShuttle
                iterNode += 1
                Y[iter] = E[iterNode].indStation

                if iter != 1
                    if zValues[(iterNode-1)*(nbNodes+1)] == 1
                        nbAjoutT += 1
                    end
                    X[iter] = (values[iter] + nbAjoutT * T)
                else
                    X[iter] = values[iter]
                end

            end

            plot(X, Y, label = "Navette $shuttle")
        end
	legend()
        return m
    end
end

function testDidactic(;T::Int = 3600, nbShuttles::Int = 2)

    nbStations = 3
    stations = [Station("Station "*string(i), rand(), rand()) for i = 1:nbStations]

    distStations = distDidact()

    return runPESP(stations = stations, distStations = distStations, T = T, nbShuttles = nbShuttles)

end

function nantesCarq(;nbStations::Union{Int, Missing} = missing, nbShuttles::Int = 2, T::Int = 3600)

	if typeof(nbStations) <: Missing
		mainPart1()
		println("Input the number of station you want to take : ")
		print(">> ")
		nbStations = parse(Int, readline())
		PyPlot.close()
	end

	allStations = getStations("stations.dat")
    allTownsites = getTownsites("townsites.dat",allStations)

	h = getDensities(allTownsites)
    x,z = MCP_model_Lucas(allStations,allTownsites,nbStations,false)

	indX = broadcast(stat -> stat.ind, x)

	distStations = calculDistLatLong(allStations)[indX, indX]

	return runPESP(stations = x, distStations = distStations, T = T, nbShuttles = nbShuttles)
end

#
# """
# Pour Jules :
#
# La fonction PESP_model te renvoie le modèle JuMP oprimisé.
#
# Fonctions utiles dessus :
#
# Récuperer le modèle :
#
# m = testDidactic()
#
# Les contraintes liées aux variables
#
# - JuMP.all_constraints(m, JuMP.VariableRef, MOI.GreaterThan{Float64})
# - JuMP.all_constraints(m, JuMP.VariableRef, MOI.LessThan{Float64})
# - JuMP.all_constraints(m, JuMP.VariableRef, MOI.ZeroOne)
#
# Les réels contraintes
#
# - JuMP.all_constraints(m, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.Interval{Float64})
#
# La fonction objectif
#
# - JuMP.objective_function(m)
#
# Le numéro d'indicage des contraintes est le numéro du sommet associé, l'ordre de construction est montré (si c'est lisible) sur le graphique que je t'ai envoyé sur Messenger.
#
# Cordialement
#
# Le problème vient des terminus sur le retour : En effet, quand on verifie les temps π_j - π_i, pour le terminus du retour on vient à faire :
#
# π_{premier de la navette} - π_{dernier de la navette} :
#
# Dans les faits on calcul : 40 <= 0 - 820 + z * 3600 <= 70
#
# Tu peux prendre z = 1 ou 0 ca ne satisfait pas la contrainte donc le problème est infaisable ! Je crois qu'on a mal compris une partie du sujet
# Donc si tu trouves l'erreur je t'appel maitre ! (Car l'erreur ne vient pas du code parce qu'il fait exactement tout ce que je veux mais plutôt de la comprehension du problème.)
# """

ext: jl

max_stars_repo_name: LucasBaussay/ProjetODC

max_issues_repo_name: LucasBaussay/ProjetODC

max_forks_repo_name: LucasBaussay/ProjetODC


text:
"""
Test script for dirBin0Rec0 model: one parameter for total number of links and one for reciprocity
"""

include("../../../add_load_paths.jl")

using StaticNets

# test sampling

nSample = 500
N=30
    θ_0 = -3.1
    η_0 = 2
    model = ErgmDirBin0Rec0()
    diadProb = diadProbFromPars(model, [θ_0, η_0])
    A_vec = [samplSingMatCan(model, diadProb, N) for i=1:nSample]

using PyPlot
hist([sum(A) for A in A_vec])
hist([sum(A'.*A) for A in A_vec])

# Test Estimates
par_est = reduce(hcat,[estimate(model, A) for A in A_vec])
hist(par_est[2,:], density=true )
 vlines(η_0, 0, 1, "r")
hist(par_est[1,:], density=true )
 vlines(θ_0, 0, 1, "r")

# Are likelihood and loglikelihoods equal??
using ErgmRcall

#


text:
function logsumexp(logx::Vector{Float64})
    maxlog = maximum(logx)
    maxlog + log(sum(exp.(logx - maxlog)))
end

function logsumexp(logx::Matrix{Float64}, dim::Int)
    maxlog = maximum(logx, dim)
    maxlog + log(sum(exp.(logx .- maxlog), dim))
end


text:
# This file specifies some saving functions for configuration
# and iteration data.
# I am not yet sure whether to make the below methods in the
# interface definitions eventually...

import FileIO;

function _jld2_filename( filename :: AbstractString )
  fn, ext = splitext(filename)
  if isempty(ext) || ext ≠ ".jld2"
    return string(fn, ".jld2")
  else
    return filename
  end
end

"""
    save_config(filename, ac :: AbstractConfig )

Save the configuration object `ac` at path `filename`.
Ensures, that the file extension is `.jld2`.
The fieldname to retrieve the database object is `database`.

Returns the save path if successful and `nothing` else.
"""
function save_config( filename :: AbstractString, ac :: AbstractConfig ) :: Union{Nothing, AbstractString}
  fn = _jld2_filename(filename)
  try
    FileIO.save(fn, Dict("abstract_config" => ac))
    return fn
  catch e
    @warn "Could not save algorithm configuration:" exception=(ex,bt)
    return nothing
  end
end

"""
    load_config(fn :: AbstractString)
Load and return the `AbstractConfig` that was saved previously
with [`save_config`](@ref).
Return `nothing` on error.
"""
function load_config( filename :: AbstractString ) :: Union{AbstractConfig, Nothing}
  try
    return FileIO.load( filename )["abstract_config"]
  catch e
    @warn "Could not load algorithm configuration:" exception=(ex,bt)
    return nothing
  end
end

"""
    save_database(filename, id :: AbstractIterData )

Save the database that is referenced by `db(id)` at path `filename`.
Ensures, that the file extension is `.jld2`.
The fieldname to retrieve the database object is `database`.

Returns the save path if successful and `nothing` else.
"""
function save_database(filename::AbstractString, id :: AbstractIterData):: Union{Nothing, AbstractString}
  return save_database(filename, db(id))
end

"""
    save_database(filename, DB :: AbstractDB )

Save the database `DB` at path `filename`.
Ensures, that the file extension is `.jld2`.
The fieldname to retrieve the database object is `database`.

Returns the save path if successful and `nothing` else.
"""
function save_database( filename :: AbstractString, DB :: AbstractDB ) :: Union{Nothing, AbstractString}
  fn = _jld2_filename(filename)
  try
    FileIO.save(fn, Dict("database" => DB))
    return fn
  catch e
    @warn "Could not save evaluation database:" exception=(ex,bt)
    return nothing
  end
end

"""
    load_database(fn :: AbstractString)
Load and return the `<:AbstractDB` that was saved previously
with [`save_database`](@ref).
Return `nothing` on error.
"""
function load_database( filename :: AbstractString ) :: Union{AbstractDB, Nothing}
  try
    return FileIO.load( filename )["database"]
  catch e
    @warn "Could not load evaluation database:" exception=(ex,bt)
    return nothing
  end
end

"""
    save_iter_data(filename, id :: AbstractIterData )

Save the whole object `id` at path `filename`.
Ensures, that the file extension is `.jld2`.
The fieldname to retrieve the database object is `database`.

Returns the save path if successful and `nothing` else.
"""
function save_iter_data( filename :: AbstractString, id :: AbstractIterData ) :: Union{Nothing, AbstractString}
  fn = _jld2_filename(filename)
  try
    FileIO.save(fn, Dict("iter_data" => id))
    return fn
  catch e
    @warn "Could not save iteration data object:" exception=(ex,bt)
    return nothing
  end
end

"""
    load_iter_data(fn :: AbstractString)
Load and return the `<:AbstractIterData` that was saved previously
with [`save_iter_data`](@ref).
Return `nothing` on error.
"""
function load_iter_data( filename :: AbstractString ) :: Union{AbstractDB, Nothing}
  try
    return FileIO.load( filename )["iter_data"]
  catch e
    @warn "Could not load iteration data object:" exception=(ex,bt)
    return nothing
  end
end


text:
function libsass_version()
    s = ccall((:libsass_version, libsass_so), Cstring, ())
    unsafe_string(s)
end


text:
# low-level wrappers of the CUSOLVER library

#helper functions
function cusolverDnCreate()
  handle = Ref{cusolverDnHandle_t}()
  @check ccall((:cusolverDnCreate, libcusolver),
               cusolverStatus_t,
               (Ptr{cusolverDnHandle_t},),
               handle)
  return handle[]
end

function cusolverDnDestroy(handle)
  @check ccall((:cusolverDnDestroy, libcusolver),
               cusolverStatus_t,
               (cusolverDnHandle_t,),
               handle)
end

function cusolverDnSetStream(handle, streamId)
  @check ccall((:cusolverDnSetStream, libcusolver),
               cusolverStatus_t,
               (cusolverDnHandle_t, CuStream_t),
               handle, streamId)
end

function cusolverDnGetStream(handle, streamId)
  @check ccall((:cusolverDnGetStream, libcusolver),
               cusolverStatus_t,
               (cusolverDnHandle_t, Ptr{CuStream_t}),
               handle, streamId)
end

function cusolverSpCreate()
  handle = Ref{cusolverSpHandle_t}()
  @check ccall((:cusolverSpCreate, libcusolver),
               cusolverStatus_t,
               (Ptr{cusolverSpHandle_t},),
               handle)
  return handle[]
end

function cusolverSpDestroy(handle)
  @check ccall((:cusolverSpDestroy, libcusolver),
               cusolverStatus_t,
               (cusolverSpHandle_t,),
               handle)
end

function cusolverSpSetStream(handle, streamId)
  @check ccall((:cusolverSpSetStream, libcusolver),
               cusolverStatus_t,
               (cusolverSpHandle_t, CuStream_t),
               handle, streamId)
end

function cusolverSpGetStream(handle, streamId)
  @check ccall((:cusolverSpGetStream, libcusolver),
               cusolverStatus_t,
               (cusolverSpHandle_t, Ptr{CuStream_t}),
               handle, streamId)
end

function cusolverSpCreateCsrqrInfo(info)
  @check ccall((:cusolverSpCreateCsrqrInfo, libcusolver),
               cusolverStatus_t,
               (Ptr{csrqrInfo_t},),
               info)
end

function cusolverSpDestroyCsrqrInfo(info)
  @check ccall((:cusolverDestroyCsrqrInfo, libcusolver),
               cusolverStatus_t,
               (csrqrInfo_t,),
               info)
end

function cusolverDnCreateGesvdjInfo(info)
  @check ccall((:cusolverDnCreateGesvdjInfo, libcusolver),
               cusolverStatus_t,
               (Ptr{gesvdjInfo_t},),
               info)
end

function cusolverDnDestroyGesvdjInfo(info)
  @check ccall((:cusolverDnDestroyGesvdjInfo, libcusolver),
               cusolverStatus_t,
               (gesvdjInfo_t,),
               info)
end

function cusolverDnXgesvdjSetTolerance(info, tolerance)
  @check ccall((:cusolverDnXgesvdjSetTolerance, libcusolver),
               cusolverStatus_t,
               (gesvdjInfo_t, Float64),
               info, Float64(tolerance))
end

function cusolverDnXgesvdjSetMaxSweeps(info, max_sweeps)
  @check ccall((:cusolverDnXgesvdjSetMaxSweeps, libcusolver),
               cusolverStatus_t,
               (gesvdjInfo_t, Cint),
               info, Cint(max_sweeps))
end

function cusolverDnCreateSyevjInfo(info)
  @check ccall((:cusolverDnCreateSyevjInfo, libcusolver),
               cusolverStatus_t,
               (Ptr{syevjInfo_t},),
               info)
end

function cusolverDnDestroySyevjInfo(info)
  @check ccall((:cusolverDnDestroySyevjInfo, libcusolver),
               cusolverStatus_t,
               (syevjInfo_t,),
               info)
end

function cusolverDnXsyevjSetTolerance(info, tolerance)
  @check ccall((:cusolverDnXsyevjSetTolerance, libcusolver),
               cusolverStatus_t,
               (syevjInfo_t, Float64),
               info, Float64(tolerance))
end

function cusolverDnXsyevjSetMaxSweeps(info, max_sweeps)
  @check ccall((:cusolverDnXsyevjSetMaxSweeps, libcusolver),
               cusolverStatus_t,
               (syevjInfo_t, Cint),
               info, Cint(max_sweeps))
end

function cusolverRfCreate(handle)
  @check ccall((:cusolverRfCreate, libcusolver),
               cusolverStatus_t,
               (Ptr{cusolverRfHandle_t},),
               handle)
end

function cusolverRfDestroy(handle)
  @check ccall((:cusolverRfDestroy, libcusolver),
               cusolverStatus_t,
               (cusolverRfHandle_t,),
               handle)
end

function cusolverRfSetStream(handle, streamId)
  @check ccall((:cusolverRfSetStream, libcusolver),
               cusolverStatus_t,
               (cusolverRfHandle_t, CuStream_t),
               handle, streamId)
end

function cusolverRfGetStream(handle, streamId)
  @check ccall((:cusolverRfGetStream, libcusolver),
               cusolverStatus_t,
               (cusolverRfHandle_t, Ptr{CuStream_t}),
               handle, streamId)
end

function cusolverGetProperty(property::CUDAapi.libraryPropertyType)
  value_ref = Ref{Cint}()
  @check ccall((:cusolverGetProperty, libcusolver),
               cusolverStatus_t,
               (Cint, Ptr{Cint}),
               property, value_ref)
  value_ref[]
end


text:
# library
# Since @rimport and @rlibrary create module objects which may be conflict with other objects,
# it is safer to place them at the end of the test.
@rimport MASS as mass
@test_approx_eq rcopy(rcall(mass.ginv, RObject([1 2; 0 4]))) [1 -0.5; 0 0.25]
@rlibrary MASS
@test_approx_eq rcopy(rcall(ginv, RObject([1 2; 0 4]))) [1 -0.5; 0 0.25]


text:
u, v, w = Vertex(), Vertex(), Vertex()
a1, a2, b1, b2, c1, c2 = [Edge() for _ in 1:6]
f1, f2 = Face(), Face()

u.edge, v.edge, w.edge = a1, b1, c1

a1.dest, b1.dest, c1.dest = u, v, w
a2.dest, b2.dest, c2.dest = w, u, v

a1.twin, a2.twin = a2, a1
b1.twin, b2.twin = b2, b1
c1.twin, c2.twin = c2, c1

a1.next, b1.next, c1.next = b1, c1, a1
a1.prev, b1.prev, c1.prev = c1, a1, b1

a2.next, b2.next, c2.next = c2, a2, b2
a2.prev, b2.prev, c2.prev = b2, c2, a2

a1.face = b1.face = c1.face = f1
a2.face = b2.face = c2.face = f2

f1.edge, f2.edge = a1, a2

@testset "Edge path" begin
    el = collect(edge_path(a1, b1))
    @test el == [a1]

    el = collect(edge_path(a1, c1))
    @test el == [a1, b1]

    el = collect(edge_path(a1, c1, reverse=true))
    @test el == [a1]

    el = collect(edge_path(a1, b1, reverse=true))
    @test el == [a1, c1]
end

@testset "Edge loop" begin
    el = collect(edge_loop(a1))
    @test el == [a1, b1, c1]

    el = collect(edge_loop(a1, reverse=true))
    @test el == [a1, c1, b1]
end

@testset "Edges - vertex" begin
    el = collect(edges(u))
    @test el == [a1, b2]

    el = collect(edges(u, reverse=true))
    @test el == [a1, b2]
end

@testset "Edges - face" begin
    el = collect(edges(f1))
    @test el == [a1, b1, c1]

    el = collect(edges(f1, reverse=true))
    @test el == [a1, c1, b1]
end

@testset "Vertices - edge" begin
    vl = collect(vertices(a1))
    @test vl == [u, w]

    vl = collect(vertices(a1, reverse=true))
    @test vl == [w, u]
end

@testset "Vertices - face" begin
    vl = collect(vertices(f1))
    @test vl == [u, v, w]

    vl = collect(vertices(f1, reverse=true))
    @test vl == [u, w, v]
end

@testset "Faces - edge" begin
    fl = collect(faces(a1))
    @test fl == [f1, f2]

    fl = collect(faces(a1, reverse=true))
    @test fl == [f2, f1]
end

@testset "Faces - vertex" begin
    fl = collect(faces(u))
    @test fl == [f1, f2]

    fl = collect(faces(u, reverse=true))
    @test fl == [f1, f2]
end

@testset "Faces - face" begin
    fl = collect(faces(f1))
    @test fl == [f2, f2, f2]

    fl = collect(faces(f1, reverse=true))
    @test fl == [f2, f2, f2]
end


text:
# literal zero for initialisation
zero(::Type{Quire8}) = Quire8(0x0000_0000)
zero(::Type{Quire16}) = Quire16(0x0000_0000_0000_0000_0000_0000_0000_0000)
zero(::Type{Quire32}) = Base.zext_int(Quire32,0x0)

# literal one
one(::Type{Quire8}) = Quire8(0x00001000)
one(::Type{Quire16}) = Quire16(0x01000000000000000000000000000000)
one(::Type{Quire32}) = fma(zero(Quire32),one(Posit32),one(Posit32))

# Fused multiply-add q+b*c
fma(q::Quire8,b::Posit8,c::Posit8) = ccall((:q8_fdp_add, SoftPositPath), Quire8, (Quire8,Posit8,Posit8),q,b,c)
fma(q::Quire16,b::Posit16,c::Posit16) = ccall((:q16_fdp_add, SoftPositPath), Quire16, (Quire16,Posit16,Posit16),q,b,c)
fma(q::Quire32,b::Posit32,c::Posit32) = ccall((:q32_fdp_add, SoftPositPath), Quire32, (Quire32,Posit32,Posit32),q,b,c)

# Fused multiply-sub q-b*c
fms(q::Quire8,a::Posit8,b::Posit8) = ccall((:q8_fdp_sub, SoftPositPath), Quire8, (Quire8,Posit8,Posit8),q,a,b)
fms(q::Quire16,a::Posit16,b::Posit16) = ccall((:q16_fdp_sub, SoftPositPath), Quire16, (Quire16,Posit16,Posit16),q,a,b)
fms(q::Quire32,a::Posit32,b::Posit32) = ccall((:q32_fdp_sub, SoftPositPath), Quire32, (Quire32,Posit32,Posit32),q,a,b)


text:
struct CartpoleEnvironment{E<:CartpoleSimulator,RNG,F,R,I,T} <: Environment
    sim::E
    rng::RNG
    actions::Vector{Float64}
    discount::Float64
    statefun::F # maps environment state to env state
    rewardfun::R # maps envionment state plus action to reward
    initfun::I
    isterminalfun::T
end

function cartpole_statefun(env::CartpoleEnvironment, simstate)
    config = configuration(simstate)
    vel = velocity(simstate)
    @inbounds begin
        cart_x = config[1]
        cart_x_dot = vel[1]
        pole_angle = config[2]
        pole_angle_dot = vel[2]
    end
    SVector(cart_x, cart_x_dot, pole_angle, pole_angle_dot)
end

function cartpole_isterminalfun(pos_limit, angle_limit)
    function (env::CartpoleEnvironment, s)
        abs(s[1]) >= pos_limit || abs(s[3]) >= angle_limit
    end
end

# use common default values
function CartpoleEnvironment(;
        M = 1., m = 0.1, l = 1., com = l/2, dt = 2e-2,
        pos_limit = 2.4, angle_limit = deg2rad(12),
        rng = Base.GLOBAL_RNG,
        actions = [-10.,0.,10.],
        discount = 0.99,
        statefun = cartpole_statefun,
        rewardfun = (env, a, simstate) -> 1,
        initfun = (env) -> (0.1 .* rand(env.rng,SVector{4}) .- 0.05),
        isterminalfun = cartpole_isterminalfun(pos_limit, angle_limit)
       )
    sim = CartpoleSimulator(M, m, l, com, (-pos_limit, pos_limit), dt)
    CartpoleEnvironment(sim, rng, convert(Vector{Float64}, actions), Float64(discount), statefun, rewardfun, initfun, isterminalfun)
end

function Base.show(io::IO, env::CartpoleEnvironment)
    println(io, typeof(env).name)
    println(io, "  time: ", env.sim.time, " s")
    print(io, "  state: ", env.statefun(env, env.sim.state))
end

POMDPs.discount(env::CartpoleEnvironment) = env.discount
POMDPs.actions(env::CartpoleEnvironment) = env.actions
POMDPs.action_index(env::CartpoleEnvironment, a) = findfirst(env.actions, a)
POMDPs.n_actions(env::CartpoleEnvironment) = length(env.actions)

function reset!(env::CartpoleEnvironment)
    simstate = reset!(env.sim, env.initfun(env))
    env.statefun(env, simstate)
end

function step!(env::CartpoleEnvironment, a)
    simstate = step!(env.sim, a)
    env.statefun(env, simstate), env.rewardfun(env, a, simstate)
end

function POMDPs.isterminal(env::CartpoleEnvironment)
    env.isterminalfun(env, env.statefun(env, env.sim.state))
end


text:
using Test, SpeedDate, BioSequences

@testset "Mutation Counting" begin
    @testset "Catching every mutation type" begin
        function generate_possibilities_tester(::Type{A}) where {A <: NucleicAcidAlphabet}
            sym = symbols(A())
            arra = Vector{eltype(A)}()
            arrb = Vector{eltype(A)}()
            for i in 1:length(sym), j in i:length(sym)
                push!(arra, sym[i])
                push!(arrb, sym[j])
            end
            return LongSequence{A}(arra), LongSequence{A}(arrb)
        end

        @testset "2 - Bit" begin
            dnaA, dnaB = generate_possibilities_tester(DNAAlphabet{2})
            @test SpeedDate.count_n_mutations(dnaA, dnaB) == SpeedDate.count_n_mutations(dnaB, dnaA) == 6
        end

        @testset "4 - Bit" begin
            dnaA, dnaB = generate_possibilities_tester(DNAAlphabet{4})
            @test SpeedDate.count_n_mutations(dnaA, dnaB) == SpeedDate.count_n_mutations(dnaB, dnaA) == (6, 10)
        end
    end

    @testset "Randomized tests" begin

    end
end


text:
# Material files
include("./data/Cu.jl")
include("./data/SiO2.jl")
include("./data/resistivity_n-Si.jl")
include("./data/resistivity_p-Si.jl")

# Create interpolations
# Cu
knots = (cu_data[:,1],)
const itp_cu_re = interpolate(knots, cu_data[:,2], Gridded(Linear()))
const itp_cu_im = interpolate(knots, cu_data[:,3], Gridded(Linear()))

const Cu = MaterialFile(:refr_ind, :wavelength, itp_cu_re,itp_cu_im)

# SiO2
knots = (sio2_data[:,1],)
const itp_sio2_re = interpolate(knots, sio2_data[:,2], Gridded(Linear()))
const itp_sio2_im = interpolate(knots, sio2_data[:,3], Gridded(Linear()))

const SiO2 = MaterialFile(:refr_ind, :wavelength, itp_sio2_re,itp_sio2_im)

# Resitivity n-Si
knots = (nSi_data[:,1],)
const itp_nSi = interpolate(knots, nSi_data[:,2], Gridded(Linear()))
const nSi_sze = ResistivityFile(itp_nSi)

# Resitivity p-Si
knots = (pSi_data[:,1],)
const itp_pSi = interpolate(knots, pSi_data[:,2], Gridded(Linear()))
const pSi_sze = ResistivityFile(itp_pSi)


text:
using FunctionalCollections
using Test
import Base.vec

function vec(r::UnitRange)
    v = PersistentVector{Int}()
    for i in r
      v = push(v, i)
    end
    v
end

function Base.Array(r::UnitRange)
    arr = Array(Int, r)
    for i in r
      arr[i] = i
    end
    arr
end

@testset "Persistent Vectors" begin

    @testset "range constructor" begin
        @test typeof(vec(1:1000)) == PersistentVector{Int}
        @test typeof(pop(vec(1:1000))) == PersistentVector{Int}
    end

    @testset "length" begin
        @test length(vec(1:32)) == 32
        @test length(vec(1:10000)) == 10000
        @test length(pop(vec(1:1000))) == 999
    end

    @testset "accessing elements" begin
        pv = vec(1:5000)
        @test length(pv) == lastindex(pv) == 5000

        @test pv[1]    == 1
        @test pv[32]   == 32
        @test pv[500]  == 500
        @test pv[2500] == 2500
        @test pv[5000] == 5000
        @test try pv[5001]; false catch e true end

        @test try vec(1:32)[33]; false catch e true end
    end

    @testset "accessing last" begin
        @test peek(vec(1:1000)) == 1000
        @test vec(1:1000)[end]  == 1000
    end

    @testset "associng" begin
        @test assoc(vec(1:1000), 500, 1)[500] == 1
    end

    @testset "appending" begin
        # inference problems in 0.5.0, seem fixed in 0.6.0-dev
        @test append(vec(1:31), 32)        == vec(1:32)
        @test append(vec(1:999), 1000)     == vec(1:1000)
        @test append(vec(1:9999), 10000)   == vec(1:10000)
        @test append(vec(1:99999), 100000) == vec(1:100000)
    end

    @testset "structural sharing" begin
        pv = vec(1:32)
        @test length(pv) == lastindex(pv) == 32
        pv2 = append(pv, 33)
        @test length(pv2) == lastindex(pv2) == 33
        @test pv2.trie[1] === pv.tail
    end

    @testset "equality" begin
        v1 = vec(1:1000)
        v2 = vec(1:1000)

        @test v1.trie !== v2.trie
        @test v1 == v2

        @test isequal(v1, v2)
    end

    @testset "iteration" begin
        arr2 = Int[]
        for i in vec(1:10000)
            push!(arr2, i)
        end
        @test collect(1:10000) == arr2
    end

    @testset "map" begin
        v1 = vec(1:5)
        @test map((x)->x+1, v1) == PersistentVector([2, 3, 4, 5, 6])
        v2 = PersistentVector{Tuple{Int,Int}}([(1,2),(4,3)])
        @test map((x)->(x[2],x[1]), v2) == PersistentVector{Tuple{Int,Int}}([(2,1),(3,4)])
    end

    @testset "filter" begin
        v1 = vec(1:5)
        @test filter(iseven, v1) == PersistentVector([2, 4])
        v2 = PersistentVector{Tuple{Int,Int}}([(1,2),(4,3)])
        @test filter((x)->x[2] > x[1], v2) == PersistentVector{Tuple{Int,Int}}([(1,2)])
    end

    @testset "hash" begin
        @test hash(vec(1:1000)) == hash(vec(1:1000))
    end

    @testset "isempty" begin
        @test isempty(PersistentVector{Int}())
        @test !isempty(PersistentVector([1]))
    end

    @testset "identity constructor" begin
        pv = vec(1:10)
        @test pv === PersistentVector(pv)
        @test pv === PersistentVector{Int}(pv)
        @test typeof(PersistentVector{Any}(pv)) == PersistentVector{Any}
    end

    @testset "convert" begin
        @test convert(PersistentVector, [1, 2, 3]) == PersistentVector([1, 2, 3])
        @test convert(PersistentVector{Int}, [1, 2, 3]) == PersistentVector([1, 2, 3])

        @test typeof(convert(PersistentVector{Any}, [1, 2, 3])) == PersistentVector{Any}
        @test convert(PersistentVector{Float64}, [1, 2, 3]) == [1.0, 2.0, 3.0]

        # Test identity conversion
        pv = vec(1:10)
        @test convert(PersistentVector, pv) === pv
        @test convert(PersistentVector{Int}, pv) === pv
        @test typeof(convert(PersistentVector{Float64}, pv)) == PersistentVector{Float64}
    end

end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export clean_for_loops

function clean_for_loops(statements)
    function walker(node, cbdata, index, top_level, read)
        if is_for(node) && length(node.args[2].args) == 1 &&
                isa(node.args[2].args[1], Expr) &&
                node.args[2].args[1].head == :block
            node.args[2].args = clean_for_loops(node.args[2].args[1].args)
            return node
        elseif isa(node, Expr) && node.head in [:loophead, :parallel_loophead, :loopend, :parallel_loopend]
            return node
        end
        ASTWALK_RECURSE
    end
    map((s) -> AstWalk(s, walker, nothing), statements)
end


text:
#=
# ADI Reduction of $\beta$ Pictoris

This example will walk through a full reduction and some analysis of $\beta$ Pictoris in order to show the basic usage of ADI.jl for high-contrast imaging. $\beta$ Pictoris is a known exoplanet host with many publications regarding the direct images of its substellar companion and circumstellar disk.

What will *not* be covered in this example are the basics of Julia, the fine details of ADI post-processing, or any reference documentation.

---

## Setup

Let's begin by importing the necessary libraries
=#
using ADI
using DataFrames
using HCIDatasets: BetaPictoris
using Plots

## set up plotting
function imshow(img, args...; kwargs...)
    ylim = extrema(axes(img, 1))
    xlim = extrema(axes(img, 2))
    heatmap(img, args...; aspect_ratio=1, xlim=xlim, ylim=ylim, kwargs...)
end;

#=
## Data Reduction

Here we load the data for $\beta$ Pictoris from NaCo at the VLT. You may be prompted to download the data; see [HCIDatasets.jl](https://github.com/JuliaHCI/HCIDatasets.jl) for more details.
=#
cube, angles = BetaPictoris[:cube, :pa];

#=
To reduce the data, we need an *algorithm*. In ADI.jl we currently have median subtraction, PCA, NMF, and fixed-point GreeDS. These algorithms are treated as "objects" in the sense that we initialize them with options and then pass them around inside the ADI.jl framework to retrieve the results we want.

The usage for fitting the speckle estimate, projecting and subtracting this estimate from the target cube, and derotating and collapsing the residual all are encompassed by calling the algorithm as a function.
=#
alg = PCA(10) # 10 components
reduced = alg(cube, angles)
imshow(reduced)

#=
To try out different algorithms, all you should have to do is change one line and re-run the remaining code. Let's briefly explore a few different algorithms
=#
algs = (alg, NMF(10), Median(), GreeDS(10))
reduced_frames = [alg(cube, angles) for alg in algs];
#-
figs = (imshow(reduced, ticks=false) for reduced in reduced_frames)
plot(
    figs...,
    title=["PCA(10)" "NMF(10)" "Median()" "GreeDS(10)"],
    layout=(2, 2),
    size=(900, 800),
    dpi=75
)

#=
You may want to mask out an interior angle since there is an inner limit for our signal to be a real planet (as opposed to systematics from the optical system or noise). We can mask out an interior circle either before processing with the algorithm or afterwards using `HCIToolbox.mask_circle`. HCIToolbox.jl is re-exported by ADI.jl, though so we have it available here
=#
mask_cube = mask_circle(cube, 10)
mask_reduced = alg(mask_cube, angles)
imshow(mask_reduced)

#=
## S/N and Significance Maps

Now that we have our reduced frame, let's look at the signal-to-noise ratio (SNR, S/N). We use the exact S/N calculation here, implemented in a fast, multi-threaded framework. In order to measure the S/N, though, we need the effective FWHM of our instrument. Normally, we would measure this from an off-axis (or non-coronagraphic) PSF, but for simplicity I'll hard-code a value.
=#
fwhm = 4.6
snrmap = detectionmap(snr, reduced, fwhm)
imshow(snrmap)

#=
If we want to get the statistical significance, we need to convert from the Student-t confidence interval derived in the S/N to the Gaussian significance. We can accomplish this by calling
=#
sigmap = detectionmap(significance, reduced, fwhm)
imshow(sigmap)

#=
Now, let's do some very basic frequentist planet detection by thresholding this significance
=#
sigmap_cutoff = @. sigmap > 5
imshow(sigmap_cutoff)

#=
looks like we've successfully pulled out the companion Beta Pictoris b from the data!

## Contrast Curve

We are also interested in analyzing how the algorithm affects our data, especially calculating the *throughput* and the *contrast curve*. These measure, effectively, how much signal is lost during the subtraction step of the algorithm and give us an idea of what the limits of our algorithm are with our data.

Before we move on, we need to create a PSF model for our data. `HCIToolbox.Kernels` includes some simple functional PSFs or we can use an empirical PSF. We will use the empirical PSF provided by HCIDatasets for our calculations
=#
psf = BetaPictoris[:psf] ./ maximum(BetaPictoris[:psf])
kern_psf = Kernels.Normal(fwhm);
#-
plot(
    imshow(psf),
    imshow(construct(kern_psf, size(psf))),
    layout=2,
    size=(500, 250),
    cbar=false,
    ticks=false
)

#=
and now we can calculate the 5σ contrast curve
=#

cc = contrast_curve(alg, cube, angles, psf; fwhm=fwhm, nbranch=6) |> DataFrame
first(filter(row -> isfinite(row.contrast_corr), cc), 5)

# and lets plot it
plot(
    cc.distance,
    [cc.contrast_corr cc.contrast],
    yscale=:log10,
    xlim=(0, NaN),
    label=["Student-t" "Gaussian"],
    ylabel="5-sigma contrast",
    xlabel="radius [px]"
)
#=
You'll notice a pretty severe bump indicating poor contrast where the original companion is! Because this companion is very bright (it can even be seen in the raw frames!) it will bias the contrast measurement.

Typically you'd like to fit the companion signal and remove it in a maximum likelihood framework. For convenience here, let's use the `:cube_empty` entry for `BetaPictoris`, which already has the companion removed.
=#
cube_empty = BetaPictoris[:cube_empty]
reduced_empty = alg(cube_empty, angles)
imshow(reduced_empty)
#-
cc_empty = contrast_curve(alg, cube_empty, angles, psf; fwhm=fwhm, nbranch=6) |> DataFrame
first(filter(row -> isfinite(row.contrast_corr), cc_empty), 5)
#-
plot(
    cc_empty.distance,
    [cc_empty.contrast_corr cc_empty.contrast],
    yscale=:log10,
    xlim=(0, NaN),
    label=["Student-t" "Gaussian"],
    ylabel="5-sigma contrast",
    xlabel="radius [px]"
)


text:
# ------------------------------------------------------------------------------------------
# # Modeling data 2
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Building a model
#
# Recall that in notebook 3, we saw that we could use a mathematical function to classify an
# image as an apple or a banana, based on the average amount of green in an image:
#
#
#
#
#
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# A common function for performing this kind of **classification** is the sigmoid that we
# saw in the last notebook, and that we will now extend by adding two **parameters**, $w$
# and $b$:
#
# $$\sigma(x; w, b) := \frac{1}{1 + \exp(-wx + b)}$$
#
# $$ x = \mathrm{data} $$
#
# \begin{align}
# \sigma(x;w,b) &\approx 0 \implies \mathrm{apple} \\
# \sigma(x;w,b) &\approx 1 \implies \mathrm{banana}
# \end{align}
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# In our mathematical notation above, the `;` in the function differentiates between the
# **data** and the **parameters**. `x` is the data and is determined from the image. The
# parameters, `w` and `b`, are numbers which we choose to make our function match the
# results it should be modeling.
#
# Note that in the code below, we don't distinguish between data and parameters - both are
# just inputs to our function, σ!
# ------------------------------------------------------------------------------------------

using Images

apple = load("d
    ata/10_100.jpg")
banana = load("data/104_100.jpg")

apple_green_amount = mean(Float64.(green.(apple)))
banana_green_amount = mean(Float64.(green.(banana)))

"Average green for apple = $apple_green_amount; " *
"Average green for banana = $banana_green_amount; "

σ(x, w, b) = 1 / (1 + exp(-w * x + b))

# ------------------------------------------------------------------------------------------
# What we want is that when we give σ as input the average green for the apple, roughly `x =
# 0.3385`, it should return as output something close to 0, meaning "apple". And when we
# give σ the input `x = 0.8808`, it should output something close to 1, meaning "banana".
#
# By changing the parameters of the function, we can change the shape of the function, and
# hence make it represent, or **fit**, the data better!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Data fitting by varying parameters
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can understand how our choice of `w` and `b` affects our model by seeing how our values
# for `w` and `b` change the plot of the $\sigma$ function.
#
# To do so, we will use the `Interact.jl` Julia package, which provides "widgets" for
# controlling parameters interactively via sliders:
# ------------------------------------------------------------------------------------------

using Plots; gr()   # GR works better for interactive manipulations
using Interact      # package for interactive manipulation

# ------------------------------------------------------------------------------------------
# Run the code in the next cell. You should see two "sliders" appear, one for `w` and one
# for `b`.
#
# **Game**:
# Move both of those sliders around until the blue curve, labeled "model", which is the
# graph of the `\sigma` function, passes through *both* of the data points at the same time.
# ------------------------------------------------------------------------------------------

@manipulate for w in -10:0.01:30, b in 0:0.1:20

    plot(x -> σ(x, w, b), xlim=(-0,1), ylim=(-0.1,1.1), label="model", legend=:topleft, lw=3)

    scatter!([apple_green_amount],  [0.0], label="apple", ms=5)   # marker size = 5
    scatter!([banana_green_amount], [1.0], label="banana", ms=5)

end

# ------------------------------------------------------------------------------------------
# Notice that the two parameters do two very different things. The **weight**, `w`,
# determines *how fast* the transition between 0 and 1 occurs. It encodes how trustworthy we
# think our data  actually is, and in what range we should be putting points between 0 and 1
# and thus calling them "unsure". The **bias**, `b`, encodes *where* on the $x$-axis the
# switch should take place. It can be seen as shifting the function left-right. We'll come
# to understand these *parameters* more in notebook 6.
#
# Here are some parameter choices that work well:
# ------------------------------------------------------------------------------------------

w = 25.58; b = 15.6

plot(x -> σ(x, w, b), xlim=(0,1), ylim=(-0.1,1.1), label="model", legend=:topleft, lw=3)

scatter!([apple_green_amount], [0.0], label="apple")
scatter!([banana_green_amount],[1.0], label="banana")

# ------------------------------------------------------------------------------------------
# (Note that in this problem there are many combinations of `w` and `b` that fit the data
# well.)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Once we have a model, we have a computational representation for how to choose between
# "apple" and "banana". So let's pull in some new images and see what our model says about
# them!
# ------------------------------------------------------------------------------------------

apple2 = load("data/107_100.jpg")

green_amount = mean(Float64.(green.(apple2)))
@show green_amount

scatter!([green_amount], [0.0], label="new apple")

# ------------------------------------------------------------------------------------------
# Our model successfully says that our new image is an apple! Pat yourself on the back:
# you've actually just trained your first neural network!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
# Load the image of a banana in `data/8_100.jpg` as `mybanana`. Edit the code below to
# calculate the amount of green in `mybanana` and to overlay data for this image with the
# existing model and data points.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # To get the desired overlay, the code we need is
#
# ```julia
# mybanana = load("data/8_100.jpg")
# mybanana_green_amount = mean(Float64.(green.(banana)))
# scatter!([mybanana_green_amount], [1.0], label="my banana")
# ```
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Closing remarks: bigger models, more data, more accuracy
#
# That last apple should start making you think: not all apples are red; some are yellow.
# "Redness" is one attribute of being an apple, but isn't the whole thing. What we need to
# do is incorporate more ideas into our model by allowing more inputs. However, more inputs
# would mean more parameters to play with. Also, we would like to have the computer start
# "learning" on its own, instead of modifying the parameters ourselves until we think it
# "looks right". How do we take the next step?
#
# The first thing to think about is, if you wanted to incorporate more data into the model,
# how would you change the sigmoid function? Play around with some ideas. But also, start
# thinking about how you chose parameters. What process did you do to finally end up at good
# parameters? These two problems (working with models with more data and automatically
# choosing parameters) are the last remaining step to understanding deep learning.
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/05. ML - Building models.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
type Distribution{FS<:FeatureSpace}
  space::FS
  features::Number
  total::Number
  smooth::Function
  smooth_data::Array
  mdata::Any
  Distribution() = new(FS(),0,0,_no_smoothing,[])
  Distribution(fv::FeatureVector) = new(fv,length(fv),get_total(fv),_no_smoothing,[])
  Distribution(c::Cluster) = new(c,length(c.vector_sum),get_total(c.vector_sum) ,_no_smoothing,[])
  Distribution(ds::DataSet) = new(ds,length(ds.vector_sum),get_total(ds.vector_sum) ,_no_smoothing,[])
  function get_total(fv::FeatureVector)
    total = 0
    for value in values(fv)
      total += value
    end
    return total
  end
end
Distribution(fv::FeatureVector) = Distribution{FeatureVector}(fv)
Distribution(c::Cluster) = Distribution{Cluster}(c)
Distribution(ds::DataSet) = Distribution{DataSet}(ds)

function getindex(d::Distribution, key)
  return d.space[key]
end

function get_num_texts_in_dist(d::Distribution{DataSet})
  texts = 0
  for clust in values(d.space.clusters)
    texts += length(clust)
  end
  return texts
end

function get_num_texts_given_feature(d::Distribution{DataSet}, feature)
  texts = 0
  for clust in values(d.space.clusters)
    for fv in values(clust.vectors)
      if haskey(fv,feature)
        texts += 1
      end
    end
  end
  return texts
end

function probability(d::Distribution, feature)
  return d.smooth(d, feature, d.smooth_data)
end

function prob_fv(d::Distribution{Cluster}, key)
  return d.space[key].total / d.total
end

function prob_fv(d::Distribution{DataSet}, clust_key, key)
  return d.space[clust_key][key].total / d.total
end

function prob_cl(d::Distribution{DataSet}, key)
  return d.space[key].vector_sum.total / d.total
end

function cond_prob_f_given_fv(d::Distribution{Cluster}, fv_key, key)
  return d.space[fv_key][key] / d.space[fv_key].total
end

function cond_prob_f_given_fv(d::Distribution{DataSet}, clust_key, fv_key, key)
  return d.space[clust_key][fv_key][key] / d.space[clust_key][fv_key].total
end

function cond_prob_f_given_fv_in_clust(d::Distribution{DataSet}, clust_key, fv_key, key)
  return d.space[clust_key][fv_key][key] / d.space[clust_key].vector_sum[key]
end

function cond_prob_f_given_clust(d::Distribution{DataSet}, clust_key, key)
  return d.space[clust_key].vector_sum[key] / d.space[clust_key].vector_sum.total
end

function cond_prob_fv_given_clust(d::Distribution{DataSet}, clust_key, key)
  return d.space[clust_key][key].total / d.space[clust_key].vector_sum.total
end

# new scannell stuff
function prob_clust_in_dataset(d::Distribution{DataSet}, clust_key)
  return length(d.space[clust_key].vectors) / get_num_texts_in_dist(d)
end

function prob_clust_given_feature(d::Distribution{DataSet}, clust_key, feature)
  num_fvs = 0
  prob = 0
  for fv in values(d.space.clusters[clust_key])
    if haskey(fv,feature)
      num_fvs += 1
    end
  end
  prob = num_fvs / get_num_texts_given_feature(d,feature)
  return prob
end

# end scannell stuff

function keys(d::Distribution)
  return keys(d.space)
end

function features(d::Distribution{FeatureVector})
  return keys(d.space)
end

function features(d::Distribution)
  return keys(d.space.vector_sum)
end

function isempty(d::Distribution)
  return isempty(d.space)
end

function entropy(d::Distribution)
  ent = 0
  for feature in features(d)
    ent -= probability(d,feature)*log2(probability(d,feature))
  end
  return ent
end

function fv_entropy(d::Distribution{Cluster})
  ent = 0
  for fv in keys(d.space)
    ent -= prob_fv(d,fv)*log2(prob_fv(d,fv))
  end
  return ent
end

function fv_entropy(d::Distribution{DataSet})
  ent = 0
  for clust in keys(d.space)
    for fv in keys(d.space[clust])
      ent -= prob_fv(d,clust,fv)*log2(prob_fv(d,clust,fv))
    end
  end
  return ent
end

function fv_entropy(d::Distribution{DataSet}, clust_key)
  ent = 0
  for fv in keys(d.space[clust_key])
    ent -= prob_fv(d,clust_key,fv)*log2(prob_fv(d,clust_key,fv))
  end
  return ent
end

function clust_entropy(d::Distribution{DataSet})
  ent = 0
  for clust in keys(d.space)
    ent -= prob_cl(d,clust)*log2(prob_cl(d,clust))
  end
  return ent
end

function feature_entropy(d::Distribution{Cluster},feature)
  ent = 0
  for fv in keys(d.space.vectors)
    x = cond_prob_f_given_fv(d,fv,feature)*log2(cond_prob_f_given_fv(d,fv,feature))
    if isnan(x)
      x = 0
    end
    ent -= x
  end
  return ent
end

function feature_entropy(d::Distribution{DataSet},feature)
  ent = 0
  for clust in keys(d.space.clusters)
    x = cond_prob_f_given_clust(d,clust,feature)*log2(cond_prob_f_given_clust(d,clust,feature))
    if isnan(x)
      x = 0
    end
    ent -= x
  end
  return ent
end

# scannell
function clust_in_dataset_entropy(d::Distribution{DataSet})
  ent = 0
  for clust in keys(d.space)
    x = prob_clust_in_dataset(d,clust)*log2(prob_clust_in_dataset(d,clust))
    if isnan(x)
      x = 0
    end
    ent -= x
  end
  return ent
end

function clust_given_feature_entropy(d::Distribution{DataSet},feature)
  ent = 0
  for clust in keys(d.space)
    x = prob_clust_given_feature(d,clust,feature)*log2(prob_clust_given_feature(d,clust,feature))
    if isnan(x)
      x = 0
    end
    ent -= x
  end
  return ent
end

function clust_info_gain(d1::Distribution{DataSet}, feature, weight::Bool=false)
  ent_of_dist = clust_in_dataset_entropy(d1)
  ent_of_word = clust_given_feature_entropy(d1, feature)
  if weight
    return (ent_of_dist - ent_of_word)/ent_of_dist
  end
  return ent_of_dist - ent_of_word
end

function ld_info_gain(d::Distribution{DataSet},feature)
  all_ig = 0
  for clust in keys(d.space.clusters)
    child_ig = 0
    f_in_clust_ent = 0
    for fv in keys(d.space[clust].vectors)
      prob = cond_prob_f_given_fv_in_clust(d,clust,fv,feature)
      x = (prob*log2(prob))
      if isnan(x)
        x = 0
      end
      f_in_clust_ent -= x
    end
    child_ig = clust_in_dataset_entropy(d) - f_in_clust_ent
    all_ig += (child_ig*prob_clust_given_feature(d,clust,feature)/log2(length(d.space[clust].vectors)))
  end
  return (clust_info_gain(d,feature)/log2(length(d.space.clusters))) - all_ig
end
# end scannell

function info_gain(d1::Distribution, d2::Distribution)
  return entropy(d1)-entropy(d2)
end

function feature_info_gain(d1::Distribution, word)
  return entropy(d1)-feature_entropy(d1, word)
end

function perplexity(d::Distribution)
  return 2^entropy(d)
end

function display(dist::Distribution)
  display(dist.space)
end

#helper function that sets the smoothing type
function set_smooth!(d::Distribution{FeatureVector}, f::Function, sd::Array)
  d.smooth = f
  d.smooth_data = sd
end

#no smoothing default
function remove_smoothing!(d::Distribution)
  set_smooth!(d,_no_smoothing,[])
end

function _no_smoothing(d::Distribution{FeatureVector}, key, data::Array)
  return d.space[key] / d.total
end

function _no_smoothing(d::Distribution, feature, data::Array)
  return d.space.vector_sum[feature] / d.total
end

#add-delta smoothing, defaults to add-one smoothing
function delta_smoothing!(d::Distribution, δ::Number=1)
  if δ <= 0
    Base.warn("δ must be greater than 0")
  end
  set_smooth!(d,_δ_smoothing,[δ,d.features,d.total])
end

function _δ_smoothing(d::Distribution{FeatureVector}, key, data::Array)
  if !haskey(d.space, key)
    return (data[1])/(data[1]*(data[2]+1)+data[3])
  end
  return (d.space[key]+data[1])/(data[1]*(data[2]+1)+data[3])
end

function _δ_smoothing(d::Distribution, key, data::Array)
  if !haskey(d.space.vector_sum, key)
    return (data[1])/(data[1]*(data[2]+1)+data[3])
  end
  return (d.space[key]+data[1])/(data[1]*(data[2]+1)+data[3])
end

#simple good-turing smoothing
function goodturing_smoothing!(d::Distribution{FeatureVector})
  freqs = FeatureVector()
  for value in values(d.space)
    freqs[value] += 1
  end
  set_smooth!(d,_gt_smoothing, [d.total, freqs])
end

function goodturing_smoothing!(d::Distribution)
  freqs = FeatureVector()
  for value in values(d.space.vector_sum)
    freqs[value] += 1
  end
  set_smooth!(d,_gt_smoothing, [d.total, freqs])
end

function _gt_smoothing(d::Distribution{FeatureVector}, key, data::Array)
  if !haskey(d.space, key)
    return data[2][1] / data[1] #num of keys that occur once / total number of keys
  end
  c = d.space[key]
  c_adjust = (c+1) * (data[2][c+1]/data[2][c])
  return c_adjust / data[1]
end

function _gt_smoothing(d::Distribution, key, data::Array)
  if !haskey(d.space.vector_sum, key)
    return data[2][1] / data[1] #num of keys that occur once / total number of keys
  end
  c = d.space.vector_sum[key]
  c_adjust = (c+1) * (data[2][c+1]/data[2][c])
  return c_adjust / data[1]
end


text:
using RecSys

function test(matfile, entryname)
    ratings = MatFile(matfile, entryname)
    rec = ALSWR(ratings, ParShmem())
    train(rec, 10, 4)
    err = rmse(rec)
    println("rmse of the model: $err")

    recommended, watched, nexcl = recommend(rec, 100)
    isempty(watched) || println("Already watched: $watched")
    (nexcl == 0) || println("Excluded $(nexcl) movies already watched")
    println("Recommended: $recommended")
    nothing
end

function test_chunks(dataset_path, model_path)
    mem = Base.Sys.free_memory()
    mem_model = mem_inputs = round(Int, mem/3)
    user_item_ratings = SparseBlobs(joinpath(dataset_path, "splits", "R_itemwise"); maxcache=mem_model)
    item_user_ratings = SparseBlobs(joinpath(dataset_path, "splits", "RT_userwise"); maxcache=mem_model)

    rec = ALSWR(user_item_ratings, item_user_ratings, ParBlob())

    train(rec, 10, 4, model_path, mem_inputs)

    err = rmse(rec)
    println("rmse of the model: $err")
    nothing
end


text:
#lowerBoundsSymmetric

#computes the maximum over the geometric grid
#sep_1(a, b, l, u) solves  min_{t in [l, u]}  at + 2bt H(t)
#sep_2(a, l, u) solves min_{t in [l, u]} a H(t)
#user must ensure sep_1, sep_2 and h are consistent.
#H(t) = 1/2t int_[-t, t] h(1+s)ds
function _rev_p_lb_symmetric(S, h, sep_1, sep_2, H, H_un, pj, numCuts, print_trace, TOL)
	m = Model(solver=GurobiSolver(OutputFlag=false))
	@variable(m, theta)
	@variable(m, lam)
	@objective(m, Min, theta)

	#add endpoints to get things started
	#t = 0
	ind = pj <= 1 ? pj : 0
	@constraint(m, lam * H(0) >= ind - theta)

	#t = abs(pj -1)
	tstar = abs(pj -1)
	@constraint(m, lam * H(tstar) >= ind - theta)
	@constraint(m, (2theta - pj) * tstar + lam * H_un(tstar) >= pj * (1 - pj))

	#t = S-1
	@constraint(m, (2theta - pj) * (S - 1) + lam * H_un(S - 1) >= pj * (1 - pj))

	for iter = 1:numCuts
		status = solve(m)
		feasible = true

		print_trace && println("Iteration: \t", iter)
		#First constraint [0, abs(pj -1)]
		rhs =  ind - getvalue(theta)
		val, tstar = sep_2(getvalue(lam), 0, abs(pj -1))
		print_trace && println("Iter: $iter Val $val RHS $rhs  : tstar $tstar")
		if val < rhs - TOL
			@constraint(m, lam * H(tstar) >= ind - theta)
			feasible = false
		end

		#2nd constraint [abs(pj-1) S-1]
		rhs = pj * (1 - pj)
		val, tstar = sep_1(2 * getvalue(theta) - pj, getvalue(lam), abs(pj-1), S-1)

		print_trace && println("Iter: $iter Val $val RHS $rhs : tstar $tstar")
		if val < rhs - TOL
			cnst = @constraint(m, (2theta - pj) * tstar + lam * H_un(tstar) >= rhs)
			feasible = false
		end

		if feasible
			if status != :Optimal
				println("Problem likely infeasible.  Dual Status:\t", status)
				return -1.0
			else #correctly solved
				return getobjectivevalue(m)
			end
		end
	end
	##Exhausted iteration count
	throw("Max Iteration Count Reached:  h(1): $(h(1))")
	return -1.
end

function vopp_lb_symmetric_MAD(mu, S, M, D; delta = S/100, numCuts=100, print_trace=false, TOL=1e-6, pj=-1, safe_fail=false)
	#standardize
	c = mu * (1 - M)
	Sc = vopp.comp_Sc(S, M)

	#round down Sc if it's too large and throw otherwise
	Sc = min(Sc, 2)
	@assert Sc > 1 "Scale $S is too small for any distribution"

	if D/M > .25
		println("Deviation $(D/M) > .25, problem likely infeasible")
		if safe_fail
			println("Safe_fail set to true. Aborting Calc")
			return -1.0
		end
	end

	#create the functions/separators closed-form
	h(t) = M/2 * abs(t-1) - D
	H(t) = M * t / 4 - D
	H_un(t) = 2t * H(t)

	#sep_1(a, b, l, u) solves  min_{t in [l, u]}  at + b * H_un(t)
	#min_{t in [l, u]} (a - 2bD) t + bM/2 * t^2
	function sep_1(a, b, l, u)
		f(t) = (a - 2b * D) * t + b * M/2 * t^2
		valstar, tstar = f(l), l
		if f(u) < valstar
			valstar, tstar = f(u), u
		end
		if b > 0
			tmid = (2b * D - a) / b / M
			valmid = Inf
			if l <= tmid <= u
				valmid = f(tmid)
			end
			if valmid < valstar
				valstar, tstar = valmid, tmid
			end
		end
		valstar, tstar
	end

	#solves min_{t in [l, u]} a H(t)
	function sep_2(a, l, u)
		if a * H(u) < a * H(l)
			return a * H(u), u
		end
		return a * H(l), l
	end

	#search over the geometric price ladder
	if pj > 0
		ps = [pj]
	else
		ps = geom_price_ladder(Sc, delta)
	end
	rstar = -Inf
	pstar = 0.
	rj = 0.

	for pj in ps
		#kick i tto the workhorse
		rj = _rev_p_lb_symmetric(Sc, h, sep_1, sep_2, H, H_un, pj, numCuts, print_trace, TOL)
		print_trace && println("\n Price Pj:\t", pj, "\t", "Rj:\t", rj, "\n")

		if rj > rstar
			rstar = rj
			pstar = pj
		end
	end

	return 1/rstar, pstar
end

function vopp_lb_symmetric_CV(mu, S, M, C; delta = S/100, numCuts=100, print_trace=false, TOL=1e-6, pj=-1, safe_fail=false)
	#standardize
	c = mu * (1 - M)
	Sc = vopp.comp_Sc(S, M)

	#round down Sc if it's too large and throw otherwise
	Sc = min(Sc, 2)
	@assert Sc > 1 "Scale $S is too small for any distribution"

	max_cv = max_std_cv_unimodal_guess(1, Sc, M, true)
	if  C / M > max_cv
		println("CV $(C/M) > $max_cv, so problem likely infeasible")
		if safe_fail
			println("Safe_fail set to true. Aborting Calc")
			return -1.0
		end
	end

	#create the functions/separators closed-form
	h(t) = M^2 * (t-1)^2 - C^2
	H(t) = M^2 / 3 * t^2 - C^2
	H_un(t) = 2t * H(t)

	#sep_1(a, b, l, u) solves
	#min_{t in [l, u]}  at + b * H_un(t)
	function sep_1(a, b, l, u)
		f(t) = a * t + 2 * b * M^2 / 3 * t^3 - 2b * t * C^2
		valstar, tstar = f(l), l
		if f(u) < valstar
			valstar, tstar = f(u), u
		end
		if (2b * C^2 - a) / (2b * M^2) > 0
			tmid = sqrt((2b * C^2 - a) / (2b * M^2))
			valmid = Inf
			if l <= tmid <= u
				valmid = f(tmid)
			end
			if valmid < valstar
				valstar, tstar = valmid, tmid
			end
		end
		valstar, tstar
	end

	#solves min_{t in [l, u]} a H(t)
	function sep_2(a, l, u)
		f(t) = a * H(t)
		tstar, valstar = l, f(l)
		if f(u) < valstar
			tstar, valstar = u, f(u)
		end
		if l <= 0. <= u
			if f(0) < valstar
				tstar, valstar = 0., f(0.)
			end
		end
		return valstar, tstar
	end

	#search over the geometric price ladder
	if pj > 0
		ps = [pj]
	else
		ps = geom_price_ladder(Sc, delta)
	end
	rstar = -Inf
	pstar = 0.
	rj = 0.

	for pj in ps
		#kick i tto the workhorse
		rj = _rev_p_lb_symmetric(Sc, h, sep_1, sep_2, H, H_un, pj, numCuts, print_trace, TOL)
		print_trace && println("\n Price Pj:\t", pj, "\t", "Rj:\t", rj, "\n")

		if rj > rstar
			rstar = rj
			pstar = pj
		end
	end

	return 1/rstar, pstar
end


text:
function ueberpruefeMatrix(Stelle,Matrixeintrag,PixelNummer,Abstand)

laengenmatrix = zeros(size(PixelNummer,2),Int(maximum(maximum(PixelNummer))));

if Stelle(4,5) == PixelNummer[4,5]
    if Matrixeintrag(5,4) ==  Abstand[5,4]
        for S = 1:size(PixelNummer,2)
            for Z = 1:size(PixelNummer,1)
                if Stelle(Z,S) > 0
                    wo = Int(Stelle(Z,S));
                    laengenmatrix[S,wo] = Matrixeintrag(Z,S);
                end
            end
        end
        println("Eure Formel für die Stelle an die der Wert geschrieben ist, ist korrekt.")
    else
        println("Schaut noch einmal eure Formel für die Bestimmung des Matrixeintrags an.")
    end
else
    if Matrixeintrag(5,4) ==  Abstand[5,4]
        println("Eure Formel für den Matrixeintrag, der in die Matrix geschrieben werden soll ist korrekt.")
        println("Schaut euch noch einmal eure Formel für die Stelle an, an die der Wert geschrieben werden soll.")
    else
        println("Beide eure Formeln sind noch nicht korrekt, überprüft eure Eingaben oder überlegt noch einmal neu.")
    end
end
return laengenmatrix;
end


text:
const CONSTANT_SYMBOLS = (:nothing, :Float64, :Float32, :Int8, :UInt8, :Int16, :UInt16, :Int32, :UInt32, :Int64, :UInt64)
function add_constant!(ls::LoopSet, var::Symbol, elementbytes::Int)
  globalconst = Base.sym_in(var, CONSTANT_SYMBOLS)
  instr = globalconst ? Instruction(GLOBALCONSTANT, var) : LOOPCONSTANT
  op = Operation(length(operations(ls)), var, elementbytes, instr, constant, NODEPENDENCY, Symbol[], NOPARENTS)
  rop = pushop!(ls, op, var)
  (!globalconst && (rop === op)) && pushpreamble!(ls, op, var)
  rop
end
# function add_constant!(ls::LoopSet, var, elementbytes::Int = 8)
#     sym = gensym(:loopconstant)
#     pushpreamble!(ls, Expr(:(=), sym, var))
#     add_constant!(ls, sym, elementbytes)
# end
function add_constant!(ls::LoopSet, var::Number, elementbytes::Int = 8, varname = gensym!(ls, "loopconstnumber"))
  op = Operation(length(operations(ls)), varname, elementbytes, LOOPCONSTANT, constant, NODEPENDENCY, Symbol[], NOPARENTS)
  ops = operations(ls)
  typ = var isa Integer ? HardInt : HardFloat
  if iszero(var)
    for (id,typ_) ∈ ls.preamble_zeros
      (instruction(ops[id]) == LOOPCONSTANT && typ == typ_) && return ops[id]
    end
    push!(ls.preamble_zeros, (identifier(op),typ))
  elseif var isa Integer
    idescript = integer_description(var)
    for (id,descript) ∈ ls.preamble_symint
      if (instruction(ops[id]) == LOOPCONSTANT) && (idescript == descript)
        return ops[id]
      end
    end
    push!(ls.preamble_symint, (identifier(op), idescript))
  else#if var isa FloatX
    for (id,fvar) ∈ ls.preamble_symfloat
      (instruction(ops[id]) == LOOPCONSTANT && fvar == var) && return ops[id]
    end
    push!(ls.preamble_symfloat, (identifier(op), var))
  end
  rop = pushop!(ls, op)
  rop === op || return rop
  pushpreamble!(ls, Expr(:(=), name(op), var))
  rop
end
function ensure_constant_lowered!(ls::LoopSet, op::Operation)
  if iscompute(op)
    call = callexpr(instruction(op))
    for opp ∈ parents(op)
      ensure_constant_lowered!(ls, opp)
      push!(call.args, name(opp))
    end
    pushpreamble!(ls, Expr(:(=), name(op), call))
  elseif isconstant(op) & !isconstantop(op)
    opid = identifier(op)
    for (id, sym) ∈ ls.preamble_symsym
      if id == opid
        pushpreamble!(ls, Expr(:(=), name(op), sym))
        return
      end
    end
    for (id,(intval,intsz,signed)) ∈ ls.preamble_symint
      if id == opid
        if intsz == 1
          pushpreamble!(ls, Expr(:(=), name(op), intval % Bool))
        elseif signed
          pushpreamble!(ls, Expr(:(=), name(op), intval))
        else
          pushpreamble!(ls, Expr(:(=), name(op), intval % UInt))
        end
        return
      end
    end
    for (id,floatval) ∈ ls.preamble_symfloat
      if id == opid
        pushpreamble!(ls, Expr(:(=), name(op), floatval))
        return
      end

    end
    for (id,typ) ∈ ls.preamble_zeros
      if id == opid
        pushpreamble!(ls, Expr(:(=), name(op), staticexpr(0)))
        return
      end
    end
    for (id,f) ∈ ls.preamble_funcofeltypes
      if id == opid
        pushpreamble!(ls, Expr(:(=), name(op), Expr(:call, reduction_zero(f), Float64)))
        return
      end
    end
  end
end
function ensure_constant_lowered!(ls::LoopSet, mpref::ArrayReferenceMetaPosition, ind::Symbol)
  length(loopdependencies(mpref)) == 0 && return
  for (id,opp) ∈ enumerate(parents(mpref))
    if name(opp) === ind
      ensure_constant_lowered!(ls, opp)
    end
  end
  return nothing
end
function add_constant!(ls::LoopSet, mpref::ArrayReferenceMetaPosition, elementbytes::Int)
  op = Operation(length(operations(ls)), varname(mpref), elementbytes, LOOPCONSTANT, constant, NODEPENDENCY, Symbol[], NOPARENTS, mpref.mref)
  add_vptr!(ls, op)
  temp = gensym!(ls, "intermediateconstref")
  vloadcall = Expr(:call, lv(:_vload), mpref.mref.ptr)
  nindices = length(getindices(op))
  # getoffsets(op) .+= 1
  if nindices > 0
    dummyloop = first(ls.loops)
    for ind ∈ getindicesonly(op)
      ensure_constant_lowered!(ls, mpref, ind)
    end
    push!(vloadcall.args, mem_offset(op, UnrollArgs(dummyloop, dummyloop, dummyloop, 0, 0, 0), fill(false,nindices), true, ls))
  end
  push!(vloadcall.args, Expr(:call, lv(:False)), staticexpr(reg_size(ls)))
  pushpreamble!(ls, Expr(:(=), temp, vloadcall))
  pushpreamble!(ls, Expr(:(=), name(op), temp))
  pushpreamble!(ls, op, temp)
  pushop!(ls, op, temp)
end
# This version has loop dependencies. var gets assigned to sym when lowering.
# value is what will get assigned within the loop.
# assignedsym will be assigned to value within the preamble
function add_constant!(
    ls::LoopSet, value::Symbol, deps::Vector{Symbol}, assignedsym::Symbol, elementbytes::Int, f::Symbol = Symbol("")
  )
  retop = get(ls.opdict, value, nothing)
  if retop === nothing
    op = Operation(length(operations(ls)), assignedsym, elementbytes, Instruction(f, value), constant, deps, NODEPENDENCY, NOPARENTS)
  else
    op = Operation(length(operations(ls)), assignedsym, elementbytes, :identity, compute, deps, reduceddependencies(retop), [retop])
  end
  pushop!(ls, op, assignedsym)
end
# function add_constant!(
#     ls::LoopSet, value, deps::Vector{Symbol}, assignedsym::Symbol, elementbytes::Int, f::Symbol = Symbol("")
# )
#     intermediary = gensym(:intermediate) # hack, passing meta info here
#     pushpreamble!(ls, Expr(:(=), intermediary, value))
#     add_constant!(ls, intermediary, deps, assignedsym, f, elementbytes)
# end
function add_constant!(
    ls::LoopSet, value::Number, deps::Vector{Symbol}, assignedsym::Symbol, elementbytes::Int
)
  op = add_constant!(ls, gensym!(ls, string(value)), deps, assignedsym, elementbytes, :numericconstant)
  pushpreamble!(ls, op, value)
  op
end


text:
using DifferentialEquations, DataFrames
using Plots
using CSV
include("odeFileExample.jl")
include("variableNames.jl")
include("scanIncludes.jl")

f=DDEFunction(ddeFile,syms=Symbol.(syms))
#f=ODEFunction(odeFile,syms=Symbol.(syms))
y0=ones(length(syms))
p=1
h(p,t)=ones(length(syms))
params=1
prob=DDEProblem(f,ones(length(syms)),h,(0.0,10.0),syms=Symbol.(syms))
#prob=ODEProblem(f,ones(length(syms)),(0.0,10.0))
print("Solving Equations \n")
@timev sol=solve(prob)
plot(sol)

modify["k_ABBinding"]=1.5

prob=DDEProblem(f,ones(length(syms)),h,(0.0,10.0),syms=Symbol.(syms))
print("Solving Equations \n")
@timev sol=solve(prob)
plot!(sol)
png("testOutScan.png")


text:
using Distributions
using EmpiricalLikelihood
using Ipopt
using JuMP
using Random
using Test
using Statistics

het_moment(X, y, beta) = [y[1] - beta[1] - X[1] * beta[2]]

@testset "Heteroskedasticity" begin
    n = 50
    Random.seed!(213)
    β = [1., 1.]
    β_test = zeros(100, 2)
    for rep in 1:100
        X = rand(LogNormal(), n)
        ϵ = rand(Normal(), n)
        u = ϵ .* @. sqrt(.1 .+ .2 * X .+ .3 * X.^2)
        y = β[1] .+ X * β[2] .+ u
        β_test[rep, :] = optimize_el(1, het_moment, reshape(X, n, 1), reshape(y, n, 1), β, .35)[1]
    end
    @test mean(β_test, dims=1)[1, :] ≈ β atol=1e-1
end

@testset "Test JuMP interface" begin
    n = 50
    Random.seed!(213)
    β = [1., 1.]
    β_test = zeros(100, 2)
    for rep in 1:100
        X = rand(LogNormal(), n)
        ϵ = rand(Normal(), n)
        u = ϵ .* @. sqrt(.1 .+ .2 * X .+ .3 * X.^2)
        y = β[1] .+ X * β[2] .+ u
        model = setup_jump_problem(1, het_moment, reshape(X, n, 1), reshape(y, n, 1), β, .35)[1]
        set_optimizer(model, optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
        optimize!(model)
        β_test[rep, :] = @. value(model[:theta])
    end
    @test mean(β_test, dims=1)[1, :] ≈ β atol=1e-1
end


text:
import DataFrames
import Dates
import TimeSeries

function verify_time_series(sys::System, num_initial_times, num_time_series, len)
    total_time_series = 0
    all_time_series = get_time_series_multiple(sys)
    for time_series in all_time_series
        if length(time_series) != len
            @error "length doesn't match" length(time_series) len
            return false
        end
        total_time_series += 1
    end

    if num_time_series != total_time_series
        @error "num_time_series doesn't match" num_time_series total_time_series
        return false
    end

    return true
end

@testset "Test read_time_series_file_metadata" begin
    for file in ["timeseries_pointers.json", "timeseries_pointers.csv"]
        filename = joinpath(RTS_GMLC_DIR, file)
        all_time_series = IS.read_time_series_file_metadata(filename)
        @test length(all_time_series) == 260

        for time_series in all_time_series
            @test isfile(time_series.data_file)
        end
    end
end

@testset "Test time_series normalization" begin
    component_name = "122_HYDRO_1"
    timeseries_file = joinpath(
        DATA_DIR,
        "forecasts",
        "RTS_GMLC_forecasts",
        "gen",
        "Hydro",
        "DAY_AHEAD_hydro.csv",
    )
    gen = ThermalStandard(nothing)
    gen.name = component_name
    resolution = Dates.Hour(1)

    # Parse the file directly in order to compare values.
    ts_base = SingleTimeSeries(component_name, timeseries_file, gen, resolution)
    timeseries = get_data(ts_base)
    max_value = maximum(TimeSeries.values(timeseries))

    file_metadata = IS.TimeSeriesFileMetadata(
        simulation = "DAY_AHEAD",
        category = "Generator",
        component_name = "122_HYDRO_1",
        name = "active_power",
        normalization_factor = 1.0,
        data_file = timeseries_file,
        percentiles = [],
        resolution = resolution,
        time_series_type_module = "InfrastructureSystems",
        time_series_type = "SingleTimeSeries",
    )

    # Test code path where no normalization occurs.
    sys = PSB.build_system(PSB.MatPowerTestSystems, "matpower_RTS_GMLC_sys")
    add_time_series!(sys, [file_metadata])
    verify_time_series(sys, 1, 1, 24)
    time_series = collect(get_time_series_multiple(sys))[1]
    @test TimeSeries.values(time_series.data) == TimeSeries.values(timeseries)

    # Test code path where timeseries is normalized by dividing by the max value.
    file_metadata.normalization_factor = "Max"
    sys = PSB.build_system(PSB.MatPowerTestSystems, "matpower_RTS_GMLC_sys")
    add_time_series!(sys, [file_metadata])
    verify_time_series(sys, 1, 1, 24)
    time_series = collect(get_time_series_multiple(sys))[1]
    @test TimeSeries.values(time_series.data) == TimeSeries.values(timeseries ./ max_value)

    # Test code path where timeseries is normalized by dividing by a custom value.
    nf = 95.0
    file_metadata.normalization_factor = nf
    sys = PSB.build_system(PSB.MatPowerTestSystems, "matpower_RTS_GMLC_sys")
    add_time_series!(sys, [file_metadata])
    verify_time_series(sys, 1, 1, 24)
    time_series = collect(get_time_series_multiple(sys))[1]
    @test TimeSeries.values(time_series.data) == TimeSeries.values(timeseries ./ nf)
end

@testset "Test single time_series addition" begin
    component_name = "122_HYDRO_1"
    name = "active_power"
    timeseries_file = joinpath(
        DATA_DIR,
        "forecasts",
        "RTS_GMLC_forecasts",
        "gen",
        "Hydro",
        "DAY_AHEAD_hydro.csv",
    )
    resolution = Dates.Hour(1)

    # Test with a filename.
    sys = PSB.build_system(PSB.MatPowerTestSystems, "matpower_RTS_GMLC_sys")
    component = get_component(HydroDispatch, sys, component_name)
    ts = SingleTimeSeries(
        name,
        timeseries_file,
        component,
        resolution;
        normalization_factor = 1.0,
    )
    ta = get_data(ts)
    add_time_series!(sys, component, ts)
    verify_time_series(sys, 1, 1, 24)
    time_series = collect(get_time_series_multiple(sys))[1]
    @test TimeSeries.timestamp(get_data(time_series)) == TimeSeries.timestamp(ta)
    @test TimeSeries.values(get_data(time_series)) == TimeSeries.values(ta)

    # Test with TimeSeries.TimeArray.
    sys = PSB.build_system(PSB.MatPowerTestSystems, "matpower_RTS_GMLC_sys")
    component = get_component(HydroDispatch, sys, component_name)
    ts = SingleTimeSeries(name, ta; normalization_factor = 1.0)
    add_time_series!(sys, component, ts)
    verify_time_series(sys, 1, 1, 24)
    time_series = collect(get_time_series_multiple(sys))[1]
    @test TimeSeries.values(get_data(time_series)) == TimeSeries.values(ta)

    # Test with DataFrames.DataFrame.
    sys = PSB.build_system(PSB.MatPowerTestSystems, "matpower_RTS_GMLC_sys")
    component = get_component(HydroDispatch, sys, component_name)
    df = DataFrames.DataFrame(ta)
    ts = SingleTimeSeries(name, df; normalization_factor = 1.0)
    add_time_series!(sys, component, ts)
    verify_time_series(sys, 1, 1, 24)
    time_series = collect(get_time_series_multiple(sys))[1]
end

@testset "TimeSeriesData data matpower" begin
    sys = PSB.build_system(PSB.MatPowerTestSystems, "matpower_case5_re_sys")
    file_metadata = joinpath(TIME_SERIES_DIR, "5bus_ts", "timeseries_pointers_da.json")
    add_time_series!(sys, file_metadata)
    @test verify_time_series(sys, 1, 5, 24)

    # Add the same files.
    # This will fail because the component-name pairs will be duplicated.
    @test_throws ArgumentError add_time_series!(sys, file_metadata)

    file_metadata = joinpath(TIME_SERIES_DIR, "5bus_ts", "timeseries_pointers_rt.json")

    ## This will fail because the resolutions are different.
    @test_throws IS.ConflictingInputsError add_time_series!(sys, file_metadata)

    ## TODO: need a dataset with same resolution but different horizon.

    # sys = System(PowerSystems.PowerModelsData(joinpath(MATPOWER_DIR, "case5_re.m")))
    sys = PSB.build_system(PSB.MatPowerTestSystems, "matpower_case5_re_sys")
    add_time_series!(sys, file_metadata)
    @test verify_time_series(sys, 1, 5, 288)
end


text:
module TrackerFlux

using Flux, Tracker

track(m) = fmap(x -> x isa AbstractArray ? Tracker.param(x) : x, m)
untrack(m) = fmap(Tracker.data, m)

function Flux.Optimise.update!(opt, xs::Tracker.Params, gs)
    for x in xs
        Flux.Optimise.update!(opt, x, gs[x])
    end
end
function Flux.Optimise.update!(opt, x, x̄)
    Tracker.update!(x, -Flux.Optimise.apply!(opt, Tracker.data(x), Tracker.data(x̄)))
end

_truncate(x::AbstractArray) = Tracker.data(x)
_truncate(x::Tuple) = _truncate.(x)
truncate!(m::Flux.Recur) = (m.state = _truncate(m.state))
truncate!(m) = foreach(truncate!, Flux.functor(m)[1])

function overload_gradient()
    @eval Flux.gradient(f, args...) = Tracker.gradient(f, args...)
end

function __init__()
    @eval Flux truncate! = $truncate!
end

end


text:
module YaoTargetQobj

export MemoryInfo, QobjQASMTarget, code_qobj, measure_ssa_uses, measure_ssa_returns

using UUIDs
using MLStyle
using IBMQClient
using YaoCompiler
using Configurations
using IBMQClient.Schema
using YaoCompiler.Intrinsics
using CompilerPluginTools
using YaoCompiler: YaoCompileTarget
using YaoLocations: plain

@option struct QobjQASMTarget <: YaoCompileTarget
    nshots::Int = 1024
    seed::Int = 1
    max_credits::Int = 3
    qobj_id::String = string(uuid1())
    schema_version::VersionNumber = v"1.3"
    memory_offset::Int = 0
    description::Maybe{String} = nothing
end

function code_qobj(op::Type{<:Operation}; kw...)
    description = get(kw, :description, string(op))
    target = QobjQASMTarget(;kw..., description)
    return YaoCompiler.compile(target, Intrinsics.apply, Tuple{AnyReg, op})
end

# TODO: move to YaoCompiler
function push_stmt!(list::Vector, stmt)
    isnothing(stmt) && return list
    if stmt isa Vector
        append!(list, stmt)
    else
        push!(list, stmt)
    end
end

include("validate.jl")
include("codegen.jl")

end