text:
"""
    cohenES(test::String = "",size::String = "")

Give the conventional effect size (small, medium, large) for the tests available in this package.
Available options for `test` are `p`, `t`, `r`, `anova`, `chisq`, and `f2`.
Available options for `size` are `small`, `medium`, and `large`.
"""
function cohenES(;test::String = "",size::String = "")

    if lower(test) in ("p","t","r","anova","chisq","f2") == false
        error("`test` must be one of `p`, `t`, `r`, `anova`, `chisq`,or `f2`")
    end

    testd = Dict("p"=>1,"t"=>2,"r"=>3,"anova"=>4,"chisq"=>5,"f2"=>6)

    if lower(size) in ("small","medium","large") == false
        error("`size` must be one of `small`, `medium`, or `large`")
    end

    sized = Dict("small" => 1, "medium" => 2, "large" => 3)

    effsize = [[.2 .5 .8],
        [.2 .5 .8],
        [.1 .3 .5],
        [.1 .25 .4],
        [.1 .3 .5],
        [.02 .15 .35]]

    return effsize[testd[test]][sized[size]]
end

max_stars_repo_name: mwsohn/pwr

max_issues_repo_name: mwsohn/pwr

max_forks_repo_name: mwsohn/pwr


text:
# Structured matrices

#####
##### `Diagonal`
#####

function rrule(::Type{<:Diagonal}, d::AbstractVector)
    function Diagonal_pullback(ȳ::AbstractMatrix)
        return (NO_FIELDS, diag(ȳ))
    end
    function Diagonal_pullback(ȳ::Composite)
        # TODO: Assert about the primal type in the Composite, It should be Diagonal
        # infact it should be exactly the type of `Diagonal(d)`
        # but right now Zygote loses primal type information so we can't use it.
        # See https://github.com/FluxML/Zygote.jl/issues/603
        return (NO_FIELDS, ȳ.diag)
    end
    return Diagonal(d), Diagonal_pullback
end

function rrule(::typeof(diag), A::AbstractMatrix)
    function diag_pullback(ȳ)
        return (NO_FIELDS, @thunk(Diagonal(ȳ)))
    end
    return diag(A), diag_pullback
end

function rrule(::typeof(*), D::Diagonal{<:Real}, V::AbstractVector{<:Real})
    function times_pullback(Ȳ)
        return (NO_FIELDS, @thunk(Diagonal(Ȳ .* V)), @thunk(D * Ȳ))
    end
    return D * V, times_pullback
end

#####
##### `Symmetric`
#####

function rrule(::Type{<:Symmetric}, A::AbstractMatrix)
    function Symmetric_pullback(ȳ)
        return (NO_FIELDS, @thunk(_symmetric_back(ȳ)))
    end
    return Symmetric(A), Symmetric_pullback
end

_symmetric_back(ΔΩ) = UpperTriangular(ΔΩ) + LowerTriangular(ΔΩ)' - Diagonal(ΔΩ)
_symmetric_back(ΔΩ::Union{Diagonal,UpperTriangular}) = ΔΩ

#####
##### `Adjoint`
#####

# ✖✖✖TODO: Deal with complex-valued arrays as well
function rrule(::Type{<:Adjoint}, A::AbstractMatrix{<:Real})
    function Adjoint_pullback(ȳ)
        return (NO_FIELDS, @thunk(adjoint(ȳ)))
    end
    return Adjoint(A), Adjoint_pullback
end

function rrule(::Type{<:Adjoint}, A::AbstractVector{<:Real})
    function Adjoint_pullback(ȳ)
        return (NO_FIELDS, @thunk(vec(adjoint(ȳ))))
    end
    return Adjoint(A), Adjoint_pullback
end

function rrule(::typeof(adjoint), A::AbstractMatrix{<:Real})
    function adjoint_pullback(ȳ)
        return (NO_FIELDS, @thunk(adjoint(ȳ)))
    end
    return adjoint(A), adjoint_pullback
end

function rrule(::typeof(adjoint), A::AbstractVector{<:Real})
    function adjoint_pullback(ȳ)
        return (NO_FIELDS, @thunk(vec(adjoint(ȳ))))
    end
    return adjoint(A), adjoint_pullback
end

#####
##### `Transpose`
#####

function rrule(::Type{<:Transpose}, A::AbstractMatrix)
    function Transpose_pullback(ȳ)
        return (NO_FIELDS, @thunk transpose(ȳ))
    end
    return Transpose(A), Transpose_pullback
end

function rrule(::Type{<:Transpose}, A::AbstractVector)
    function Transpose_pullback(ȳ)
        return (NO_FIELDS, @thunk vec(transpose(ȳ)))
    end
    return Transpose(A), Transpose_pullback
end

function rrule(::typeof(transpose), A::AbstractMatrix)
    function transpose_pullback(ȳ)
        return (NO_FIELDS, @thunk transpose(ȳ))
    end
    return transpose(A), transpose_pullback
end

function rrule(::typeof(transpose), A::AbstractVector)
    function transpose_pullback(ȳ)
        return (NO_FIELDS, @thunk vec(transpose(ȳ)))
    end
    return transpose(A), transpose_pullback
end

#####
##### Triangular matrices
#####

function rrule(::Type{<:UpperTriangular}, A::AbstractMatrix)
    function UpperTriangular_pullback(ȳ)
        return (NO_FIELDS, @thunk Matrix(ȳ))
    end
    return UpperTriangular(A), UpperTriangular_pullback
end

function rrule(::Type{<:LowerTriangular}, A::AbstractMatrix)
    function LowerTriangular_pullback(ȳ)
        return (NO_FIELDS, @thunk Matrix(ȳ))
    end
    return LowerTriangular(A), LowerTriangular_pullback
end

ext: jl


text:
function solve_modelST_TR(H, g, nlp_stop, PData::PDataST, δ::T; cgtol::T = 0.1) where {T}
    # cas particulier Steihaug-Toint
    # ϵ = sqrt(eps(T)) # * 100.0 # old
    # cgtol = max(ϵ, min(cgtol, 9 * cgtol / 10, 0.01 * norm(g)^(1.0 + PData.ζ))) # old

    ζ, ξ, maxtol, mintol = PData.ζ, PData.ξ, PData.maxtol, PData.mintol
    n = length(g)
    gNorm2 = norm(g)
    # precision = max(1e-12, min(0.5, (gNorm2^ζ)))
    # Tolerance used in Assumption 2.6b in the paper ( ξ > 0, 0 < ζ ≤ 1 )
    cgatol = min(maxtol, ξ * gNorm2^(1 + ζ))
    cgatol = max(mintol, cgatol) # add some feasible limit
    cgrtol = min(maxtol, ξ * gNorm2^ζ)
    cgrtol = max(mintol, cgrtol) # add some feasible limit

    solver = PData.solver
    cg!(
        solver,
        H,
        -g,
        atol = cgatol,
        rtol = cgrtol,
        radius = δ,
        itmax = max(2 * n, 50),
        verbose = 0,
    )

    PData.d .= solver.x
    PData.OK = solver.stats.solved

    return PData.d, PData.λ
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""
    ObjsChannel(;
        accepted_user=nothing,
        created=nothing,
        creator=nothing,
        id=nothing,
        is_archived=nothing,
        is_channel=nothing,
        is_frozen=nothing,
        is_general=nothing,
        is_member=nothing,
        is_moved=nothing,
        is_mpim=nothing,
        is_non_threadable=nothing,
        is_org_shared=nothing,
        is_pending_ext_shared=nothing,
        is_private=nothing,
        is_read_only=nothing,
        is_shared=nothing,
        is_thread_only=nothing,
        last_read=nothing,
        latest=nothing,
        members=nothing,
        name=nothing,
        name_normalized=nothing,
        num_members=nothing,
        pending_shared=nothing,
        previous_names=nothing,
        priority=nothing,
        purpose=nothing,
        topic=nothing,
        unlinked=nothing,
        unread_count=nothing,
        unread_count_display=nothing,
    )

    - accepted_user::DefsUserId
    - created::Int32
    - creator::DefsUserId
    - id::DefsChannelId
    - is_archived::Bool
    - is_channel::Bool
    - is_frozen::Bool
    - is_general::Bool
    - is_member::Bool
    - is_moved::Int32
    - is_mpim::Bool
    - is_non_threadable::Bool
    - is_org_shared::Bool
    - is_pending_ext_shared::Bool
    - is_private::Bool
    - is_read_only::Bool
    - is_shared::Bool
    - is_thread_only::Bool
    - last_read::DefsTs
    - latest::Any
    - members::Vector{DefsUserId}
    - name::String
    - name_normalized::String
    - num_members::Int32
    - pending_shared::Vector{DefsTeam}
    - previous_names::Vector{DefsChannelName}
    - priority::Float32
    - purpose::ObjsChannelPurpose
    - topic::ObjsChannelPurpose
    - unlinked::Int32
    - unread_count::Int32
    - unread_count_display::Int32
