text:
struct Page
    session::Session
    child_sessions::Dict{String, Session}

    # Inlines all CSS & JS dependencies, and doesn't use websocket connection
    # to load initial data. This makes JSServe a bit slower, but allows you do
    # export the Page to static HTML e.g. in Pluto
    exportable::Bool
    # Doesn't even try to connect to a julia process
    offline::Bool
end

const CURRENT_PAGE = Ref{Page}()

"""
    Page(;
        session=nothing,
        exportable=false,
        offline=false,
        server_config...
    )

A Page should be used for anything that displays multiple outputs, like Pluto/IJulia/Documenter.
It activates a special html mime show mode, which is more efficient in that scenario.

A Page creates a single entry point, to connect to the Julia process and load dependencies.
For Documenter, the page needs to be set to `exportable=true, offline=true`.
Exportable has the effect of inlining all data & js dependencies, so that everything can be loaded in a single HTML object.
`offline=true` will make the Page not even try to connect to a running Julia
process, which makes sense for the kind of static export we do in Documenter.

For convenience, one can also pass additional server configurations, which will directly
get put into `configure_server!(;server_config...)`.
Have a look at the docs for `configure_server!` to see the parameters.
"""
function Page(;
        session=nothing,
        exportable=false,
        offline=false,
        server_config...
    )
    configure_server!(;server_config...)
    if session === nothing
        serializer = UrlSerializer(inline_all=exportable)
        session = Session(url_serializer=serializer)
    end

    p = Page(
        session, Dict{String, Session}(),
        exportable,
        offline,
    )
    CURRENT_PAGE[] = p
    return p
end

function Base.close(page::Page)
    for (k, s) in page.child_sessions
        close(s)
    end
    empty!(page.child_sessions)
    close(page.session)
    empty!(page.session.unique_object_cache)
end

function include_all_assets(session, assets)
    # protect against require being loaded by someone else
    # (e.g. Documenter.jl)
    return [
        jsrender(session, js"""
            window.__define = window.define;
            window.__require = window.require;
            window.define = undefined;
            window.require = undefined;
        """),
        include_asset.(assets, (session.url_serializer,))...,
        jsrender(session, js"""
            window.define = window.__define;
            window.require = window.__require;
            window.__define = undefined;
            window.__require = undefined;
        """)
    ]
end

function Base.show(io::IO, ::MIME"text/html", page::Page)
    if !(page.offline)
        server = get_server()
        insert_session!(server, page.session)
    end
    serializer = page.session.url_serializer
    websocket_url = JSSERVE_CONFIGURATION.external_url[]
    delete_session = Observable("")

    push!(page.session, page.session.js_comm)

    on(delete_session) do session_id
        if !haskey(page.child_sessions, session_id)
            error("Session to delete not found ($(session_id)), please open an issue with JSServe.jl")
        end

        @debug("deleting child session: $(session_id)")
        child = page.child_sessions[session_id]
        # Set connection to nothing,
        # so that no clean up happens over the active websocket connection
        # (e.g. close will try to close the ws connection)
        child.on_close[] = true
        child.connection[] = nothing
        close(child)
        delete!(page.child_sessions, session_id)
    end

    init = if page.offline
        js"JSServe.setup_connection({offline: true})"
    else
        js"""
            const proxy_url = $(websocket_url)
            const session_id = $(page.session.id)
            // track if our child session doms get removed from the document (dom)
            JSServe.track_deleted_sessions($(delete_session))
            JSServe.setup_connection({proxy_url, session_id})
            JSServe.sent_done_loading()
        """
    end
    deps = [
        MsgPackLib,
        PakoLib,
        Base64Lib,
        JSServeLib,
    ]
    page_init_dom = DOM.div(
        include_all_assets(page.session, deps)...,
        init
    )
    println(io, node_html(page.session, page_init_dom))
end