"""
mutable struct ObjsChannel <: SwaggerModel
    accepted_user::Any # spec type: Union{ Nothing, DefsUserId } # spec name: accepted_user
    created::Any # spec type: Union{ Nothing, Int32 } # spec name: created
    creator::Any # spec type: Union{ Nothing, DefsUserId } # spec name: creator
    id::Any # spec type: Union{ Nothing, DefsChannelId } # spec name: id
    is_archived::Any # spec type: Union{ Nothing, Bool } # spec name: is_archived
    is_channel::Any # spec type: Union{ Nothing, Bool } # spec name: is_channel
    is_frozen::Any # spec type: Union{ Nothing, Bool } # spec name: is_frozen
    is_general::Any # spec type: Union{ Nothing, Bool } # spec name: is_general
    is_member::Any # spec type: Union{ Nothing, Bool } # spec name: is_member
    is_moved::Any # spec type: Union{ Nothing, Int32 } # spec name: is_moved
    is_mpim::Any # spec type: Union{ Nothing, Bool } # spec name: is_mpim
    is_non_threadable::Any # spec type: Union{ Nothing, Bool } # spec name: is_non_threadable
    is_org_shared::Any # spec type: Union{ Nothing, Bool } # spec name: is_org_shared
    is_pending_ext_shared::Any # spec type: Union{ Nothing, Bool } # spec name: is_pending_ext_shared
    is_private::Any # spec type: Union{ Nothing, Bool } # spec name: is_private
    is_read_only::Any # spec type: Union{ Nothing, Bool } # spec name: is_read_only
    is_shared::Any # spec type: Union{ Nothing, Bool } # spec name: is_shared
    is_thread_only::Any # spec type: Union{ Nothing, Bool } # spec name: is_thread_only
    last_read::Any # spec type: Union{ Nothing, DefsTs } # spec name: last_read
    latest::Any # spec type: Union{ Nothing, Any } # spec name: latest
    members::Any # spec type: Union{ Nothing, Vector{DefsUserId} } # spec name: members
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    name_normalized::Any # spec type: Union{ Nothing, String } # spec name: name_normalized
    num_members::Any # spec type: Union{ Nothing, Int32 } # spec name: num_members
    pending_shared::Any # spec type: Union{ Nothing, Vector{DefsTeam} } # spec name: pending_shared
    previous_names::Any # spec type: Union{ Nothing, Vector{DefsChannelName} } # spec name: previous_names
    priority::Any # spec type: Union{ Nothing, Float32 } # spec name: priority
    purpose::Any # spec type: Union{ Nothing, ObjsChannelPurpose } # spec name: purpose
    topic::Any # spec type: Union{ Nothing, ObjsChannelPurpose } # spec name: topic
    unlinked::Any # spec type: Union{ Nothing, Int32 } # spec name: unlinked
    unread_count::Any # spec type: Union{ Nothing, Int32 } # spec name: unread_count
    unread_count_display::Any # spec type: Union{ Nothing, Int32 } # spec name: unread_count_display

    function ObjsChannel(;accepted_user=nothing, created=nothing, creator=nothing, id=nothing, is_archived=nothing, is_channel=nothing, is_frozen=nothing, is_general=nothing, is_member=nothing, is_moved=nothing, is_mpim=nothing, is_non_threadable=nothing, is_org_shared=nothing, is_pending_ext_shared=nothing, is_private=nothing, is_read_only=nothing, is_shared=nothing, is_thread_only=nothing, last_read=nothing, latest=nothing, members=nothing, name=nothing, name_normalized=nothing, num_members=nothing, pending_shared=nothing, previous_names=nothing, priority=nothing, purpose=nothing, topic=nothing, unlinked=nothing, unread_count=nothing, unread_count_display=nothing)
        o = new()
        validate_property(ObjsChannel, Symbol("accepted_user"), accepted_user)
        setfield!(o, Symbol("accepted_user"), accepted_user)
        validate_property(ObjsChannel, Symbol("created"), created)
        setfield!(o, Symbol("created"), created)
        validate_property(ObjsChannel, Symbol("creator"), creator)
        setfield!(o, Symbol("creator"), creator)
        validate_property(ObjsChannel, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(ObjsChannel, Symbol("is_archived"), is_archived)
        setfield!(o, Symbol("is_archived"), is_archived)
        validate_property(ObjsChannel, Symbol("is_channel"), is_channel)
        setfield!(o, Symbol("is_channel"), is_channel)
        validate_property(ObjsChannel, Symbol("is_frozen"), is_frozen)
        setfield!(o, Symbol("is_frozen"), is_frozen)
        validate_property(ObjsChannel, Symbol("is_general"), is_general)
        setfield!(o, Symbol("is_general"), is_general)
        validate_property(ObjsChannel, Symbol("is_member"), is_member)
        setfield!(o, Symbol("is_member"), is_member)
        validate_property(ObjsChannel, Symbol("is_moved"), is_moved)
        setfield!(o, Symbol("is_moved"), is_moved)
        validate_property(ObjsChannel, Symbol("is_mpim"), is_mpim)
        setfield!(o, Symbol("is_mpim"), is_mpim)
        validate_property(ObjsChannel, Symbol("is_non_threadable"), is_non_threadable)
        setfield!(o, Symbol("is_non_threadable"), is_non_threadable)
        validate_property(ObjsChannel, Symbol("is_org_shared"), is_org_shared)
        setfield!(o, Symbol("is_org_shared"), is_org_shared)
        validate_property(ObjsChannel, Symbol("is_pending_ext_shared"), is_pending_ext_shared)
        setfield!(o, Symbol("is_pending_ext_shared"), is_pending_ext_shared)
        validate_property(ObjsChannel, Symbol("is_private"), is_private)
        setfield!(o, Symbol("is_private"), is_private)
        validate_property(ObjsChannel, Symbol("is_read_only"), is_read_only)
        setfield!(o, Symbol("is_read_only"), is_read_only)
        validate_property(ObjsChannel, Symbol("is_shared"), is_shared)
        setfield!(o, Symbol("is_shared"), is_shared)
        validate_property(ObjsChannel, Symbol("is_thread_only"), is_thread_only)
        setfield!(o, Symbol("is_thread_only"), is_thread_only)
        validate_property(ObjsChannel, Symbol("last_read"), last_read)
        setfield!(o, Symbol("last_read"), last_read)
        validate_property(ObjsChannel, Symbol("latest"), latest)
        setfield!(o, Symbol("latest"), latest)
        validate_property(ObjsChannel, Symbol("members"), members)
        setfield!(o, Symbol("members"), members)
        validate_property(ObjsChannel, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(ObjsChannel, Symbol("name_normalized"), name_normalized)
        setfield!(o, Symbol("name_normalized"), name_normalized)
        validate_property(ObjsChannel, Symbol("num_members"), num_members)
        setfield!(o, Symbol("num_members"), num_members)
        validate_property(ObjsChannel, Symbol("pending_shared"), pending_shared)
        setfield!(o, Symbol("pending_shared"), pending_shared)
        validate_property(ObjsChannel, Symbol("previous_names"), previous_names)
        setfield!(o, Symbol("previous_names"), previous_names)
        validate_property(ObjsChannel, Symbol("priority"), priority)
        setfield!(o, Symbol("priority"), priority)
        validate_property(ObjsChannel, Symbol("purpose"), purpose)
        setfield!(o, Symbol("purpose"), purpose)
        validate_property(ObjsChannel, Symbol("topic"), topic)
        setfield!(o, Symbol("topic"), topic)
        validate_property(ObjsChannel, Symbol("unlinked"), unlinked)
        setfield!(o, Symbol("unlinked"), unlinked)
        validate_property(ObjsChannel, Symbol("unread_count"), unread_count)
        setfield!(o, Symbol("unread_count"), unread_count)
        validate_property(ObjsChannel, Symbol("unread_count_display"), unread_count_display)
        setfield!(o, Symbol("unread_count_display"), unread_count_display)
        o
    end
end # type ObjsChannel

const _property_map_ObjsChannel = Dict{Symbol,Symbol}(Symbol("accepted_user")=>Symbol("accepted_user"), Symbol("created")=>Symbol("created"), Symbol("creator")=>Symbol("creator"), Symbol("id")=>Symbol("id"), Symbol("is_archived")=>Symbol("is_archived"), Symbol("is_channel")=>Symbol("is_channel"), Symbol("is_frozen")=>Symbol("is_frozen"), Symbol("is_general")=>Symbol("is_general"), Symbol("is_member")=>Symbol("is_member"), Symbol("is_moved")=>Symbol("is_moved"), Symbol("is_mpim")=>Symbol("is_mpim"), Symbol("is_non_threadable")=>Symbol("is_non_threadable"), Symbol("is_org_shared")=>Symbol("is_org_shared"), Symbol("is_pending_ext_shared")=>Symbol("is_pending_ext_shared"), Symbol("is_private")=>Symbol("is_private"), Symbol("is_read_only")=>Symbol("is_read_only"), Symbol("is_shared")=>Symbol("is_shared"), Symbol("is_thread_only")=>Symbol("is_thread_only"), Symbol("last_read")=>Symbol("last_read"), Symbol("latest")=>Symbol("latest"), Symbol("members")=>Symbol("members"), Symbol("name")=>Symbol("name"), Symbol("name_normalized")=>Symbol("name_normalized"), Symbol("num_members")=>Symbol("num_members"), Symbol("pending_shared")=>Symbol("pending_shared"), Symbol("previous_names")=>Symbol("previous_names"), Symbol("priority")=>Symbol("priority"), Symbol("purpose")=>Symbol("purpose"), Symbol("topic")=>Symbol("topic"), Symbol("unlinked")=>Symbol("unlinked"), Symbol("unread_count")=>Symbol("unread_count"), Symbol("unread_count_display")=>Symbol("unread_count_display"))
const _property_types_ObjsChannel = Dict{Symbol,String}(Symbol("accepted_user")=>"DefsUserId", Symbol("created")=>"Int32", Symbol("creator")=>"DefsUserId", Symbol("id")=>"DefsChannelId", Symbol("is_archived")=>"Bool", Symbol("is_channel")=>"Bool", Symbol("is_frozen")=>"Bool", Symbol("is_general")=>"Bool", Symbol("is_member")=>"Bool", Symbol("is_moved")=>"Int32", Symbol("is_mpim")=>"Bool", Symbol("is_non_threadable")=>"Bool", Symbol("is_org_shared")=>"Bool", Symbol("is_pending_ext_shared")=>"Bool", Symbol("is_private")=>"Bool", Symbol("is_read_only")=>"Bool", Symbol("is_shared")=>"Bool", Symbol("is_thread_only")=>"Bool", Symbol("last_read")=>"DefsTs", Symbol("latest")=>"Any", Symbol("members")=>"Vector{DefsUserId}", Symbol("name")=>"String", Symbol("name_normalized")=>"String", Symbol("num_members")=>"Int32", Symbol("pending_shared")=>"Vector{DefsTeam}", Symbol("previous_names")=>"Vector{DefsChannelName}", Symbol("priority")=>"Float32", Symbol("purpose")=>"ObjsChannelPurpose", Symbol("topic")=>"ObjsChannelPurpose", Symbol("unlinked")=>"Int32", Symbol("unread_count")=>"Int32", Symbol("unread_count_display")=>"Int32")
Base.propertynames(::Type{ ObjsChannel }) = collect(keys(_property_map_ObjsChannel))
Swagger.property_type(::Type{ ObjsChannel }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsChannel[name]))}
Swagger.field_name(::Type{ ObjsChannel }, property_name::Symbol) =  _property_map_ObjsChannel[property_name]

function check_required(o::ObjsChannel)
    (getproperty(o, Symbol("created")) === nothing) && (return false)
    (getproperty(o, Symbol("creator")) === nothing) && (return false)
    (getproperty(o, Symbol("id")) === nothing) && (return false)
    (getproperty(o, Symbol("is_channel")) === nothing) && (return false)
    (getproperty(o, Symbol("is_mpim")) === nothing) && (return false)
    (getproperty(o, Symbol("is_org_shared")) === nothing) && (return false)
    (getproperty(o, Symbol("is_private")) === nothing) && (return false)
    (getproperty(o, Symbol("is_shared")) === nothing) && (return false)
    (getproperty(o, Symbol("members")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("name_normalized")) === nothing) && (return false)
    (getproperty(o, Symbol("purpose")) === nothing) && (return false)
    (getproperty(o, Symbol("topic")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ObjsChannel }, name::Symbol, val)
end


text:
using GridGames

@gridgame(
    name=:Tac,
    grid=(3,3),
    players=[:X,:O],
    actions=[
        (3,1), (3,2), (3,3),
        (2,1), (2,2), (2,3),
        (1,1), (1,2), (1,3),
    ],
    pieces=[:X,:O],
)

function Tac()
    board = emptyboard(Tac)
    player = getplayer(Tac, :X)
    states = [TacState(board,player)]
    Tac([],[],states)
end

function GridGames.getnextstate(s::TacState, a::TacAction)
    next_state = TacState(
        copy(s.board),
        getnexttoplay(s)
    )
    next_state[a.key...] = getpieceidx(Tac, s.to_play.name)
    next_state
end

function GridGames.getlegalactions(s::TacState)
    if isgamewon(s)
        return []
    end
    map(
        cell->getaction(s,(cell.row,cell.col)),
        emptycells(s)
    )
end

function GridGames.getwinner(s::TacState)
    function getwinner(cells)
        piece = getmatchingpiece(s, cells)
        if !isnothing(piece)
            return getplayer(s, piece)
        end
        nothing
    end
    tests = []
    append!(tests, [s[row, 1:3] for row in 1:3])
    append!(tests, [s[1:3, col] for col in 1:3])
    push!(tests, [s[1,1], s[2,2], s[3,3]])
    push!(tests, [s[3,1], s[2,2], s[1,3]])
    for test in tests
        player = getwinner(test)
        if !isnothing(player)
            return player
        end
    end
    nothing
end


text:
@inline function countlines_filename(filename::AbstractString)
    num_lines = open(filename, "r") do io
        return countlines(io)
    end
    return num_lines
end


text:
module gs2ls
using NetCDF

export createInputFile, GS2species, GS2resolution, GS2params, GS2geometry, timeavg

"""
GS2species: an encapsulation of data related to each species in a gs2 input file
"""
type GS2species
   z::Float64
   mass::Float64
   dens::Float64
   temp::Float64
   fprim::Float64
   tprim::Float64
   vnewk::Float64
   electron::Bool
end
GS2species(;mass=1.0,z=1.0,dens=1.0,temp=1.0,vnewk=0.0,fprim=0.0,tprim=0.0,electron=false) =GS2species(z,mass,dens,temp,fprim,tprim,vnewk,electron)

"""
GS2resolution(ntheta,negrid,ngauss,delt,linear,kx,ky,nperiod,nx,ny,y0)
"""
type GS2resolution
   # Resolution parameters
   ntheta::Int64
   negrid::Int64
   ngauss::Int64
   delt::Float64
   linear::Bool
   kx::Float64
   ky::Float64
   nperiod::Int64
   nx::Int64
   ny::Int64
   y0::Float64
end
# Set defaults
GS2resolution(;ntheta=24,negrid=16,ngauss=8,delt=0.1,linear=true,kx=0.0,ky=0.5,nperiod=1,nx=-1,ny=64,y0=10.0) = GS2resolution(ntheta,negrid,ngauss,delt,linear,kx,ky,nperiod,nx,ny,y0)

"""
GS2params(runname,cpu_time,fapar,fbpar,beta,tite,fluxe,fluxemu,writephi)
"""
type GS2params
   # Simulation and diagnostic parameters
   runname::String
   cpu_time::Float64
   fapar::Float64
   fbpar::Float64
   beta::Float64
   tite::Float64
   fluxe::Bool
   fluxemu::Bool
   writephi::Bool
end
GS2params(;runname="testcase",cpu_time=3600.0,fapar=0.0,fbpar=0.0,beta=-0.01,tite=1.0,fluxe=true,fluxemu=false,writephi=false) = GS2params(runname,cpu_time,fapar,fbpar,beta,tite,fluxe,fluxemu,writephi)

"""
GS2geometry(equilibrium_option,rhoc,irho,qinp,shat,akappa,kappri,tri,tripri,rmaj,r_geo)
"""
type GS2geometry
   # Geometry parameters
   equilibrium_option::String
   rhoc::Float64
   irho::Int64
   qinp::Float64
   shat::Float64
   shift::Float64
   akappa::Float64
   akappri::Float64
   tri::Float64
   tripri::Float64
   rmaj::Float64
   r_geo::Float64
end
GS2geometry(;equilibrium_option="eik",rhoc="0.5",irho=2,qinp=1.0,shat=1.0,shift=0.0,akappa=0.0,akappri=0.0,tri=0.0,tripri=0.0,rmaj=3.0,r_geo=-1.0)=GS2geometry(equilibrium_option,rhoc,irho,qinp,shat,shift,akappa,akappri,tri,tripri,rmaj,r_geo)

function timeavg(time,y,tavgfrac=-1.0)

   tavgmin = time[1] + (1.0-tavgfrac)*(time[end]-time[1])
   itmin = indmin(abs(time-tavgmin))
   return  sum( 0.5*(time[itmin+1:end]-time[itmin:end-1]).*(y[itmin+1:end]+y[itmin:end-1]))/(time[end]-time[itmin])
end

   # internal logic: linear

   # required always: ntheta, p.rhoc, equilibrium_option,shat, qinp, ngauss, negrid, fapar, fbpar, delt, cpu_time, spec, runname, fluxe, fluxemu,writephi

   # optional parameters: beta, tite, ky, kx, nx,ny, y0, nperiod, irho, rmaj, akappa, akappri,tri,tripri

"""
writeNamelist(file::IOStream,namelist::String,dict::Dict)

Takes a julia dictionary and string, and writes a fortran namelist to an already opened file.
"""
function writeNamelist(file::IOStream,namelist::String,dict::Dict)
   write(file,"\n\n&"*namelist*"\n")
   for key in keys(dict)
      if typeof(dict[key]) == String
         line = " "*key*"=\""*dict[key]*"\""
      elseif typeof(dict[key]) == Bool
         line = " "*key*"=."*string(dict[key])*"."
      else
         line = " "*key*"="*string(dict[key])
      end
      write(file,line*"\n")
   end
   write(file,"/")
end

"""
createInputFile(filename::String,p::GS2params,r::GS2resolution,g::GS2geometry,spec::Array{GS2species,1})

Writes an individual GS2 input file base on provided data and internal logic
"""
function createInputFile(filename::String,p::GS2params,r::GS2resolution,g::GS2geometry,spec::Array{GS2species,1})

   # Check for defaults
   if g.r_geo < 0.0
      g.r_geo = g.rmaj
   end
   if r.nx < 0
      r.nx = r.ny
   end
   if abs(p.fapar) > 0.0 && p.beta < 0.0
      error("Must specify beta when fapar>0")
   end
   if g.equilibrium_option == "s-alpha" && p.beta < 0.0
      warn("shift set equal to zero because beta=0.0")
      beta = 0.0
   end
   if !r.linear
      # TODO: Add logic to guess nx based on given nperiod for nonlinear runs
      r.nperiod = 1
   end

   # Open file
   f = open(filename,"w")
   write(f,"! GS2 input file generated by gs2ls.jl\n")

   # Calculate derived quantities
   nspec = length(spec)

   maxmass = 0.0
   for is in 1:nspec
      maxmass = max(maxmass,spec[is].mass)
   end
   for is in 1:nspec
      if (spec[is].mass < maxmass*1.0e-3) && (spec[is].z == -1.0)
         spec[is].electron = true
         ne = spec[is].dens
      end
   end
   zeff = 0.0
   ne = 1.0
   for is in 1:nspec
      if spec[is].electron
         ne = spec[is].dens
      else
         zeff = zeff + spec[is].z^2*spec[is].dens
      end
   end
   zeff = zeff/ne

   if r.linear
      grid_option = "single"
      boundary_option="default"
   else
      grid_option = "box"
      boundary_option="linked"
   end

   eps = g.rhoc/g.rmaj
   epsl = 2.0/g.rmaj
   pk = epsl/g.qinp

   # Calculate pressure gradient and shift (for salpha)
   prestot = 0.0
   prestot_prime = 0.0
   totcharge = 0.0
   totcharge_prim = 0.0
   for is in 1:nspec
      prestot = prestot + spec[is].dens*spec[is].temp
      prestot_prime = prestot_prime - spec[is].dens*spec[is].temp.*(spec[is].fprim + spec[is].tprim)
      totcharge = totcharge + spec[is].dens.*spec[is].z
      totcharge_prim = totcharge_prim + spec[is].dens.*spec[is].z.*spec[is].fprim
   end
   beta_prime = p.beta * prestot_prime/prestot

   if g.equilibrium_option == "s-alpha"
      g.shift = -g.qinp^2 * rmaj * beta_prime
   end

   # Check quasineutrality
   abs(totcharge) < 1.0e-3 || warn("Equilibrium quasineutrality not satisfied")
   abs(totcharge_prim) < 1.0e-3 || warn("Equilibrium quasineutrality gradient not satisfied")

   # Enforce electronhood

   ########################################
   # parameters
   writeNamelist(f, "parameters", Dict("beta"=>p.beta,"zeff"=>zeff,"tite"=>p.tite))

   ########################################
   # kt_grids_knobs
   writeNamelist(f, "kt_grids_knobs", Dict("grid_option"=>grid_option))

   ########################################
   # kt_grids_single_parameters
   # kt_grids_box_parameters
   if r.linear
      dict = Dict("aky"=>r.ky,"akx"=>r.kx,"theta0"=>0.0)
      writeNamelist(f, "kt_grids_single_parameters", dict)
   else
      jtwist = round(Int64,2.0*pi*g.shat)
      dict = Dict("nx"=>r.nx,"ny"=>r.ny,"jtwist"=>jtwist,"y0"=>r.y0,"x0"=>r.y0)
      writeNamelist(f, "kt_grids_box_parameters", dict)
   end

   ########################################
   # theta_grid_knobs
   writeNamelist(f, "theta_grid_knobs", Dict("equilibrium_option"=>g.equilibrium_option))

   ########################################
   # theta_grid_parameters
   # theta_grid_salpha_knobs
   # theta_grid_eik_knobs
   if g.equilibrium_option == "s-alpha"
      dict = Dict("ntheta"=>r.ntheta,"nperiod"=>r.nperiod,"shat"=>g.shat,"shift",g.shift,"eps"=>eps,"epsl"=>epsl,"pk",pk)
      writeNamelist(f, "theta_grid_parameters", dict)

      writeNamelist(f, "theta_grid_salpha_knobs", Dict("model_option"=>"s-alpha"))
   elseif g.equilibrium_option == "eik"

      dict = Dict("ntheta"=>r.ntheta,"nperiod"=>r.nperiod,"rhoc"=>g.rhoc,"qinp"=>g.qinp,"shat"=>g.shat,"shift"=>g.shift,"rmaj"=>g.rmaj,"r_geo"=>g.r_geo,"akappa"=>g.akappa,"akappri"=>g.akappri,"tri"=>g.tri,"tripri"=>g.tripri)
      writeNamelist(f, "theta_grid_parameters", dict)

      dict = Dict("iflux"=>0,"itor"=>1,"irho"=>g.irho,"local_eq"=>true,"bishop"=>4,"delrho"=>0.001,"s_hat_input"=>g.shat,"beta_prime_input"=>beta_prime)
      writeNamelist(f, "theta_grid_eik_knobs", dict)
   else
      error("equilibrium_option invalid")
      # TODO: Add numerical geom options
   end

   ########################################
   # le_grids_knobs
   writeNamelist(f, "le_grids_knobs", Dict("ngauss"=>r.ngauss,"negrid"=>r.negrid))

   ########################################
   # dist_fn_knobs
   writeNamelist(f, "dist_fn_knobs", Dict("boundary_option"=>boundary_option,"nonad_zero"=>true))

   ########################################
   # fields_knobs
   writeNamelist(f, "fields_knobs", Dict("field_option"=>"implicit"))

   ########################################
   # knobs
   # TODO: Guess timestep and cpu_time
   dict = Dict("fphi"=>1.0,"fapar"=>p.fapar,"fbpar"=>p.fbpar,"delt"=>r.delt,"nstep"=>1000000,"avail_cpu_time"=>p.cpu_time)
   writeNamelist(f, "knobs", dict)

   ########################################
   # reinit_knobs
   writeNamelist(f, "reinit_knobs", Dict("delt_adj"=>4.0,"delt_minimum"=>1.0e-8,"abort_rapid_time_step_change"=>true))

   ########################################
   # layouts_knobs
   writeNamelist(f, "layouts_knobs", Dict("layout"=>"xyles"))

   ########################################
   # collisions_knobs
   writeNamelist(f, "collisions_knobs", Dict("collision_model"=>"default"))

   ########################################
   # nonlinear_terms_knobs
   if r.linear
      writeNamelist(f, "nonlinear_terms_knobs", Dict("nonlinear_mode"=>"off","flow_mode"=>"off","cfl"=>0.5))
   else
      writeNamelist(f, "nonlinear_terms_knobs", Dict("nonlinear_mode"=>"on","flow_mode"=>"off","cfl"=>0.5))
   end

   ########################################
   # species_knobs
   # dist_fn_species_knobs
   writeNamelist(f, "species_knobs", Dict("nspec"=>nspec))

   for is in 1:nspec
      if spec[is].electron
         spectype = "electron"
      else
         spectype = "ion"
      end

      dict = Dict("z"=>spec[is].z,"mass"=>spec[is].mass,"dens"=>spec[is].dens,"temp"=>spec[is].temp,"fprim"=>spec[is].fprim,"tprim"=>spec[is].tprim,"vnewk"=>spec[is].vnewk,"type"=>spectype)
      writeNamelist(f, "species_parameters_"*string(is), dict)

      writeNamelist(f, "dist_fn_species_knobs_"*string(is), Dict("fexpr"=>0.45,"bakdif"=>0.05))
   end

   ########################################
   # init_g_knobs
   dict = Dict("phiinit"=>0.001,"restart_file"=>p.runname*".nc","ginit_option"=>"noise","restart_dir"=>"restartdir")
   writeNamelist(f, "init_g_knobs",dict)

   ########################################
   # gs2_diagnostics_knobs
   dict = Dict{String,Any}()
   dict["write_flux_emu"] = p.fluxemu
   dict["print_flux_line"] = false
   dict["write_flux_e"] = p.fluxe
   dict["write_eflux"] = p.fluxe

   if r.linear
      dict["write_nl_flux"] = false
      dict["print_line"] = true
      dict["write_line"] = true
      dict["write_omega"] = true
      dict["write_omavg"] = true
      dict["omegatol"] = 0.0001
      dict["save_for_restart"] = false
   else
      dict["write_nl_flux"] = true
      dict["print_line"] = false
      dict["write_line"] = false
      dict["write_omega"] = false
      dict["write_omavg"] = false
      dict["omegatol"] = -0.0001
      dict["save_for_restart"] = true
   end
   dict["write_hrate"] = false
   dict["write_eigenfunc"] = true
   dict["write_final_fields"] = true
   dict["write_final_moments"] = true
   dict["nsave"] = 500
   dict["nwrite"] = 50
   dict["navg"] = 50
   dict["omegatinst"] = 500.0
   dict["write_phi_over_time"] = p.writephi

   writeNamelist(f,"gs2_diagnostics_knobs",dict)

   # Close file
   close(f)

end

end

max_stars_repo_name: gjwilkie/PlasmaTools

max_issues_repo_name: gjwilkie/PlasmaTools

max_forks_repo_name: gjwilkie/PlasmaTools


text:
using ModelingToolkit, OrdinaryDiffEq
using AnalogMTK
using AnalogMTK.EasyPlot #set, cons
import Printf: @sprintf
j=im

#==Circuit parameters
===============================================================================#
𝑓 = 60 #Hz
Tsim = 1 #sec

#@named Vsrc = VSource(1.0)
@named Vsrc = VSource(Sine(A=1, f=𝑓))
@named ZS = Resistor(R=100)
@named RL = Resistor(R=100)
@named CL = Capacitor(C=10e-9) #10e-6, 10e-9, 47e-9

#==Basic calculations & input validation
===============================================================================#
#=
𝜔 = 2π*𝑓
ZC = 1/(j*𝜔*C)

println()
@info("Simulation parameters:")
@show V, 𝑓, 𝜔
@show R, C, abs(ZC)

@show 𝜏 = R*C
=#

#==Circuit definition
===============================================================================#
#@named D1 = Diode(Ideal())
@named D1 = Diode(Shockley(), IS=1e-3)
#@named D1 = Resistor(R=1)
@named GND = Ground()

rc_eqs = [
	connect(Vsrc.p, ZS.p)
	connect(ZS.n, D1.p)
	connect(D1.n, RL.p, CL.p)
	connect(GND.g, Vsrc.n, RL.n, CL.n)
]

systems = [GND, Vsrc, ZS, D1, RL, CL]
@named rc_model = ODESystem(rc_eqs, AnalogMTK.t; systems)

println()
@info("Initial ODESystem (model definition):"); flush(stdout); flush(stderr)
display(rc_model)

println()
@info("Computing simplified model (ODESystem)..."); flush(stdout); flush(stderr)
#sys = rc_model
sys = structural_simplify(rc_model)
display(sys)
println("\tdone.")

#==Circuit simulation
===============================================================================#
println()
@info("Running simulation..."); flush(stdout); flush(stderr)
u0 = [
      ZS.v => 0.0,
      CL.v => 0.0,
      CL.p.i => 0.0,
      D1.v => 0.0,
     ]
prob = ODEProblem(sys, u0, (0, Float64(Tsim)))
sol = solve(prob, Trapezoid())
	#Trapezoid(), ImplicitEuler()
	#Rodas4(), RadauIIA5(), Tsit5()
@show length(sol.t)
println("\tdone.")

#==Plot results
===============================================================================#
println()
@info("Plotting results..."); flush(stdout); flush(stderr)
ylabels = ["Voltage [V]", "Current [A]"]
plot = plot_transient(sol.t, sol[Vsrc.p.v], id="Vsrc"; ylabels)
plot_transient!(plot, sol.t, sol[D1.v], id="D1.v")
plot_transient!(plot, sol.t, sol[CL.p.v], id="Vout")
plot_transient!(plot, sol.t, sol[RL.p.i], id="Iout", strip=2)
plot_transient!(plot, sol.t, sol[D1.i], id="D1.i", strip=2)
plot_transient!(plot, sol.t, sol[Vsrc.i], id="Vsrc.i", strip=2)
set(plot, ystrip1=set(min=-2, max=2))
displaygui(plot)
println("\tdone.")

:TEST_DONE


text:
using PyTensorFlow
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

if Sys.islinux()
py"""
import tensorflow as tf
libUpwpsOp = tf.load_op_library('build/libUpwpsOp.so')
@tf.custom_gradient
def upwps_op(permi,mobi,src,funcref,h,rhograv,index):
    pres = libUpwpsOp.upwps_op(permi,mobi,src,funcref,h,rhograv,index)
    def grad(dy):
        return libUpwpsOp.upwps_op_grad(dy, pres, permi,mobi,src,funcref,h,rhograv,index)
    return pres, grad