function assure_ready(page::Page)
    # Nothing to wait for when offline
    page.offline && return
    session = page.session
    filter!(page.child_sessions) do (id, session)
        never_opened = isnothing(session.connection[])
        if never_opened
            @warn("Removing unopened Session: $(id)")
            close(session)
        end
        return !never_opened
    end
    tstart = time()
    warned = false
    while time() - tstart < 10
        children_loaded = all(isready, values(page.child_sessions))
        isready(session) && children_loaded && break
        if (time() - tstart) > 1 && !warned
            warned = true
            @warn("Waiting for page sessions to load.
                This can happen for the first cells to run, or is indicative of faulty state")
        end
    end
    send(session, fused_messages!(session))
end

function render_sub_session(parent_session, html_dom)

    session = Session(parent_session;
        connection=Base.RefValue{Union{Nothing, WebSocket, IOBuffer}}(nothing),
        dependencies=Set{Asset}(),
        id=string(uuid4()),
        observables=Dict{String, Tuple{Bool, Observable}}(),
        on_close=Observable(false),
        message_queue=Dict{Symbol, Any}[],
        js_fully_loaded=Channel{Bool}(1),
        deregister_callbacks=Observables.ObserverFunction[]
    )

    on_init = Observable(false)
    # Manually register on_init - this is a bit fragile, but
    # we need to normally register on_init with `session`, BUT
    # it already needs to be registered upfront with JSServe in the browser
    # so that it can properly trigger the observable early in the loaded html/js
    send(parent_session, payload=on_init[], id=on_init.id, msg_type=RegisterObservable)
    # Render the app and register all the resources with the session
    # Note, since we set the connection to nothing, nothing gets sent yet
    # This is important, since we can only sent the messages after the HTML has been rendered
    # since we must assume the messages / js contains references to HTML elements
    js_dom = jsrender(session, html_dom)
    register_resource!(session, js_dom)

    new_deps = setdiff(session.dependencies, parent_session.dependencies)
    union!(parent_session.dependencies, new_deps)

    on(session, on_init) do is_init
        if !is_init
            error("The html didn't initialize correctly")
        end
        init_session(session)
    end

    init = js"""
        // register this session so it gets deleted when it gets removed from dom
        JSServe.register_sub_session($(session.id))
        console.log("Initializing session!!!")
        JSServe.update_obs($(on_init), true)
    """

    final_dom = DOM.span(
        js_dom,
        include_all_assets(parent_session, new_deps)...,
        jsrender(session, init),
        id=session.id,
    )

    obs_shared_with_parent = intersect(keys(session.observables), keys(parent_session.observables))
    # Those obs are managed by page session, so we don't need to have them in here!
    # This is important when we clean them up, since a child session shouldn't
    # delete any observable already managed by the parent session
    for obs in obs_shared_with_parent
        delete!(session.observables, obs)
    end
    # but in the end, all observables need to be registered
    # with the page session, since that's where javascript will sent all the events
    merge!(parent_session.observables, session.observables)

    on(session, session.on_close) do closed
        # remove the sub-session specific js resources
        if closed
            obs_ids = collect(keys(session.observables))
            JSServeLib.delete_observables(parent_session, obs_ids)
            # since we deleted all obs shared with the parent session,
            # we can delete savely delete all of them from there
            # Note, that we previously added them all here!
            for (k, o) in session.observables
                delete!(parent_session.observables, k)
            end
        end
    end

    # finally give page a connection! :)
    session.connection[] = parent_session.connection[]

    return final_dom
end

function show_in_page(page::Page, app::App)
    page_session = page.session
    is_offline = page.offline
    # The page session must be ready to display anything!
    # Since the page display is rendered async in the browser and we have no
    # idea when it's done on the Julia side, so we need to wait here!
    assure_ready(page)

    # Create a child session, to track the per app resources
    session = Session(page_session;
        connection=Base.RefValue{Union{Nothing, WebSocket, IOBuffer}}(nothing),
        dependencies=Set{Asset}(),
        id=string(uuid4()),
        observables=Dict{String, Tuple{Bool, Observable}}(),
        on_close=Observable(false),
        message_queue=Dict{Symbol, Any}[],
        js_fully_loaded=Channel{Bool}(1),
        deregister_callbacks=Observables.ObserverFunction[]
    )
    # register with page session for proper clean up!
    page.child_sessions[session.id] = session

    on_init = Observable(false)
    # Manually register on_init - this is a bit fragile, but
    # we need to normally register on_init with `session`, BUT
    # it already needs to be registered upfront with JSServe in the browser
    # so that it can properly trigger the observable early in the loaded html/js
    send(page_session, payload=on_init[], id=on_init.id, msg_type=RegisterObservable)
    # Render the app and register all the resources with the session
    # Note, since we set the connection to nothing, nothing gets sent yet
    # This is important, since we can only sent the messages after the HTML has been rendered
    # since we must assume the messages / js contains references to HTML elements
    html_dom = Base.invokelatest(app.handler, session, (; show="/show_inline"))
    js_dom = jsrender(session, html_dom)
    register_resource!(session, js_dom)

    new_deps = setdiff(session.dependencies, page_session.dependencies)
    union!(page_session.dependencies, new_deps)

    exportable = page.exportable

    on(session, on_init) do is_init
        if !is_init
            error("The html didn't initialize correctly")
        end
        init_session(session)
    end

    init = if exportable
        # We take all messages and serialize them directly into the init js
        messages = fused_messages!(session)
        js"""(()=> {
            JSServe.register_sub_session($(session.id))
            const init_data_b64 = $(serialize_string(session, messages))
            JSServe.init_from_b64(init_data_b64)
            if (!$(is_offline)){
                JSServe.update_obs($(on_init), true)
            }
        })()"""
    else
        js"""
            // register this session so it gets deleted when it gets removed from dom
            JSServe.register_sub_session($(session.id))
            JSServe.update_obs($(on_init), true)
        """
    end

    final_dom = DOM.div(
        js_dom,
        include_all_assets(page_session, new_deps)...,
        jsrender(session, init),
        id=session.id,
    )

    obs_shared_with_parent = intersect(keys(session.observables), keys(page_session.observables))
    # Those obs are managed by page session, so we don't need to have them in here!
    # This is important when we clean them up, since a child session shouldn't
    # delete any observable already managed by the parent session
    for obs in obs_shared_with_parent
        delete!(session.observables, obs)
    end
    # but in the end, all observables need to be registered
    # with the page session, since that's where javascript will sent all the events
    merge!(page_session.observables, session.observables)

    on(session, session.on_close) do closed
        # remove the sub-session specific js resources
        if closed
            obs_ids = collect(keys(session.observables))
            JSServeLib.delete_observables(page_session, obs_ids)
            # since we deleted all obs shared with the parent session,
            # we can delete savely delete all of them from there
            # Note, that we previously added them all here!
            for (k, o) in session.observables
                delete!(page_session.observables, k)
            end
        end
    end

    # finally give page a connection! :)
    session.connection[] = page_session.connection[]

    return repr(MIME"text/html"(), Hyperscript.Pretty(final_dom))
end

function show_in_iframe(server, session, app)
    session_route = "/$(session.id)"
    # Our default is to display the app in an IFrame, which is a bit complicated
    # we need to resize the iframe based on its content, which is a bit complicated
    # because we can't directly access it. So we sent a message here
    # to the parent iframe, for which we register an
    # event handler via resize_iframe_parent, which then
    # resizes the parent iframe accordingly
    route!(server, session_route => app)
    return jsrender(session, iframe_html(server, session, session_route))
end

function Base.show(io::IO, m::Union{MIME"text/html", MIME"application/prs.juno.plotpane+html"}, app::App)
    if isassigned(CURRENT_PAGE)
        # We are in Page rendering mode!
        page = CURRENT_PAGE[]
        println(io, show_in_page(page, app))
    else
        server = get_server()
        session = Session()
        println(io, show_in_iframe(server, session, app))
    end
end

function iframe_html(server::Server, session::Session, route::String)
    # Display the route we just added in an iframe inline:
    url = online_url(server, route)
    remote_origin = online_url(server, "")
    style = "position: absolute; width: 100%; height: 100%; padding: 0; overflow: hidden; border: none"
    return DOM.iframe(src=url, id=session.id, style=style, scrolling="no")
end

function node_html(session::Session, node::Hyperscript.Node)
    js_dom = DOM.div(jsrender(session, node), id="application-dom")
    # register resources (e.g. observables, assets)
    register_resource!(session, js_dom)
    return repr(MIME"text/html"(), Hyperscript.Pretty(js_dom))
end

"""
    page_html(session::Session, html_body)
Embedds the html_body in a standalone html document!
"""
function page_html(session::Session, html)
    proxy_url = JSSERVE_CONFIGURATION.external_url[]
    serializer = session.url_serializer
    rendered = jsrender(session, html)
    register_resource!(session, rendered)
    session_deps = include_asset.(session.dependencies, (serializer,))
    html_body = DOM.html(
        DOM.head(
            DOM.meta(charset="UTF-8"),
            include_asset(PakoLib, serializer),
            include_asset(MsgPackLib, serializer),
            include_asset(JSServeLib, serializer),
            include_asset(Base64Lib, serializer),
            session_deps...
        ),
        DOM.body(
            DOM.div(rendered, id="application-dom"),
            onload=DontEscape("""
                const proxy_url = '$(proxy_url)'
                const session_id = '$(session.id)'
                JSServe.setup_connection({proxy_url, session_id})
                JSServe.sent_done_loading()
            """)
        )
    )
    return repr(MIME"text/html"(), Hyperscript.Pretty(html_body))
end

function Base.show(io::IOContext, m::MIME"application/vnd.jsserve.application+html", dom::App)
    if get(io, :use_offline_mode, false)
        export_folder = get(io, :export_folder, "/")
        absolute_urls = get(io, :absolute_urls, false)
        content_delivery_url = get(io, :content_delivery_url, "")
        html, session = export_standalone(dom.handler, export_folder;
            absolute_urls=absolute_urls,
            content_delivery_url=content_delivery_url,
            write_index_html=false)
        # We prepare for being offline, but we still start a server while things
        # are online!
        application = get_server()
        application.sessions[session.id] = session
        println(io, html)
    else
        show(io.io, MIME"text/html"(), dom)
    end
end

function Base.show(io::IO, m::MIME"application/vnd.jsserve.application+html", app::App)
    show(IOContext(io), m, app)
end

function Base.show(io::IO, ::MIME"juliavscode/html", app::App)
    show(IOContext(io), MIME"text/html"(), app)
end


text:
# MIT License
# Copyright (c) 2017: Xavier Gandibleux, Anthony Przybylski, Gauthier Soleilhac, and contributors.
__precompile__()
module vOptGeneric
using Combinatorics, Suppressor

using JuMP, MathProgBase

export vModel,
    getvOptData,
    solve,
    writeMOP,
    parseMOP,
    printX_E,
    getvalue,
    getY_N,
    @addobjective,
    @variable,
    @constraint

include("MOP.jl")
include("algorithms.jl")

mutable struct vOptData
    objs::Vector{QuadExpr}               #Objectives
	objSenses::Vector{Symbol}            #Objective Senses
    Y_N::Vector{Vector{Float64}}        #Objective values for each point
    X_E::Vector{Vector{Float64}}        #Variable values for each point
    isacopy::Bool
end

Base.copy(vd::vOptData) = vOptData(deepcopy(vd.objs), deepcopy(vd.objSenses), [], [], true)

function getvOptData(m::Model)
    !haskey(m.ext, :vOpt) && error("This model wasn't created with vOptGeneric")
    return m.ext[:vOpt]::vOptData
end

function vModel(;solver=JuMP.UnsetSolver())
    m = Model(solver=solver)
    m.solvehook = solvehook
    m.printhook = printhook
    m.ext[:vOpt] = vOptData(Vector{QuadExpr}(), #objs
                              Vector{Symbol}(), #objSenses
                              Vector{Vector{Float64}}(), #Y_N
                              Vector{Vector{Float64}}(), #X_E
                              false) #isacopy
    return m
end

function solvehook(m::Model; relax=false, method=nothing, step = 0.5, round_results = false, verbose = true, args...)::Symbol

    vd = getvOptData(m)
    if vd.isacopy
        vd.objs .= copy.(vd.objs, m)
    end

    if method == :epsilon
        return solve_eps(m, step, round_results, verbose ; relaxation=relax, args...)
    elseif method == :dicho || method == :dichotomy
        return solve_dicho(m, round_results ; relaxation=relax, args...)
    elseif method == :Chalmet || method == :chalmet
        return solve_Chalmet(m, step ; relaxation=relax, args...)
    elseif method == :lex || method == :lexico
        return solve_lexico(m, verbose ; relaxation=relax, args...)
    else
        warn("use solve(m, method = :(epsilon | dichotomy | chalmet | lexico) )")
        return :Unsolved
    end

end

function printhook(io::IO, m::Model)
    vd = getvOptData(m)
    for i = 1:length(vd.objs)
        println(vd.objSenses[i] == :Min ? "Min " : "Max ", vd.objs[i])
    end
    str = JuMP.model_str(JuMP.REPLMode, m)
    index = findfirst("Subject to", str)
    index !== nothing && print(str[first(index):end])
end

macro addobjective(m, args...)
    if length(args) != 2
        error("in @addobjective: needs three arguments: model,
                objective sense (Max or Min) and linear expression.")
    end
    m = esc(m)
    sense = args[1]
    if sense == :Min || sense == :Max
        sense = Expr(:quote,sense)
    end
    expr = esc(args[2])
    return quote
        f = @expression($m, $expr)
        !isa(f, JuMP.GenericAffExpr) && error("in @addobjective : vOptGeneric only supports linear objectives")
        vd = $m.ext[:vOpt]
        push!(vd.objSenses, $(esc(sense)))
        push!(vd.objs, convert(QuadExpr, f))
    end
end

# Returns coefficients for the affine part of an objective
function prepAffObjective(m, objaff::JuMP.GenericQuadExpr)
    # Check that no coefficients are NaN/Inf
    JuMP.assert_isfinite(objaff)
    if !JuMP.verify_ownership(m, objaff.aff.vars)
        error("Variable not owned by model present in objective")
    end
    f = zeros(m.numCols)
    @inbounds for ind in 1:length(objaff.aff.vars)
        f[objaff.aff.vars[ind].col] += objaff.aff.coeffs[ind]
    end
    return f
end

@deprecate print_X_E printX_E
function printX_E(m::Model)
    vd = getvOptData(m)
    for i = 1:length(vd.Y_N)
        print(vd.Y_N[i]," : ")
        for j = 1:m.numCols
            if vd.X_E[i][j] != 0
                if m.colCat[j] == :Int || m.colCat[j] == :Bin
                    print(getname(m,j),"=",round(Int,vd.X_E[i][j])," ")
                else
                    print(getname(m,j),"=",vd.X_E[i][j]," ")
                end
            end
        end
        println()
   end
end

function getvalue(v::Variable, i::Int)
    vd = getvOptData(v.m)
    return vd.X_E[i][v.col]
end

function getvalue(arr::Array{Variable}, i::Int)

    ret = similar(arr,Float64)
    isempty(ret) && return ret

    for I in eachindex(arr)
        v = arr[I]
        value = getvalue(v, i)
        ret[I] = value
    end

    return ret
end

getvalue(x::JuMP.JuMPContainer, i::Int) = JuMP._map(v -> getvalue(v, i), x)

function getY_N(m::Model)
    return getvOptData(m).Y_N
end

end#module

ext: jl


text:
module OnlinePCA

using HDF5:
    HDF5File, HDF5Group, h5open
using SparseArrays:
    SparseArrays, SparseMatrixCSC, sortSparseMatrixCSC!
using DelimitedFiles:
    writedlm, readdlm
using Statistics:
    mean, var, median
using LinearAlgebra:
    Diagonal, lu!, qr!, svd, svd!, dot, norm, eigvecs, tr
using Random:
    randperm
using ProgressMeter:
	Progress, next!
using ArgParse:
    ArgParseSettings, parse_args, @add_arg_table
using StatsBase:
	percentile
using DataFrames:
	DataFrame
using GLM:
	glm, coef, IdentityLink, @formula
using Distributions:
	Gamma, ccdf, Chisq
using CodecZstd:
	ZstdCompressorStream, ZstdDecompressorStream

export output, common_parse_commandline, csv2bin, sumr, tenxsumr, filtering, hvg, sgd, oja, ccipca, gd, rsgd, svrg, rsvrg, orthiter, arnoldi, lanczos, halko, algorithm971, rbkiter, singlepass, singlepass2, tenxpca

include("Utils.jl")
include("csv2bin.jl")
include("sumr.jl")
include("tenxsumr.jl")
include("filtering.jl")
include("hvg.jl")
include("sgd.jl")
include("oja.jl")
include("ccipca.jl")
include("gd.jl")
include("rsgd.jl")
include("svrg.jl")
include("rsvrg.jl")
include("orthiter.jl")
include("arnoldi.jl")
include("lanczos.jl")
include("halko.jl")
include("algorithm971.jl")
include("rbkiter.jl")
include("singlepass.jl")
include("singlepass2.jl")
include("tenxpca.jl")

end


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################

isdefined(Base, :__precompile__) && __precompile__()

module JuMP

importall Base.Operators

import MathProgBase

using Calculus
using ReverseDiffSparse
using ForwardDiff
using Compat

import Compat: UTF8String, view

@compat import Base.show

function __init__()
    ENABLE_NLP_RESOLVE[1] = false
end

export
# Objects
    Model, Variable, Norm, AffExpr, QuadExpr, SOCExpr,
    LinearConstraint, QuadConstraint, SDConstraint, SOCConstraint,
    NonlinearConstraint,
    ConstraintRef,
# Functions
    # Model related
    getobjectivevalue, getobjective,
    getobjectivesense, setobjectivesense, setsolver,
    writeLP, writeMPS,
    addSOS1, addSOS2, solve,
    internalmodel,
    # Variable
    setname, getname, setlowerbound, setupperbound, getlowerbound, getupperbound,
    getvalue, setvalue, getdual, setcategory, getcategory,
    getvariable, getconstraint,
    linearindex,
    # Expressions and constraints
    linearterms,

# Macros and support functions
    @LinearConstraint, @LinearConstraints, @QuadConstraint, @QuadConstraints,
    @SOCConstraint, @SOCConstraints,
    @expression, @expressions, @NLexpression, @NLexpressions,
    @variable, @variables, @constraint, @constraints,
    @NLconstraint, @NLconstraints,
    @SDconstraint, @SDconstraints,
    @objective, @NLobjective,
    @NLparameter, @constraintref

include("JuMPContainer.jl")
include("utils.jl")

###############################################################################
# Model class
# Keeps track of all model and column info
abstract AbstractModel
type Model <: AbstractModel
    obj#::QuadExpr
    objSense::Symbol

    linconstr#::Vector{LinearConstraint}
    quadconstr
    sosconstr
    socconstr
    sdpconstr

    # Column data
    numCols::Int
    colNames::Vector{UTF8String}
    colNamesIJulia::Vector{UTF8String}
    colLower::Vector{Float64}
    colUpper::Vector{Float64}
    colCat::Vector{Symbol}

    customNames::Vector

    # Variable cones of the form, e.g. (:SDP, 1:9)
    varCones::Vector{Tuple{Symbol,Any}}

    # Solution data
    objVal
    colVal::Vector{Float64}
    redCosts::Vector{Float64}
    linconstrDuals::Vector{Float64}
    conicconstrDuals::Vector{Float64}
    constrDualMap::Vector{Vector{Int}}
    # internal solver model object
    internalModel
    # Solver+option object from MathProgBase
    solver::MathProgBase.AbstractMathProgSolver
    internalModelLoaded::Bool
    # callbacks
    callbacks
    # lazycallback
    # cutcallback
    # heurcallback

    # hook into a solve call...function of the form f(m::Model; kwargs...),
    # where kwargs get passed along to subsequent solve calls
    solvehook
    # ditto for a print hook
    printhook

    # List of JuMPContainer{Variables} associated with model
    dictList::Vector

    # storage vector for merging duplicate terms
    indexedVector::IndexedVector{Float64}

    nlpdata#::NLPData
    simplify_nonlinear_expressions::Bool

    varDict::Dict{Symbol,Any} # dictionary from variable names to variable objects
    conDict::Dict{Symbol,Any} # dictionary from constraint names to constraint objects
    varData::ObjectIdDict

    getvalue_counter::Int # number of times we call getvalue on a JuMPContainer, so that we can print out a warning
    operator_counter::Int # number of times we add large expressions

    # Extension dictionary - e.g. for robust
    # Extensions should define a type to hold information particular to
    # their functionality, and store an instance of the type in this
    # dictionary keyed on an extension-specific symbol
    ext::Dict{Symbol,Any}
end

# dummy solver
type UnsetSolver <: MathProgBase.AbstractMathProgSolver
end

# Default constructor
function Model(;solver=UnsetSolver(), simplify_nonlinear_expressions::Bool=false)
    if !isa(solver,MathProgBase.AbstractMathProgSolver)
        error("solver argument ($solver) must be an AbstractMathProgSolver")
    end
    Model(zero(QuadExpr),              # obj
          :Min,                        # objSense
          LinearConstraint[],          # linconstr
          QuadConstraint[],            # quadconstr
          SOSConstraint[],             # sosconstr
          SOCConstraint[],             # socconstr
          SDConstraint[],             # sdpconstr
          0,                           # numCols
          UTF8String[],                # colNames
          UTF8String[],                # colNamesIJulia
          Float64[],                   # colLower
          Float64[],                   # colUpper
          Symbol[],                    # colCat
          Variable[],                  # customNames
          Vector{Tuple{Symbol,Any}}[], # varCones
          0,                           # objVal
          Float64[],                   # colVal
          Float64[],                   # redCosts
          Float64[],                   # linconstrDuals
          Float64[],                   # conicconstrDuals
          Any[],                       # constrDualsMap
          nothing,                     # internalModel
          solver,                      # solver
          false,                       # internalModelLoaded
          Any[],                       # callbacks
          nothing,                     # solvehook
          nothing,                     # printhook
          Any[],                       # dictList
          IndexedVector(Float64,0),    # indexedVector
          nothing,                     # nlpdata
          simplify_nonlinear_expressions, # ...
          Dict{Symbol,Any}(),          # varDict
          Dict{Symbol,Any}(),          # conDict
          ObjectIdDict(),              # varData
          0,                           # getvalue_counter
          0,                           # operator_counter
          Dict{Symbol,Any}(),          # ext
    )
end

# Getters/setters
MathProgBase.numvar(m::Model) = m.numCols
MathProgBase.numlinconstr(m::Model) = length(m.linconstr)
MathProgBase.numquadconstr(m::Model) = length(m.quadconstr)
function MathProgBase.numconstr(m::Model)
    c = length(m.linconstr) + length(m.quadconstr) + length(m.sosconstr)
    if m.nlpdata !== nothing
        c += length(m.nlpdata.nlconstr)
    end
    return c
end

for f in MathProgBase.SolverInterface.methods_by_tag[:rewrap]
    eval(Expr(:import,:MathProgBase,f))
    @eval function $f(m::Model)
        # check internal model exists
        if !m.internalModelLoaded
            error("Model not solved")
        else
            return $f(internalmodel(m))
        end
    end
    eval(Expr(:export,f))
end

function getobjective(m::Model)
    traits = ProblemTraits(m)
    if traits.nlp
        error("getobjective() not supported for nonlinear models")
    end
    return m.obj
end

getobjectivevalue(m::Model) = m.objVal
getobjectivesense(m::Model) = m.objSense
function setobjectivesense(m::Model, newSense::Symbol)
    if (newSense != :Max && newSense != :Min)
        error("Model sense must be :Max or :Min")
    end
    m.objSense = newSense
end
setobjective(m::Model, something::Any) =
    error("in setobjective: needs three arguments: model, objective sense (:Max or :Min), and expression.")

setobjective(::Model, ::Symbol, x::Array) =
    error("in setobjective: array of size $(size(x)) passed as objective; only scalar objectives are allowed")

function setsolver(m::Model, solver::MathProgBase.AbstractMathProgSolver)
    m.solver = solver
    m.internalModel = nothing
    m.internalModelLoaded = false
    nothing
end
# Deep copy the model
function Base.copy(source::Model)

    dest = Model()
    dest.solver = source.solver  # The two models are linked by this

    # Objective
    dest.obj = copy(source.obj, dest)
    dest.objSense = source.objSense

    # Constraints
    dest.linconstr  = map(c->copy(c, dest), source.linconstr)
    dest.quadconstr = map(c->copy(c, dest), source.quadconstr)
    dest.sosconstr  = map(c->copy(c, dest), source.sosconstr)
    dest.sdpconstr  = map(c->copy(c, dest), source.sdpconstr)

    # Variables
    dest.numCols = source.numCols
    dest.colNames = source.colNames[:]
    dest.colNamesIJulia = source.colNamesIJulia[:]
    dest.colLower = source.colLower[:]
    dest.colUpper = source.colUpper[:]
    dest.colCat = source.colCat[:]

    # varCones
    dest.varCones = copy(source.varCones)

    # callbacks and hooks
    if !isempty(source.callbacks)
        error("Copying callbacks is not supported")
    end
    if source.solvehook !== nothing
        dest.solvehook = source.solvehook
    end
    if source.printhook !== nothing
        dest.printhook = source.printhook
    end

    # variable/extension dicts
    if !isempty(source.ext)
        dest.ext = similar(source.ext)
        for (key, val) in source.ext
            dest.ext[key] = try
                copy(source.ext[key])
            catch
                error("Error copying extension dictionary. Is `copy` defined for all your user types?")
            end
        end
    end
    dest.varDict = Dict{Symbol,Any}()
    for (symb,v) in source.varDict
        dest.varDict[symb] = copy(v, dest)
    end

    dest.conDict = Dict{Symbol,Any}()
    # TODO: implement constraint copying
    # for (symb,v) in source.conDict
    #     dest.conDict[symb] = copy(v, dest)
    # end

    # varData---possibly shouldn't copy

    if source.nlpdata !== nothing
        dest.nlpdata = copy(source.nlpdata)
    end

    return dest
end

internalmodel(m::Model) = m.internalModel

setsolvehook(m::Model, f) = (m.solvehook = f)
setprinthook(m::Model, f) = (m.printhook = f)

#############################################################################
# AbstractConstraint
# Abstract base type for all constraint types
abstract AbstractConstraint
# Abstract base type for all scalar types
# In JuMP, used only for Variable. Useful primarily for extensions
abstract AbstractJuMPScalar

Base.start(::AbstractJuMPScalar) = false
Base.next(x::AbstractJuMPScalar, state) = (x, true)
Base.done(::AbstractJuMPScalar, state) = state
Base.isempty(::AbstractJuMPScalar) = false

#############################################################################
# Variable class
# Doesn't actually do much, just a pointer back to the model
immutable Variable <: AbstractJuMPScalar
    m::Model
    col::Int
end

linearindex(x::Variable) = x.col
Base.isequal(x::Variable,y::Variable) = (x.col == y.col) && (x.m === y.m)

Variable(m::Model, lower, upper, cat::Symbol, name::AbstractString="", value::Number=NaN) =
    error("Attempt to create scalar Variable with lower bound of type $(typeof(lower)) and upper bound of type $(typeof(upper)). Bounds must be scalars in Variable constructor.")

function Variable(m::Model,lower::Number,upper::Number,cat::Symbol,name::AbstractString="",value::Number=NaN)
    m.numCols += 1
    push!(m.colNames, name)
    push!(m.colNamesIJulia, name)
    push!(m.colLower, convert(Float64,lower))
    push!(m.colUpper, convert(Float64,upper))
    push!(m.colCat, cat)
    push!(m.colVal,value)
    if cat == :Fixed
        @assert lower == upper
        m.colVal[end] = lower
    end
    if m.internalModelLoaded
        if method_exists(MathProgBase.addvar!, (typeof(m.internalModel),Vector{Int},Vector{Float64},Float64,Float64,Float64))
            MathProgBase.addvar!(m.internalModel,float(lower),float(upper),0.0)
        else
            Base.warn_once("Solver does not appear to support adding variables to an existing model. Hot-start is disabled.")
            m.internalModelLoaded = false
        end
    end
    return Variable(m, m.numCols)
end

# Name setter/getters
function setname(v::Variable,n::AbstractString)
    push!(v.m.customNames, v)
    v.m.colNames[v.col] = n
    v.m.colNamesIJulia[v.col] = n
end
getname(m::Model, col) = var_str(REPLMode, m, col)
getname(v::Variable) = var_str(REPLMode, v.m, v.col)

# Bound setter/getters
function setlowerbound(v::Variable,lower::Number)
    v.m.colCat[v.col] == :Fixed && error("use setvalue for changing the value of a fixed variable")
    v.m.colLower[v.col] = lower
end
function setupperbound(v::Variable,upper::Number)
    v.m.colCat[v.col] == :Fixed && error("use setvalue for changing the value of a fixed variable")
    v.m.colUpper[v.col] = upper
end
getlowerbound(v::Variable) = v.m.colLower[v.col]
getupperbound(v::Variable) = v.m.colUpper[v.col]

# Value setter/getter
function setvalue(v::Variable, val::Number)
    v.m.colVal[v.col] = val
    if v.m.colCat[v.col] == :Fixed
        v.m.colLower[v.col] = val
        v.m.colUpper[v.col] = val
    end
end

# internal method that doesn't print a warning if the value is NaN
_getValue(v::Variable) = v.m.colVal[v.col]

function getvalue(v::Variable)
    ret = _getValue(v)
    if isnan(ret)
        Base.warn("Variable value not defined for $(getname(v)). Check that the model was properly solved.")
    end
    ret
end

function getvalue(arr::Array{Variable})
    ret = similar(arr, Float64)
    # return immediately for empty array
    if isempty(ret)
        return ret
    end
    # warnedyet is set to true if we've already warned for a component of a JuMPContainer
    warnedyet = false
    m = first(arr).m
    # whether this was constructed via @variable, essentially
    registered = haskey(m.varData, arr)
    for I in eachindex(arr)
        v = arr[I]
        value = _getValue(v)
        ret[I] = value
        if !warnedyet && isnan(value)
            if registered
                Base.warn("Variable value not defined for component of $(m.varData[arr].name). Check that the model was properly solved.")
                warnedyet = true
            else
                Base.warn("Variable value not defined for $(m.colNames[v.col]). Check that the model was properly solved.")
            end
        end
    end
    # Copy printing data from @defVar for Array{Variable} to corresponding Array{Float64} of values
    if registered
        m.varData[ret] = m.varData[arr]
    end
    ret
end

# Dual value (reduced cost) getter
function getdual(v::Variable)
    if length(v.m.redCosts) < MathProgBase.numvar(v.m)
        error("Variable bound duals (reduced costs) not available. Check that the model was properly solved and no integer variables are present.")
    end
    return v.m.redCosts[v.col]
end

const var_cats = [:Cont, :Int, :Bin, :SemiCont, :SemiInt]
function setcategory(v::Variable, cat::Symbol)
    cat in var_cats || error("Unrecognized variable category $cat. Should be one of:\n    $var_cats")
    v.m.colCat[v.col] = cat
end

getcategory(v::Variable) = v.m.colCat[v.col]

Base.zero(::Type{Variable}) = AffExpr(Variable[],Float64[],0.0)
Base.zero(::Variable) = zero(Variable)
Base.one(::Type{Variable}) = AffExpr(Variable[],Float64[],1.0)
Base.one(::Variable) = one(Variable)

function verify_ownership(m::Model, vec::Vector{Variable})
    n = length(vec)
    @inbounds for i in 1:n
        vec[i].m !== m && return false
    end
    return true
end

Base.copy(v::Variable, new_model::Model) = Variable(new_model, v.col)
Base.copy(x::Void, new_model::Model) = nothing
function Base.copy(v::Array{Variable}, new_model::Model)
    ret = similar(v, Variable, size(v))
    for I in eachindex(v)
        ret[I] = Variable(new_model, v[I].col)
    end
    ret
end

# Copy methods for variable containers
Base.copy(d::JuMPContainer) = map(copy, d)
Base.copy(d::JuMPContainer, new_model::Model) = map(x -> copy(x, new_model), d)

###############################################################################
# GenericAffineExpression, AffExpr
# GenericRangeConstraint, LinearConstraint
include("affexpr.jl")

###############################################################################
# GenericQuadExpr, QuadExpr
# GenericQuadConstraint, QuadConstraint
include("quadexpr.jl")

##########################################################################
# GenericNorm, Norm
# GenericNormExpr. GenericSOCExpr, SOCExpr
# GenericSOCConstraint, SOCConstraint
include("norms.jl")

##########################################################################
# SOSConstraint  (special ordered set constraints)
include("sos.jl")

##########################################################################
# SDConstraint is a (dual) semidefinite constraint of the form
# ∑ cᵢ Xᵢ ≥ D, where D is a n×n symmetric data matrix, cᵢ are
# scalars, and Xᵢ are n×n symmetric variable matrices. The inequality
# is taken w.r.t. the psd partial order.
type SDConstraint <: AbstractConstraint
    terms
end

# Special-case X ≥ 0, which is often convenient
function SDConstraint(lhs::Matrix, rhs::Number)
    rhs == 0 || error("Cannot construct a semidefinite constraint with nonzero scalar bound $rhs")
    SDConstraint(lhs)
end

function addconstraint(m::Model, c::SDConstraint)
    push!(m.sdpconstr,c)
    m.internalModelLoaded = false
    ConstraintRef{Model,SDConstraint}(m,length(m.sdpconstr))
end

# helper method for mapping going on below
Base.copy(x::Number, new_model::Model) = copy(x)

Base.copy(c::SDConstraint, new_model::Model) =
    SDConstraint(map(t -> copy(t, new_model), c.terms))

##########################################################################
# ConstraintRef
# Reference to a constraint for retrieving solution info
immutable ConstraintRef{M<:AbstractModel,T<:AbstractConstraint}
    m::M
    idx::Int
end

typealias LinConstrRef ConstraintRef{Model,LinearConstraint}

LinearConstraint(ref::LinConstrRef) = ref.m.linconstr[ref.idx]::LinearConstraint

linearindex(x::ConstraintRef) = x.idx

function getdual(c::LinConstrRef)
    if length(c.m.linconstrDuals) != MathProgBase.numlinconstr(c.m)
        error("Dual solution not available. Check that the model was properly solved and no integer variables are present.")
    end
    return c.m.linconstrDuals[c.idx]
end

# Returns the number of non-infinity and nonzero bounds on variables
function getNumBndRows(m::Model)
    numBounds = 0
    for i in 1:m.numCols
        seen = false
        lb, ub = m.colLower[i], m.colUpper[i]
        for (_,cone) in m.varCones
            if i in cone
                seen = true
                @assert lb == -Inf && ub == Inf
                break
            end
        end

        if !seen
            if lb != -Inf && lb != 0
                numBounds += 1
            end
            if ub != Inf && ub != 0
                numBounds += 1
            end
        end
    end
    return numBounds
end

# Returns the number of second-order cone constraints
function getNumSOCRows(m::Model)
    numSOCRows = 0
    for con in m.socconstr
        numSOCRows += length(con.normexpr.norm.terms) + 1
    end
    return numSOCRows
end

function getdual(c::ConstraintRef{Model,SOCConstraint})
    numBndRows = getNumBndRows(c.m)
    numSOCRows = getNumSOCRows(c.m)
    if length(c.m.conicconstrDuals) != (MathProgBase.numlinconstr(c.m) + numBndRows + numSOCRows)
        error("Dual solution not available. Check that the model was properly solved and no integer variables are present.")
    end
    return c.m.conicconstrDuals[
        c.m.constrDualMap[MathProgBase.numlinconstr(c.m) + numBndRows + c.idx]]
end

function setRHS(c::LinConstrRef, rhs::Number)
    constr = c.m.linconstr[c.idx]
    sen = sense(constr)
    if sen == :range
        error("Modifying range constraints is currently unsupported.")
    elseif sen == :(==)
        constr.lb = float(rhs)
        constr.ub = float(rhs)
    elseif sen == :>=
        constr.lb = float(rhs)
    else
        @assert sen == :<=
        constr.ub = float(rhs)
    end
end

Variable(m::Model,lower::Number,upper::Number,cat::Symbol,objcoef::Number,
    constraints::JuMPArray,coefficients::Vector{Float64}, name::AbstractString="", value::Number=NaN) =
    Variable(m, lower, upper, cat, objcoef, constraints.innerArray, coefficients, name, value)

# add variable to existing constraints
function Variable(m::Model,lower::Number,upper::Number,cat::Symbol,objcoef::Number,
    constraints::Vector,coefficients::Vector{Float64}, name::AbstractString="", value::Number=NaN)
    for c in constraints
        if !isa(c, LinConstrRef)
            error("Unexpected constraint of type $(typeof(c)). Column-wise modeling only supported for linear constraints")
        end
    end
    @assert cat != :Fixed || (lower == upper)
    m.numCols += 1
    push!(m.colNames, name)
    push!(m.colNamesIJulia, name)
    push!(m.colLower, convert(Float64,lower))
    push!(m.colUpper, convert(Float64,upper))
    push!(m.colCat, cat)
    push!(m.colVal, value)
    if cat == :Fixed
        @assert lower == upper
        m.colVal[end] = lower
    end
    v = Variable(m,m.numCols)
    # add to existing constraints
    @assert length(constraints) == length(coefficients)
    for i in 1:length(constraints)
        c::LinearConstraint = m.linconstr[constraints[i].idx]
        coef = coefficients[i]
        push!(c.terms.vars,v)
        push!(c.terms.coeffs,coef)
    end
    push!(m.obj.aff.vars, v)
    push!(m.obj.aff.coeffs,objcoef)

    if m.internalModelLoaded
        if method_exists(MathProgBase.addvar!, (typeof(m.internalModel),Vector{Int},Vector{Float64},Float64,Float64,Float64))
            MathProgBase.addvar!(m.internalModel,Int[c.idx for c in constraints],coefficients,float(lower),float(upper),float(objcoef))
        else
            Base.warn_once("Solver does not appear to support adding variables to an existing model. Hot-start is disabled.")
            m.internalModelLoaded = false
        end
    end

    return v
end

# handle dictionary of variables
function registervar(m::Model, varname::Symbol, value)
    if haskey(m.varDict, varname)
        m.varDict[varname] = nothing # indicate duplicate variable
    else
        m.varDict[varname] = value
    end
    return value
end
registervar(m::Model, varname, value) = value # variable name isn't a simple symbol, ignore

function registercon(m::Model, conname::Symbol, value)
    if haskey(m.conDict, conname)
        m.conDict[conname] = nothing # indicate duplicate variable
    else
        m.conDict[conname] = value
    end
    return value
end
registercon(m::Model, conname, value) = value # constraint name isn't a simple symbol, ignore

function getvariable(m::Model, varname::Symbol)
    if !haskey(m.varDict, varname)
        error("No variable with name $varname")
    elseif m.varDict[varname] === nothing
        error("Multiple variables with name $varname")
    else
        return m.varDict[varname]
    end
end

function getconstraint(m::Model, conname::Symbol)
    if !haskey(m.conDict, conname)
        error("No constraint with name $conname")
    elseif m.conDict[conname] === nothing
        error("Multiple constraints with name $conname")
    else
        return m.conDict[conname]
    end
end

# usage warnings
function getvalue_warn(x::JuMPContainer{Variable})
    isempty(x) && return
    v = first(values(x))
    m = v.m
    m.getvalue_counter += 1
    if m.getvalue_counter > 400
        Base.warn_once("getvalue has been called on a collection of variables a large number of times. For performance reasons, this should be avoided. Instead of getvalue(x)[a,b,c], use getvalue(x[a,b,c]) to avoid temporary allocations.")
    end
    return
end
getvalue_warn(x::JuMPContainer) = nothing

function operator_warn(lhs::AffExpr,rhs::AffExpr)
    if length(lhs.vars) > 50 || length(rhs.vars) > 50
        if length(lhs.vars) > 1
            m = lhs.vars[1].m
            m.operator_counter += 1
            if m.operator_counter > 20000
                Base.warn_once("The addition operator has been used on JuMP expressions a large number of times. This warning is safe to ignore but may indicate that model generation is slower than necessary. For performance reasons, you should not add expressions in a loop. Instead of x += y, use append!(x,y) to modify x in place. If y is a single variable, you may also use push!(x, coef, y) in place of x += coef*y.")
            end
        end
    end
    return
end
operator_warn(lhs,rhs) = nothing

##########################################################################
# Types used in the nonlinear code
immutable NonlinearExpression
    m::Model
    index::Int
end

immutable NonlinearParameter <: AbstractJuMPScalar
    m::Model
    index::Int
end

##########################################################################
# Behavior that's uniform across all JuMP "scalar" objects

typealias JuMPTypes Union{AbstractJuMPScalar,
                          NonlinearExpression,
                          Norm,
                          GenericAffExpr,
                          QuadExpr,
                          SOCExpr}
typealias JuMPScalars Union{Number,JuMPTypes}

# would really want to do this on ::Type{T}, but doesn't work on v0.4
Base.eltype{T<:JuMPTypes}(::T) = T
Base.size(::JuMPTypes) = ()
Base.size(x::JuMPTypes,d::Int) = 1
Base.ndims(::JuMPTypes) = 0

##########################################################################
# Operator overloads
include("operators.jl")
# Writers - we support MPS (MILP + QuadObj), LP (MILP)
include("writers.jl")
# Macros - @defVar, sum{}, etc.
include("macros.jl")
# Solvers
include("solvers.jl")
# Callbacks - lazy, cuts, ...
include("callbacks.jl")
# Nonlinear-specific code
include("nlp.jl")
# Pretty-printing of JuMP-defined types.
include("print.jl")
# Deprecations
include("deprecated.jl")

getvalue{T<:JuMPTypes}(arr::Array{T}) = map(getvalue, arr)

function setvalue{T<:AbstractJuMPScalar}(set::Array{T}, val::Array)
    promote_shape(size(set), size(val)) # Check dimensions match
    for I in eachindex(set)
        setvalue(set[I], val[I])
    end
    nothing
end

##########################################################################
end


text:
#########################################################################
# A wrapper around an error that happened in a Julia callback

struct PyJlError <: Exception
    err
    traceback

    # Eventually hold a python backtrace
    pytraceback
end

PyJlError(err, trace) = PyJlError(err, trace, nothing)
PyJlError(e::PyJlError, pytrace) = PyJlError(e.err, e.traceback, pytrace)

function show(io::IO, e::PyJlError)
    println(io, "(in a Julia function called from Python):")

    er = e.err
    bt = e.traceback
    io = io
    bt = scrub_pycall_backtrace(bt)
    lvl = get(io, :level, 0)

    if er isa PyJlError
        i = length(bt)
        i > length(er.trace) || deleteat!(er.trace, i:length(er.trace))
        io = IOContext(io, :level=>lvl+1)
        println(io, "----------------------")
        printstyled(io, "ERROR: "; bold=true, color=Base.error_color())
    else
        printstyled(io, "JULIA: "; bold=true, color=Base.error_color())
    end

    showerror(IOContext(io, :limit => true), er, bt, backtrace = bt!==nothing)

    if e.pytraceback !== nothing && !ispynull(e.pytraceback)
        o = pycall(format_traceback, PyObject, e.pytraceback)
        if !ispynull(o)
            println(io)
            println(io)
            println(io, "Stacktrace (Python):")

            for s in PyVector{AbstractString}(o)
                print(io, s)
            end
        end
    end

    if lvl != 0
        println(io)
        print(io, "----------------------")
    end
end

function scrub_pycall_backtrace(bt)
    if bt !== nothing && !(bt isa Vector{Any}) # ignore our sentinel value types
        bt = stacktrace(bt, true)
        # remove REPL-related frames from interactive printing
        eval_ind = findlast(frame -> !frame.from_c && frame.func === :_pyjlwrap_call, bt)
        eval_ind === nothing || deleteat!(bt, eval_ind+1:length(bt))
    end
    return bt
end

#########################################################################
# Wrapper around Python exceptions

struct PyError <: Exception
    msg::String # message string from Julia context, or "" if none

    # info returned by PyErr_Fetch/PyErr_Normalize
    T::PyObject
    val::PyObject
    traceback::PyObject

    # generate a PyError object.  Should normally only be called when
    # PyErr_Occurred returns non-NULL, and clears the Python error
    # indicator.
end

function show(io::IO, e::PyError)
    #print(io, "PyError",
    #      isempty(e.msg) ? e.msg : string(" (",e.msg,")"),
    #      " ")
    print(io, "PyError--")

    if ispynull(e.T)
        println(io, "None")
    else
        println(io, pystring(e.T), "\n", pystring(e.val))
    end
    if !ispynull(e.traceback)
        o = pycall(format_traceback, PyObject, e.traceback)
        if !ispynull(o)
            for s in PyVector{AbstractString}(o)
                print(io, s)
            end
        end
    end
end

#########################################################################
# Conversion of Python exceptions into Julia exceptions

# whether a Python exception has occurred
pyerr_occurred() = ccall((@pysym :PyErr_Occurred), PyPtr, ()) != C_NULL

# call to discard Python exceptions
pyerr_clear() = ccall((@pysym :PyErr_Clear), Cvoid, ())

function pyerr_check(msg::AbstractString, val::Any)
    pyerr_occurred() && throw(pyerror(msg))
    val # the val argument is there just to pass through to the return value
end

pyerr_check(msg::AbstractString) = pyerr_check(msg, nothing)
pyerr_check() = pyerr_check("")

# extract function name from ccall((@pysym :Foo)...) or ccall(:Foo,...) exprs
callsym(s::Symbol) = s
callsym(s::QuoteNode) = s.value
import Base.Meta.isexpr
callsym(ex::Expr) = isexpr(ex,:macrocall,2) ? callsym(ex.args[2]) : isexpr(ex,:ccall) ? callsym(ex.args[1]) : ex

"""
    _handle_error(msg)

Throw a PyError if available, otherwise throw ErrorException.
This is a hack to manually do the optimization described in
https://github.com/JuliaLang/julia/issues/29688
"""
@noinline function _handle_error(msg)
    pyerr_check(msg)
    error(msg, " failed")
end

# Macros for common pyerr_check("Foo", ccall((@pysym :Foo), ...)) pattern.
macro pycheck(ex)
    :(pyerr_check($(string(callsym(ex))), $(esc(ex))))
end

# Macros to check that ccall((@pysym :Foo), ...) returns value != bad
macro pycheckv(ex, bad)
    quote
        val = $(esc(ex))
        if val == $(esc(bad))
            _handle_error($(string(callsym(ex))))
        end
        val
    end
end
macro pycheckn(ex)
    :(@pycheckv $(esc(ex)) C_NULL)
end
macro pycheckz(ex)
    :(@pycheckv $(esc(ex)) -1)
end

function pyerror(msg::AbstractString)
    ptype, pvalue, ptraceback = Ref{PyPtr}(), Ref{PyPtr}(), Ref{PyPtr}()
    # equivalent of passing C pointers &exc[1], &exc[2], &exc[3]:
    ccall((@pysym :PyErr_Fetch), Cvoid, (Ref{PyPtr},Ref{PyPtr},Ref{PyPtr}), ptype, pvalue, ptraceback)
    ccall((@pysym :PyErr_NormalizeException), Cvoid, (Ref{PyPtr},Ref{PyPtr},Ref{PyPtr}), ptype, pvalue, ptraceback)
    pyerror(msg, PyObject(ptype[]), PyObject(pvalue[]), PyObject(ptraceback[]))
end

function pyerror(msg::AbstractString, e::PyError)
    pyerror(msg, e.T, e.val, e.traceback)
end

function pyerror(msg::AbstractString, ptype::PyObject, pvalue::PyObject, ptraceback::PyObject)
    pargs = _getproperty(pvalue, "args")

    # If the value of the error is a PyJlError, it was generated in a pyjlwrap callback, and
    # we forward it.
    if pargs != C_NULL
        args = PyObject(pargs)
        if length(args) > 0
            arg = PyObject(@pycheckn ccall((@pysym :PySequence_GetItem), PyPtr, (PyPtr,Int), args, 0))
            if is_pyjlwrap(arg)
                jarg = unsafe_pyjlwrap_to_objref(arg)
                if jarg isa PyJlError
                    return PyJlError(jarg, ptraceback)
                end
            end
        end
    end

    return PyError(msg, ptype, pvalue, ptraceback)
end

#########################################################################
# Mapping of Julia Exception types to Python exceptions

const pyexc = Dict{DataType, PyPtr}()
mutable struct PyIOError <: Exception end

function pyexc_initialize()
    pyexc[Exception] = @pyglobalobjptr :PyExc_RuntimeError
    pyexc[ErrorException] = @pyglobalobjptr :PyExc_RuntimeError
    pyexc[SystemError] = @pyglobalobjptr :PyExc_SystemError
    pyexc[TypeError] = @pyglobalobjptr :PyExc_TypeError
    pyexc[Meta.ParseError] = @pyglobalobjptr :PyExc_SyntaxError
    pyexc[ArgumentError] = @pyglobalobjptr :PyExc_ValueError
    pyexc[KeyError] = @pyglobalobjptr :PyExc_KeyError
    pyexc[LoadError] = @pyglobalobjptr :PyExc_ImportError
    pyexc[MethodError] = @pyglobalobjptr :PyExc_RuntimeError
    pyexc[EOFError] = @pyglobalobjptr :PyExc_EOFError
    pyexc[BoundsError] = @pyglobalobjptr :PyExc_IndexError
    pyexc[DivideError] = @pyglobalobjptr :PyExc_ZeroDivisionError
    pyexc[DomainError] = @pyglobalobjptr :PyExc_RuntimeError
    pyexc[OverflowError] = @pyglobalobjptr :PyExc_OverflowError
    pyexc[InexactError] = @pyglobalobjptr :PyExc_ArithmeticError
    pyexc[OutOfMemoryError] = @pyglobalobjptr :PyExc_MemoryError
    pyexc[StackOverflowError] = @pyglobalobjptr :PyExc_MemoryError
    pyexc[UndefRefError] = @pyglobalobjptr :PyExc_RuntimeError
    pyexc[InterruptException] = @pyglobalobjptr :PyExc_KeyboardInterrupt
    pyexc[PyIOError] = @pyglobalobjptr :PyExc_IOError
end

_showerror_string(io::IO, e, ::Nothing) = showerror(io, e)
_showerror_string(io::IO, e, bt) = showerror(io, e, bt)

# bt argument defaults to nothing, to delay dispatching on the presence of a
# backtrace until after the try-catch block
"""
    showerror_string(e) :: String

Convert output of `showerror` to a `String`.  Since this function may
be called via Python C-API, it tries to not throw at all cost.
"""
function showerror_string(e::T, bt = nothing) where {T}
    try
        io = IOBuffer()
        _showerror_string(io, e, bt)
        return String(take!(io))
    catch
        try
            return """
                   $e
                   ERROR: showerror(::IO, ::$T) failed!"""
        catch
            try
                return """
                       $T
                       ERROR: showerror(::IO, ::$T) failed!
                       ERROR: string(::$T) failed!"""
            catch
                name = try
                    io = IOBuffer()
                    Base.show_datatype(io, T)
                    String(take!(io))
                catch
                    "_UNKNOWN_TYPE_"
                end
                return """
                       Unprintable error.
                       ERROR: showerror(::IO, ::$name) failed!
                       ERROR: string(::$name) failed!
                       ERROR: string($name) failed!"""
            end
        end
    end
end

function pyraise(e, bt = nothing)
    eT = typeof(e)
    pyeT = haskey(pyexc::Dict, eT) ? pyexc[eT] : pyexc[Exception]
    err = PyJlError(e, bt)
    ccall((@pysym :PyErr_SetObject), Cvoid, (PyPtr, PyPtr),
          pyeT, PyObject(err))
end

# Second argument allows for backtraces passed to `pyraise` to be ignored.
function pyraise(e::PyError, ::Vector = [])
    ccall((@pysym :PyErr_Restore), Cvoid, (PyPtr, PyPtr, PyPtr),
          e.T, e.val, e.traceback)
    # refs were stolen
    setfield!(e.T, :o, PyPtr_NULL)
    setfield!(e.val, :o, PyPtr_NULL)
    setfield!(e.traceback, :o, PyPtr_NULL)
end

"""
    @pyraise e

Throw the exception `e` to Python, attaching a backtrace.  This macro should only be
used in a `catch` block so that `catch_backtrace()` is valid.
"""
macro pyraise(e)
    :(pyraise($(esc(e)), catch_backtrace()))
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -10.727396450032298 <= q <= 10.341307641992039)

                     add_NL_constraint(m, :((-0.7594921333542755 + 0.29083593057294754*$(x[1]) + -0.853464816888847*$(x[2]) + 0.355548149730335*$(x[3]) + -0.9799338681626466*$(x[4]))*(1 + erf((-0.7594921333542755 + 0.29083593057294754*$(x[1]) + -0.853464816888847*$(x[2]) + 0.355548149730335*$(x[3]) + -0.9799338681626466*$(x[4]))/sqrt(2)))/2 + (0.7967948470186679 + 0.4125187490566202*$(x[1]) + -0.07888951328395777*$(x[2]) + 0.5307209485984372*$(x[3]) + -0.6243253014424854*$(x[4]))*(1 + erf((0.7967948470186679 + 0.4125187490566202*$(x[1]) + -0.07888951328395777*$(x[2]) + 0.5307209485984372*$(x[3]) + -0.6243253014424854*$(x[4]))/sqrt(2)))/2 + (0.5596303002912069 + 0.036981069983024195*$(x[1]) + -0.36792399865103764*$(x[2]) + 0.8775927268629178*$(x[3]) + -0.988765970566825*$(x[4]))*(1 + erf((0.5596303002912069 + 0.036981069983024195*$(x[1]) + -0.36792399865103764*$(x[2]) + 0.8775927268629178*$(x[3]) + -0.988765970566825*$(x[4]))/sqrt(2)))/2 + (-0.4093924155176669 + 0.2071206814621389*$(x[1]) + -0.48217712463878026*$(x[2]) + -0.6673158871950284*$(x[3]) + -0.9659564229118476*$(x[4]))*(1 + erf((-0.4093924155176669 + 0.2071206814621389*$(x[1]) + -0.48217712463878026*$(x[2]) + -0.6673158871950284*$(x[3]) + -0.9659564229118476*$(x[4]))/sqrt(2)))/2 + (-0.38058500245806215 + 0.5074069251520332*$(x[1]) + 0.613305475914625*$(x[2]) + -0.0940610867820495*$(x[3]) + 0.5995073981555823*$(x[4]))*(1 + erf((-0.38058500245806215 + 0.5074069251520332*$(x[1]) + 0.613305475914625*$(x[2]) + -0.0940610867820495*$(x[3]) + 0.5995073981555823*$(x[4]))/sqrt(2)))/2 - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Example: Compartmental models in epidemiology
#=
- https://github.com/epirecipes/sir-julia/blob/master/markdown/function_map/function_map.md
- https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#Deterministic_versus_stochastic_epidemic_models
=#
using ModelingToolkit, Test

@inline function rate_to_proportion(r,t)
    1-exp(-r*t)
end;

# Independent and dependent variables and parameters
@parameters t c nsteps δt β γ
D = Difference(t; dt=0.1)
@variables S(t) I(t) R(t) next_S(t) next_I(t) next_R(t)

infection = rate_to_proportion(β*c*I/(S+I+R),δt)*S
recovery = rate_to_proportion(γ,δt)*I

# Equations
eqs = [D(S) ~ S-infection,
       D(I) ~ I+infection-recovery,
       D(R) ~ R+recovery]

# System
@named sys = DiscreteSystem(eqs,t,[S,I,R],[c,nsteps,δt,β,γ]; controls = [β, γ])

# Problem
u0 = [S => 990.0, I => 10.0, R => 0.0]
p = [β => 0.05, c => 10.0, γ => 0.25, δt => 0.1, nsteps => 400]
tspan = (0.0,ModelingToolkit.value(substitute(nsteps,p))) # value function (from Symbolics) is used to convert a Num to Float64
prob_map = DiscreteProblem(sys,u0,tspan,p)

# Solution
using OrdinaryDiffEq
sol_map = solve(prob_map,FunctionMap());

# Direct Implementation

function sir_map!(u_diff,u,p,t)
    (S,I,R) = u
    (β,c,γ,δt) = p
    N = S+I+R
    infection = rate_to_proportion(β*c*I/N,δt)*S
    recovery = rate_to_proportion(γ,δt)*I
    @inbounds begin
        u_diff[1] = S-infection
        u_diff[2] = I+infection-recovery
        u_diff[3] = R+recovery
    end
    nothing
end;
u0 = [990.0,10.0,0.0];
p = [0.05,10.0,0.25,0.1];
prob_map = DiscreteProblem(sir_map!,u0,tspan,p);
sol_map2 = solve(prob_map,FunctionMap());

@test Array(sol_map) ≈ Array(sol_map2)

# Delayed difference equation
@parameters t
@variables x(..) y(..) z(t)
D1 = Difference(t; dt=1.5)
D2 = Difference(t; dt=2)

@test ModelingToolkit.is_delay_var(Symbolics.value(t), Symbolics.value(x(t-2)))
@test ModelingToolkit.is_delay_var(Symbolics.value(t), Symbolics.value(y(t-1)))
@test !ModelingToolkit.is_delay_var(Symbolics.value(t), Symbolics.value(z))
@test_throws ErrorException ModelingToolkit.get_delay_val(Symbolics.value(t), Symbolics.arguments(Symbolics.value(x(t+2)))[1])
@test_throws ErrorException z(t)

# Equations
eqs = [
    D1(x(t)) ~ 0.4x(t) + 0.3x(t-1.5) + 0.1x(t-3),
    D2(y(t)) ~ 0.3y(t) + 0.7y(t-2) + 0.1z,
]

# System
@named sys = DiscreteSystem(eqs,t,[x(t),x(t-1.5),x(t-3),y(t),y(t-2),z],[])

eqs2, max_delay = ModelingToolkit.linearize_eqs(sys; return_max_delay=true)

@test max_delay[Symbolics.operation(Symbolics.value(x(t)))] ≈ 3
@test max_delay[Symbolics.operation(Symbolics.value(y(t)))] ≈ 2

linearized_eqs = [
    eqs
    x(t - 3.0) ~ x(t - 1.5)
    x(t - 1.5) ~ x(t)
    y(t - 2.0) ~ y(t)
]
@test all(eqs2 .== linearized_eqs)


text:
module QuCmp

<<<<<<< HEAD
<<<<<<< HEAD
=======
>>>>>>> roger-work
abstract QuComput

include("const.jl")
include("utils/LogicExpr.jl")
include("utils/math.jl")
include("Adiabatic/Adiabatic.jl")

<<<<<<< HEAD
=======
# package code goes here
>>>>>>> ad742772d3e82260d55f64ba87c31ee32eb94575
=======
>>>>>>> roger-work

end # module

max_stars_repo_name: Roger-luo/QuCmp

max_issues_repo_name: Roger-luo/QuCmp

max_forks_repo_name: Roger-luo/QuCmp


text:
n = parse(Int64, readline())

if n % 2 == 1
    println("Weird")
elseif n % 2 == 0
    if n >= 2 && n <= 5
        println("Not Weird")
    elseif n >= 6 && n <= 20
        println("Weird")
    else
        println("Not Weird")
    end
end

max_stars_repo_path: Hackerrank/30 Days of Code/Julia/day 03.jl

max_stars_repo_name: Next-Gen-UI/Code-Dynamics

max_issues_repo_name: Next-Gen-UI/Code-Dynamics

max_forks_repo_name: Next-Gen-UI/Code-Dynamics


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function mphase(jd::AbstractFloat)
    ram, decm, dism = moonpos(jd, radians=true)
    ras, decs = sunpos(jd, radians=true)
    # phi: geocentric elongation of the Moon from the Sun
    # inc: selenocentric (Moon centered) elongation of the Earth from the Sun
    sin_decs, cos_decs = sincos(decs)
    sin_decm, cos_decm = sincos(decm)
    phi = acos(sin_decs * sin_decm + cos_decs * cos_decm * cos(ras - ram))
    # "dism" is in kilometers, AU in meters
    sin_phi, cos_phi = sincos(phi)
    inc = atan(AU * sin_phi, dism * 1000 - AU * cos_phi)
    return (1 + cos(inc))/2
end

"""
    mphase(jd) -> k

### Purpose ###

Return the illuminated fraction of the Moon at given Julian date(s).

### Arguments ###

* `jd`: the Julian ephemeris date.

### Output ###

The illuminated fraction \$k\$ of Moon's disk, with \$0 \\leq k \\leq 1\$. \$k
= 0\$ indicates a new moon, while \$k = 1\$ stands for a full moon.

### Method ###

Algorithm from Chapter 46 of "Astronomical Algorithms" by Jean Meeus
(Willmann-Bell, Richmond) 1991.  `sunpos` and `moonpos` are used to get
positions of the Sun and the Moon, and the Moon distance.  The selenocentric
elongation of the Earth from the Sun (phase angle) is then computed, and used to
determine the illuminated fraction.

### Example ###

Plot the illuminated fraction of the Moon for every day in January 2018 with a
hourly sampling.  Use [PyPlot.jl](https://github.com/JuliaPlots/Plots.jl/) for
plotting

```julia
using PyPlot
points = DateTime(2018,01,01):Dates.Hour(1):DateTime(2018,01,31,23,59,59);
plot(points, mphase.(jdcnv.(points)))
```

Note that in this calendar month there are two full moons, this event is called
[blue moon](https://en.wikipedia.org/wiki/Blue_moon).

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
mphase(jd::Real) = mphase(float(jd))


text:
"""
    Vec3D{T} <: FieldVector{3,T}

A custom type for mutable static arrays. It completely resembles a 3D point
object in some geometric space.

# Fields
- `x::T`: first field of `Vec3D`
- `y::T`: second field of `Vec3D`
- `z::T`: third field of `Vec3D`
"""
mutable struct Vec3D{T} <: FieldVector{3,T}
    x::T
    y::T
    z::T
end

"Define `similar_type` to preserve array type through array operations."
StaticArrays.similar_type(::Type{<:Vec3D}, ::Type{T}, s::Size{(3,)}) where {T} = Vec3D{T}

"""
    System{V<:Real}

Holds all relevant information for the simulation system.

# Fields
- `xpos::VT`: stores the array for the positions of the particles.
- `density::T`: the density of the system.
- `temperature::T`: the temperature of the system.
- `box::CellListMap.Box{UnitCellType,N,T,M}`: simulation box that will get updated as the
    simulation runs.
- `cutoff::T`: cutoff radius for the molecular interactions
- `rng::Random.AbstractRNG`: holds the RNG object for the system
- `npart::Int`: total number of particles in the system
"""
mutable struct System{B,T,VT,N,R}
    xpos::VT
    density::T
    temperature::T
    box::B
    rng::Random.AbstractRNG
    npart::N
    npartrange::R
end

function System(
    density::T, temp::T, particles::I, cutoff::T; dims=3, lcell=2
) where {T<:Real,I<:Int}
    box_size = cbrt(T(particles) / density)
    box = CellListMap.Box(fill(box_size, dims), cutoff; lcell=lcell)
    rng = Xorshifts.Xoroshiro128Plus()
    xpos = _initialize_positions(box_size, rng, particles)
    syst = System(xpos, density, temp, box, rng, particles, collect(1:particles))

    return syst
end

function _initialize_positions(box_size::T, rng, particles) where {T<:Real}
    range = (zero(T), box_size)
    xpos = [random_vec(Vec3D{T}, range; rng=rng) for _ in 1:particles]

    return xpos
end

"""
    Simulation

A type to hold all relevant information for the simulation.

# Fields
- `system::S`: ideally, a `System` type that has all the information about the positions
    and physical information on the system.
- `ensemble::E`: this is the _statistical ensemble_ used for doing the simulation, ideally
    of type `Ensemble`.
- `potential::P`: the interaction potential between the particles, ideally an object of
"""
struct Simulation{S,E,P}
    system::S
    ensemble::E
    potential::P
end


text:
@testset "5.2 Inverse cosine" begin
    include("5.2.2 (d x)^m (a+b arccos(c x))^n.jl")
    include("5.2.4 (f x)^m (d+e x^2)^p (a+b arccos(c x))^n.jl")
    include("5.2.5 Inverse cosine functions.jl")
end

max_stars_repo_path: test/5 Inverse trig functions/5.2 Inverse cosine/.jl


text:
# ------------------------------------------------------------------------------------------
# # Multiple dispatch
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# In this notebook we'll explore **multiple dispatch**, which is a key feature of Julia.
#
# Multiple dispatch makes software *generic* and *fast*!
#
# #### Starting with the familiar
#
# To understand multiple dispatch in Julia, let's start with what we've already seen.
#
# We can declare functions in Julia without giving Julia any information about the types of
# the input arguments that function will receive:
# ------------------------------------------------------------------------------------------

f(x) = x^2

# ------------------------------------------------------------------------------------------
# and then Julia will determine on its own which input argument types make sense and which
# do not:
# ------------------------------------------------------------------------------------------

f(10)

f([1, 2, 3])

# ------------------------------------------------------------------------------------------
# #### Specifying the types of our input arguments
#
# However, we also have the *option* to tell Julia explicitly what types our input arguments
# are allowed to have.
#
# For example, let's write a function `foo` that only takes strings as inputs.
# ------------------------------------------------------------------------------------------

foo(x::String, y::String) = println("My inputs x and y are both strings!")

# ------------------------------------------------------------------------------------------
# We see here that in order to restrict the type of `x` and `y` to `String`s, we just follow
# the input argument name by a double colon and the keyword `String`.
#
# Now we'll see that `foo` works on `String`s and doesn't work on other input argument
# types.
# ------------------------------------------------------------------------------------------

foo("hello", "hi!")

foo(3, 4)

# ------------------------------------------------------------------------------------------
# To get `foo` to work on integer (`Int`) inputs, let's tack `::Int` onto our input
# arguments when we declare `foo`.
# ------------------------------------------------------------------------------------------

foo(x::Int, y::Int) = println("My inputs x and y are both integers!")

foo(3, 4)

# ------------------------------------------------------------------------------------------
# Now `foo` works on integers! But look, `foo` also still works when `x` and `y` are
# strings!
# ------------------------------------------------------------------------------------------

foo("hello", "hi!")

# ------------------------------------------------------------------------------------------
# This is starting to get to the heart of multiple dispatch. When we declared
#
# ```julia
# foo(x::Int, y::Int) = println("My inputs x and y are both integers!")
# ```
# we didn't overwrite or replace
# ```julia
# foo(y::String, y::String)```
#
# Instead, we just added an additional ***method*** to the ***generic function*** called
# `foo`.
#
# A ***generic function*** is the abstract concept associated with a particular operation.
#
# For example, the generic function `+` represents the concept of addition.
#
# A ***method*** is a specific implementation of a generic function for *particular argument
# types*.
#
# For example, `+` has methods that accept floating point numbers, integers, matrices, etc.
#
# We can use the `methods` to see how many methods there are for `foo`.
# ------------------------------------------------------------------------------------------

methods(foo)

methods(+)

# ------------------------------------------------------------------------------------------
# So, we now can call `foo` on integers or strings. When you call `foo` on a particular set
# of arguments, Julia will infer the types of the inputs and dispatch the appropriate
# method. *This* is multiple dispatch.
#
# Multiple dispatch makes our code generic and fast. Our code can be generic and flexible
# because we can write code in terms of abstract operations such as addition and
# multiplication, rather than in terms of specific implementations. At the same time, our
# code runs quickly because Julia is able to call efficient methods for the relevant types.
#
# To see which method is being dispatched when we call a generic function, we can use the
# @which macro:
# ------------------------------------------------------------------------------------------

@which foo(3, 4)

@which 3.0 + 3.0

# ------------------------------------------------------------------------------------------
# Given that a method written specifically for floating point numbers is dispatched on `3.0
# + 3.0`, the LLVM code generated is extremely terse:
# ------------------------------------------------------------------------------------------

@code_llvm 3.0 + 3.0

# ------------------------------------------------------------------------------------------
# Note that Julia is fast even when we write generic function definitions because, at the
# end of the day, specific/tailored methods are called under the hood.
#
# For example, note that we can declare the adding function `myadd` without providing any
# type annotations -
# ------------------------------------------------------------------------------------------

myadd(x, y) = x + y

# ------------------------------------------------------------------------------------------
# and though we haven't constrained the types of `x` and `y`, we'll see that the LLVM code
# generated for `myadd(3.0, 3.0)` looks like that of `3.0 + 3.0`
# ------------------------------------------------------------------------------------------

@code_llvm myadd(3.0, 3.0)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
function perm(A)
  m, n = size(A)
  if m != n; throw(ArgumentError("permanent is for square matrices only")); end
  sum(σ -> prod(i -> A[i,σ[i]], 1:n), permutations(1:n))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# The license for this code is available at https://github.com/avik-pal/FastStyleTransfer.jl/blob/master/LICENSE.md

#-------------------Instance Normalization-----------------------------------

# NOTE: The Instance Normalization code is slow and can act as a huge bottleneck.
# Hence until this issue is fixed we shall be using BatchNorm

struct InstanceNorm
    β
    γ
end

Flux.treelike(InstanceNorm)

InstanceNorm(chs::Int; initβ = zeros, initγ = ones) = InstanceNorm(param(initβ(chs)), param(initγ(chs)))

function (IN::InstanceNorm)(x)
    local chs = length(IN.β.data)
    reshape(IN.γ, (1,1,chs,1)) .* ((x .- mean(x, [1,2,3])) ./ std(x, [1,2,3])) .+ reshape(IN.β, (1,1,chs,1))
end

#---------------------------Residual Block-----------------------------------

struct ResidualBlock
    conv_layers
    norm_layers
end

Flux.treelike(ResidualBlock)

ResidualBlock(chs::Int) =
   ResidualBlock((Conv((3,3), chs=>chs, pad = (1,1)), Conv((3,3), chs=>chs, pad = (1,1))), (BatchNorm(chs), BatchNorm(chs)))

function (r::ResidualBlock)(x)
    value = relu.(r.norm_layers[1](r.conv_layers[1](x)))
    r.norm_layers[2](r.conv_layers[2](value)) + x
end

#--------------------------Reflection Pad-------------------------------------

# Paper suggests using Reflection Padding. However normal padding is being used until this layer is implemented

struct ReflectionPad
    dim::Int
end

Flux.treelike(ReflectionPad)

#----------------------Convolution Block--------------------------------------

ConvBlock(chs::Pair{<:Int,<:Int}, kernel::Tuple{Int,Int}, stride::Tuple{Int,Int} = (1,1), pad::Tuple{Int,Int} = (0,0)) =
    Chain(Conv(kernel, chs, stride = stride, pad = pad), ReflectionPad(kernel[1]÷2))

#-------------------------Upsample--------------------------------------------

Upsample(x) = repeat(x, inner = (2,2,1,1))

#----------------------Upsampling BLock---------------------------------------

# TODO: Use reflection padding instead of zero padding once its implemented

UpsamplingBlock(chs::Pair{<:Int,<:Int}, kernel::Tuple{Int,Int}, stride::Tuple{Int,Int}, upsample::Int, pad::Tuple{Int,Int} = (0,0)) =
    Chain(Conv(kernel, chs, stride = stride, pad = (kernel[1]÷2, kernel[2]÷2)), x -> Upsample(x))

#---------------------Convolution Transpose-----------------------------------

function out_size(stride, pad, dilation, kernel, xdims, output_pad)
    dims = []
    for i in zip(stride, pad, dilation, kernel, xdims, output_pad)
        push!(dims, i[1] * (i[5] - 1) + (i[4] - 1) * i[3] - 2 * i[2] + 1 + i[6])
    end
    dims
end

function convtranspose(x, w; stride = 1, pad = 0, dilation = 1, output_pad = 0)
    stride, pad, dilation, output_pad = NNlib.padtuple(x, stride), NNlib.padtuple(x, pad), NNlib.padtuple(x, dilation), NNlib.padtuple(x, output_pad)
    y = similar(x, out_size(stride, pad, dilation, size(w)[1:end-2], size(x)[1:end-2], output_pad)...,size(w)[end-1],size(x)[end])
    NNlib.∇conv_data(x, y, w, stride = stride, pad = pad, dilation = dilation)
end

convtranspose(x::TrackedArray, w::TrackedArray; kw...) = track(convtranspose, x, w; kw...)
convtranspose(x::AbstractArray, w::TrackedArray; kw...) = track(convtranspose, x, w; kw...)
convtranspose(x::TrackedArray, w::AbstractArray; kw...) = track(convtranspose, x, w; kw...)

@grad convtranspose(x, w; kw...) =
    convtranspose(data.((x, w))...; kw...), Δ -> nobacksies(:convtranspose, (NNlib.conv(data.((Δ, w))...; kw[1:end-1]...), NNlib.∇conv_filter(data.((x, Δ, w))...; kw[1:end-1]...)))

struct ConvTranspose{N,F,A,V}
    σ::F
    weight::A
    bias::V
    stride::NTuple{N,Int}
    pad::NTuple{N,Int}
    dilation::NTuple{N,Int}
    output_pad::NTuple{N,Int}
end

ConvTranspose(w::AbstractArray{T,N}, b::AbstractVector{T}, σ = identity;
    stride = 1, pad = 0, dilation = 1, output_pad = 0) where {T,N} =
    ConvTranspose(σ, w, b, expand.(sub2(Val{N}), (stride, pad, dilation, output_pad))...)

ConvTranspose(k::NTuple{N,Integer}, ch::Pair{<:Integer,<:Integer}, σ = identity; init = Flux.initn,
    stride = 1, pad = 0, dilation = 1, output_pad = 0) where N =
    ConvTranspose(param(init(k..., ch[2], ch[1])), param(zeros(ch[2])), σ,
                  stride = stride, pad = pad, dilation = dilation, output_pad = output_pad)

Flux.treelike(ConvTranspose)

function (c::ConvTranspose)(x)
    σ, b = c.σ, reshape(c.bias, map(_->1, c.stride)..., :, 1)
    σ.(convtranspose(x, c.weight, stride = c.stride, pad = c.pad, dilation = c.dilation, output_pad = c.output_pad) .+ b)
end

function Base.show(io::IO, l::ConvTranspose)
    print(io, "ConvTranspose(", size(l.weight)[1:ndims(l.weight)-2])
    print(io, ", ", size(l.weight, ndims(l.weight)-1), "=>", size(l.weight, ndims(l.weight)))
    l.σ == identity || print(io, ", ", l.σ)
    print(io, ")")
end