"""
elseif Sys.isapple()
py"""
import tensorflow as tf
libUpwpsOp = tf.load_op_library('build/libUpwpsOp.dylib')
@tf.custom_gradient
def upwps_op(permi,mobi,src,funcref,h,rhograv,index):
    pres = libUpwpsOp.upwps_op(permi,mobi,src,funcref,h,rhograv,index)
    def grad(dy):
        return libUpwpsOp.upwps_op_grad(dy, pres, permi,mobi,src,funcref,h,rhograv,index)
    return pres, grad
"""
elseif Sys.iswindows()
py"""
import tensorflow as tf
libUpwpsOp = tf.load_op_library('build/libUpwpsOp.dll')
@tf.custom_gradient
def upwps_op(permi,mobi,src,funcref,h,rhograv,index):
    pres = libUpwpsOp.upwps_op(permi,mobi,src,funcref,h,rhograv,index)
    def grad(dy):
        return libUpwpsOp.upwps_op_grad(dy, pres, permi,mobi,src,funcref,h,rhograv,index)
    return pres, grad
"""
end

upwps_op = py"upwps_op"

len_z = 4
len_x = 4

rho = 1000.0
G = 9.8

nScale = 5
tf_g = Array{Any}(undef, nScale)
tf_coef = Array{Any}(undef, nScale)
tf_h = Array{Any}(undef, nScale)
tf_p = Array{Any}(undef, nScale)
p_true_array = Array{Any}(undef, nScale)
p_inv_array = Array{Any}(undef, nScale)

h_array = @. 1 / 2^(1:nScale)

for iScale = 1:nScale
    h = h_array[iScale]
    nz = Int(len_z/h + 1)
    nx = Int(len_x/h + 1)
    g = zeros(nz, nx)
    coef = zeros(nz, nx)
    p_true = zeros(nz, nx)
    for i = 1:nz
        for j = 1:nx
            g[i,j] = -sin(2*pi/len_z*(i-1)*h) * sin(2*pi/len_x*(j-1)*h)
            coef[i,j] = 1.0 - cos(2*pi/len_z*(i-1)*h) * sin(2*pi/len_x*(j-1)*h) * len_z / (2*pi*rho*G)

            # g[i,j] = 2.0*(i-1)*h*exp(-(((i-1)*h)^2) -(((j-1)*h)^2)) * rho * G
            # coef[i,j] = 1.0 + exp(-(((i-1)*h)^2) -(((j-1)*h)^2))

            p_true[i,j] = rho*G*(i-1)*h
        end
    end
    p_true_array[iScale] = p_true .- mean(p_true)
    # p_true_array[iScale] = p_true
    tf_g[iScale] = constant(g)
    tf_coef[iScale] = constant(coef)
    tf_h[iScale] = constant(h)
    tf_ones = constant(ones(nz, nx))
    tf_p[iScale] = upwps_op(tf_coef[iScale], tf_ones, tf_g[iScale], tf_ones, tf_h[iScale], constant(rho*G), constant(0))
end

sess = Session()
init(sess)
p_inv_array = run(sess, tf_p)
for iScale = 1:nScale
    p_inv_array[iScale] = p_inv_array[iScale] .- mean(p_inv_array[iScale])
end

function l2_error(p_true, p_inv, iScale)
    l2_error = 0.0
    l2_norm = 0.0
    h = h_array[iScale]
    nz = size(p_true)[1]
    nx = size(p_true)[2]
    for i = 1:nz
        for j = 1:nx
            l2_error += (p_true[i,j]-p_inv[i,j])^2 * h^2
            l2_norm += (p_true[i,j])^2 * h^2
        end
    end
    return sqrt(l2_error)/sqrt(l2_norm)
end

Error_array = Array{Any}(undef, nScale)
for iScale = 1:nScale
    Error_array[iScale] = l2_error(p_true_array[iScale], p_inv_array[iScale], iScale)
end

loglog(h_array, Error_array, "*-", label="MMS convergence")
loglog(h_array, h_array.^2 * 0.5*Error_array[1]/h_array[1]^2, "--",label="\$\\mathcal{O}(h^2)\$")
loglog(h_array, h_array * 0.5*Error_array[1]/h_array[1], "-",label="\$\\mathcal{O}(h)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$h\$")
ylabel("Error")

# imshow(p)


text:
__precompile__(true)
module YourSubmodule

export your_function
"""
Write something about the purpose of the function up top.

Arguments:
- `arg_1::Integer` Desribute first function argument. [input_units]
- `arg_2::Integer` Desribute second function argument. [input_units]

Returns:
- `ret_1::Float64` Describe first return value. [input_units]
- `ret_2::Float64` Describe second return value. [input_units]

Notes:
1. Make specific comments about corner cases or things to know about usage here.

References:
1. List citations to references used in writing the function here
"""
function your_function(arg_1::Integer, arg_2::Real)
    if arg_2 < 0
        return 0, 0
    end

    return arg_1 + arg_2, arg_1 - arg_2
end

end # End of module YourSubmodule (This comment is nice to distinguish from an end)


text:
using LinearAlgebra

function polar_force(force,angle)
    return force.*[sin(deg2rad(angle)), cos(deg2rad(angle))]
end

function moments(forces, locations)
    sum = zeros(3)
    for force_loc in zip(forces, locations)
        force, location = force_loc
        force = vcat(force,[0])
        location = vcat(location,[0])
        res = cross(force, location)
        print(res)
        sum+=res
    end
    sum
end

#forces = [polar_force(718.4, 180 - 30), polar_force(879.54, 45), polar_force(100, -90)]
forces = [[0, -2000], [0, -1200], [0, 15600], [0, -12400]]
location = [[0, 0], [6, 0], [10, 0], [12, 0]]
print(moments(forces, location))


text:
export ft_getdeviceinfo

include("FT_DEVICE.jl")

function ft_getdeviceinfo(ft_handle::Culong)
  ft_device = Ref{Cuint}() # see ftdevicetypedict
  id = Ref{Cuint}()
  serialnumber = Array(UInt8,17)
  description = Array(UInt8,65)
  serialnumber[17] = 0x00
  description[65] = 0x00
  ft_status = ccall((:FT_GetDeviceInfo, d2xx),
                     Cuint,
                     (Culong, Ref{Cuint}, Ref{Cuint}, Ptr{UInt8}, Ptr{UInt8}, Ptr{Void}),
                     ft_handle, ft_device, id, serialnumber, description, C_NULL)
  sn = convert(ASCIIString,serialnumber[1:findfirst(serialnumber,0)-1])
  d = convert(ASCIIString,description[1:findfirst(description,0)-1])
  checkstatus(ft_status)
  return (ft_device[], id[], sn, d)
end


text:
struct ChainDict{K, V}
    cur::Dict{K, V}
    init::Ref{ChainDict{K, V}}
end

ChainDict(cur::Dict{K, V}, init::ChainDict{K, V}) where {K, V} = ChainDict(cur, Ref(init))
ChainDict(cur::Dict{K, V}) where {K, V} = ChainDict(cur, Ref{ChainDict{K, V}}())
ChainDict{K, V}() where {K, V} = ChainDict(Dict{K, V}())

function Base.get!(f::Function, c::ChainDict{K, V}, k::K)::V where {K, V}
    get!(c.cur, k) do
        f()
    end
end

function Base.get(c::ChainDict{K, V}, k::K)::V where {K, V}
    get(c.cur, k) do
        isassigned(c.init) ? get(c.init[], k) : throw(KeyError(k))
    end
end

function Base.get(f::Function, c::ChainDict{K, V}, k::K) where {K, V}
    get(c.cur, k) do
        isassigned(c.init) ? get(f, c.init[], k) : f()
    end
end

function for_chaindict(f::Function, d::ChainDict{K, V}) where {K, V}
    keys = Set{K}()
    while true
        for (k, v) in d.cur
            if k in keys
                continue
            else
                push!(keys, k)
            end
            f(k, v)
        end
        if isassigned(d.init)
            d = d.init[]
        else
            return
        end
    end
end

function for_chaindict_dup(f::Function, d::ChainDict{K, V}) where {K, V}
    while true
        for (k, v) in d.cur
            f(k, v)
        end
        if isassigned(d.init)
            d = d.init[]
        else
            return
        end
    end
end

Base.getindex(c::ChainDict, k) = Base.get(c, k)

function Base.setindex!(c::ChainDict{K, V}, value::V, k::K) where {K, V}
    c.cur[k] = value
end

function child(c::ChainDict{K, V}) where {K, V}
    ChainDict(Dict{K, V}(), c)
end

function parent(c::ChainDict{K, V}) where {K, V}
    c.init[]
end


text:
export multi_fit

function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: DefaultRecipe
    Nspec = length(source.domain)

    elapsed = time()
    mzer = GFit.cmpfit()
    mzer.Δfitstat_theshold = 1.e-5

    # Initialize components and guess initial values
    println(logio(source), "\nFit continuum components...")

    multi = MultiModel()
    for id in 1:Nspec
        λ = source.domain[id][:]
        model = Model(source.domain[id],
                      :qso_cont => QSFit.qso_cont_component(TRecipe),
                      :Continuum => SumReducer([:qso_cont]))

        # TODO if source.options[:instr_broadening]
        # TODO     GFit.set_instr_response!(model, (l, f) -> instrumental_broadening(l, f, source.spectra[id].resolution))
        # TODO end

        c = model[:qso_cont]
        c.x0.val = median(λ)
        c.norm.val = Spline1D(λ, source.data[id].val, k=1, bc="error")(c.x0.val)
        push!(multi, model)
    end

    for id in 1:Nspec
        model = multi[id]
        λ = source.domain[id][:]

        # Host galaxy template
        if source.options[:use_host_template]   &&
            (minimum(λ) .< 5500 .< maximum(λ))
            model[:galaxy] = QSFit.hostgalaxy(source.options[:host_template])
            model[:Continuum] = SumReducer([:qso_cont, :galaxy])

            # Split total flux between continuum and host galaxy
            vv = Spline1D(λ, source.data[id].val, k=1, bc="error")(5500.)
            model[:galaxy].norm.val  = 1/2 * vv
            model[:qso_cont].x0.val *= 1/2 * vv / Spline1D(λ, model(:qso_cont), k=1, bc="error")(5500.)

            if id != ref_id
                @patch! multi[id][:galaxy].norm = multi[ref_id][:galaxy].norm
            end
        end

        # Balmer continuum and pseudo-continuum
        if source.options[:use_balmer]
            tmp = [:qso_cont, :balmer]
            (:galaxy in keys(model))  &&  push!(tmp, :galaxy)
            model[:balmer] = QSFit.balmercont(0.1, 0.5)
            model[:Continuum] = SumReducer(tmp)
            c = model[:balmer]
            c.norm.val  = 0.1
            c.norm.fixed = false
            c.norm.high = 0.5
            c.ratio.val = 0.5
            c.ratio.fixed = false
            c.ratio.low  = 0.1
            c.ratio.high = 1
            @patch! model[:balmer].norm *= model[:qso_cont].norm
        end
    end
    bestfit = fit!(multi, source.data, minimizer=mzer);  show(logio(source), bestfit)

    # QSO continuum renormalization
    for id in 1:Nspec
        model = multi[id]
        freeze(model, :qso_cont)
        c = model[:qso_cont]
        initialnorm = c.norm.val
        if c.norm.val > 0
            println(logio(source), "$id: Cont. norm. (before): ", c.norm.val)
            while true
                residuals = (model() - source.data[id].val) ./ source.data[id].unc
                (count(residuals .< 0) / length(residuals) > 0.9)  &&  break
                (c.norm.val < initialnorm / 5)  &&  break # give up
                c.norm.val *= 0.99
                evaluate!(model)
            end
            println(logio(source), "$id : Cont. norm. (after) : ", c.norm.val)
        else
            println(logio(source), "$id: Skipping cont. renormalization")
        end
        freeze(model, :qso_cont)
        (:galaxy in keys(model))  &&  freeze(model, :galaxy)
        (:balmer in keys(model))  &&  freeze(model, :balmer)
    end
    evaluate!(multi)

    # Fit iron templates
    println(logio(source), "\nFit iron templates...")
    for id in 1:Nspec
        model = multi[id]
        λ = source.domain[id][:]

        iron_components = Vector{Symbol}()
        if source.options[:use_ironuv]
            fwhm = 3000.
            source.options[:instr_broadening]  ||  (fwhm = sqrt(fwhm^2 + (source.spectra[id].resolution * 2.355)^2))
            comp = QSFit.ironuv(fwhm)
            (_1, _2, coverage) = spectral_coverage(λ, source.spectra[id].resolution, comp)
            threshold = get(source.options[:min_spectral_coverage], :ironuv, source.options[:min_spectral_coverage][:default])
            if coverage >= threshold
                model[:ironuv] = comp
                model[:ironuv].norm.val = 0.5
                push!(iron_components, :ironuv)
            else
                println(logio(source), "Ignoring ironuv component on prediction $id (threshold: $threshold)")
            end
        end

        if source.options[:use_ironopt]
            fwhm = 3000.
            source.options[:instr_broadening]  ||  (fwhm = sqrt(fwhm^2 + (source.spectra[id].resolution * 2.355)^2))
            comp = QSFit.ironopt_broad(fwhm)
            (_1, _2, coverage) = spectral_coverage(λ, source.spectra[id].resolution, comp)
            threshold = get(source.options[:min_spectral_coverage], :ironopt, source.options[:min_spectral_coverage][:default])
            if coverage >= threshold
                fwhm = 500.
                source.options[:instr_broadening]  ||  (fwhm = sqrt(fwhm^2 + (source.spectra[id].resolution * 2.355)^2))
                model[:ironoptbr] = comp
                model[:ironoptna] = QSFit.ironopt_narrow(fwhm)
                model[:ironoptbr].norm.val = 0.1  # TODO: guess a sensible value
                model[:ironoptna].norm.val = 0.0
                freeze(model, :ironoptna)  # will be freed during last run
                push!(iron_components, :ironoptbr, :ironoptna)
            else
                println(logio(source), "Ignoring ironopt component on prediction $id (threshold: $threshold)")
            end
        end
        if length(iron_components) > 0
            model[:Iron] = SumReducer(iron_components)
            model[:main] = SumReducer([:Continuum, :Iron])
            evaluate!(model)
            bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)
        else
            model[:Iron] = @expr m -> [0.]
            model[:main] = SumReducer([:Continuum, :Iron])
        end
        (:ironuv    in keys(model))  &&  freeze(model, :ironuv)
        (:ironoptbr in keys(model))  &&  freeze(model, :ironoptbr)
        (:ironoptna in keys(model))  &&  freeze(model, :ironoptna)
    end
    evaluate!(multi)

    # Add emission lines
    line_names = [collect(keys(source.line_names[id])) for id in 1:Nspec]
    line_groups = [unique(collect(values(source.line_names[id]))) for id in 1:Nspec]
    println(logio(source), "\nFit known emission lines...")
    for id in 1:Nspec
        model = multi[id]
        λ = source.domain[id][:]
        resid = source.data[id].val - model()  # will be used to guess line normalization
        for (cname, comp) in source.line_comps[id]
            model[cname] = comp
        end
        for (group, lnames) in QSFit.invert_dictionary(source.line_names[id])
            model[group] = SumReducer(lnames)
        end
        model[:main] = SumReducer([:Continuum, :Iron, line_groups[id]...])

        if haskey(model, :MgII_2798)
            model[:MgII_2798].voff.low  = -1000
            model[:MgII_2798].voff.high =  1000
        end
        if haskey(model, :OIII_5007_bw)
            model[:OIII_5007_bw].fwhm.val  = 500
            model[:OIII_5007_bw].fwhm.low  = 1e2
            model[:OIII_5007_bw].fwhm.high = 1e3
            model[:OIII_5007_bw].voff.low  = 0
            model[:OIII_5007_bw].voff.high = 2e3
        end
        for cname in line_names[id]
            model[cname].norm_integrated = source.options[:norm_integrated]
        end

        # Guess values
        evaluate!(model)
        for cname in line_names[id]
            c = model[cname]
            resid_at_line = Spline1D(λ, resid, k=1, bc="nearest")(c.center.val)
            c.norm.val *= abs(resid_at_line) / maximum(model(cname))

            # If instrumental broadening is not used and the line profile
            # is a Gaussian one take spectral resolution into account.
            # This is significantly faster than convolving with an
            # instrument response but has some limitations:
            # - works only with Gaussian profiles;
            # - all components must be additive (i.e. no absorptions)
            # - further narrow components (besides known emission lines)
            #   will not be corrected for instrumental resolution
            if !source.options[:instr_broadening]
                if isa(c, QSFit.SpecLineGauss)
                    c.spec_res_kms = source.spectra[id].resolution
                else
                    println(logio(source), "Line $cname is not a Gaussian profile: Can't take spectral resolution into account")
                end
            end
        end

        # Patch parameters
        @patch! begin
         # model[:OIII_4959].norm = model[:OIII_5007].norm / 3
            model[:OIII_4959].voff = model[:OIII_5007].voff
        end
        @patch! begin
            model[:OIII_5007_bw].voff += model[:OIII_5007].voff
            model[:OIII_5007_bw].fwhm += model[:OIII_5007].fwhm
        end
        @patch! begin
            # model[:OI_6300].norm = model[:OI_6364].norm / 3
            model[:OI_6300].voff = model[:OI_6364].voff
        end
        @patch! begin
            # model[:NII_6549].norm = model[:NII_6583].norm / 3
            model[:NII_6549].voff = model[:NII_6583].voff
        end
        @patch! begin
            # model[:SII_6716].norm = model[:SII_6731].norm / 1.5
            model[:SII_6716].voff = model[:SII_6731].voff
        end

        @patch! model[:na_Hb].voff = model[:na_Ha].voff

        # The following are required to avoid degeneracy with iron
        # template
        @patch! begin
            model[:Hg].voff = model[:br_Hb].voff
            model[:Hg].fwhm = model[:br_Hb].fwhm
        end
        @patch! begin
            model[:br_Hg].voff = model[:br_Hb].voff
            model[:br_Hg].fwhm = model[:br_Hb].fwhm
        end
        @patch! begin
            model[:na_Hg].voff = model[:na_Hb].voff
            model[:na_Hg].fwhm = model[:na_Hb].fwhm
        end

        # Ensure luminosity at peak of the broad base component is
        # smaller than the associated broad component:
        if  haskey(model, :br_Hb)  &&
            haskey(model, :bb_Hb)
            model[:bb_Hb].norm.high = 1
            model[:bb_Hb].norm.val  = 0.5
            @patch! model[:bb_Hb].norm *= model[:br_Hb].norm / model[:br_Hb].fwhm * model[:bb_Hb].fwhm
        end
        if  haskey(model, :br_Ha)  &&
            haskey(model, :bb_Ha)
            model[:bb_Ha].norm.high = 1
            model[:bb_Ha].norm.val  = 0.5
            @patch! model[:bb_Ha].norm *= model[:br_Ha].norm / model[:br_Ha].fwhm * model[:bb_Ha].fwhm
        end

        bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)

        for lname in line_names[id]
            freeze(model, lname)
        end
    end

    # Add unknown lines
    println(logio(source), "\nFit unknown emission lines...")
    if source.options[:n_unk] > 0
        for id in 1:Nspec
            model = multi[id]
            tmp = OrderedDict{Symbol, GFit.AbstractComponent}()
            for j in 1:source.options[:n_unk]
                tmp[Symbol(:unk, j)] = line_component(TRecipe, QSFit.UnkLine(5e3))
                tmp[Symbol(:unk, j)].norm_integrated = source.options[:norm_integrated]
            end
            for (cname, comp) in tmp
                model[cname] = comp
            end
            model[:UnkLines] = SumReducer(collect(keys(tmp)))
            model[:main] = SumReducer([:Continuum, :Iron, line_groups[id]..., :UnkLines])
            evaluate!(model)
            for j in 1:source.options[:n_unk]
                freeze(model, Symbol(:unk, j))
            end
        end
    else
        # Here we need a :UnkLines reducer, even when n_unk is 0
        for id in 1:Nspec
            model = multi[id]
            model[:UnkLines] = @expr m -> [0.]
            model[:main] = SumReducer([:Continuum, :Iron, line_groups[id]...])
        end
    end
    evaluate!(multi)

    # Set "unknown" line center wavelength where there is a maximum in
    # the fit residuals, and re-run a fit.
    for id in 1:Nspec
        model = multi[id]
        λ = source.domain[id][:]
        λunk = Vector{Float64}()
        while true
            (length(λunk) >= source.options[:n_unk])  &&  break
            evaluate!(model)
            Δ = (source.data[id].val - model()) ./ source.data[id].unc

            # Avoid considering again the same region (within 1A) TODO: within resolution
            for l in λunk
                Δ[findall(abs.(l .- λ) .< 1)] .= 0.
            end

            # Avoidance regions
            for rr in source.options[:unk_avoid]
                Δ[findall(rr[1] .< λ .< rr[2])] .= 0.
            end

            # Do not add lines close to from the edges since these may
            # affect qso_cont fitting
            Δ[findall((λ .< minimum(λ)*1.02)  .|
                      (λ .> maximum(λ)*0.98))] .= 0.
            iadd = argmax(Δ)
            (Δ[iadd] <= 0)  &&  break  # No residual is greater than 0, skip further residuals....
            push!(λunk, λ[iadd])

            cname = Symbol(:unk, length(λunk))
            model[cname].norm.val = 1.
            model[cname].center.val  = λ[iadd]
            model[cname].center.low  = λ[iadd] - λ[iadd]/10. # allow to shift 10%
            model[cname].center.high = λ[iadd] + λ[iadd]/10.

            thaw(model, cname)
            bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)
            freeze(model, cname)
        end
    end
    evaluate!(multi)

    # ----------------------------------------------------------------
    # Last run with all parameters free
    println(logio(source), "\nLast run with all parameters free...")
    for id in 1:Nspec
        model = multi[id]
        thaw(model, :qso_cont)
        (:galaxy in keys(model))        &&  thaw(model, :galaxy)
        (:balmer in keys(model))        &&  thaw(model, :balmer)
        (:ironuv    in keys(model))     &&  thaw(model, :ironuv)
        (:ironoptbr in keys(model))     &&  thaw(model, :ironoptbr)
        (:ironoptna in keys(model))     &&  thaw(model, :ironoptna)

        for lname in line_names[id]
            thaw(model, lname)
        end
        for j in 1:source.options[:n_unk]
            cname = Symbol(:unk, j)
            if model[cname].norm.val > 0
                thaw(model, cname)
            else
                freeze(model, cname)
            end
        end
    end
    bestfit = fit!(multi, source.data, minimizer=mzer)

    # Disable "unknown" lines whose normalization uncertainty is larger
    # than 3 times the normalization
    needs_fitting = false
    for id in 1:Nspec
        model = multi[id]
        for ii in 1:source.options[:n_unk]
            cname = Symbol(:unk, ii)
            isfixed(model, cname)  &&  continue
            if bestfit[id][cname].norm.val == 0.
                freeze(model, cname)
                needs_fitting = true
                println(logio(source), "Disabling $cname (norm. = 0)")
            elseif bestfit[id][cname].norm.unc / bestfit[id][cname].norm.val > 3
                model[cname].norm.val = 0.
                freeze(model, cname)
                needs_fitting = true
                println(logio(source), "Disabling $cname (unc. / norm. > 3)")
            end
        end
    end
    if needs_fitting
        println(logio(source), "\nRe-run fit...")
        bestfit = fit!(multi, source.data, minimizer=mzer)
    end

    println(logio(source))
    show(logio(source), bestfit)

    out = QSFit.QSFitMultiResults(source, multi, bestfit)
    elapsed = time() - elapsed
    println(logio(source), "\nElapsed time: $elapsed s")
    QSFit.close_logio(source)
    return out
end


text:
# type Dyadic
#    num::Int64
#    log2den::Int64
#    function Dyadic(num::Integer, den::Integer)
#        tmp = log2(abs(den))
#        if isinteger(tmp)
#            return new( sign(den)*num, Int(tmp) )
#        else
#            error( DomainError() )
#        end
#    end
#end
#convert(Float64, x::Dyadic) = x.num / 2^(x.log2den)


text:
# Max function

export Maximum

"""
**Maximum coefficient**

    Maximum(λ=1)

For a nonnegative parameter `λ ⩾ 0`, returns the function
```math
f(x) = \\lambda \\cdot \\max \\{x_i : i = 1,\\ldots, n \\}.
```
"""
Maximum(lambda::R=1) where {R <: Real} = SumLargest(one(Int32), lambda)


text:
using DataFrames
using DynamicsPlots
using Plots, RecipesBase, Test

df = DataFrame(x = 1:10, y = rand(10))
timeevolution(df, legend = :left)

include("../src/pointcloud.jl")
P = randn(2,100)
pointcloud(P)
P = [randn(2) for _ in 1:100]
pointcloud!(P, color = :red)

include("../src/geometry.jl")
circle([-0,0],2,color = :red, fill = (0,0.1,:red))


text:
#=
Copyright 2017 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

# Documentation to be completed.
# See references[^1][^2][^3][^4]
function computeFlux(v :: Vessel, A :: Array{Float64, 1},
                     Q :: Array{Float64, 1}, Flux :: Array{Float64, 2})

  for i in 1:v.M+2
    Flux[1,i] = Q[i]
    Flux[2,i] = Q[i]*Q[i]/A[i] + v.gamma_ghost[i] * A[i]*sqrt(A[i])
  end

  return Flux
end

function maxmod(a :: Float64, b :: Float64)
  if a > b
    return a

  else
    return b

  end
end

function minmod(a :: Float64, b :: Float64)
  if a <= 0. || b <= 0.
    return 0.

  elseif a < b
    return a

  else
    return b

  end
end

function superBee(v :: Vessel, dU :: Array{Float64, 2}, slopes :: Array{Float64, 1})

  for i in 1:v.M+2
    s1 = minmod(dU[1,i], 2*dU[2,i])
    s2 = minmod(2*dU[1,i], dU[2,i])
    slopes[i] = maxmod(s1, s2)
  end

  return slopes
end

function computeLimiter(v :: Vessel, U :: Array{Float64, 1},
                        invDx :: Float64, dU :: Array{Float64, 2},
                        slopes :: Array{Float64, 1})

  for i = 2:v.M+2
    dU[1, i]   = (U[i] - U[i-1])*invDx
    dU[2, i-1] = (U[i] - U[i-1])*invDx
  end

  return superBee(v, dU, slopes)

end

function computeLimiter(v :: Vessel, U :: Array{Float64, 2}, idx :: Int64,
                        invDx :: Float64, dU :: Array{Float64, 2},
                        slopes :: Array{Float64, 1})
  U = U[idx,:]
  for i = 2:v.M+2
    dU[1, i]   = (U[i] - U[i-1])*invDx
    dU[2, i-1] = (U[i] - U[i-1])*invDx
  end

  return superBee(v, dU, slopes)

end

function MUSCL(v :: Vessel, dt :: Float64, b :: Blood)

  v.vA[1] = v.U00A
  v.vA[end] = v.UM1A

  v.vQ[1] = v.U00Q
  v.vQ[end] = v.UM1Q

  for i = 2:v.M+1
    v.vA[i] = v.A[i-1]
    v.vQ[i] = v.Q[i-1]
  end

  v.slopesA = computeLimiter(v, v.vA, v.invDx, v.dU, v.slopesA)
  v.slopesQ = computeLimiter(v, v.vQ, v.invDx, v.dU, v.slopesQ)

  for i = 1:v.M+2
    v.Al[i] = v.vA[i] + v.slopesA[i]*v.halfDx
    v.Ar[i] = v.vA[i] - v.slopesA[i]*v.halfDx

    v.Ql[i] = v.vQ[i] + v.slopesQ[i]*v.halfDx
    v.Qr[i] = v.vQ[i] - v.slopesQ[i]*v.halfDx
  end

  v.Fl = computeFlux(v, v.Al, v.Ql, v.Fl)
  v.Fr = computeFlux(v, v.Ar, v.Qr, v.Fr)

  dxDt = v.dx/dt
  invDxDt = 1./dxDt

  for i = 1:v.M+1
    v.flux[1, i] = 0.5 * ( v.Fr[1, i+1] + v.Fl[1, i] - dxDt *
                        (v.Ar[i+1]    - v.Al[i]) )
    v.flux[2, i] = 0.5 * ( v.Fr[2, i+1] + v.Fl[2, i] - dxDt *
                        (v.Qr[i+1]    - v.Ql[i]) )
  end

  for i = 2:v.M+1
    v.uStar[1, i] = v.vA[i] + invDxDt * (v.flux[1, i-1] - v.flux[1, i])
    v.uStar[2, i] = v.vQ[i] + invDxDt * (v.flux[2, i-1] - v.flux[2, i])
  end

  v.uStar[1,1] = v.uStar[1,2]
  v.uStar[2,1] = v.uStar[2,2]
  v.uStar[1,end] = v.uStar[1,end-1]
  v.uStar[2,end] = v.uStar[2,end-1]

  # v.slopesA = computeLimiter(v, vec(v.uStar[1,:]), v.invDx, v.dU, v.slopesA)
  # v.slopesQ = computeLimiter(v, vec(v.uStar[2,:]), v.invDx, v.dU, v.slopesQ)
  v.slopesA = computeLimiter(v, v.uStar, 1, v.invDx, v.dU, v.slopesA)
  v.slopesQ = computeLimiter(v, v.uStar, 2, v.invDx, v.dU, v.slopesQ)

  for i = 1:v.M+2
    v.Al[i] = v.uStar[1,i] + v.slopesA[i]*v.halfDx
    v.Ar[i] = v.uStar[1,i] - v.slopesA[i]*v.halfDx

    v.Ql[i] = v.uStar[2,i] + v.slopesQ[i]*v.halfDx
    v.Qr[i] = v.uStar[2,i] - v.slopesQ[i]*v.halfDx
  end

  v.Fl = computeFlux(v, v.Al, v.Ql, v.Fl)
  v.Fr = computeFlux(v, v.Ar, v.Qr, v.Fr)

  for i = 1:v.M+1
    v.flux[1, i] = 0.5 * ( v.Fr[1, i+1] + v.Fl[1, i] - dxDt *
                        (v.Ar[i+1]    - v.Al[i]) )
    v.flux[2, i] = 0.5 * ( v.Fr[2, i+1] + v.Fl[2, i] - dxDt *
                        (v.Qr[i+1]    - v.Ql[i]) )
  end

  for i = 2:v.M+1
    v.A[i-1] = 0.5*(v.A[i-1] + v.uStar[1, i] + invDxDt *
                    (v.flux[1, i-1] - v.flux[1, i]) )
    v.Q[i-1] = 0.5*(v.Q[i-1] + v.uStar[2, i] + invDxDt *
                    (v.flux[2, i-1] - v.flux[2, i]))
  end

  #source term
  rho_inv = 1./b.rho
  viscT = 2*(b.gamma_profile + 2)*pi*b.mu
  for i = 1:v.M
    s_A0_inv = sqrt(1./v.A0[i])
    s_A = sqrt(v.A[i])
    s_A_inv = 1./s_A
    # v.Q[i] -= 2*(b.gamma_profile+2)*pi*b.mu*v.Q[i]/(v.A[i])*rho_inv*dt   #viscosity
    # # v.Q[i] += dt*0.5*v.beta[i]*v.A[i]^1.5/(v.A0[i]*b.rho)*v.dA0dx[i]   #dP/dA0
    # v.Q[i] += dt*0.5*v.beta[i]*v.A[i]^1.5*A0_inv*rho_inv*v.dA0dx[i]   #dP/dA0
    # v.Q[i] -= dt*(v.A[i]*rho_inv)*(sqrt(v.A[i]*A0_inv)-1.)*v.dTaudx[i]  #dP/dh0
    # v.Q[i] += dt*v.A[i]*s_A*rho_inv *
    #           (0.5*v.beta[i]*s_A0_inv*s_A0_inv*v.dA0dx[i] -
    #           (s_A0_inv - s_A_inv)*v.dTaudx[i] -
    #            2*(b.gamma_profile + 2)*pi*b.mu*v.Q[i]*
    #            s_A_inv*s_A_inv*s_A_inv*s_A_inv*s_A_inv )

    v.Q[i] += dt*rho_inv*( -viscT*v.Q[i]*s_A_inv*s_A_inv +
                        v.A[i]*(v.beta[i]*0.5*v.dA0dx[i] -
                                (s_A*s_A0_inv-1.)*v.dTaudx[i]))

    v.P[i] = pressure(v.A[i], v.A0[i], v.beta[i], v.Pext)
    v.u[i] = v.Q[i]/v.A[i]
    v.c[i] = waveSpeed(v.A[i], v.gamma[i])
  end

end
#
# ### References
#
# [^1]: LeVeque RJ. Finite volume methods for hyperbolic problems. Cambridge university press; 2002.
# [^2]: LeVeque RJ. Numerical methods for conservation laws. Birkhauser; 1992.
# [^3]: Toro, Eleuterio F. Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Springer Science and Business Media, 2009.
# [^4]: Van Leer B. Towards the ultimate conservative difference scheme. Journal of Computational Physics. 1997;135(2):229-48.


text:
function TransformedDomain{T,S}(
    domain::TransformedDomain{T,S},
    parent::ArbField,
) where {T,S}
    # TODO: Might need to recompute values here
    return TransformedDomain(
        S(domain.original, parent),
        domain.rotation,
        domain.scaling,
        domain.translation,
    )
end

function TransformedDomain(domain, rotation::T, scaling, translation) where {T}
    if T == fmpq || T <: Rational
        rotation = fmpq(rotation)
    else
        rotation = domain.parent(rotation)
    end
    scaling = domain.parent(scaling)
    translation = SVector{2,arb}(domain.parent.(translation))
    return TransformedDomain(domain, rotation, scaling, translation)
end

function Base.show(io::IO, domain::TransformedDomain{T}) where {T}
    println(io, "Transformed Domain")
    if T == arb
        println(io, "Rotation: $(domain.rotation)")
    else
        println(
            io,
            "Rotation: $(numerator(domain.rotation))π/$(denominator(domain.rotation))",
        )
    end
    println(io, "Scaling: $(domain.scaling)")
    println(io, "Translation: $(domain.translation)")
    print(io, "Domain: $(domain.original)")
end

function Base.getproperty(domain::TransformedDomain{T}, name::Symbol) where {T}
    if name == :parent
        return domain.original.parent
    elseif name == :map
        if T == arb
            s, c = sincos(domain.rotation)
        elseif T == fmpq
            s, c = sincospi(domain.rotation, domain.parent)
        end
        M = SMatrix{2,2}(c, s, -s, c)
        return xy -> domain.scaling .* M * xy + domain.translation
    elseif name == :invmap
        if T == arb
            s, c = sincos(-domain.rotation)
        elseif T == fmpq
            s, c = sincospi(-domain.rotation, domain.parent)
        end
        M = SMatrix{2,2}(c, s, -s, c)
        return xy -> M * (xy - domain.translation) ./ domain.scaling
    else
        return getfield(domain, name)
    end
end

boundaries(domain::TransformedDomain) = boundaries(domain.original)

angle(domain::TransformedDomain, i::Integer) = angle(domain.original, i)
angledivπ(domain::TransformedDomain, i::Integer) = angledivπ(domain.original, i)
angles(domain::TransformedDomain) = angles(domain.original)
angledivπ(domain::TransformedDomain) = angledivπ(domain.original)

vertex(domain::TransformedDomain, i::Integer) = domain.map(vertex(domain.original, i))
vertices(domain::TransformedDomain) = domain.map.(vertices(domain.original))

center(domain::TransformedDomain) = domain.map(center(domain.original))

area(domain::TransformedDomain) = domain.scaling * area(domain.original)

Base.in(xy, domain::TransformedDomain) = domain.invmap(xy) ∈ domain.original

boundary_parameterization(t, domain::TransformedDomain, i::Integer) =
    domain.map(boundary_parameterization(t, domain.original, i))

function boundary_points(
    domain::TransformedDomain,
    i::Integer,
    n::Integer;
    distribution = :chebyshev,
)
    pts, idx = boundary_points(domain.original, i, n; distribution)
    return domain.map.(pts), idx
end

interior_points(domain::TransformedDomain, n::Integer; rng = MersenneTwister(42)) =
    domain.map.(interior_points(domain.original, n, rng = rng))


text:
module BlockDiagonals

using Base: @propagate_inbounds
using ChainRulesCore
using FillArrays: Zeros
using FiniteDifferences
using LinearAlgebra

import ChainRulesCore.ProjectTo

export BlockDiagonal, blocks
export blocksize, blocksizes, nblocks

include("blockdiagonal.jl")
include("base_maths.jl")
include("chainrules.jl")
include("linalg.jl")

end  # end module


text:
using JuMP
using CPLEX

m = read_from_file("academic_timetables/academictimetablesmall.mps")

set_optimizer(m, CPLEX.Optimizer)
optimize!(m)

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
using YaoSym, YaoBlocks, YaoArrayRegister, SymEngine
using YaoSym: simplify_expi
using SymEngine
using Test

@testset "imag and conj" begin
    @vars a b c x
    @test imag(Basic(2 + im * a)) == a
    @test imag(Basic((2 + im) * (im * a + 1))) == 2 * a + 1
    @test real(sin(a)) == sin(a)
    @test imag(sin(a)) == 0
    @test imag(x^2) == 0
    @test real(x^2) == x^2
    @test real(exp(im * a)) == cos(a)
    @test imag(exp(im * a)) == sin(a)
    @test abs(exp(im * a)) == 1
    @test abs(sin(a)) == sin(a)
    @test abs(sin(a) + 2) == sin(a) + 2
end

@testset "mat" begin
    @vars θ γ η
    for G in [X, Y, Z, ConstGate.T, H]
        @test Matrix(mat(Basic, G)) ≈ Matrix(mat(G))
        @test Matrix(mat(Basic, control(4, 3, 2 => G))) ≈ Matrix(control(4, 3, 2 => G))
    end

    for GC in [Rx, shift, phase]
        G = GC(θ)
        m = mat(Basic, G)
        m = subs.(m, θ, 0.5)
        @test Matrix(mat(GC(0.5))) ≈ Matrix(m)
        @test subs.(Matrix(mat(Basic, put(4, 2 => G))), θ, 0.5) ≈ Matrix(put(4, 2 => GC(0.5)))
    end

    G = Rz(θ) * Rx(γ) * Rz(θ)
    m = mat(Basic, G)
    @test mat(G) == m
    @test_throws ArgumentError mat(Float64, G)
    m = subs.(m, Ref(θ => Basic(π) / 2), Ref(γ => Basic(π) / 6))
    @test Matrix(mat(Rz(π / 2) * Rx(π / 6) * Rz(π / 2))) ≈ Matrix(m)

    A = randn(ComplexF64, 4, 4)
    mb = matblock(Basic.(A))
    @test Matrix(mb) ≈ A

    @test pswap(4, 2, 1, θ) == pswap(2, 1, θ)(4)
    @test pswap(4, 2, 1, θ) isa PSwap{4,Basic}
end

@testset "sub" begin
    @vars θ ϕ
    @test subs(Float64, Rx(θ), θ => 0.5) == Rx(0.5)
    @test subs(Float64, shift(θ), θ => 0.5) == shift(0.5)
    @test subs(Float64, phase(θ), θ => 0.5) == phase(0.5)
    @test subs(
        Float64,
        chain(control(2, 1, 2 => shift(θ)), put(2, 1 => chain(Rx(θ), kron(Rx(θ))))),
        θ => 0.5,
    ) == chain(control(2, 1, 2 => shift(0.5)), put(2, 1 => chain(Rx(0.5), kron(Rx(0.5)))))
    @test subs(Float64, time_evolve(X, θ), θ => 0.5) == time_evolve(X, 0.5)
    @test subs(Rx(θ), θ => ϕ) == Rx(ϕ)
end

@testset "Sym Engine Conj" begin
    @vars a b
    @test cos(exp(a + im * b))' == cos(exp(a - im * b))
    @test cos(exp(a + im * b))' == cos(exp(a - im * b))
    @test mat(Rx(a)) == [
        cos(a / 2) -im*sin(a / 2)
        -im*sin(a / 2) cos(a / 2)
    ]
    @test mat(Rx(a)') == [
        cos(a / 2) im*sin(a / 2)
        im*sin(a / 2) cos(a / 2)
    ]
end

@testset "expect" begin
    @vars a b
    reg = ket"111" |> control(3, 1, 2 => Rx(a))
    op = put(3, 2 => Z)
    reg2 = ArrayReg(bit"111") |> control(3, 1, 2 => Rx(0.5))
    ex = expect(op, reg)
    ex = subs(ex, a => 0.5)
    @test ComplexF64(ex) ≈ expect(op, reg2)
end

@testset "grad" begin
    @vars a b
    reg = ArrayReg(Basic, bit"111") => control(3, 1, 2 => Rx(a))
    op = put(3, 2 => Z)
    reg2 = ArrayReg(bit"111") => control(3, 1, 2 => Rx(0.5))
    ex = expect'(op, reg)[2]
    ex = subs(ex[], a => 0.5)
    @test ComplexF64(ex) ≈ expect'(op, reg2)[2][]
end

@testset "grad" begin
    @vars a b
    reg = ArrayReg(Basic, bit"000") => put(3, 2 => Rx(a))
    op = put(3, 2 => Z)
    reg2 = ArrayReg(bit"000") => put(3, 2 => Rx(0.9))
    ex = expect(op, reg)
    ex = subs(ex[], a => 0.9)
    @test ComplexF64(ex) ≈ expect(op, reg2)
    ex = expect'(op, reg)[2]
    ex = subs(ex[], a => 0.9)
    @test ComplexF64(ex) ≈ expect'(op, reg2)[2][]
end

function check_dumpload(gate::AbstractBlock{N}) where {N}
    gate2 = eval(YaoBlocks.parse_ex(YaoBlocks.dump_gate(gate), N))
    gate2 == gate || mat(gate2) ≈ mat(gate)
end

@testset "dumpload" begin
    # basic types
    @vars θ
    @test check_dumpload(phase(θ))
    @test check_dumpload(shift(θ))
    @test check_dumpload(time_evolve(X, θ))
    @test check_dumpload(rot(X, θ))
end


text:
# Power flow module. The implementation is a modification of
# MATPOWER's code. We attach the original MATPOWER's license in makeYbus.m:
#
# MATPOWER
# Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC)
# by Ray Zimmerman, PSERC Cornell
#
# Covered by the 3-clause BSD License.
#
# Expression of Jacobians and Hessians from MATPOWER

#
function matpower_residual_jacobian(V, Ybus)
    n = size(V, 1)
    Ibus = Ybus*V
    diagV       = sparse(1:n, 1:n, V, n, n)
    diagIbus    = sparse(1:n, 1:n, Ibus, n, n)
    diagVnorm   = sparse(1:n, 1:n, V./abs.(V), n, n)

    dSbus_dVm = diagV * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm
    dSbus_dVa = 1im * diagV * conj(diagIbus - Ybus * diagV)
    return (dSbus_dVm, dSbus_dVa)
end

function _matpower_lineflow_jacobian(V, branches::Branches)
    nb = size(V, 1)
    f = branches.from_buses
    t = branches.to_buses
    nl = length(f)

    diagV     = sparse(1:nb, 1:nb, V, nb, nb)
    diagVnorm = sparse(1:nb, 1:nb, V./abs.(V), nb, nb)

    # Connection matrices
    Cf = sparse(1:nl, f, ones(nl), nl, nb)
    Ct = sparse(1:nl, t, ones(nl), nl, nb)

    i = [1:nl; 1:nl]
    Yf = sparse(i, [f; t], [branches.Yff; branches.Yft], nl, nb)
    Yt = sparse(i, [f; t], [branches.Ytf; branches.Ytt], nl, nb)

    If = Yf * V
    It = Yt * V

    diagCfV = sparse(1:nl, 1:nl, Cf * V, nl, nl)
    diagCtV = sparse(1:nl, 1:nl, Ct * V, nl, nl)

    Sf = diagCfV * conj(If)
    St = diagCtV * conj(It)

    diagIf = sparse(1:nl, 1:nl, If, nl, nl)
    diagIt = sparse(1:nl, 1:nl, It, nl, nl)

    dSf_dVm = conj(diagIf) * Cf * diagVnorm + diagCfV * conj(Yf * diagVnorm)
    dSf_dVa = 1im * (conj(diagIf) * Cf * diagV - diagCfV * conj(Yf * diagV))

    dSt_dVm = conj(diagIt) * Ct * diagVnorm + diagCtV * conj(Yt * diagVnorm)
    dSt_dVa = 1im * (conj(diagIt) * Ct * diagV - diagCtV * conj(Yt * diagV))

    return (Sf, St, dSf_dVm, dSf_dVa, dSt_dVm, dSt_dVa)
end

function matpower_lineflow_power_jacobian(V, branches::Branches)
    nl = length(branches.from_buses)
    (Sf, St, dSf_dVm, dSf_dVa, dSt_dVm, dSt_dVa) = _matpower_lineflow_jacobian(V, branches)

    dSf = sparse(1:nl, 1:nl, Sf, nl, nl)
    dSt = sparse(1:nl, 1:nl, St, nl, nl)

    dHf_dVm = 2 * (real(dSf) * real(dSf_dVm) + imag(dSf) * imag(dSf_dVm))
    dHf_dVa = 2 * (real(dSf) * real(dSf_dVa) + imag(dSf) * imag(dSf_dVa))

    dHt_dVm = 2 * (real(dSt) * real(dSt_dVm) + imag(dSt) * imag(dSt_dVm))
    dHt_dVa = 2 * (real(dSt) * real(dSt_dVa) + imag(dSt) * imag(dSt_dVa))

    dH_dVm = [dHf_dVm; dHt_dVm]
    dH_dVa = [dHf_dVa; dHt_dVa]
    return dH_dVm, dH_dVa
end

# Hessian vector-product H*λ, H = Jₓₓ (from Matpower)
# Suppose ordering is correct
function _matpower_hessian(V, Ybus, λ)
    n = size(V, 1)
    # build up auxiliary matrices
    Ibus = Ybus * V
    diagV       = sparse(1:n, 1:n, V, n, n)
    diagIbus    = sparse(1:n, 1:n, Ibus, n, n)
    d_invV      = sparse(1:n, 1:n, 1 ./ abs.(V), n, n)
    d_lambda    = sparse(1:n, 1:n, λ, n, n)

    A = sparse(1:n, 1:n, λ .* V, n, n)
    B = Ybus * diagV
    C = A * conj(B)
    D = Ybus' * diagV
    Dl = sparse(1:n, 1:n, D * λ, n, n)
    E = conj(diagV)  * (D * d_lambda - Dl)
    F = C - A * sparse(1:n, 1:n, conj(Ibus), n, n)
    G = d_invV

    # θθ
    G11 = E + F
    # vθ
    G21 = 1im .* G * (E - F)
    # vv
    G22 = G * (C + transpose(C)) * G
    return (G11, transpose(G21), G22)
end

function residual_hessian(V, Ybus, λ, pv, pq, ref)
    # decompose vector
    n = length(V)
    λp = zeros(n) ; λq = zeros(n)
    npv = length(pv)
    npq = length(pq)
    nref = length(ref)
    λp[pv] = λ[1:npv]
    λp[pq] = λ[npv+1:npv+npq]
    λq[pq] = λ[npv+npq+1:end]

    Gp11, Gp12, Gp22 = _matpower_hessian(V, Ybus, λp)
    Pθθ = real.(Gp11)
    Pvθ = real.(Gp12)
    Pvv = real.(Gp22)

    Gq11, Gq12, Gq22 = _matpower_hessian(V, Ybus, λq)
    Qθθ = imag.(Gq11)
    Qvθ = imag.(Gq12)
    Qvv = imag.(Gq22)

    # w.r.t. xx
    H11 = Pθθ[pv, pv] + Qθθ[pv, pv]
    H12 = Pθθ[pv, pq] + Qθθ[pv, pq]
    H13 = Pvθ[pv, pq] + Qvθ[pv, pq]
    H22 = Pθθ[pq, pq] + Qθθ[pq, pq]
    H23 = Pvθ[pq, pq] + Qvθ[pq, pq]
    H33 = Pvv[pq, pq] + Qvv[pq, pq]

    Hxx = [
        H11  H12  H13
        H12' H22  H23
        H13' H23' H33
    ]::SparseMatrixCSC{Float64, Int}

    # w.r.t. uu
    H11 = Pvv[ref, ref] + Qvv[ref, ref]
    H12 = Pvv[ref,  pv] + Qvv[ref,  pv]
    H22 = Pvv[pv,   pv] + Qvv[pv,   pv]

    Huu = [
         H11  H12 spzeros(nref, npv)
         H12' H22 spzeros(npv, npv)
         spzeros(npv, nref + 2 * npv)
    ]::SparseMatrixCSC{Float64, Int}

    # w.r.t. xu
    Pvθ = real.(transpose(Gp12))
    Qvθ = imag.(transpose(Gq12))
    H11 = Pvθ[ref, pv] + Qvθ[ref, pv]
    H12 = Pvθ[ref, pq] + Qvθ[ref, pq]
    H13 = Pvv[ref, pq] + Qvv[ref, pq]
    H21 = Pvθ[pv,  pv] + Qvθ[pv,  pv]
    H22 = Pvθ[pv,  pq] + Qvθ[pv,  pq]
    H23 = Pvv[pv,  pq] + Qvv[pv,  pq]

    Hxu = [
        H11  H12  H13
        H21  H22  H23
        spzeros(npv, npv + 2 * npq)
    ]::SparseMatrixCSC{Float64, Int}

    return (Hxx, Hxu, Huu)
end

# ∂²pg_ref / ∂²x
function active_power_hessian(V, Ybus, pv, pq, ref)
    # decompose vector
    n = length(V)
    npv = length(pv)
    npq = length(pq)
    nref = length(ref)

    λp = zeros(n)
    # Pick only components wrt ref nodes
    λp[ref] .= 1.0

    G11, G12, G22 = _matpower_hessian(V, Ybus, λp)
    Pθθ = real.(G11)
    Pvθ = real.(G12)
    Pvv = real.(G22)

    H11 = Pθθ[pv, pv]
    H12 = Pθθ[pv, pq]
    H13 = Pvθ[pv, pq]
    H22 = Pθθ[pq, pq]
    H23 = Pvθ[pq, pq]
    H33 = Pvv[pq, pq]

    # w.r.t. xx
    Hxx = [
        H11  H12  H13
        H12' H22  H23
        H13' H23' H33
    ]::SparseMatrixCSC{Float64, Int}

    # w.r.t. uu
    H11 = Pvv[ref, ref]
    H12 = Pvv[ref, pv]
    H22 = Pvv[pv, pv]

    Huu = [
         H11  H12 spzeros(nref, npv)
         H12' H22 spzeros(npv, npv)
         spzeros(npv, nref + 2 * npv)
    ]::SparseMatrixCSC{Float64, Int}

    # w.r.t. xu
    Pvθ = real.(transpose(G12))
    Qvθ = imag.(transpose(G12))
    H11 = Pvθ[ref, pv]
    H12 = Pvθ[ref, pq]
    H13 = Pvv[ref, pq]
    H21 = Pvθ[pv, pv]
    H22 = Pvθ[pv, pq]
    H23 = Pvv[pv, pq]

    Hxu = [
        H11  H12  H13
        H21  H22  H23
        spzeros(npv, npv + 2 * npq)
    ]::SparseMatrixCSC{Float64, Int}

    return (Hxx, Hxu, Huu)
end

# ∂²qg / ∂²x * λ
function reactive_power_hessian(V, Ybus, λ, pv, pq, ref)
    n = length(V)
    npv = length(pv)
    npq = length(pq)
    nref = length(ref)

    G11, G12, G22 = _matpower_hessian(V, Ybus, λ)
    Qθθ = imag.(G11)
    Qvθ = imag.(G12)
    Qvv = imag.(G22)

    H11 = Qθθ[pv, pv]
    H12 = Qθθ[pv, pq]
    H13 = Qvθ[pv, pq]
    H22 = Qθθ[pq, pq]
    H23 = Qvθ[pq, pq]
    H33 = Qvv[pq, pq]

    Hxx = [
        H11  H12  H13
        H12' H22  H23
        H13' H23' H33
    ]::SparseMatrixCSC{Float64, Int}

    H11 = Qvv[ref, ref]
    H12 = Qvv[ref, pv]
    H22 = Qvv[pv, pv]

    Huu = [
         H11  H12 spzeros(nref, npv)
         H12' H22 spzeros(npv, npv)
         spzeros(npv, nref + 2* npv)
    ]::SparseMatrixCSC{Float64, Int}

    Pvθ = real.(transpose(G12))
    Qvθ = imag.(transpose(G12))
    H11 = Qvθ[ref, pv]
    H12 = Qvθ[ref, pq]
    H13 = Qvv[ref, pq]
    H21 = Qvθ[pv, pv]
    H22 = Qvθ[pv, pq]
    H23 = Qvv[pv, pq]

    Hxu = [
        H11  H12  H13
        H21  H22  H23
        spzeros(npv, 2*npq+npv)
    ]::SparseMatrixCSC{Float64, Int}
    return (Hxx, Hxu, Huu)
end


text:
using GLVisualize, GeometryTypes, Colors
using GLAbstraction
if !isdefined(:runtests)
	window = glscreen()
	timesignal = bounce(linspace(0f0,1f0,360))
end

description = """
Animated 2D particles in 3D space with glow.
"""

function sin_torus(radius, thickness, N, t)
    pirange = linspace(0, 2pi, N^2)
    points = zeros(Point3f0, N^2)
    p,q = 6,9
    for (i,x)=enumerate(pirange)
        r = cos(q*(x+t))+thickness
        points[i] = Point3f0(r*cos(p*x),r*sin(p*x), -sin(q*x))/radius
    end
    points
end

t = const_lift(*, 1000f0, timesignal)
ps1 = const_lift(sin_torus, 2, 6, 10, t)
ps2 = const_lift(sin_torus, 5, 4, 10, t)
ps3 = const_lift(sin_torus, 10, 2, 10, t)
ps = map(vcat, ps1, ps2, ps3)
l1, l2, l3 = map(x->length(value(x)), (ps1, ps2, ps3))
colors = [
    fill(RGBA(0.9f0, 0f0, 0.7f0, 1f0), l1);
    fill(RGBA(0.7f0, 0.2f0, 0.8f0, 1f0), l2);
    fill(RGBA(1.0f0, 0.4f0, 0.95f0, 1f0), l3);
]
scales = [
    fill(Vec2f0(0.1), l3);
    fill(Vec2f0(0.05), l2);
    fill(Vec2f0(0.04), l1);
]

_view(visualize(
    (Circle, ps),
    glow_color = RGBA(0.8f0, 0.6f0, 0.95f0, 0.8f0),
    glow_width = const_lift(/, timesignal, 20f0),
    scale = scales, color = colors
), camera=:perspective)

if !isdefined(:runtests)
    renderloop(window)
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    isdefined(COSMO, Symbol("#19#20")) && precompile(Tuple{getfield(COSMO, Symbol("#19#20")),Tuple{SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int}},true},COSMO.ExponentialCone{Float64}}})
    isdefined(COSMO, Symbol("#19#20")) && precompile(Tuple{getfield(COSMO, Symbol("#19#20")),Tuple{SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int}},true},COSMO.Nonnegatives{Float64}}})
    isdefined(COSMO, Symbol("#19#20")) && precompile(Tuple{getfield(COSMO, Symbol("#19#20")),Tuple{SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int}},true},COSMO.PsdConeTriangle{Float64}}})
    isdefined(COSMO, Symbol("#19#20")) && precompile(Tuple{getfield(COSMO, Symbol("#19#20")),Tuple{SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int}},true},COSMO.SecondOrderCone{Float64}}})
    isdefined(COSMO, Symbol("#19#20")) && precompile(Tuple{getfield(COSMO, Symbol("#19#20")),Tuple{SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int}},true},COSMO.ZeroSet{Float64}}})
    isdefined(COSMO, Symbol("#25#26")) && precompile(Tuple{getfield(COSMO, Symbol("#25#26")),Tuple{SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int}},true},COSMO.SecondOrderCone{Float64}}})
    precompile(Tuple{Type{CholmodKKTSolver},SparseArrays.SparseMatrixCSC{Float64,Int},SparseArrays.SparseMatrixCSC{Float64,Int},Float64,Array{Float64,1}})
    precompile(Tuple{typeof(Base.Broadcast._broadcast_getindex_evalf),typeof(MathOptInterface.add_constraint),MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}},MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.ExponentialCone})
    precompile(Tuple{typeof(Base.Broadcast._broadcast_getindex_evalf),typeof(MathOptInterface.add_constraint),MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}},MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.Nonnegatives})
    precompile(Tuple{typeof(Base.Broadcast._broadcast_getindex_evalf),typeof(MathOptInterface.add_constraint),MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}},MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.PositiveSemidefiniteConeTriangle})
    precompile(Tuple{typeof(Base.Broadcast._broadcast_getindex_evalf),typeof(MathOptInterface.add_constraint),MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}},MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.SecondOrderCone})
    precompile(Tuple{typeof(Base.Broadcast._broadcast_getindex_evalf),typeof(MathOptInterface.add_constraint),MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}},MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.Zeros})

    precompile(Tuple{typeof(COSMO.admm_step!),Array{Float64,1},COSMO.SplitVector{Float64},Array{Float64,1},SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int}},true},SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int}},true},Array{Float64,1},Array{Float64,1},Array{Float64,1},CholmodKKTSolver{Float64,Int},Array{Float64,1},Array{Float64,1},Array{Float64,1},Float64,Float64,Int,Int,COSMO.CompositeConvexSet{Float64}})
    precompile(Tuple{typeof(COSMO.apply_sense!),MathOptInterface.OptimizationSense,SparseArrays.SparseMatrixCSC{Float64,Int},Array{Float64,1},Float64})
    precompile(Tuple{typeof(COSMO.dual_kkt_condition!),Array{Float64,1},Array{Float64,1},SparseArrays.SparseMatrixCSC{Float64,Int},Array{Float64,1},Array{Float64,1},SparseArrays.SparseMatrixCSC{Float64,Int},Array{Float64,1}})
    precompile(Tuple{typeof(COSMO.kkt_col_norms!),SparseArrays.SparseMatrixCSC{Float64,Int},SparseArrays.SparseMatrixCSC{Float64,Int},Array{Float64,1},Array{Float64,1}})
    precompile(Tuple{typeof(COSMO.optimize!),COSMO.Workspace{Float64}})
    precompile(Tuple{typeof(COSMO.pass_attributes!),Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}},MathOptInterface.Utilities.IndexMap,Function})
    precompile(Tuple{typeof(COSMO.pre_allocate_variables!),COSMO.Workspace{Float64}})
    precompile(Tuple{typeof(COSMO.primal_kkt_condition!),Array{Float64,1},SparseArrays.SparseMatrixCSC{Float64,Int},Array{Float64,1},COSMO.SplitVector{Float64},Array{Float64,1}})
    precompile(Tuple{typeof(COSMO.processconstraints!),Tuple{Array{Int,1},Array{Int,1},Array{Float64,1}},Array{Float64,1},Array{COSMO.AbstractConvexSet{Float64},1},Array{Float64,1},Array{Float64,1},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}},MathOptInterface.Utilities.IndexMap,Dict{Int,UnitRange{Int}},Type{MathOptInterface.VectorAffineFunction{Float64}},Type{MathOptInterface.ExponentialCone}})
    precompile(Tuple{typeof(COSMO.processconstraints!),Tuple{Array{Int,1},Array{Int,1},Array{Float64,1}},Array{Float64,1},Array{COSMO.AbstractConvexSet{Float64},1},Array{Float64,1},Array{Float64,1},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}},MathOptInterface.Utilities.IndexMap,Dict{Int,UnitRange{Int}},Type{MathOptInterface.VectorAffineFunction{Float64}},Type{MathOptInterface.Nonnegatives}})
    precompile(Tuple{typeof(COSMO.processconstraints!),Tuple{Array{Int,1},Array{Int,1},Array{Float64,1}},Array{Float64,1},Array{COSMO.AbstractConvexSet{Float64},1},Array{Float64,1},Array{Float64,1},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}},MathOptInterface.Utilities.IndexMap,Dict{Int,UnitRange{Int}},Type{MathOptInterface.VectorAffineFunction{Float64}},Type{MathOptInterface.PositiveSemidefiniteConeTriangle}})
    precompile(Tuple{typeof(COSMO.processconstraints!),Tuple{Array{Int,1},Array{Int,1},Array{Float64,1}},Array{Float64,1},Array{COSMO.AbstractConvexSet{Float64},1},Array{Float64,1},Array{Float64,1},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}},MathOptInterface.Utilities.IndexMap,Dict{Int,UnitRange{Int}},Type{MathOptInterface.VectorAffineFunction{Float64}},Type{MathOptInterface.SecondOrderCone}})
    precompile(Tuple{typeof(COSMO.processconstraints!),Tuple{Array{Int,1},Array{Int,1},Array{Float64,1}},Array{Float64,1},Array{COSMO.AbstractConvexSet{Float64},1},Array{Float64,1},Array{Float64,1},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}},MathOptInterface.Utilities.IndexMap,Dict{Int,UnitRange{Int}},Type{MathOptInterface.VectorAffineFunction{Float64}},Type{MathOptInterface.Zeros}})
    precompile(Tuple{typeof(COSMO.rectify_scaling!),SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int}},true},SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int}},true},COSMO.SecondOrderCone{Float64}})
    precompile(Tuple{typeof(COSMO.scale_data!),SparseArrays.SparseMatrixCSC{Float64,Int},SparseArrays.SparseMatrixCSC{Float64,Int},Array{Float64,1},Array{Float64,1},Diagonal{Float64,Array{Float64,1}},Diagonal{Float64,Array{Float64,1}},Float64})
    precompile(Tuple{typeof(COSMO.scale_data!),SparseArrays.SparseMatrixCSC{Float64,Int},SparseArrays.SparseMatrixCSC{Float64,Int},Array{Float64,1},Array{Float64,1},LinearAlgebra.UniformScaling{Bool},Diagonal{Float64,Array{Float64,1}},Float64})
    precompile(Tuple{typeof(MathOptInterface.Bridges.Objective.add_all_bridges),MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}},Type})
    precompile(Tuple{typeof(MathOptInterface.Bridges.Variable.add_all_bridges),MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}},Type})
    precompile(Tuple{typeof(MathOptInterface.Utilities._pass_attributes),Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}},Bool,MathOptInterface.Utilities.IndexMap,Array{MathOptInterface.AbstractConstraintAttribute,1},Tuple{DataType},Tuple{Array{MathOptInterface.ConstraintIndex{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.ExponentialCone},1}},Tuple{Array{MathOptInterface.ConstraintIndex{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.ExponentialCone},1}},typeof(MathOptInterface.Utilities.load)})
    precompile(Tuple{typeof(MathOptInterface.Utilities._pass_attributes),Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}},Bool,MathOptInterface.Utilities.IndexMap,Array{MathOptInterface.AbstractConstraintAttribute,1},Tuple{DataType},Tuple{Array{MathOptInterface.ConstraintIndex{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.Nonnegatives},1}},Tuple{Array{MathOptInterface.ConstraintIndex{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.Nonnegatives},1}},typeof(MathOptInterface.Utilities.load)})
    precompile(Tuple{typeof(MathOptInterface.Utilities._pass_attributes),Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}},Bool,MathOptInterface.Utilities.IndexMap,Array{MathOptInterface.AbstractConstraintAttribute,1},Tuple{DataType},Tuple{Array{MathOptInterface.ConstraintIndex{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.PositiveSemidefiniteConeTriangle},1}},Tuple{Array{MathOptInterface.ConstraintIndex{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.PositiveSemidefiniteConeTriangle},1}},typeof(MathOptInterface.Utilities.load)})
    precompile(Tuple{typeof(MathOptInterface.Utilities._pass_attributes),Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}},Bool,MathOptInterface.Utilities.IndexMap,Array{MathOptInterface.AbstractConstraintAttribute,1},Tuple{DataType},Tuple{Array{MathOptInterface.ConstraintIndex{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.SecondOrderCone},1}},Tuple{Array{MathOptInterface.ConstraintIndex{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.SecondOrderCone},1}},typeof(MathOptInterface.Utilities.load)})
    precompile(Tuple{typeof(MathOptInterface.Utilities._pass_attributes),Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}},Bool,MathOptInterface.Utilities.IndexMap,Array{MathOptInterface.AbstractConstraintAttribute,1},Tuple{DataType},Tuple{Array{MathOptInterface.ConstraintIndex{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.Zeros},1}},Tuple{Array{MathOptInterface.ConstraintIndex{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.Zeros},1}},typeof(MathOptInterface.Utilities.load)})
    precompile(Tuple{typeof(MathOptInterface.Utilities.load),Optimizer,MathOptInterface.ObjectiveFunction{MathOptInterface.SingleVariable},MathOptInterface.SingleVariable})
    precompile(Tuple{typeof(MathOptInterface.add_variables),MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}},Int})
    precompile(Tuple{typeof(MathOptInterface.get),MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}},MathOptInterface.ConstraintDual,MathOptInterface.ConstraintIndex{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.Nonnegatives}})
    precompile(Tuple{typeof(MathOptInterface.get),MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}},MathOptInterface.ConstraintDual,MathOptInterface.ConstraintIndex{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.Zeros}})
    precompile(Tuple{typeof(MathOptInterface.get),MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}},MathOptInterface.ObjectiveValue})
    precompile(Tuple{typeof(copy),Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Tuple{Base.OneTo{Int}},typeof(MathOptInterface.add_constraint),Tuple{Base.RefValue{MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}}},Array{Union{MathOptInterface.SingleVariable, MathOptInterface.VectorAffineFunction{Float64}},1},Array{Any,1}}}})
    precompile(Tuple{typeof(empty_model!),COSMO.Workspace{Float64}})
    precompile(Tuple{typeof(setproperty!),COSMO.ProblemData{Float64},Symbol,SparseArrays.SparseMatrixCSC{Float64,Int}})
end